Merge branch 'features/r2_11_2/JAL-3821_reinstate_patch' into develop
[jalview.git] / 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 }