src/javajs/util/Measure.java

   1 /* $RCSfile$\r
   2  * $Author: egonw $\r
   3  * $Date: 2005-11-10 09:52:44 -0600 (Thu, 10 Nov 2005) $\r
   4  * $Revision: 4255 $\r
   5  *\r
   6  * Copyright (C) 2003-2005  The Jmol Development Team\r
   7  *\r
   8  * Contact: jmol-developers@lists.sf.net\r
   9  *\r
  10  *  This library is free software; you can redistribute it and/or\r
  11  *  modify it under the terms of the GNU Lesser General Public\r
  12  *  License as published by the Free Software Foundation; either\r
  13  *  version 2.1 of the License, or (at your option) any later version.\r
  14  *\r
  15  *  This library is distributed in the hope that it will be useful,\r
  16  *  but WITHOUT ANY WARRANTY; without even the implied warranty of\r
  17  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU\r
  18  *  Lesser General Public License for more details.\r
  19  *\r
  20  *  You should have received a copy of the GNU Lesser General Public\r
  21  *  License along with this library; if not, write to the Free Software\r
  22  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.\r
  23  */\r
  24 package javajs.util;\r
  25 \r
  26 import javajs.api.EigenInterface;\r
  27 \r
  28 import javajs.api.Interface;\r
  29 \r
  30 \r
  31 \r
  32 \r
  33 //import org.jmol.script.T;\r
  34 \r
  35 final public class Measure {\r
  36 \r
  37   public final static float radiansPerDegree = (float) (2 * Math.PI / 360);\r
  38   \r
  39   public static float computeAngle(T3 pointA, T3 pointB, T3 pointC, V3 vectorBA, V3 vectorBC, boolean asDegrees) {\r
  40     vectorBA.sub2(pointA, pointB);\r
  41     vectorBC.sub2(pointC, pointB);\r
  42     float angle = vectorBA.angle(vectorBC);\r
  43     return (asDegrees ? angle / radiansPerDegree : angle);\r
  44   }\r
  45 \r
  46   public static float computeAngleABC(T3 pointA, T3 pointB, T3 pointC, boolean asDegrees) {\r
  47     V3 vectorBA = new V3();\r
  48     V3 vectorBC = new V3();        \r
  49     return computeAngle(pointA, pointB, pointC, vectorBA, vectorBC, asDegrees);\r
  50   }\r
  51 \r
  52   public static float computeTorsion(T3 p1, T3 p2, T3 p3, T3 p4, boolean asDegrees) {\r
  53   \r
  54     float ijx = p1.x - p2.x;\r
  55     float ijy = p1.y - p2.y;\r
  56     float ijz = p1.z - p2.z;\r
  57   \r
  58     float kjx = p3.x - p2.x;\r
  59     float kjy = p3.y - p2.y;\r
  60     float kjz = p3.z - p2.z;\r
  61   \r
  62     float klx = p3.x - p4.x;\r
  63     float kly = p3.y - p4.y;\r
  64     float klz = p3.z - p4.z;\r
  65   \r
  66     float ax = ijy * kjz - ijz * kjy;\r
  67     float ay = ijz * kjx - ijx * kjz;\r
  68     float az = ijx * kjy - ijy * kjx;\r
  69     float cx = kjy * klz - kjz * kly;\r
  70     float cy = kjz * klx - kjx * klz;\r
  71     float cz = kjx * kly - kjy * klx;\r
  72   \r
  73     float ai2 = 1f / (ax * ax + ay * ay + az * az);\r
  74     float ci2 = 1f / (cx * cx + cy * cy + cz * cz);\r
  75   \r
  76     float ai = (float) Math.sqrt(ai2);\r
  77     float ci = (float) Math.sqrt(ci2);\r
  78     float denom = ai * ci;\r
  79     float cross = ax * cx + ay * cy + az * cz;\r
  80     float cosang = cross * denom;\r
  81     if (cosang > 1) {\r
  82       cosang = 1;\r
  83     }\r
  84     if (cosang < -1) {\r
  85       cosang = -1;\r
  86     }\r
  87   \r
  88     float torsion = (float) Math.acos(cosang);\r
  89     float dot = ijx * cx + ijy * cy + ijz * cz;\r
  90     float absDot = Math.abs(dot);\r
  91     torsion = (dot / absDot > 0) ? torsion : -torsion;\r
  92     return (asDegrees ? torsion / radiansPerDegree : torsion);\r
  93   }\r
  94 \r
  95   /**\r
  96    * This method calculates measures relating to two points in space \r
  97    * with related quaternion frame difference. It is used in Jmol for\r
  98    * calculating straightness and many other helical quantities.\r
  99    * \r
 100    * @param a\r
 101    * @param b\r
 102    * @param dq\r
 103    * @return  new T3[] { pt_a_prime, n, r, P3.new3(theta, pitch, residuesPerTurn), pt_b_prime };\r
 104    */\r
 105   public static T3[] computeHelicalAxis(P3 a, P3 b, Quat dq) {\r
 106     \r
 107     //                b\r
 108     //           |   /|\r
 109     //           |  / |\r
 110     //           | /  |\r
 111     //           |/   c\r
 112     //         b'+   / \\r
 113     //           |  /   \      Vcb = Vab . n\r
 114     //         n | /     \d    Vda = (Vcb - Vab) / 2\r
 115     //           |/theta  \\r
 116     //         a'+---------a\r
 117     //                r \r
 118 \r
 119     V3 vab = new V3();\r
 120     vab.sub2(b, a);\r
 121     /*\r
 122      * testing here to see if directing the normal makes any difference -- oddly\r
 123      * enough, it does not. When n = -n and theta = -theta vab.n is reversed,\r
 124      * and that magnitude is multiplied by n in generating the A'-B' vector.\r
 125      * \r
 126      * a negative angle implies a left-handed axis (sheets)\r
 127      */\r
 128     float theta = dq.getTheta();\r
 129     V3 n = dq.getNormal();\r
 130     float v_dot_n = vab.dot(n);\r
 131     if (Math.abs(v_dot_n) < 0.0001f)\r
 132       v_dot_n = 0;\r
 133     V3 va_prime_d = new V3();\r
 134     va_prime_d.cross(vab, n);\r
 135     if (va_prime_d.dot(va_prime_d) != 0)\r
 136       va_prime_d.normalize();\r
 137     V3 vda = new V3();\r
 138     V3 vcb = V3.newV(n);\r
 139     if (v_dot_n == 0)\r
 140       v_dot_n = PT.FLOAT_MIN_SAFE; // allow for perpendicular axis to vab\r
 141     vcb.scale(v_dot_n);\r
 142     vda.sub2(vcb, vab);\r
 143     vda.scale(0.5f);\r
 144     va_prime_d.scale(theta == 0 ? 0 : (float) (vda.length() / Math.tan(theta\r
 145         / 2 / 180 * Math.PI)));\r
 146     V3 r = V3.newV(va_prime_d);\r
 147     if (theta != 0)\r
 148       r.add(vda);\r
 149     P3 pt_a_prime = P3.newP(a);\r
 150     pt_a_prime.sub(r);\r
 151     // already done this. ??\r
 152     if (v_dot_n != PT.FLOAT_MIN_SAFE)\r
 153       n.scale(v_dot_n);\r
 154     // must calculate directed angle:\r
 155     P3 pt_b_prime = P3.newP(pt_a_prime);\r
 156     pt_b_prime.add(n);\r
 157     theta = computeTorsion(a, pt_a_prime, pt_b_prime, b, true);\r
 158     if (Float.isNaN(theta) || r.length() < 0.0001f)\r
 159       theta = dq.getThetaDirectedV(n); // allow for r = 0\r
 160     // anything else is an array\r
 161     float residuesPerTurn = Math.abs(theta == 0 ? 0 : 360f / theta);\r
 162     float pitch = Math.abs(v_dot_n == PT.FLOAT_MIN_SAFE ? 0 : n.length()\r
 163         * (theta == 0 ? 1 : 360f / theta));\r
 164     return new T3[] { pt_a_prime, n, r, P3.new3(theta, pitch, residuesPerTurn), pt_b_prime };\r
 165   }\r
 166 \r
 167   public static P4 getPlaneThroughPoints(T3 pointA,\r
 168                                               T3 pointB,\r
 169                                               T3 pointC, V3 vNorm,\r
 170                                               V3 vAB, P4 plane) {\r
 171     float w = getNormalThroughPoints(pointA, pointB, pointC, vNorm, vAB);\r
 172     plane.set4(vNorm.x, vNorm.y, vNorm.z, w);\r
 173     return plane;\r
 174   }\r
 175   \r
 176   public static void getPlaneThroughPoint(T3 pt, V3 normal, P4 plane) {\r
 177     plane.set4(normal.x, normal.y, normal.z, -normal.dot(pt));\r
 178   }\r
 179   \r
 180   public static float distanceToPlane(P4 plane, T3 pt) {\r
 181     return (plane == null ? Float.NaN \r
 182         : (plane.dot(pt) + plane.w) / (float) Math.sqrt(plane.dot(plane)));\r
 183   }\r
 184 \r
 185   public static float directedDistanceToPlane(P3 pt, P4 plane, P3 ptref) {\r
 186     float f = plane.dot(pt) + plane.w;\r
 187     float f1 = plane.dot(ptref) + plane.w;\r
 188     return Math.signum(f1) * f /  (float) Math.sqrt(plane.dot(plane));\r
 189   }\r
 190 \r
 191   public static float distanceToPlaneD(P4 plane, float d, P3 pt) {\r
 192     return (plane == null ? Float.NaN : (plane.dot(pt) + plane.w) / d);\r
 193   }\r
 194 \r
 195   public static float distanceToPlaneV(V3 norm, float w, P3 pt) {\r
 196     return (norm == null ? Float.NaN \r
 197         : (norm.dot(pt) + w)  / (float) Math.sqrt(norm.dot(norm)));\r
 198   }\r
 199 \r
 200   /**\r
 201    * note that if vAB or vAC is dispensible, vNormNorm can be one of them\r
 202    * @param pointA\r
 203    * @param pointB\r
 204    * @param pointC\r
 205    * @param vNormNorm\r
 206    * @param vAB\r
 207    */\r
 208   public static void calcNormalizedNormal(T3 pointA, T3 pointB,\r
 209          T3 pointC, V3 vNormNorm, V3 vAB) {\r
 210     vAB.sub2(pointB, pointA);\r
 211     vNormNorm.sub2(pointC, pointA);\r
 212     vNormNorm.cross(vAB, vNormNorm);\r
 213     vNormNorm.normalize();\r
 214   }\r
 215 \r
 216   public static float getDirectedNormalThroughPoints(T3 pointA, \r
 217          T3 pointB, T3 pointC, T3 ptRef, V3 vNorm, \r
 218          V3 vAB) {\r
 219     // for x = plane({atomno=1}, {atomno=2}, {atomno=3}, {atomno=4})\r
 220     float nd = getNormalThroughPoints(pointA, pointB, pointC, vNorm, vAB);\r
 221     if (ptRef != null) {\r
 222       P3 pt0 = P3.newP(pointA);\r
 223       pt0.add(vNorm);\r
 224       float d = pt0.distance(ptRef);\r
 225       pt0.sub2(pointA, vNorm);\r
 226       if (d > pt0.distance(ptRef)) {\r
 227         vNorm.scale(-1);\r
 228         nd = -nd;\r
 229       }\r
 230     }\r
 231     return nd;\r
 232   }\r
 233   \r
 234   /**\r
 235    * if vAC is dispensible vNorm can be vAC\r
 236    * @param pointA\r
 237    * @param pointB\r
 238    * @param pointC\r
 239    * @param vNorm\r
 240    * @param vTemp\r
 241    * @return w\r
 242    */\r
 243   public static float getNormalThroughPoints(T3 pointA, T3 pointB,\r
 244                                    T3 pointC, V3 vNorm, V3 vTemp) {\r
 245     // for Polyhedra\r
 246     calcNormalizedNormal(pointA, pointB, pointC, vNorm, vTemp);\r
 247     // ax + by + cz + d = 0\r
 248     // so if a point is in the plane, then N dot X = -d\r
 249     vTemp.setT(pointA);\r
 250     return -vTemp.dot(vNorm);\r
 251   }\r
 252 \r
 253   public static void getPlaneProjection(P3 pt, P4 plane, P3 ptProj, V3 vNorm) {\r
 254     float dist = distanceToPlane(plane, pt);\r
 255     vNorm.set(plane.x, plane.y, plane.z);\r
 256     vNorm.normalize();\r
 257     vNorm.scale(-dist);\r
 258     ptProj.add2(pt, vNorm);\r
 259   }\r
 260 \r
 261   public final static V3 axisY = V3.new3(0, 1, 0);\r
 262   \r
 263   public static void getNormalToLine(P3 pointA, P3 pointB,\r
 264                                    V3 vNormNorm) {\r
 265     // vector in xy plane perpendicular to a line between two points RMH\r
 266     vNormNorm.sub2(pointA, pointB);\r
 267     vNormNorm.cross(vNormNorm, axisY);\r
 268     vNormNorm.normalize();\r
 269     if (Float.isNaN(vNormNorm.x))\r
 270       vNormNorm.set(1, 0, 0);\r
 271   }\r
 272   \r
 273   public static void getBisectingPlane(P3 pointA, V3 vAB,\r
 274                                                  T3 ptTemp, V3 vTemp, P4 plane) {\r
 275     ptTemp.scaleAdd2(0.5f, vAB, pointA);\r
 276     vTemp.setT(vAB);\r
 277     vTemp.normalize();\r
 278     getPlaneThroughPoint(ptTemp, vTemp, plane);\r
 279     }\r
 280     \r
 281   public static void projectOntoAxis(P3 point, P3 axisA,\r
 282                                      V3 axisUnitVector,\r
 283                                      V3 vectorProjection) {\r
 284     vectorProjection.sub2(point, axisA);\r
 285     float projectedLength = vectorProjection.dot(axisUnitVector);\r
 286     point.scaleAdd2(projectedLength, axisUnitVector, axisA);\r
 287     vectorProjection.sub2(point, axisA);\r
 288   }\r
 289   \r
 290   public static void calcBestAxisThroughPoints(P3[] points, P3 axisA,\r
 291                                                V3 axisUnitVector,\r
 292                                                V3 vectorProjection,\r
 293                                                int nTriesMax) {\r
 294     // just a crude starting point.\r
 295 \r
 296     int nPoints = points.length;\r
 297     axisA.setT(points[0]);\r
 298     axisUnitVector.sub2(points[nPoints - 1], axisA);\r
 299     axisUnitVector.normalize();\r
 300 \r
 301     /*\r
 302      * We now calculate the least-squares 3D axis\r
 303      * through the helix alpha carbons starting with Vo\r
 304      * as a first approximation.\r
 305      * \r
 306      * This uses the simple 0-centered least squares fit:\r
 307      * \r
 308      * Y = M cross Xi\r
 309      * \r
 310      * minimizing R^2 = SUM(|Y - Yi|^2) \r
 311      * \r
 312      * where Yi is the vector PERPENDICULAR of the point onto axis Vo\r
 313      * and Xi is the vector PROJECTION of the point onto axis Vo\r
 314      * and M is a vector adjustment \r
 315      * \r
 316      * M = SUM_(Xi cross Yi) / sum(|Xi|^2)\r
 317      * \r
 318      * from which we arrive at:\r
 319      * \r
 320      * V = Vo + (M cross Vo)\r
 321      * \r
 322      * Basically, this is just a 3D version of a \r
 323      * standard 2D least squares fit to a line, where we would say:\r
 324      * \r
 325      * y = m xi + b\r
 326      * \r
 327      * D = n (sum xi^2) - (sum xi)^2\r
 328      * \r
 329      * m = [(n sum xiyi) - (sum xi)(sum yi)] / D\r
 330      * b = [(sum yi) (sum xi^2) - (sum xi)(sum xiyi)] / D\r
 331      * \r
 332      * but here we demand that the line go through the center, so we\r
 333      * require (sum xi) = (sum yi) = 0, so b = 0 and\r
 334      * \r
 335      * m = (sum xiyi) / (sum xi^2)\r
 336      * \r
 337      * In 3D we do the same but \r
 338      * instead of x we have Vo,\r
 339      * instead of multiplication we use cross products\r
 340      * \r
 341      * A bit of iteration is necessary.\r
 342      * \r
 343      * Bob Hanson 11/2006\r
 344      * \r
 345      */\r
 346 \r
 347     calcAveragePointN(points, nPoints, axisA);\r
 348 \r
 349     int nTries = 0;\r
 350     while (nTries++ < nTriesMax\r
 351         && findAxis(points, nPoints, axisA, axisUnitVector, vectorProjection) > 0.001) {\r
 352     }\r
 353 \r
 354     /*\r
 355      * Iteration here gets the job done.\r
 356      * We now find the projections of the endpoints onto the axis\r
 357      * \r
 358      */\r
 359 \r
 360     P3 tempA = P3.newP(points[0]);\r
 361     projectOntoAxis(tempA, axisA, axisUnitVector, vectorProjection);\r
 362     axisA.setT(tempA);\r
 363   }\r
 364 \r
 365   public static float findAxis(P3[] points, int nPoints, P3 axisA,\r
 366                         V3 axisUnitVector, V3 vectorProjection) {\r
 367     V3 sumXiYi = new V3();\r
 368     V3 vTemp = new V3();\r
 369     P3 pt = new P3();\r
 370     P3 ptProj = new P3();\r
 371     V3 a = V3.newV(axisUnitVector);\r
 372 \r
 373     float sum_Xi2 = 0;\r
 374     for (int i = nPoints; --i >= 0;) {\r
 375       pt.setT(points[i]);\r
 376       ptProj.setT(pt);\r
 377       projectOntoAxis(ptProj, axisA, axisUnitVector,\r
 378           vectorProjection);\r
 379       vTemp.sub2(pt, ptProj);\r
 380       //sum_Yi2 += vTemp.lengthSquared();\r
 381       vTemp.cross(vectorProjection, vTemp);\r
 382       sumXiYi.add(vTemp);\r
 383       sum_Xi2 += vectorProjection.lengthSquared();\r
 384     }\r
 385     V3 m = V3.newV(sumXiYi);\r
 386     m.scale(1 / sum_Xi2);\r
 387     vTemp.cross(m, axisUnitVector);\r
 388     axisUnitVector.add(vTemp);\r
 389     axisUnitVector.normalize();  \r
 390     //check for change in direction by measuring vector difference length\r
 391     vTemp.sub2(axisUnitVector, a);\r
 392     return vTemp.length();\r
 393   }\r
 394   \r
 395   \r
 396   public static void calcAveragePoint(P3 pointA, P3 pointB,\r
 397                                       P3 pointC) {\r
 398     pointC.set((pointA.x + pointB.x) / 2, (pointA.y + pointB.y) / 2,\r
 399         (pointA.z + pointB.z) / 2);\r
 400   }\r
 401   \r
 402   public static void calcAveragePointN(P3[] points, int nPoints,\r
 403                                 P3 averagePoint) {\r
 404     averagePoint.setT(points[0]);\r
 405     for (int i = 1; i < nPoints; i++)\r
 406       averagePoint.add(points[i]);\r
 407     averagePoint.scale(1f / nPoints);\r
 408   }\r
 409 \r
 410   public static Lst<P3> transformPoints(Lst<P3> vPts, M4 m4, P3 center) {\r
 411     Lst<P3> v = new  Lst<P3>();\r
 412     for (int i = 0; i < vPts.size(); i++) {\r
 413       P3 pt = P3.newP(vPts.get(i));\r
 414       pt.sub(center);\r
 415       m4.rotTrans(pt);\r
 416       pt.add(center);\r
 417       v.addLast(pt);\r
 418     }\r
 419     return v;\r
 420   }\r
 421 \r
 422   public static boolean isInTetrahedron(P3 pt, P3 ptA, P3 ptB,\r
 423                                         P3 ptC, P3 ptD,\r
 424                                         P4 plane, V3 vTemp,\r
 425                                         V3 vTemp2, boolean fullyEnclosed) {\r
 426     boolean b = (distanceToPlane(getPlaneThroughPoints(ptC, ptD, ptA, vTemp, vTemp2, plane), pt) >= 0);\r
 427     if (b != (distanceToPlane(getPlaneThroughPoints(ptA, ptD, ptB, vTemp, vTemp2, plane), pt) >= 0))\r
 428       return false;\r
 429     if (b != (distanceToPlane(getPlaneThroughPoints(ptB, ptD, ptC, vTemp, vTemp2, plane), pt) >= 0))\r
 430       return false;\r
 431     float d = distanceToPlane(getPlaneThroughPoints(ptA, ptB, ptC, vTemp, vTemp2, plane), pt);\r
 432     if (fullyEnclosed)\r
 433       return (b == (d >= 0));\r
 434     float d1 = distanceToPlane(plane, ptD);\r
 435     return d1 * d <= 0 || Math.abs(d1) > Math.abs(d);\r
 436   }\r
 437 \r
 438 \r
 439   /**\r
 440    * \r
 441    * @param plane1\r
 442    * @param plane2\r
 443    * @return       [ point, vector ] or []\r
 444    */\r
 445   public static Lst<Object> getIntersectionPP(P4 plane1, P4 plane2) {\r
 446     float a1 = plane1.x;\r
 447     float b1 = plane1.y;\r
 448     float c1 = plane1.z;\r
 449     float d1 = plane1.w;\r
 450     float a2 = plane2.x;\r
 451     float b2 = plane2.y;\r
 452     float c2 = plane2.z;\r
 453     float d2 = plane2.w;\r
 454     V3 norm1 = V3.new3(a1, b1, c1);\r
 455     V3 norm2 = V3.new3(a2, b2, c2);\r
 456     V3 nxn = new V3();\r
 457     nxn.cross(norm1, norm2);\r
 458     float ax = Math.abs(nxn.x);\r
 459     float ay = Math.abs(nxn.y);\r
 460     float az = Math.abs(nxn.z);\r
 461     float x, y, z, diff;\r
 462     int type = (ax > ay ? (ax > az ? 1 : 3) : ay > az ? 2 : 3);\r
 463     switch(type) {\r
 464     case 1:\r
 465       x = 0;\r
 466       diff = (b1 * c2 - b2 * c1);\r
 467       if (Math.abs(diff) < 0.01) return null;\r
 468       y = (c1 * d2 - c2 * d1) / diff;\r
 469       z = (b2 * d1 - d2 * b1) / diff;\r
 470       break;\r
 471     case 2:\r
 472       diff = (a1 * c2 - a2 * c1);\r
 473       if (Math.abs(diff) < 0.01) return null;\r
 474       x = (c1 * d2 - c2 * d1) / diff;\r
 475       y = 0;\r
 476       z = (a2 * d1 - d2 * a1) / diff;\r
 477       break;\r
 478     case 3:\r
 479     default:\r
 480       diff = (a1 * b2 - a2 * b1);\r
 481       if (Math.abs(diff) < 0.01) return null;\r
 482       x = (b1 * d2 - b2 * d1) / diff;\r
 483       y = (a2 * d1 - d2 * a1) / diff;\r
 484       z = 0;\r
 485     }\r
 486     Lst<Object>list = new  Lst<Object>();\r
 487     list.addLast(P3.new3(x, y, z));\r
 488     nxn.normalize();\r
 489     list.addLast(nxn);\r
 490     return list;\r
 491   }\r
 492 \r
 493   /**\r
 494    * \r
 495    * @param pt1  point on line\r
 496    * @param v    unit vector of line\r
 497    * @param plane \r
 498    * @param ptRet  point of intersection of line with plane\r
 499    * @param tempNorm \r
 500    * @param vTemp \r
 501    * @return       ptRtet\r
 502    */\r
 503   public static P3 getIntersection(P3 pt1, V3 v,\r
 504                                                P4 plane, P3 ptRet, V3 tempNorm, V3 vTemp) {\r
 505     getPlaneProjection(pt1, plane, ptRet, tempNorm);\r
 506     tempNorm.set(plane.x, plane.y, plane.z);\r
 507     tempNorm.normalize();\r
 508     if (v == null)\r
 509       v = V3.newV(tempNorm);\r
 510     float l_dot_n = v.dot(tempNorm);\r
 511     if (Math.abs(l_dot_n) < 0.01) return null;\r
 512     vTemp.sub2(ptRet, pt1);\r
 513     ptRet.scaleAdd2(vTemp.dot(tempNorm) / l_dot_n, v, pt1);\r
 514     return ptRet;\r
 515   }\r
 516 \r
 517   /*\r
 518     public static Point3f getTriangleIntersection(Point3f a1, Point3f a2,\r
 519                                                  Point3f a3, Point4f plane,\r
 520                                                  Point3f b1,\r
 521                                                  Point3f b2, Point3f b3,\r
 522                                                  Vector3f vNorm, Vector3f vTemp, \r
 523                                                  Point3f ptRet, Point3f ptTemp, Vector3f vTemp2, Point4f pTemp, Vector3f vTemp3) {\r
 524       \r
 525       if (getTriangleIntersection(b1, b2, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))\r
 526         return ptRet;\r
 527       if (getTriangleIntersection(b2, b3, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))\r
 528         return ptRet;\r
 529       if (getTriangleIntersection(b3, b1, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))\r
 530         return ptRet;\r
 531       return null;\r
 532     }\r
 533   */\r
 534   /*  \r
 535     public static boolean getTriangleIntersection(Point3f b1, Point3f b2,\r
 536                                                   Point3f a1, Point3f a2,\r
 537                                                   Point3f a3, Vector3f vTemp,\r
 538                                                   Point4f plane, Vector3f vNorm,\r
 539                                                   Vector3f vTemp2, Vector3f vTemp3,\r
 540                                                   Point3f ptRet,\r
 541                                                   Point3f ptTemp) {\r
 542       if (distanceToPlane(plane, b1) * distanceToPlane(plane, b2) >= 0)\r
 543         return false;\r
 544       vTemp.sub(b2, b1);\r
 545       vTemp.normalize();\r
 546       if (getIntersection(b1, vTemp, plane, ptRet, vNorm, vTemp2) != null) {\r
 547         if (isInTriangle(ptRet, a1, a2, a3, vTemp, vTemp2, vTemp3))\r
 548           return true;\r
 549       }\r
 550       return false;\r
 551     }\r
 552     private static boolean isInTriangle(Point3f p, Point3f a, Point3f b,\r
 553                                         Point3f c, Vector3f v0, Vector3f v1,\r
 554                                         Vector3f v2) {\r
 555       // from http://www.blackpawn.com/texts/pointinpoly/default.html\r
 556       // Compute barycentric coordinates\r
 557       v0.sub(c, a);\r
 558       v1.sub(b, a);\r
 559       v2.sub(p, a);\r
 560       float dot00 = v0.dot(v0);\r
 561       float dot01 = v0.dot(v1);\r
 562       float dot02 = v0.dot(v2);\r
 563       float dot11 = v1.dot(v1);\r
 564       float dot12 = v1.dot(v2);\r
 565       float invDenom = 1 / (dot00 * dot11 - dot01 * dot01);\r
 566       float u = (dot11 * dot02 - dot01 * dot12) * invDenom;\r
 567       float v = (dot00 * dot12 - dot01 * dot02) * invDenom;\r
 568       return (u > 0 && v > 0 && u + v < 1);\r
 569     }\r
 570   */\r
 571 \r
 572   /**\r
 573    * Closed-form solution of absolute orientation requiring 1:1 mapping of\r
 574    * positions.\r
 575    * \r
 576    * @param centerAndPoints\r
 577    * @param retStddev\r
 578    * @return unit quaternion representation rotation\r
 579    * \r
 580    * @author hansonr Bob Hanson\r
 581    * \r
 582    */\r
 583   public static Quat calculateQuaternionRotation(P3[][] centerAndPoints,\r
 584                                                  float[] retStddev) {\r
 585 \r
 586     retStddev[1] = Float.NaN;\r
 587     Quat q = new Quat();\r
 588     if (centerAndPoints[0].length == 1\r
 589         || centerAndPoints[0].length != centerAndPoints[1].length)\r
 590       return q;\r
 591 \r
 592     /*\r
 593      * see Berthold K. P. Horn,\r
 594      * "Closed-form solution of absolute orientation using unit quaternions" J.\r
 595      * Opt. Soc. Amer. A, 1987, Vol. 4, pp. 629-642\r
 596      * http://www.opticsinfobase.org/viewmedia.cfm?uri=josaa-4-4-629&seq=0\r
 597      * \r
 598      * \r
 599      * A similar treatment was developed independently (and later!) \r
 600      * by G. Kramer, in G. R. Kramer,\r
 601      * "Superposition of Molecular Structures Using Quaternions"\r
 602      * Molecular Simulation, 1991, Vol. 7, pp. 113-119. \r
 603      * \r
 604      *  In that treatment there is a lot of unnecessary calculation \r
 605      *  along the trace of matrix M (eqn 20). \r
 606      *  I'm not sure why the extra x^2 + y^2 + z^2 + x'^2 + y'^2 + z'^2\r
 607      *  is in there, but they are unnecessary and only contribute to larger\r
 608      *  numerical averaging errors and additional processing time, as far as\r
 609      *  I can tell. Adding aI, where a is a scalar and I is the 4x4 identity\r
 610      *  just offsets the eigenvalues but doesn't change the eigenvectors.\r
 611      * \r
 612      * and Lydia E. Kavraki, "Molecular Distance Measures"\r
 613      * http://cnx.org/content/m11608/latest/\r
 614      * \r
 615      */\r
 616 \r
 617     int n = centerAndPoints[0].length - 1;\r
 618     if (n < 2)\r
 619       return q;\r
 620 \r
 621     double Sxx = 0, Sxy = 0, Sxz = 0, Syx = 0, Syy = 0, Syz = 0, Szx = 0, Szy = 0, Szz = 0;\r
 622     P3 ptA = new P3();\r
 623     P3 ptB = new P3();\r
 624     for (int i = n + 1; --i >= 1;) {\r
 625       P3 aij = centerAndPoints[0][i];\r
 626       P3 bij = centerAndPoints[1][i];\r
 627       ptA.sub2(aij, centerAndPoints[0][0]);\r
 628       ptB.sub2(bij, centerAndPoints[0][1]);\r
 629       Sxx += (double) ptA.x * (double) ptB.x;\r
 630       Sxy += (double) ptA.x * (double) ptB.y;\r
 631       Sxz += (double) ptA.x * (double) ptB.z;\r
 632       Syx += (double) ptA.y * (double) ptB.x;\r
 633       Syy += (double) ptA.y * (double) ptB.y;\r
 634       Syz += (double) ptA.y * (double) ptB.z;\r
 635       Szx += (double) ptA.z * (double) ptB.x;\r
 636       Szy += (double) ptA.z * (double) ptB.y;\r
 637       Szz += (double) ptA.z * (double) ptB.z;\r
 638     }\r
 639     retStddev[0] = getRmsd(centerAndPoints, q);\r
 640     double[][] N = new double[4][4];\r
 641     N[0][0] = Sxx + Syy + Szz;\r
 642     N[0][1] = N[1][0] = Syz - Szy;\r
 643     N[0][2] = N[2][0] = Szx - Sxz;\r
 644     N[0][3] = N[3][0] = Sxy - Syx;\r
 645 \r
 646     N[1][1] = Sxx - Syy - Szz;\r
 647     N[1][2] = N[2][1] = Sxy + Syx;\r
 648     N[1][3] = N[3][1] = Szx + Sxz;\r
 649 \r
 650     N[2][2] = -Sxx + Syy - Szz;\r
 651     N[2][3] = N[3][2] = Syz + Szy;\r
 652 \r
 653     N[3][3] = -Sxx - Syy + Szz;\r
 654 \r
 655     //this construction prevents JavaScript from requiring preloading of Eigen\r
 656     \r
 657     float[] v = ((EigenInterface) Interface.getInterface("javajs.util.Eigen"))\r
 658         .setM(N).getEigenvectorsFloatTransposed()[3];\r
 659     q = Quat.newP4(P4.new4(v[1], v[2], v[3], v[0]));\r
 660     retStddev[1] = getRmsd(centerAndPoints, q);\r
 661     return q;\r
 662   }\r
 663 \r
 664   /**\r
 665    * Fills a 4x4 matrix with rotation-translation of mapped points A to B.\r
 666    * If centerA is null, this is a standard 4x4 rotation-translation matrix;\r
 667    * otherwise, this 4x4 matrix is a rotation around a vector through the center of ptsA,\r
 668    * and centerA is filled with that center; \r
 669    * Prior to Jmol 14.3.12_2014.02.14, when used from the JmolScript compare() function,\r
 670    * this method returned the second of these options instead of the first.\r
 671    * \r
 672    * @param ptsA\r
 673    * @param ptsB\r
 674    * @param m  4x4 matrix to be returned \r
 675    * @param centerA return center of rotation; if null, then standard 4x4 matrix is returned\r
 676    * @return stdDev\r
 677    */\r
 678   public static float getTransformMatrix4(Lst<P3> ptsA, Lst<P3> ptsB, M4 m,\r
 679                                           P3 centerA) {\r
 680     P3[] cptsA = getCenterAndPoints(ptsA);\r
 681     P3[] cptsB = getCenterAndPoints(ptsB);\r
 682     float[] retStddev = new float[2];\r
 683     Quat q = calculateQuaternionRotation(new P3[][] { cptsA, cptsB },\r
 684         retStddev);\r
 685     M3 r = q.getMatrix();\r
 686     if (centerA == null)\r
 687       r.rotate(cptsA[0]);\r
 688     else\r
 689       centerA.setT(cptsA[0]);\r
 690     V3 t = V3.newVsub(cptsB[0], cptsA[0]);\r
 691     m.setMV(r, t);\r
 692     return retStddev[1];\r
 693   }\r
 694 \r
 695   /**\r
 696    * from a list of points, create an array that includes the center\r
 697    * point as the first point. This array is used as a starting point for\r
 698    * a quaternion analysis of superposition.\r
 699    * \r
 700    * @param vPts\r
 701    * @return  array of points with first point center\r
 702    */\r
 703         public static P3[] getCenterAndPoints(Lst<P3> vPts) {\r
 704           int n = vPts.size();\r
 705           P3[] pts = new P3[n + 1];\r
 706           pts[0] = new P3();\r
 707           if (n > 0) {\r
 708             for (int i = 0; i < n; i++) {\r
 709               pts[0].add(pts[i + 1] = vPts.get(i));\r
 710             }\r
 711             pts[0].scale(1f / n);\r
 712           }\r
 713           return pts;\r
 714         }\r
 715 \r
 716   public static float getRmsd(P3[][] centerAndPoints, Quat q) {\r
 717     double sum2 = 0;\r
 718     P3[] ptsA = centerAndPoints[0];\r
 719     P3[] ptsB = centerAndPoints[1];\r
 720     P3 cA = ptsA[0];\r
 721     P3 cB = ptsB[0];\r
 722     int n = ptsA.length - 1;\r
 723     P3 ptAnew = new P3();\r
 724     \r
 725     for (int i = n + 1; --i >= 1;) {\r
 726       ptAnew.sub2(ptsA[i], cA);\r
 727       q.transform2(ptAnew, ptAnew).add(cB);\r
 728       sum2 += ptAnew.distanceSquared(ptsB[i]);\r
 729     }\r
 730     return (float) Math.sqrt(sum2 / n);\r
 731   }\r
 732 \r
 733 }\r