unused/javajs/util/Measure.java

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