-/* $RCSfile$\r
- * $Author: egonw $\r
- * $Date: 2005-11-10 09:52:44 -0600 (Thu, 10 Nov 2005) $\r
- * $Revision: 4255 $\r
- *\r
- * Copyright (C) 2003-2005 The Jmol Development Team\r
- *\r
- * Contact: jmol-developers@lists.sf.net\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it under the terms of the GNU Lesser General Public\r
- * License as published by the Free Software Foundation; either\r
- * version 2.1 of the License, or (at your option) any later version.\r
- *\r
- * This library is distributed in the hope that it will be useful,\r
- * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU\r
- * Lesser General Public License for more details.\r
- *\r
- * You should have received a copy of the GNU Lesser General Public\r
- * License along with this library; if not, write to the Free Software\r
- * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.\r
- */\r
-package javajs.util;\r
-\r
-import javajs.api.EigenInterface;\r
-\r
-import javajs.api.Interface;\r
-\r
-\r
-\r
-\r
-//import org.jmol.script.T;\r
-\r
-final public class Measure {\r
-\r
- public final static float radiansPerDegree = (float) (2 * Math.PI / 360);\r
- \r
- public static float computeAngle(T3 pointA, T3 pointB, T3 pointC, V3 vectorBA, V3 vectorBC, boolean asDegrees) {\r
- vectorBA.sub2(pointA, pointB);\r
- vectorBC.sub2(pointC, pointB);\r
- float angle = vectorBA.angle(vectorBC);\r
- return (asDegrees ? angle / radiansPerDegree : angle);\r
- }\r
-\r
- public static float computeAngleABC(T3 pointA, T3 pointB, T3 pointC, boolean asDegrees) {\r
- V3 vectorBA = new V3();\r
- V3 vectorBC = new V3(); \r
- return computeAngle(pointA, pointB, pointC, vectorBA, vectorBC, asDegrees);\r
- }\r
-\r
- public static float computeTorsion(T3 p1, T3 p2, T3 p3, T3 p4, boolean asDegrees) {\r
- \r
- float ijx = p1.x - p2.x;\r
- float ijy = p1.y - p2.y;\r
- float ijz = p1.z - p2.z;\r
- \r
- float kjx = p3.x - p2.x;\r
- float kjy = p3.y - p2.y;\r
- float kjz = p3.z - p2.z;\r
- \r
- float klx = p3.x - p4.x;\r
- float kly = p3.y - p4.y;\r
- float klz = p3.z - p4.z;\r
- \r
- float ax = ijy * kjz - ijz * kjy;\r
- float ay = ijz * kjx - ijx * kjz;\r
- float az = ijx * kjy - ijy * kjx;\r
- float cx = kjy * klz - kjz * kly;\r
- float cy = kjz * klx - kjx * klz;\r
- float cz = kjx * kly - kjy * klx;\r
- \r
- float ai2 = 1f / (ax * ax + ay * ay + az * az);\r
- float ci2 = 1f / (cx * cx + cy * cy + cz * cz);\r
- \r
- float ai = (float) Math.sqrt(ai2);\r
- float ci = (float) Math.sqrt(ci2);\r
- float denom = ai * ci;\r
- float cross = ax * cx + ay * cy + az * cz;\r
- float cosang = cross * denom;\r
- if (cosang > 1) {\r
- cosang = 1;\r
- }\r
- if (cosang < -1) {\r
- cosang = -1;\r
- }\r
- \r
- float torsion = (float) Math.acos(cosang);\r
- float dot = ijx * cx + ijy * cy + ijz * cz;\r
- float absDot = Math.abs(dot);\r
- torsion = (dot / absDot > 0) ? torsion : -torsion;\r
- return (asDegrees ? torsion / radiansPerDegree : torsion);\r
- }\r
-\r
- /**\r
- * This method calculates measures relating to two points in space \r
- * with related quaternion frame difference. It is used in Jmol for\r
- * calculating straightness and many other helical quantities.\r
- * \r
- * @param a\r
- * @param b\r
- * @param dq\r
- * @return new T3[] { pt_a_prime, n, r, P3.new3(theta, pitch, residuesPerTurn), pt_b_prime };\r
- */\r
- public static T3[] computeHelicalAxis(P3 a, P3 b, Quat dq) {\r
- \r
- // b\r
- // | /|\r
- // | / |\r
- // | / |\r
- // |/ c\r
- // b'+ / \\r
- // | / \ Vcb = Vab . n\r
- // n | / \d Vda = (Vcb - Vab) / 2\r
- // |/theta \\r
- // a'+---------a\r
- // r \r
-\r
- V3 vab = new V3();\r
- vab.sub2(b, a);\r
- /*\r
- * testing here to see if directing the normal makes any difference -- oddly\r
- * enough, it does not. When n = -n and theta = -theta vab.n is reversed,\r
- * and that magnitude is multiplied by n in generating the A'-B' vector.\r
- * \r
- * a negative angle implies a left-handed axis (sheets)\r
- */\r
- float theta = dq.getTheta();\r
- V3 n = dq.getNormal();\r
- float v_dot_n = vab.dot(n);\r
- if (Math.abs(v_dot_n) < 0.0001f)\r
- v_dot_n = 0;\r
- V3 va_prime_d = new V3();\r
- va_prime_d.cross(vab, n);\r
- if (va_prime_d.dot(va_prime_d) != 0)\r
- va_prime_d.normalize();\r
- V3 vda = new V3();\r
- V3 vcb = V3.newV(n);\r
- if (v_dot_n == 0)\r
- v_dot_n = PT.FLOAT_MIN_SAFE; // allow for perpendicular axis to vab\r
- vcb.scale(v_dot_n);\r
- vda.sub2(vcb, vab);\r
- vda.scale(0.5f);\r
- va_prime_d.scale(theta == 0 ? 0 : (float) (vda.length() / Math.tan(theta\r
- / 2 / 180 * Math.PI)));\r
- V3 r = V3.newV(va_prime_d);\r
- if (theta != 0)\r
- r.add(vda);\r
- P3 pt_a_prime = P3.newP(a);\r
- pt_a_prime.sub(r);\r
- // already done this. ??\r
- if (v_dot_n != PT.FLOAT_MIN_SAFE)\r
- n.scale(v_dot_n);\r
- // must calculate directed angle:\r
- P3 pt_b_prime = P3.newP(pt_a_prime);\r
- pt_b_prime.add(n);\r
- theta = computeTorsion(a, pt_a_prime, pt_b_prime, b, true);\r
- if (Float.isNaN(theta) || r.length() < 0.0001f)\r
- theta = dq.getThetaDirectedV(n); // allow for r = 0\r
- // anything else is an array\r
- float residuesPerTurn = Math.abs(theta == 0 ? 0 : 360f / theta);\r
- float pitch = Math.abs(v_dot_n == PT.FLOAT_MIN_SAFE ? 0 : n.length()\r
- * (theta == 0 ? 1 : 360f / theta));\r
- return new T3[] { pt_a_prime, n, r, P3.new3(theta, pitch, residuesPerTurn), pt_b_prime };\r
- }\r
-\r
- public static P4 getPlaneThroughPoints(T3 pointA,\r
- T3 pointB,\r
- T3 pointC, V3 vNorm,\r
- V3 vAB, P4 plane) {\r
- float w = getNormalThroughPoints(pointA, pointB, pointC, vNorm, vAB);\r
- plane.set4(vNorm.x, vNorm.y, vNorm.z, w);\r
- return plane;\r
- }\r
- \r
- public static void getPlaneThroughPoint(T3 pt, V3 normal, P4 plane) {\r
- plane.set4(normal.x, normal.y, normal.z, -normal.dot(pt));\r
- }\r
- \r
- public static float distanceToPlane(P4 plane, T3 pt) {\r
- return (plane == null ? Float.NaN \r
- : (plane.dot(pt) + plane.w) / (float) Math.sqrt(plane.dot(plane)));\r
- }\r
-\r
- public static float directedDistanceToPlane(P3 pt, P4 plane, P3 ptref) {\r
- float f = plane.dot(pt) + plane.w;\r
- float f1 = plane.dot(ptref) + plane.w;\r
- return Math.signum(f1) * f / (float) Math.sqrt(plane.dot(plane));\r
- }\r
-\r
- public static float distanceToPlaneD(P4 plane, float d, P3 pt) {\r
- return (plane == null ? Float.NaN : (plane.dot(pt) + plane.w) / d);\r
- }\r
-\r
- public static float distanceToPlaneV(V3 norm, float w, P3 pt) {\r
- return (norm == null ? Float.NaN \r
- : (norm.dot(pt) + w) / (float) Math.sqrt(norm.dot(norm)));\r
- }\r
-\r
- /**\r
- * note that if vAB or vAC is dispensible, vNormNorm can be one of them\r
- * @param pointA\r
- * @param pointB\r
- * @param pointC\r
- * @param vNormNorm\r
- * @param vAB\r
- */\r
- public static void calcNormalizedNormal(T3 pointA, T3 pointB,\r
- T3 pointC, V3 vNormNorm, V3 vAB) {\r
- vAB.sub2(pointB, pointA);\r
- vNormNorm.sub2(pointC, pointA);\r
- vNormNorm.cross(vAB, vNormNorm);\r
- vNormNorm.normalize();\r
- }\r
-\r
- public static float getDirectedNormalThroughPoints(T3 pointA, \r
- T3 pointB, T3 pointC, T3 ptRef, V3 vNorm, \r
- V3 vAB) {\r
- // for x = plane({atomno=1}, {atomno=2}, {atomno=3}, {atomno=4})\r
- float nd = getNormalThroughPoints(pointA, pointB, pointC, vNorm, vAB);\r
- if (ptRef != null) {\r
- P3 pt0 = P3.newP(pointA);\r
- pt0.add(vNorm);\r
- float d = pt0.distance(ptRef);\r
- pt0.sub2(pointA, vNorm);\r
- if (d > pt0.distance(ptRef)) {\r
- vNorm.scale(-1);\r
- nd = -nd;\r
- }\r
- }\r
- return nd;\r
- }\r
- \r
- /**\r
- * if vAC is dispensible vNorm can be vAC\r
- * @param pointA\r
- * @param pointB\r
- * @param pointC\r
- * @param vNorm\r
- * @param vTemp\r
- * @return w\r
- */\r
- public static float getNormalThroughPoints(T3 pointA, T3 pointB,\r
- T3 pointC, V3 vNorm, V3 vTemp) {\r
- // for Polyhedra\r
- calcNormalizedNormal(pointA, pointB, pointC, vNorm, vTemp);\r
- // ax + by + cz + d = 0\r
- // so if a point is in the plane, then N dot X = -d\r
- vTemp.setT(pointA);\r
- return -vTemp.dot(vNorm);\r
- }\r
-\r
- public static void getPlaneProjection(P3 pt, P4 plane, P3 ptProj, V3 vNorm) {\r
- float dist = distanceToPlane(plane, pt);\r
- vNorm.set(plane.x, plane.y, plane.z);\r
- vNorm.normalize();\r
- vNorm.scale(-dist);\r
- ptProj.add2(pt, vNorm);\r
- }\r
-\r
- public final static V3 axisY = V3.new3(0, 1, 0);\r
- \r
- public static void getNormalToLine(P3 pointA, P3 pointB,\r
- V3 vNormNorm) {\r
- // vector in xy plane perpendicular to a line between two points RMH\r
- vNormNorm.sub2(pointA, pointB);\r
- vNormNorm.cross(vNormNorm, axisY);\r
- vNormNorm.normalize();\r
- if (Float.isNaN(vNormNorm.x))\r
- vNormNorm.set(1, 0, 0);\r
- }\r
- \r
- public static void getBisectingPlane(P3 pointA, V3 vAB,\r
- T3 ptTemp, V3 vTemp, P4 plane) {\r
- ptTemp.scaleAdd2(0.5f, vAB, pointA);\r
- vTemp.setT(vAB);\r
- vTemp.normalize();\r
- getPlaneThroughPoint(ptTemp, vTemp, plane);\r
- }\r
- \r
- public static void projectOntoAxis(P3 point, P3 axisA,\r
- V3 axisUnitVector,\r
- V3 vectorProjection) {\r
- vectorProjection.sub2(point, axisA);\r
- float projectedLength = vectorProjection.dot(axisUnitVector);\r
- point.scaleAdd2(projectedLength, axisUnitVector, axisA);\r
- vectorProjection.sub2(point, axisA);\r
- }\r
- \r
- public static void calcBestAxisThroughPoints(P3[] points, P3 axisA,\r
- V3 axisUnitVector,\r
- V3 vectorProjection,\r
- int nTriesMax) {\r
- // just a crude starting point.\r
-\r
- int nPoints = points.length;\r
- axisA.setT(points[0]);\r
- axisUnitVector.sub2(points[nPoints - 1], axisA);\r
- axisUnitVector.normalize();\r
-\r
- /*\r
- * We now calculate the least-squares 3D axis\r
- * through the helix alpha carbons starting with Vo\r
- * as a first approximation.\r
- * \r
- * This uses the simple 0-centered least squares fit:\r
- * \r
- * Y = M cross Xi\r
- * \r
- * minimizing R^2 = SUM(|Y - Yi|^2) \r
- * \r
- * where Yi is the vector PERPENDICULAR of the point onto axis Vo\r
- * and Xi is the vector PROJECTION of the point onto axis Vo\r
- * and M is a vector adjustment \r
- * \r
- * M = SUM_(Xi cross Yi) / sum(|Xi|^2)\r
- * \r
- * from which we arrive at:\r
- * \r
- * V = Vo + (M cross Vo)\r
- * \r
- * Basically, this is just a 3D version of a \r
- * standard 2D least squares fit to a line, where we would say:\r
- * \r
- * y = m xi + b\r
- * \r
- * D = n (sum xi^2) - (sum xi)^2\r
- * \r
- * m = [(n sum xiyi) - (sum xi)(sum yi)] / D\r
- * b = [(sum yi) (sum xi^2) - (sum xi)(sum xiyi)] / D\r
- * \r
- * but here we demand that the line go through the center, so we\r
- * require (sum xi) = (sum yi) = 0, so b = 0 and\r
- * \r
- * m = (sum xiyi) / (sum xi^2)\r
- * \r
- * In 3D we do the same but \r
- * instead of x we have Vo,\r
- * instead of multiplication we use cross products\r
- * \r
- * A bit of iteration is necessary.\r
- * \r
- * Bob Hanson 11/2006\r
- * \r
- */\r
-\r
- calcAveragePointN(points, nPoints, axisA);\r
-\r
- int nTries = 0;\r
- while (nTries++ < nTriesMax\r
- && findAxis(points, nPoints, axisA, axisUnitVector, vectorProjection) > 0.001) {\r
- }\r
-\r
- /*\r
- * Iteration here gets the job done.\r
- * We now find the projections of the endpoints onto the axis\r
- * \r
- */\r
-\r
- P3 tempA = P3.newP(points[0]);\r
- projectOntoAxis(tempA, axisA, axisUnitVector, vectorProjection);\r
- axisA.setT(tempA);\r
- }\r
-\r
- public static float findAxis(P3[] points, int nPoints, P3 axisA,\r
- V3 axisUnitVector, V3 vectorProjection) {\r
- V3 sumXiYi = new V3();\r
- V3 vTemp = new V3();\r
- P3 pt = new P3();\r
- P3 ptProj = new P3();\r
- V3 a = V3.newV(axisUnitVector);\r
-\r
- float sum_Xi2 = 0;\r
- for (int i = nPoints; --i >= 0;) {\r
- pt.setT(points[i]);\r
- ptProj.setT(pt);\r
- projectOntoAxis(ptProj, axisA, axisUnitVector,\r
- vectorProjection);\r
- vTemp.sub2(pt, ptProj);\r
- //sum_Yi2 += vTemp.lengthSquared();\r
- vTemp.cross(vectorProjection, vTemp);\r
- sumXiYi.add(vTemp);\r
- sum_Xi2 += vectorProjection.lengthSquared();\r
- }\r
- V3 m = V3.newV(sumXiYi);\r
- m.scale(1 / sum_Xi2);\r
- vTemp.cross(m, axisUnitVector);\r
- axisUnitVector.add(vTemp);\r
- axisUnitVector.normalize(); \r
- //check for change in direction by measuring vector difference length\r
- vTemp.sub2(axisUnitVector, a);\r
- return vTemp.length();\r
- }\r
- \r
- \r
- public static void calcAveragePoint(P3 pointA, P3 pointB,\r
- P3 pointC) {\r
- pointC.set((pointA.x + pointB.x) / 2, (pointA.y + pointB.y) / 2,\r
- (pointA.z + pointB.z) / 2);\r
- }\r
- \r
- public static void calcAveragePointN(P3[] points, int nPoints,\r
- P3 averagePoint) {\r
- averagePoint.setT(points[0]);\r
- for (int i = 1; i < nPoints; i++)\r
- averagePoint.add(points[i]);\r
- averagePoint.scale(1f / nPoints);\r
- }\r
-\r
- public static Lst<P3> transformPoints(Lst<P3> vPts, M4 m4, P3 center) {\r
- Lst<P3> v = new Lst<P3>();\r
- for (int i = 0; i < vPts.size(); i++) {\r
- P3 pt = P3.newP(vPts.get(i));\r
- pt.sub(center);\r
- m4.rotTrans(pt);\r
- pt.add(center);\r
- v.addLast(pt);\r
- }\r
- return v;\r
- }\r
-\r
- public static boolean isInTetrahedron(P3 pt, P3 ptA, P3 ptB,\r
- P3 ptC, P3 ptD,\r
- P4 plane, V3 vTemp,\r
- V3 vTemp2, boolean fullyEnclosed) {\r
- boolean b = (distanceToPlane(getPlaneThroughPoints(ptC, ptD, ptA, vTemp, vTemp2, plane), pt) >= 0);\r
- if (b != (distanceToPlane(getPlaneThroughPoints(ptA, ptD, ptB, vTemp, vTemp2, plane), pt) >= 0))\r
- return false;\r
- if (b != (distanceToPlane(getPlaneThroughPoints(ptB, ptD, ptC, vTemp, vTemp2, plane), pt) >= 0))\r
- return false;\r
- float d = distanceToPlane(getPlaneThroughPoints(ptA, ptB, ptC, vTemp, vTemp2, plane), pt);\r
- if (fullyEnclosed)\r
- return (b == (d >= 0));\r
- float d1 = distanceToPlane(plane, ptD);\r
- return d1 * d <= 0 || Math.abs(d1) > Math.abs(d);\r
- }\r
-\r
-\r
- /**\r
- * \r
- * @param plane1\r
- * @param plane2\r
- * @return [ point, vector ] or []\r
- */\r
- public static Lst<Object> getIntersectionPP(P4 plane1, P4 plane2) {\r
- float a1 = plane1.x;\r
- float b1 = plane1.y;\r
- float c1 = plane1.z;\r
- float d1 = plane1.w;\r
- float a2 = plane2.x;\r
- float b2 = plane2.y;\r
- float c2 = plane2.z;\r
- float d2 = plane2.w;\r
- V3 norm1 = V3.new3(a1, b1, c1);\r
- V3 norm2 = V3.new3(a2, b2, c2);\r
- V3 nxn = new V3();\r
- nxn.cross(norm1, norm2);\r
- float ax = Math.abs(nxn.x);\r
- float ay = Math.abs(nxn.y);\r
- float az = Math.abs(nxn.z);\r
- float x, y, z, diff;\r
- int type = (ax > ay ? (ax > az ? 1 : 3) : ay > az ? 2 : 3);\r
- switch(type) {\r
- case 1:\r
- x = 0;\r
- diff = (b1 * c2 - b2 * c1);\r
- if (Math.abs(diff) < 0.01) return null;\r
- y = (c1 * d2 - c2 * d1) / diff;\r
- z = (b2 * d1 - d2 * b1) / diff;\r
- break;\r
- case 2:\r
- diff = (a1 * c2 - a2 * c1);\r
- if (Math.abs(diff) < 0.01) return null;\r
- x = (c1 * d2 - c2 * d1) / diff;\r
- y = 0;\r
- z = (a2 * d1 - d2 * a1) / diff;\r
- break;\r
- case 3:\r
- default:\r
- diff = (a1 * b2 - a2 * b1);\r
- if (Math.abs(diff) < 0.01) return null;\r
- x = (b1 * d2 - b2 * d1) / diff;\r
- y = (a2 * d1 - d2 * a1) / diff;\r
- z = 0;\r
- }\r
- Lst<Object>list = new Lst<Object>();\r
- list.addLast(P3.new3(x, y, z));\r
- nxn.normalize();\r
- list.addLast(nxn);\r
- return list;\r
- }\r
-\r
- /**\r
- * \r
- * @param pt1 point on line\r
- * @param v unit vector of line\r
- * @param plane \r
- * @param ptRet point of intersection of line with plane\r
- * @param tempNorm \r
- * @param vTemp \r
- * @return ptRtet\r
- */\r
- public static P3 getIntersection(P3 pt1, V3 v,\r
- P4 plane, P3 ptRet, V3 tempNorm, V3 vTemp) {\r
- getPlaneProjection(pt1, plane, ptRet, tempNorm);\r
- tempNorm.set(plane.x, plane.y, plane.z);\r
- tempNorm.normalize();\r
- if (v == null)\r
- v = V3.newV(tempNorm);\r
- float l_dot_n = v.dot(tempNorm);\r
- if (Math.abs(l_dot_n) < 0.01) return null;\r
- vTemp.sub2(ptRet, pt1);\r
- ptRet.scaleAdd2(vTemp.dot(tempNorm) / l_dot_n, v, pt1);\r
- return ptRet;\r
- }\r
-\r
- /*\r
- public static Point3f getTriangleIntersection(Point3f a1, Point3f a2,\r
- Point3f a3, Point4f plane,\r
- Point3f b1,\r
- Point3f b2, Point3f b3,\r
- Vector3f vNorm, Vector3f vTemp, \r
- Point3f ptRet, Point3f ptTemp, Vector3f vTemp2, Point4f pTemp, Vector3f vTemp3) {\r
- \r
- if (getTriangleIntersection(b1, b2, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))\r
- return ptRet;\r
- if (getTriangleIntersection(b2, b3, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))\r
- return ptRet;\r
- if (getTriangleIntersection(b3, b1, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))\r
- return ptRet;\r
- return null;\r
- }\r
- */\r
- /* \r
- public static boolean getTriangleIntersection(Point3f b1, Point3f b2,\r
- Point3f a1, Point3f a2,\r
- Point3f a3, Vector3f vTemp,\r
- Point4f plane, Vector3f vNorm,\r
- Vector3f vTemp2, Vector3f vTemp3,\r
- Point3f ptRet,\r
- Point3f ptTemp) {\r
- if (distanceToPlane(plane, b1) * distanceToPlane(plane, b2) >= 0)\r
- return false;\r
- vTemp.sub(b2, b1);\r
- vTemp.normalize();\r
- if (getIntersection(b1, vTemp, plane, ptRet, vNorm, vTemp2) != null) {\r
- if (isInTriangle(ptRet, a1, a2, a3, vTemp, vTemp2, vTemp3))\r
- return true;\r
- }\r
- return false;\r
- }\r
- private static boolean isInTriangle(Point3f p, Point3f a, Point3f b,\r
- Point3f c, Vector3f v0, Vector3f v1,\r
- Vector3f v2) {\r
- // from http://www.blackpawn.com/texts/pointinpoly/default.html\r
- // Compute barycentric coordinates\r
- v0.sub(c, a);\r
- v1.sub(b, a);\r
- v2.sub(p, a);\r
- float dot00 = v0.dot(v0);\r
- float dot01 = v0.dot(v1);\r
- float dot02 = v0.dot(v2);\r
- float dot11 = v1.dot(v1);\r
- float dot12 = v1.dot(v2);\r
- float invDenom = 1 / (dot00 * dot11 - dot01 * dot01);\r
- float u = (dot11 * dot02 - dot01 * dot12) * invDenom;\r
- float v = (dot00 * dot12 - dot01 * dot02) * invDenom;\r
- return (u > 0 && v > 0 && u + v < 1);\r
- }\r
- */\r
-\r
- /**\r
- * Closed-form solution of absolute orientation requiring 1:1 mapping of\r
- * positions.\r
- * \r
- * @param centerAndPoints\r
- * @param retStddev\r
- * @return unit quaternion representation rotation\r
- * \r
- * @author hansonr Bob Hanson\r
- * \r
- */\r
- public static Quat calculateQuaternionRotation(P3[][] centerAndPoints,\r
- float[] retStddev) {\r
-\r
- retStddev[1] = Float.NaN;\r
- Quat q = new Quat();\r
- if (centerAndPoints[0].length == 1\r
- || centerAndPoints[0].length != centerAndPoints[1].length)\r
- return q;\r
-\r
- /*\r
- * see Berthold K. P. Horn,\r
- * "Closed-form solution of absolute orientation using unit quaternions" J.\r
- * Opt. Soc. Amer. A, 1987, Vol. 4, pp. 629-642\r
- * http://www.opticsinfobase.org/viewmedia.cfm?uri=josaa-4-4-629&seq=0\r
- * \r
- * \r
- * A similar treatment was developed independently (and later!) \r
- * by G. Kramer, in G. R. Kramer,\r
- * "Superposition of Molecular Structures Using Quaternions"\r
- * Molecular Simulation, 1991, Vol. 7, pp. 113-119. \r
- * \r
- * In that treatment there is a lot of unnecessary calculation \r
- * along the trace of matrix M (eqn 20). \r
- * I'm not sure why the extra x^2 + y^2 + z^2 + x'^2 + y'^2 + z'^2\r
- * is in there, but they are unnecessary and only contribute to larger\r
- * numerical averaging errors and additional processing time, as far as\r
- * I can tell. Adding aI, where a is a scalar and I is the 4x4 identity\r
- * just offsets the eigenvalues but doesn't change the eigenvectors.\r
- * \r
- * and Lydia E. Kavraki, "Molecular Distance Measures"\r
- * http://cnx.org/content/m11608/latest/\r
- * \r
- */\r
-\r
- int n = centerAndPoints[0].length - 1;\r
- if (n < 2)\r
- return q;\r
-\r
- double Sxx = 0, Sxy = 0, Sxz = 0, Syx = 0, Syy = 0, Syz = 0, Szx = 0, Szy = 0, Szz = 0;\r
- P3 ptA = new P3();\r
- P3 ptB = new P3();\r
- for (int i = n + 1; --i >= 1;) {\r
- P3 aij = centerAndPoints[0][i];\r
- P3 bij = centerAndPoints[1][i];\r
- ptA.sub2(aij, centerAndPoints[0][0]);\r
- ptB.sub2(bij, centerAndPoints[0][1]);\r
- Sxx += (double) ptA.x * (double) ptB.x;\r
- Sxy += (double) ptA.x * (double) ptB.y;\r
- Sxz += (double) ptA.x * (double) ptB.z;\r
- Syx += (double) ptA.y * (double) ptB.x;\r
- Syy += (double) ptA.y * (double) ptB.y;\r
- Syz += (double) ptA.y * (double) ptB.z;\r
- Szx += (double) ptA.z * (double) ptB.x;\r
- Szy += (double) ptA.z * (double) ptB.y;\r
- Szz += (double) ptA.z * (double) ptB.z;\r
- }\r
- retStddev[0] = getRmsd(centerAndPoints, q);\r
- double[][] N = new double[4][4];\r
- N[0][0] = Sxx + Syy + Szz;\r
- N[0][1] = N[1][0] = Syz - Szy;\r
- N[0][2] = N[2][0] = Szx - Sxz;\r
- N[0][3] = N[3][0] = Sxy - Syx;\r
-\r
- N[1][1] = Sxx - Syy - Szz;\r
- N[1][2] = N[2][1] = Sxy + Syx;\r
- N[1][3] = N[3][1] = Szx + Sxz;\r
-\r
- N[2][2] = -Sxx + Syy - Szz;\r
- N[2][3] = N[3][2] = Syz + Szy;\r
-\r
- N[3][3] = -Sxx - Syy + Szz;\r
-\r
- //this construction prevents JavaScript from requiring preloading of Eigen\r
- \r
- float[] v = ((EigenInterface) Interface.getInterface("javajs.util.Eigen"))\r
- .setM(N).getEigenvectorsFloatTransposed()[3];\r
- q = Quat.newP4(P4.new4(v[1], v[2], v[3], v[0]));\r
- retStddev[1] = getRmsd(centerAndPoints, q);\r
- return q;\r
- }\r
-\r
- /**\r
- * Fills a 4x4 matrix with rotation-translation of mapped points A to B.\r
- * If centerA is null, this is a standard 4x4 rotation-translation matrix;\r
- * otherwise, this 4x4 matrix is a rotation around a vector through the center of ptsA,\r
- * and centerA is filled with that center; \r
- * Prior to Jmol 14.3.12_2014.02.14, when used from the JmolScript compare() function,\r
- * this method returned the second of these options instead of the first.\r
- * \r
- * @param ptsA\r
- * @param ptsB\r
- * @param m 4x4 matrix to be returned \r
- * @param centerA return center of rotation; if null, then standard 4x4 matrix is returned\r
- * @return stdDev\r
- */\r
- public static float getTransformMatrix4(Lst<P3> ptsA, Lst<P3> ptsB, M4 m,\r
- P3 centerA) {\r
- P3[] cptsA = getCenterAndPoints(ptsA);\r
- P3[] cptsB = getCenterAndPoints(ptsB);\r
- float[] retStddev = new float[2];\r
- Quat q = calculateQuaternionRotation(new P3[][] { cptsA, cptsB },\r
- retStddev);\r
- M3 r = q.getMatrix();\r
- if (centerA == null)\r
- r.rotate(cptsA[0]);\r
- else\r
- centerA.setT(cptsA[0]);\r
- V3 t = V3.newVsub(cptsB[0], cptsA[0]);\r
- m.setMV(r, t);\r
- return retStddev[1];\r
- }\r
-\r
- /**\r
- * from a list of points, create an array that includes the center\r
- * point as the first point. This array is used as a starting point for\r
- * a quaternion analysis of superposition.\r
- * \r
- * @param vPts\r
- * @return array of points with first point center\r
- */\r
- public static P3[] getCenterAndPoints(Lst<P3> vPts) {\r
- int n = vPts.size();\r
- P3[] pts = new P3[n + 1];\r
- pts[0] = new P3();\r
- if (n > 0) {\r
- for (int i = 0; i < n; i++) {\r
- pts[0].add(pts[i + 1] = vPts.get(i));\r
- }\r
- pts[0].scale(1f / n);\r
- }\r
- return pts;\r
- }\r
-\r
- public static float getRmsd(P3[][] centerAndPoints, Quat q) {\r
- double sum2 = 0;\r
- P3[] ptsA = centerAndPoints[0];\r
- P3[] ptsB = centerAndPoints[1];\r
- P3 cA = ptsA[0];\r
- P3 cB = ptsB[0];\r
- int n = ptsA.length - 1;\r
- P3 ptAnew = new P3();\r
- \r
- for (int i = n + 1; --i >= 1;) {\r
- ptAnew.sub2(ptsA[i], cA);\r
- q.transform2(ptAnew, ptAnew).add(cB);\r
- sum2 += ptAnew.distanceSquared(ptsB[i]);\r
- }\r
- return (float) Math.sqrt(sum2 / n);\r
- }\r
-\r
-}\r
+/* $RCSfile$
+ * $Author: egonw $
+ * $Date: 2005-11-10 09:52:44 -0600 (Thu, 10 Nov 2005) $
+ * $Revision: 4255 $
+ *
+ * Copyright (C) 2003-2005 The Jmol Development Team
+ *
+ * Contact: jmol-developers@lists.sf.net
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
+ */
+package javajs.util;
+
+import javajs.api.EigenInterface;
+
+import javajs.api.Interface;
+
+
+
+
+//import org.jmol.script.T;
+
+final public class Measure {
+
+ public final static float radiansPerDegree = (float) (2 * Math.PI / 360);
+
+ public static float computeAngle(T3 pointA, T3 pointB, T3 pointC, V3 vectorBA, V3 vectorBC, boolean asDegrees) {
+ vectorBA.sub2(pointA, pointB);
+ vectorBC.sub2(pointC, pointB);
+ float angle = vectorBA.angle(vectorBC);
+ return (asDegrees ? angle / radiansPerDegree : angle);
+ }
+
+ public static float computeAngleABC(T3 pointA, T3 pointB, T3 pointC, boolean asDegrees) {
+ V3 vectorBA = new V3();
+ V3 vectorBC = new V3();
+ return computeAngle(pointA, pointB, pointC, vectorBA, vectorBC, asDegrees);
+ }
+
+ public static float computeTorsion(T3 p1, T3 p2, T3 p3, T3 p4, boolean asDegrees) {
+
+ float ijx = p1.x - p2.x;
+ float ijy = p1.y - p2.y;
+ float ijz = p1.z - p2.z;
+
+ float kjx = p3.x - p2.x;
+ float kjy = p3.y - p2.y;
+ float kjz = p3.z - p2.z;
+
+ float klx = p3.x - p4.x;
+ float kly = p3.y - p4.y;
+ float klz = p3.z - p4.z;
+
+ float ax = ijy * kjz - ijz * kjy;
+ float ay = ijz * kjx - ijx * kjz;
+ float az = ijx * kjy - ijy * kjx;
+ float cx = kjy * klz - kjz * kly;
+ float cy = kjz * klx - kjx * klz;
+ float cz = kjx * kly - kjy * klx;
+
+ float ai2 = 1f / (ax * ax + ay * ay + az * az);
+ float ci2 = 1f / (cx * cx + cy * cy + cz * cz);
+
+ float ai = (float) Math.sqrt(ai2);
+ float ci = (float) Math.sqrt(ci2);
+ float denom = ai * ci;
+ float cross = ax * cx + ay * cy + az * cz;
+ float cosang = cross * denom;
+ if (cosang > 1) {
+ cosang = 1;
+ }
+ if (cosang < -1) {
+ cosang = -1;
+ }
+
+ float torsion = (float) Math.acos(cosang);
+ float dot = ijx * cx + ijy * cy + ijz * cz;
+ float absDot = Math.abs(dot);
+ torsion = (dot / absDot > 0) ? torsion : -torsion;
+ return (asDegrees ? torsion / radiansPerDegree : torsion);
+ }
+
+ /**
+ * This method calculates measures relating to two points in space
+ * with related quaternion frame difference. It is used in Jmol for
+ * calculating straightness and many other helical quantities.
+ *
+ * @param a
+ * @param b
+ * @param dq
+ * @return new T3[] { pt_a_prime, n, r, P3.new3(theta, pitch, residuesPerTurn), pt_b_prime };
+ */
+ public static T3[] computeHelicalAxis(P3 a, P3 b, Quat dq) {
+
+ // b
+ // | /|
+ // | / |
+ // | / |
+ // |/ c
+ // b'+ / \
+ // | / \ Vcb = Vab . n
+ // n | / \d Vda = (Vcb - Vab) / 2
+ // |/theta \
+ // a'+---------a
+ // r
+
+ V3 vab = new V3();
+ vab.sub2(b, a);
+ /*
+ * testing here to see if directing the normal makes any difference -- oddly
+ * enough, it does not. When n = -n and theta = -theta vab.n is reversed,
+ * and that magnitude is multiplied by n in generating the A'-B' vector.
+ *
+ * a negative angle implies a left-handed axis (sheets)
+ */
+ float theta = dq.getTheta();
+ V3 n = dq.getNormal();
+ float v_dot_n = vab.dot(n);
+ if (Math.abs(v_dot_n) < 0.0001f)
+ v_dot_n = 0;
+ V3 va_prime_d = new V3();
+ va_prime_d.cross(vab, n);
+ if (va_prime_d.dot(va_prime_d) != 0)
+ va_prime_d.normalize();
+ V3 vda = new V3();
+ V3 vcb = V3.newV(n);
+ if (v_dot_n == 0)
+ v_dot_n = PT.FLOAT_MIN_SAFE; // allow for perpendicular axis to vab
+ vcb.scale(v_dot_n);
+ vda.sub2(vcb, vab);
+ vda.scale(0.5f);
+ va_prime_d.scale(theta == 0 ? 0 : (float) (vda.length() / Math.tan(theta
+ / 2 / 180 * Math.PI)));
+ V3 r = V3.newV(va_prime_d);
+ if (theta != 0)
+ r.add(vda);
+ P3 pt_a_prime = P3.newP(a);
+ pt_a_prime.sub(r);
+ // already done this. ??
+ if (v_dot_n != PT.FLOAT_MIN_SAFE)
+ n.scale(v_dot_n);
+ // must calculate directed angle:
+ P3 pt_b_prime = P3.newP(pt_a_prime);
+ pt_b_prime.add(n);
+ theta = computeTorsion(a, pt_a_prime, pt_b_prime, b, true);
+ if (Float.isNaN(theta) || r.length() < 0.0001f)
+ theta = dq.getThetaDirectedV(n); // allow for r = 0
+ // anything else is an array
+ float residuesPerTurn = Math.abs(theta == 0 ? 0 : 360f / theta);
+ float pitch = Math.abs(v_dot_n == PT.FLOAT_MIN_SAFE ? 0 : n.length()
+ * (theta == 0 ? 1 : 360f / theta));
+ return new T3[] { pt_a_prime, n, r, P3.new3(theta, pitch, residuesPerTurn), pt_b_prime };
+ }
+
+ public static P4 getPlaneThroughPoints(T3 pointA,
+ T3 pointB,
+ T3 pointC, V3 vNorm,
+ V3 vAB, P4 plane) {
+ float w = getNormalThroughPoints(pointA, pointB, pointC, vNorm, vAB);
+ plane.set4(vNorm.x, vNorm.y, vNorm.z, w);
+ return plane;
+ }
+
+ public static void getPlaneThroughPoint(T3 pt, V3 normal, P4 plane) {
+ plane.set4(normal.x, normal.y, normal.z, -normal.dot(pt));
+ }
+
+ public static float distanceToPlane(P4 plane, T3 pt) {
+ return (plane == null ? Float.NaN
+ : (plane.dot(pt) + plane.w) / (float) Math.sqrt(plane.dot(plane)));
+ }
+
+ public static float directedDistanceToPlane(P3 pt, P4 plane, P3 ptref) {
+ float f = plane.dot(pt) + plane.w;
+ float f1 = plane.dot(ptref) + plane.w;
+ return Math.signum(f1) * f / (float) Math.sqrt(plane.dot(plane));
+ }
+
+ public static float distanceToPlaneD(P4 plane, float d, P3 pt) {
+ return (plane == null ? Float.NaN : (plane.dot(pt) + plane.w) / d);
+ }
+
+ public static float distanceToPlaneV(V3 norm, float w, P3 pt) {
+ return (norm == null ? Float.NaN
+ : (norm.dot(pt) + w) / (float) Math.sqrt(norm.dot(norm)));
+ }
+
+ /**
+ * note that if vAB or vAC is dispensible, vNormNorm can be one of them
+ * @param pointA
+ * @param pointB
+ * @param pointC
+ * @param vNormNorm
+ * @param vAB
+ */
+ public static void calcNormalizedNormal(T3 pointA, T3 pointB,
+ T3 pointC, V3 vNormNorm, V3 vAB) {
+ vAB.sub2(pointB, pointA);
+ vNormNorm.sub2(pointC, pointA);
+ vNormNorm.cross(vAB, vNormNorm);
+ vNormNorm.normalize();
+ }
+
+ public static float getDirectedNormalThroughPoints(T3 pointA,
+ T3 pointB, T3 pointC, T3 ptRef, V3 vNorm,
+ V3 vAB) {
+ // for x = plane({atomno=1}, {atomno=2}, {atomno=3}, {atomno=4})
+ float nd = getNormalThroughPoints(pointA, pointB, pointC, vNorm, vAB);
+ if (ptRef != null) {
+ P3 pt0 = P3.newP(pointA);
+ pt0.add(vNorm);
+ float d = pt0.distance(ptRef);
+ pt0.sub2(pointA, vNorm);
+ if (d > pt0.distance(ptRef)) {
+ vNorm.scale(-1);
+ nd = -nd;
+ }
+ }
+ return nd;
+ }
+
+ /**
+ * if vAC is dispensible vNorm can be vAC
+ * @param pointA
+ * @param pointB
+ * @param pointC
+ * @param vNorm
+ * @param vTemp
+ * @return w
+ */
+ public static float getNormalThroughPoints(T3 pointA, T3 pointB,
+ T3 pointC, V3 vNorm, V3 vTemp) {
+ // for Polyhedra
+ calcNormalizedNormal(pointA, pointB, pointC, vNorm, vTemp);
+ // ax + by + cz + d = 0
+ // so if a point is in the plane, then N dot X = -d
+ vTemp.setT(pointA);
+ return -vTemp.dot(vNorm);
+ }
+
+ public static void getPlaneProjection(P3 pt, P4 plane, P3 ptProj, V3 vNorm) {
+ float dist = distanceToPlane(plane, pt);
+ vNorm.set(plane.x, plane.y, plane.z);
+ vNorm.normalize();
+ vNorm.scale(-dist);
+ ptProj.add2(pt, vNorm);
+ }
+
+ public final static V3 axisY = V3.new3(0, 1, 0);
+
+ public static void getNormalToLine(P3 pointA, P3 pointB,
+ V3 vNormNorm) {
+ // vector in xy plane perpendicular to a line between two points RMH
+ vNormNorm.sub2(pointA, pointB);
+ vNormNorm.cross(vNormNorm, axisY);
+ vNormNorm.normalize();
+ if (Float.isNaN(vNormNorm.x))
+ vNormNorm.set(1, 0, 0);
+ }
+
+ public static void getBisectingPlane(P3 pointA, V3 vAB,
+ T3 ptTemp, V3 vTemp, P4 plane) {
+ ptTemp.scaleAdd2(0.5f, vAB, pointA);
+ vTemp.setT(vAB);
+ vTemp.normalize();
+ getPlaneThroughPoint(ptTemp, vTemp, plane);
+ }
+
+ public static void projectOntoAxis(P3 point, P3 axisA,
+ V3 axisUnitVector,
+ V3 vectorProjection) {
+ vectorProjection.sub2(point, axisA);
+ float projectedLength = vectorProjection.dot(axisUnitVector);
+ point.scaleAdd2(projectedLength, axisUnitVector, axisA);
+ vectorProjection.sub2(point, axisA);
+ }
+
+ public static void calcBestAxisThroughPoints(P3[] points, P3 axisA,
+ V3 axisUnitVector,
+ V3 vectorProjection,
+ int nTriesMax) {
+ // just a crude starting point.
+
+ int nPoints = points.length;
+ axisA.setT(points[0]);
+ axisUnitVector.sub2(points[nPoints - 1], axisA);
+ axisUnitVector.normalize();
+
+ /*
+ * We now calculate the least-squares 3D axis
+ * through the helix alpha carbons starting with Vo
+ * as a first approximation.
+ *
+ * This uses the simple 0-centered least squares fit:
+ *
+ * Y = M cross Xi
+ *
+ * minimizing R^2 = SUM(|Y - Yi|^2)
+ *
+ * where Yi is the vector PERPENDICULAR of the point onto axis Vo
+ * and Xi is the vector PROJECTION of the point onto axis Vo
+ * and M is a vector adjustment
+ *
+ * M = SUM_(Xi cross Yi) / sum(|Xi|^2)
+ *
+ * from which we arrive at:
+ *
+ * V = Vo + (M cross Vo)
+ *
+ * Basically, this is just a 3D version of a
+ * standard 2D least squares fit to a line, where we would say:
+ *
+ * y = m xi + b
+ *
+ * D = n (sum xi^2) - (sum xi)^2
+ *
+ * m = [(n sum xiyi) - (sum xi)(sum yi)] / D
+ * b = [(sum yi) (sum xi^2) - (sum xi)(sum xiyi)] / D
+ *
+ * but here we demand that the line go through the center, so we
+ * require (sum xi) = (sum yi) = 0, so b = 0 and
+ *
+ * m = (sum xiyi) / (sum xi^2)
+ *
+ * In 3D we do the same but
+ * instead of x we have Vo,
+ * instead of multiplication we use cross products
+ *
+ * A bit of iteration is necessary.
+ *
+ * Bob Hanson 11/2006
+ *
+ */
+
+ calcAveragePointN(points, nPoints, axisA);
+
+ int nTries = 0;
+ while (nTries++ < nTriesMax
+ && findAxis(points, nPoints, axisA, axisUnitVector, vectorProjection) > 0.001) {
+ }
+
+ /*
+ * Iteration here gets the job done.
+ * We now find the projections of the endpoints onto the axis
+ *
+ */
+
+ P3 tempA = P3.newP(points[0]);
+ projectOntoAxis(tempA, axisA, axisUnitVector, vectorProjection);
+ axisA.setT(tempA);
+ }
+
+ public static float findAxis(P3[] points, int nPoints, P3 axisA,
+ V3 axisUnitVector, V3 vectorProjection) {
+ V3 sumXiYi = new V3();
+ V3 vTemp = new V3();
+ P3 pt = new P3();
+ P3 ptProj = new P3();
+ V3 a = V3.newV(axisUnitVector);
+
+ float sum_Xi2 = 0;
+ for (int i = nPoints; --i >= 0;) {
+ pt.setT(points[i]);
+ ptProj.setT(pt);
+ projectOntoAxis(ptProj, axisA, axisUnitVector,
+ vectorProjection);
+ vTemp.sub2(pt, ptProj);
+ //sum_Yi2 += vTemp.lengthSquared();
+ vTemp.cross(vectorProjection, vTemp);
+ sumXiYi.add(vTemp);
+ sum_Xi2 += vectorProjection.lengthSquared();
+ }
+ V3 m = V3.newV(sumXiYi);
+ m.scale(1 / sum_Xi2);
+ vTemp.cross(m, axisUnitVector);
+ axisUnitVector.add(vTemp);
+ axisUnitVector.normalize();
+ //check for change in direction by measuring vector difference length
+ vTemp.sub2(axisUnitVector, a);
+ return vTemp.length();
+ }
+
+
+ public static void calcAveragePoint(P3 pointA, P3 pointB,
+ P3 pointC) {
+ pointC.set((pointA.x + pointB.x) / 2, (pointA.y + pointB.y) / 2,
+ (pointA.z + pointB.z) / 2);
+ }
+
+ public static void calcAveragePointN(P3[] points, int nPoints,
+ P3 averagePoint) {
+ averagePoint.setT(points[0]);
+ for (int i = 1; i < nPoints; i++)
+ averagePoint.add(points[i]);
+ averagePoint.scale(1f / nPoints);
+ }
+
+ public static Lst<P3> transformPoints(Lst<P3> vPts, M4 m4, P3 center) {
+ Lst<P3> v = new Lst<P3>();
+ for (int i = 0; i < vPts.size(); i++) {
+ P3 pt = P3.newP(vPts.get(i));
+ pt.sub(center);
+ m4.rotTrans(pt);
+ pt.add(center);
+ v.addLast(pt);
+ }
+ return v;
+ }
+
+ public static boolean isInTetrahedron(P3 pt, P3 ptA, P3 ptB,
+ P3 ptC, P3 ptD,
+ P4 plane, V3 vTemp,
+ V3 vTemp2, boolean fullyEnclosed) {
+ boolean b = (distanceToPlane(getPlaneThroughPoints(ptC, ptD, ptA, vTemp, vTemp2, plane), pt) >= 0);
+ if (b != (distanceToPlane(getPlaneThroughPoints(ptA, ptD, ptB, vTemp, vTemp2, plane), pt) >= 0))
+ return false;
+ if (b != (distanceToPlane(getPlaneThroughPoints(ptB, ptD, ptC, vTemp, vTemp2, plane), pt) >= 0))
+ return false;
+ float d = distanceToPlane(getPlaneThroughPoints(ptA, ptB, ptC, vTemp, vTemp2, plane), pt);
+ if (fullyEnclosed)
+ return (b == (d >= 0));
+ float d1 = distanceToPlane(plane, ptD);
+ return d1 * d <= 0 || Math.abs(d1) > Math.abs(d);
+ }
+
+
+ /**
+ *
+ * @param plane1
+ * @param plane2
+ * @return [ point, vector ] or []
+ */
+ public static Lst<Object> getIntersectionPP(P4 plane1, P4 plane2) {
+ float a1 = plane1.x;
+ float b1 = plane1.y;
+ float c1 = plane1.z;
+ float d1 = plane1.w;
+ float a2 = plane2.x;
+ float b2 = plane2.y;
+ float c2 = plane2.z;
+ float d2 = plane2.w;
+ V3 norm1 = V3.new3(a1, b1, c1);
+ V3 norm2 = V3.new3(a2, b2, c2);
+ V3 nxn = new V3();
+ nxn.cross(norm1, norm2);
+ float ax = Math.abs(nxn.x);
+ float ay = Math.abs(nxn.y);
+ float az = Math.abs(nxn.z);
+ float x, y, z, diff;
+ int type = (ax > ay ? (ax > az ? 1 : 3) : ay > az ? 2 : 3);
+ switch(type) {
+ case 1:
+ x = 0;
+ diff = (b1 * c2 - b2 * c1);
+ if (Math.abs(diff) < 0.01) return null;
+ y = (c1 * d2 - c2 * d1) / diff;
+ z = (b2 * d1 - d2 * b1) / diff;
+ break;
+ case 2:
+ diff = (a1 * c2 - a2 * c1);
+ if (Math.abs(diff) < 0.01) return null;
+ x = (c1 * d2 - c2 * d1) / diff;
+ y = 0;
+ z = (a2 * d1 - d2 * a1) / diff;
+ break;
+ case 3:
+ default:
+ diff = (a1 * b2 - a2 * b1);
+ if (Math.abs(diff) < 0.01) return null;
+ x = (b1 * d2 - b2 * d1) / diff;
+ y = (a2 * d1 - d2 * a1) / diff;
+ z = 0;
+ }
+ Lst<Object>list = new Lst<Object>();
+ list.addLast(P3.new3(x, y, z));
+ nxn.normalize();
+ list.addLast(nxn);
+ return list;
+ }
+
+ /**
+ *
+ * @param pt1 point on line
+ * @param v unit vector of line
+ * @param plane
+ * @param ptRet point of intersection of line with plane
+ * @param tempNorm
+ * @param vTemp
+ * @return ptRtet
+ */
+ public static P3 getIntersection(P3 pt1, V3 v,
+ P4 plane, P3 ptRet, V3 tempNorm, V3 vTemp) {
+ getPlaneProjection(pt1, plane, ptRet, tempNorm);
+ tempNorm.set(plane.x, plane.y, plane.z);
+ tempNorm.normalize();
+ if (v == null)
+ v = V3.newV(tempNorm);
+ float l_dot_n = v.dot(tempNorm);
+ if (Math.abs(l_dot_n) < 0.01) return null;
+ vTemp.sub2(ptRet, pt1);
+ ptRet.scaleAdd2(vTemp.dot(tempNorm) / l_dot_n, v, pt1);
+ return ptRet;
+ }
+
+ /*
+ public static Point3f getTriangleIntersection(Point3f a1, Point3f a2,
+ Point3f a3, Point4f plane,
+ Point3f b1,
+ Point3f b2, Point3f b3,
+ Vector3f vNorm, Vector3f vTemp,
+ Point3f ptRet, Point3f ptTemp, Vector3f vTemp2, Point4f pTemp, Vector3f vTemp3) {
+
+ if (getTriangleIntersection(b1, b2, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))
+ return ptRet;
+ if (getTriangleIntersection(b2, b3, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))
+ return ptRet;
+ if (getTriangleIntersection(b3, b1, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))
+ return ptRet;
+ return null;
+ }
+ */
+ /*
+ public static boolean getTriangleIntersection(Point3f b1, Point3f b2,
+ Point3f a1, Point3f a2,
+ Point3f a3, Vector3f vTemp,
+ Point4f plane, Vector3f vNorm,
+ Vector3f vTemp2, Vector3f vTemp3,
+ Point3f ptRet,
+ Point3f ptTemp) {
+ if (distanceToPlane(plane, b1) * distanceToPlane(plane, b2) >= 0)
+ return false;
+ vTemp.sub(b2, b1);
+ vTemp.normalize();
+ if (getIntersection(b1, vTemp, plane, ptRet, vNorm, vTemp2) != null) {
+ if (isInTriangle(ptRet, a1, a2, a3, vTemp, vTemp2, vTemp3))
+ return true;
+ }
+ return false;
+ }
+ private static boolean isInTriangle(Point3f p, Point3f a, Point3f b,
+ Point3f c, Vector3f v0, Vector3f v1,
+ Vector3f v2) {
+ // from http://www.blackpawn.com/texts/pointinpoly/default.html
+ // Compute barycentric coordinates
+ v0.sub(c, a);
+ v1.sub(b, a);
+ v2.sub(p, a);
+ float dot00 = v0.dot(v0);
+ float dot01 = v0.dot(v1);
+ float dot02 = v0.dot(v2);
+ float dot11 = v1.dot(v1);
+ float dot12 = v1.dot(v2);
+ float invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
+ float u = (dot11 * dot02 - dot01 * dot12) * invDenom;
+ float v = (dot00 * dot12 - dot01 * dot02) * invDenom;
+ return (u > 0 && v > 0 && u + v < 1);
+ }
+ */
+
+ /**
+ * Closed-form solution of absolute orientation requiring 1:1 mapping of
+ * positions.
+ *
+ * @param centerAndPoints
+ * @param retStddev
+ * @return unit quaternion representation rotation
+ *
+ * @author hansonr Bob Hanson
+ *
+ */
+ public static Quat calculateQuaternionRotation(P3[][] centerAndPoints,
+ float[] retStddev) {
+
+ retStddev[1] = Float.NaN;
+ Quat q = new Quat();
+ if (centerAndPoints[0].length == 1
+ || centerAndPoints[0].length != centerAndPoints[1].length)
+ return q;
+
+ /*
+ * see Berthold K. P. Horn,
+ * "Closed-form solution of absolute orientation using unit quaternions" J.
+ * Opt. Soc. Amer. A, 1987, Vol. 4, pp. 629-642
+ * http://www.opticsinfobase.org/viewmedia.cfm?uri=josaa-4-4-629&seq=0
+ *
+ *
+ * A similar treatment was developed independently (and later!)
+ * by G. Kramer, in G. R. Kramer,
+ * "Superposition of Molecular Structures Using Quaternions"
+ * Molecular Simulation, 1991, Vol. 7, pp. 113-119.
+ *
+ * In that treatment there is a lot of unnecessary calculation
+ * along the trace of matrix M (eqn 20).
+ * I'm not sure why the extra x^2 + y^2 + z^2 + x'^2 + y'^2 + z'^2
+ * is in there, but they are unnecessary and only contribute to larger
+ * numerical averaging errors and additional processing time, as far as
+ * I can tell. Adding aI, where a is a scalar and I is the 4x4 identity
+ * just offsets the eigenvalues but doesn't change the eigenvectors.
+ *
+ * and Lydia E. Kavraki, "Molecular Distance Measures"
+ * http://cnx.org/content/m11608/latest/
+ *
+ */
+
+ int n = centerAndPoints[0].length - 1;
+ if (n < 2)
+ return q;
+
+ double Sxx = 0, Sxy = 0, Sxz = 0, Syx = 0, Syy = 0, Syz = 0, Szx = 0, Szy = 0, Szz = 0;
+ P3 ptA = new P3();
+ P3 ptB = new P3();
+ for (int i = n + 1; --i >= 1;) {
+ P3 aij = centerAndPoints[0][i];
+ P3 bij = centerAndPoints[1][i];
+ ptA.sub2(aij, centerAndPoints[0][0]);
+ ptB.sub2(bij, centerAndPoints[0][1]);
+ Sxx += (double) ptA.x * (double) ptB.x;
+ Sxy += (double) ptA.x * (double) ptB.y;
+ Sxz += (double) ptA.x * (double) ptB.z;
+ Syx += (double) ptA.y * (double) ptB.x;
+ Syy += (double) ptA.y * (double) ptB.y;
+ Syz += (double) ptA.y * (double) ptB.z;
+ Szx += (double) ptA.z * (double) ptB.x;
+ Szy += (double) ptA.z * (double) ptB.y;
+ Szz += (double) ptA.z * (double) ptB.z;
+ }
+ retStddev[0] = getRmsd(centerAndPoints, q);
+ double[][] N = new double[4][4];
+ N[0][0] = Sxx + Syy + Szz;
+ N[0][1] = N[1][0] = Syz - Szy;
+ N[0][2] = N[2][0] = Szx - Sxz;
+ N[0][3] = N[3][0] = Sxy - Syx;
+
+ N[1][1] = Sxx - Syy - Szz;
+ N[1][2] = N[2][1] = Sxy + Syx;
+ N[1][3] = N[3][1] = Szx + Sxz;
+
+ N[2][2] = -Sxx + Syy - Szz;
+ N[2][3] = N[3][2] = Syz + Szy;
+
+ N[3][3] = -Sxx - Syy + Szz;
+
+ //this construction prevents JavaScript from requiring preloading of Eigen
+
+ float[] v = ((EigenInterface) Interface.getInterface("javajs.util.Eigen"))
+ .setM(N).getEigenvectorsFloatTransposed()[3];
+ q = Quat.newP4(P4.new4(v[1], v[2], v[3], v[0]));
+ retStddev[1] = getRmsd(centerAndPoints, q);
+ return q;
+ }
+
+ /**
+ * Fills a 4x4 matrix with rotation-translation of mapped points A to B.
+ * If centerA is null, this is a standard 4x4 rotation-translation matrix;
+ * otherwise, this 4x4 matrix is a rotation around a vector through the center of ptsA,
+ * and centerA is filled with that center;
+ * Prior to Jmol 14.3.12_2014.02.14, when used from the JmolScript compare() function,
+ * this method returned the second of these options instead of the first.
+ *
+ * @param ptsA
+ * @param ptsB
+ * @param m 4x4 matrix to be returned
+ * @param centerA return center of rotation; if null, then standard 4x4 matrix is returned
+ * @return stdDev
+ */
+ public static float getTransformMatrix4(Lst<P3> ptsA, Lst<P3> ptsB, M4 m,
+ P3 centerA) {
+ P3[] cptsA = getCenterAndPoints(ptsA);
+ P3[] cptsB = getCenterAndPoints(ptsB);
+ float[] retStddev = new float[2];
+ Quat q = calculateQuaternionRotation(new P3[][] { cptsA, cptsB },
+ retStddev);
+ M3 r = q.getMatrix();
+ if (centerA == null)
+ r.rotate(cptsA[0]);
+ else
+ centerA.setT(cptsA[0]);
+ V3 t = V3.newVsub(cptsB[0], cptsA[0]);
+ m.setMV(r, t);
+ return retStddev[1];
+ }
+
+ /**
+ * from a list of points, create an array that includes the center
+ * point as the first point. This array is used as a starting point for
+ * a quaternion analysis of superposition.
+ *
+ * @param vPts
+ * @return array of points with first point center
+ */
+ public static P3[] getCenterAndPoints(Lst<P3> vPts) {
+ int n = vPts.size();
+ P3[] pts = new P3[n + 1];
+ pts[0] = new P3();
+ if (n > 0) {
+ for (int i = 0; i < n; i++) {
+ pts[0].add(pts[i + 1] = vPts.get(i));
+ }
+ pts[0].scale(1f / n);
+ }
+ return pts;
+ }
+
+ public static float getRmsd(P3[][] centerAndPoints, Quat q) {
+ double sum2 = 0;
+ P3[] ptsA = centerAndPoints[0];
+ P3[] ptsB = centerAndPoints[1];
+ P3 cA = ptsA[0];
+ P3 cB = ptsB[0];
+ int n = ptsA.length - 1;
+ P3 ptAnew = new P3();
+
+ for (int i = n + 1; --i >= 1;) {
+ ptAnew.sub2(ptsA[i], cA);
+ q.transform2(ptAnew, ptAnew).add(cB);
+ sum2 += ptAnew.distanceSquared(ptsB[i]);
+ }
+ return (float) Math.sqrt(sum2 / n);
+ }
+
+}
git://source.jalview.org / jalviewjs.git / blobdiff