3 * $Date: 2005-11-10 09:52:44 -0600 (Thu, 10 Nov 2005) $
6 * Copyright (C) 2003-2005 The Jmol Development Team
8 * Contact: jmol-developers@lists.sf.net
10 * This library is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public
12 * License as published by the Free Software Foundation; either
13 * version 2.1 of the License, or (at your option) any later version.
15 * This library is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with this library; if not, write to the Free Software
22 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
26 import javajs.api.EigenInterface;
28 import javajs.api.Interface;
33 //import org.jmol.script.T;
35 final public class Measure {
37 public final static float radiansPerDegree = (float) (2 * Math.PI / 360);
39 public static float computeAngle(T3 pointA, T3 pointB, T3 pointC, V3 vectorBA, V3 vectorBC, boolean asDegrees) {
40 vectorBA.sub2(pointA, pointB);
41 vectorBC.sub2(pointC, pointB);
42 float angle = vectorBA.angle(vectorBC);
43 return (asDegrees ? angle / radiansPerDegree : angle);
46 public static float computeAngleABC(T3 pointA, T3 pointB, T3 pointC, boolean asDegrees) {
47 V3 vectorBA = new V3();
48 V3 vectorBC = new V3();
49 return computeAngle(pointA, pointB, pointC, vectorBA, vectorBC, asDegrees);
52 public static float computeTorsion(T3 p1, T3 p2, T3 p3, T3 p4, boolean asDegrees) {
54 float ijx = p1.x - p2.x;
55 float ijy = p1.y - p2.y;
56 float ijz = p1.z - p2.z;
58 float kjx = p3.x - p2.x;
59 float kjy = p3.y - p2.y;
60 float kjz = p3.z - p2.z;
62 float klx = p3.x - p4.x;
63 float kly = p3.y - p4.y;
64 float klz = p3.z - p4.z;
66 float ax = ijy * kjz - ijz * kjy;
67 float ay = ijz * kjx - ijx * kjz;
68 float az = ijx * kjy - ijy * kjx;
69 float cx = kjy * klz - kjz * kly;
70 float cy = kjz * klx - kjx * klz;
71 float cz = kjx * kly - kjy * klx;
73 float ai2 = 1f / (ax * ax + ay * ay + az * az);
74 float ci2 = 1f / (cx * cx + cy * cy + cz * cz);
76 float ai = (float) Math.sqrt(ai2);
77 float ci = (float) Math.sqrt(ci2);
78 float denom = ai * ci;
79 float cross = ax * cx + ay * cy + az * cz;
80 float cosang = cross * denom;
88 float torsion = (float) Math.acos(cosang);
89 float dot = ijx * cx + ijy * cy + ijz * cz;
90 float absDot = Math.abs(dot);
91 torsion = (dot / absDot > 0) ? torsion : -torsion;
92 return (asDegrees ? torsion / radiansPerDegree : torsion);
96 * This method calculates measures relating to two points in space
97 * with related quaternion frame difference. It is used in Jmol for
98 * calculating straightness and many other helical quantities.
103 * @return new T3[] { pt_a_prime, n, r, P3.new3(theta, pitch, residuesPerTurn), pt_b_prime };
105 public static T3[] computeHelicalAxis(P3 a, P3 b, Quat dq) {
113 // | / \ Vcb = Vab . n
114 // n | / \d Vda = (Vcb - Vab) / 2
122 * testing here to see if directing the normal makes any difference -- oddly
123 * enough, it does not. When n = -n and theta = -theta vab.n is reversed,
124 * and that magnitude is multiplied by n in generating the A'-B' vector.
126 * a negative angle implies a left-handed axis (sheets)
128 float theta = dq.getTheta();
129 V3 n = dq.getNormal();
130 float v_dot_n = vab.dot(n);
131 if (Math.abs(v_dot_n) < 0.0001f)
133 V3 va_prime_d = new V3();
134 va_prime_d.cross(vab, n);
135 if (va_prime_d.dot(va_prime_d) != 0)
136 va_prime_d.normalize();
140 v_dot_n = PT.FLOAT_MIN_SAFE; // allow for perpendicular axis to vab
144 va_prime_d.scale(theta == 0 ? 0 : (float) (vda.length() / Math.tan(theta
145 / 2 / 180 * Math.PI)));
146 V3 r = V3.newV(va_prime_d);
149 P3 pt_a_prime = P3.newP(a);
151 // already done this. ??
152 if (v_dot_n != PT.FLOAT_MIN_SAFE)
154 // must calculate directed angle:
155 P3 pt_b_prime = P3.newP(pt_a_prime);
157 theta = computeTorsion(a, pt_a_prime, pt_b_prime, b, true);
158 if (Float.isNaN(theta) || r.length() < 0.0001f)
159 theta = dq.getThetaDirectedV(n); // allow for r = 0
160 // anything else is an array
161 float residuesPerTurn = Math.abs(theta == 0 ? 0 : 360f / theta);
162 float pitch = Math.abs(v_dot_n == PT.FLOAT_MIN_SAFE ? 0 : n.length()
163 * (theta == 0 ? 1 : 360f / theta));
164 return new T3[] { pt_a_prime, n, r, P3.new3(theta, pitch, residuesPerTurn), pt_b_prime };
167 public static P4 getPlaneThroughPoints(T3 pointA,
171 float w = getNormalThroughPoints(pointA, pointB, pointC, vNorm, vAB);
172 plane.set4(vNorm.x, vNorm.y, vNorm.z, w);
176 public static void getPlaneThroughPoint(T3 pt, V3 normal, P4 plane) {
177 plane.set4(normal.x, normal.y, normal.z, -normal.dot(pt));
180 public static float distanceToPlane(P4 plane, T3 pt) {
181 return (plane == null ? Float.NaN
182 : (plane.dot(pt) + plane.w) / (float) Math.sqrt(plane.dot(plane)));
185 public static float directedDistanceToPlane(P3 pt, P4 plane, P3 ptref) {
186 float f = plane.dot(pt) + plane.w;
187 float f1 = plane.dot(ptref) + plane.w;
188 return Math.signum(f1) * f / (float) Math.sqrt(plane.dot(plane));
191 public static float distanceToPlaneD(P4 plane, float d, P3 pt) {
192 return (plane == null ? Float.NaN : (plane.dot(pt) + plane.w) / d);
195 public static float distanceToPlaneV(V3 norm, float w, P3 pt) {
196 return (norm == null ? Float.NaN
197 : (norm.dot(pt) + w) / (float) Math.sqrt(norm.dot(norm)));
201 * note that if vAB or vAC is dispensible, vNormNorm can be one of them
208 public static void calcNormalizedNormal(T3 pointA, T3 pointB,
209 T3 pointC, V3 vNormNorm, V3 vAB) {
210 vAB.sub2(pointB, pointA);
211 vNormNorm.sub2(pointC, pointA);
212 vNormNorm.cross(vAB, vNormNorm);
213 vNormNorm.normalize();
216 public static float getDirectedNormalThroughPoints(T3 pointA,
217 T3 pointB, T3 pointC, T3 ptRef, V3 vNorm,
219 // for x = plane({atomno=1}, {atomno=2}, {atomno=3}, {atomno=4})
220 float nd = getNormalThroughPoints(pointA, pointB, pointC, vNorm, vAB);
222 P3 pt0 = P3.newP(pointA);
224 float d = pt0.distance(ptRef);
225 pt0.sub2(pointA, vNorm);
226 if (d > pt0.distance(ptRef)) {
235 * if vAC is dispensible vNorm can be vAC
243 public static float getNormalThroughPoints(T3 pointA, T3 pointB,
244 T3 pointC, V3 vNorm, V3 vTemp) {
246 calcNormalizedNormal(pointA, pointB, pointC, vNorm, vTemp);
247 // ax + by + cz + d = 0
248 // so if a point is in the plane, then N dot X = -d
250 return -vTemp.dot(vNorm);
253 public static void getPlaneProjection(P3 pt, P4 plane, P3 ptProj, V3 vNorm) {
254 float dist = distanceToPlane(plane, pt);
255 vNorm.set(plane.x, plane.y, plane.z);
258 ptProj.add2(pt, vNorm);
261 public final static V3 axisY = V3.new3(0, 1, 0);
263 public static void getNormalToLine(P3 pointA, P3 pointB,
265 // vector in xy plane perpendicular to a line between two points RMH
266 vNormNorm.sub2(pointA, pointB);
267 vNormNorm.cross(vNormNorm, axisY);
268 vNormNorm.normalize();
269 if (Float.isNaN(vNormNorm.x))
270 vNormNorm.set(1, 0, 0);
273 public static void getBisectingPlane(P3 pointA, V3 vAB,
274 T3 ptTemp, V3 vTemp, P4 plane) {
275 ptTemp.scaleAdd2(0.5f, vAB, pointA);
278 getPlaneThroughPoint(ptTemp, vTemp, plane);
281 public static void projectOntoAxis(P3 point, P3 axisA,
283 V3 vectorProjection) {
284 vectorProjection.sub2(point, axisA);
285 float projectedLength = vectorProjection.dot(axisUnitVector);
286 point.scaleAdd2(projectedLength, axisUnitVector, axisA);
287 vectorProjection.sub2(point, axisA);
290 public static void calcBestAxisThroughPoints(P3[] points, P3 axisA,
294 // just a crude starting point.
296 int nPoints = points.length;
297 axisA.setT(points[0]);
298 axisUnitVector.sub2(points[nPoints - 1], axisA);
299 axisUnitVector.normalize();
302 * We now calculate the least-squares 3D axis
303 * through the helix alpha carbons starting with Vo
304 * as a first approximation.
306 * This uses the simple 0-centered least squares fit:
310 * minimizing R^2 = SUM(|Y - Yi|^2)
312 * where Yi is the vector PERPENDICULAR of the point onto axis Vo
313 * and Xi is the vector PROJECTION of the point onto axis Vo
314 * and M is a vector adjustment
316 * M = SUM_(Xi cross Yi) / sum(|Xi|^2)
318 * from which we arrive at:
320 * V = Vo + (M cross Vo)
322 * Basically, this is just a 3D version of a
323 * standard 2D least squares fit to a line, where we would say:
327 * D = n (sum xi^2) - (sum xi)^2
329 * m = [(n sum xiyi) - (sum xi)(sum yi)] / D
330 * b = [(sum yi) (sum xi^2) - (sum xi)(sum xiyi)] / D
332 * but here we demand that the line go through the center, so we
333 * require (sum xi) = (sum yi) = 0, so b = 0 and
335 * m = (sum xiyi) / (sum xi^2)
337 * In 3D we do the same but
338 * instead of x we have Vo,
339 * instead of multiplication we use cross products
341 * A bit of iteration is necessary.
347 calcAveragePointN(points, nPoints, axisA);
350 while (nTries++ < nTriesMax
351 && findAxis(points, nPoints, axisA, axisUnitVector, vectorProjection) > 0.001) {
355 * Iteration here gets the job done.
356 * We now find the projections of the endpoints onto the axis
360 P3 tempA = P3.newP(points[0]);
361 projectOntoAxis(tempA, axisA, axisUnitVector, vectorProjection);
365 public static float findAxis(P3[] points, int nPoints, P3 axisA,
366 V3 axisUnitVector, V3 vectorProjection) {
367 V3 sumXiYi = new V3();
370 P3 ptProj = new P3();
371 V3 a = V3.newV(axisUnitVector);
374 for (int i = nPoints; --i >= 0;) {
377 projectOntoAxis(ptProj, axisA, axisUnitVector,
379 vTemp.sub2(pt, ptProj);
380 //sum_Yi2 += vTemp.lengthSquared();
381 vTemp.cross(vectorProjection, vTemp);
383 sum_Xi2 += vectorProjection.lengthSquared();
385 V3 m = V3.newV(sumXiYi);
386 m.scale(1 / sum_Xi2);
387 vTemp.cross(m, axisUnitVector);
388 axisUnitVector.add(vTemp);
389 axisUnitVector.normalize();
390 //check for change in direction by measuring vector difference length
391 vTemp.sub2(axisUnitVector, a);
392 return vTemp.length();
396 public static void calcAveragePoint(P3 pointA, P3 pointB,
398 pointC.set((pointA.x + pointB.x) / 2, (pointA.y + pointB.y) / 2,
399 (pointA.z + pointB.z) / 2);
402 public static void calcAveragePointN(P3[] points, int nPoints,
404 averagePoint.setT(points[0]);
405 for (int i = 1; i < nPoints; i++)
406 averagePoint.add(points[i]);
407 averagePoint.scale(1f / nPoints);
410 public static Lst<P3> transformPoints(Lst<P3> vPts, M4 m4, P3 center) {
411 Lst<P3> v = new Lst<P3>();
412 for (int i = 0; i < vPts.size(); i++) {
413 P3 pt = P3.newP(vPts.get(i));
422 public static boolean isInTetrahedron(P3 pt, P3 ptA, P3 ptB,
425 V3 vTemp2, boolean fullyEnclosed) {
426 boolean b = (distanceToPlane(getPlaneThroughPoints(ptC, ptD, ptA, vTemp, vTemp2, plane), pt) >= 0);
427 if (b != (distanceToPlane(getPlaneThroughPoints(ptA, ptD, ptB, vTemp, vTemp2, plane), pt) >= 0))
429 if (b != (distanceToPlane(getPlaneThroughPoints(ptB, ptD, ptC, vTemp, vTemp2, plane), pt) >= 0))
431 float d = distanceToPlane(getPlaneThroughPoints(ptA, ptB, ptC, vTemp, vTemp2, plane), pt);
433 return (b == (d >= 0));
434 float d1 = distanceToPlane(plane, ptD);
435 return d1 * d <= 0 || Math.abs(d1) > Math.abs(d);
443 * @return [ point, vector ] or []
445 public static Lst<Object> getIntersectionPP(P4 plane1, P4 plane2) {
454 V3 norm1 = V3.new3(a1, b1, c1);
455 V3 norm2 = V3.new3(a2, b2, c2);
457 nxn.cross(norm1, norm2);
458 float ax = Math.abs(nxn.x);
459 float ay = Math.abs(nxn.y);
460 float az = Math.abs(nxn.z);
462 int type = (ax > ay ? (ax > az ? 1 : 3) : ay > az ? 2 : 3);
466 diff = (b1 * c2 - b2 * c1);
467 if (Math.abs(diff) < 0.01) return null;
468 y = (c1 * d2 - c2 * d1) / diff;
469 z = (b2 * d1 - d2 * b1) / diff;
472 diff = (a1 * c2 - a2 * c1);
473 if (Math.abs(diff) < 0.01) return null;
474 x = (c1 * d2 - c2 * d1) / diff;
476 z = (a2 * d1 - d2 * a1) / diff;
480 diff = (a1 * b2 - a2 * b1);
481 if (Math.abs(diff) < 0.01) return null;
482 x = (b1 * d2 - b2 * d1) / diff;
483 y = (a2 * d1 - d2 * a1) / diff;
486 Lst<Object>list = new Lst<Object>();
487 list.addLast(P3.new3(x, y, z));
495 * @param pt1 point on line
496 * @param v unit vector of line
498 * @param ptRet point of intersection of line with plane
503 public static P3 getIntersection(P3 pt1, V3 v,
504 P4 plane, P3 ptRet, V3 tempNorm, V3 vTemp) {
505 getPlaneProjection(pt1, plane, ptRet, tempNorm);
506 tempNorm.set(plane.x, plane.y, plane.z);
507 tempNorm.normalize();
509 v = V3.newV(tempNorm);
510 float l_dot_n = v.dot(tempNorm);
511 if (Math.abs(l_dot_n) < 0.01) return null;
512 vTemp.sub2(ptRet, pt1);
513 ptRet.scaleAdd2(vTemp.dot(tempNorm) / l_dot_n, v, pt1);
518 public static Point3f getTriangleIntersection(Point3f a1, Point3f a2,
519 Point3f a3, Point4f plane,
521 Point3f b2, Point3f b3,
522 Vector3f vNorm, Vector3f vTemp,
523 Point3f ptRet, Point3f ptTemp, Vector3f vTemp2, Point4f pTemp, Vector3f vTemp3) {
525 if (getTriangleIntersection(b1, b2, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))
527 if (getTriangleIntersection(b2, b3, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))
529 if (getTriangleIntersection(b3, b1, a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp))
535 public static boolean getTriangleIntersection(Point3f b1, Point3f b2,
536 Point3f a1, Point3f a2,
537 Point3f a3, Vector3f vTemp,
538 Point4f plane, Vector3f vNorm,
539 Vector3f vTemp2, Vector3f vTemp3,
542 if (distanceToPlane(plane, b1) * distanceToPlane(plane, b2) >= 0)
546 if (getIntersection(b1, vTemp, plane, ptRet, vNorm, vTemp2) != null) {
547 if (isInTriangle(ptRet, a1, a2, a3, vTemp, vTemp2, vTemp3))
552 private static boolean isInTriangle(Point3f p, Point3f a, Point3f b,
553 Point3f c, Vector3f v0, Vector3f v1,
555 // from http://www.blackpawn.com/texts/pointinpoly/default.html
556 // Compute barycentric coordinates
560 float dot00 = v0.dot(v0);
561 float dot01 = v0.dot(v1);
562 float dot02 = v0.dot(v2);
563 float dot11 = v1.dot(v1);
564 float dot12 = v1.dot(v2);
565 float invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
566 float u = (dot11 * dot02 - dot01 * dot12) * invDenom;
567 float v = (dot00 * dot12 - dot01 * dot02) * invDenom;
568 return (u > 0 && v > 0 && u + v < 1);
573 * Closed-form solution of absolute orientation requiring 1:1 mapping of
576 * @param centerAndPoints
578 * @return unit quaternion representation rotation
580 * @author hansonr Bob Hanson
583 public static Quat calculateQuaternionRotation(P3[][] centerAndPoints,
586 retStddev[1] = Float.NaN;
588 if (centerAndPoints[0].length == 1
589 || centerAndPoints[0].length != centerAndPoints[1].length)
593 * see Berthold K. P. Horn,
594 * "Closed-form solution of absolute orientation using unit quaternions" J.
595 * Opt. Soc. Amer. A, 1987, Vol. 4, pp. 629-642
596 * http://www.opticsinfobase.org/viewmedia.cfm?uri=josaa-4-4-629&seq=0
599 * A similar treatment was developed independently (and later!)
600 * by G. Kramer, in G. R. Kramer,
601 * "Superposition of Molecular Structures Using Quaternions"
602 * Molecular Simulation, 1991, Vol. 7, pp. 113-119.
604 * In that treatment there is a lot of unnecessary calculation
605 * along the trace of matrix M (eqn 20).
606 * I'm not sure why the extra x^2 + y^2 + z^2 + x'^2 + y'^2 + z'^2
607 * is in there, but they are unnecessary and only contribute to larger
608 * numerical averaging errors and additional processing time, as far as
609 * I can tell. Adding aI, where a is a scalar and I is the 4x4 identity
610 * just offsets the eigenvalues but doesn't change the eigenvectors.
612 * and Lydia E. Kavraki, "Molecular Distance Measures"
613 * http://cnx.org/content/m11608/latest/
617 int n = centerAndPoints[0].length - 1;
621 double Sxx = 0, Sxy = 0, Sxz = 0, Syx = 0, Syy = 0, Syz = 0, Szx = 0, Szy = 0, Szz = 0;
624 for (int i = n + 1; --i >= 1;) {
625 P3 aij = centerAndPoints[0][i];
626 P3 bij = centerAndPoints[1][i];
627 ptA.sub2(aij, centerAndPoints[0][0]);
628 ptB.sub2(bij, centerAndPoints[0][1]);
629 Sxx += (double) ptA.x * (double) ptB.x;
630 Sxy += (double) ptA.x * (double) ptB.y;
631 Sxz += (double) ptA.x * (double) ptB.z;
632 Syx += (double) ptA.y * (double) ptB.x;
633 Syy += (double) ptA.y * (double) ptB.y;
634 Syz += (double) ptA.y * (double) ptB.z;
635 Szx += (double) ptA.z * (double) ptB.x;
636 Szy += (double) ptA.z * (double) ptB.y;
637 Szz += (double) ptA.z * (double) ptB.z;
639 retStddev[0] = getRmsd(centerAndPoints, q);
640 double[][] N = new double[4][4];
641 N[0][0] = Sxx + Syy + Szz;
642 N[0][1] = N[1][0] = Syz - Szy;
643 N[0][2] = N[2][0] = Szx - Sxz;
644 N[0][3] = N[3][0] = Sxy - Syx;
646 N[1][1] = Sxx - Syy - Szz;
647 N[1][2] = N[2][1] = Sxy + Syx;
648 N[1][3] = N[3][1] = Szx + Sxz;
650 N[2][2] = -Sxx + Syy - Szz;
651 N[2][3] = N[3][2] = Syz + Szy;
653 N[3][3] = -Sxx - Syy + Szz;
655 //this construction prevents JavaScript from requiring preloading of Eigen
657 float[] v = ((EigenInterface) Interface.getInterface("javajs.util.Eigen"))
658 .setM(N).getEigenvectorsFloatTransposed()[3];
659 q = Quat.newP4(P4.new4(v[1], v[2], v[3], v[0]));
660 retStddev[1] = getRmsd(centerAndPoints, q);
665 * Fills a 4x4 matrix with rotation-translation of mapped points A to B.
666 * If centerA is null, this is a standard 4x4 rotation-translation matrix;
667 * otherwise, this 4x4 matrix is a rotation around a vector through the center of ptsA,
668 * and centerA is filled with that center;
669 * Prior to Jmol 14.3.12_2014.02.14, when used from the JmolScript compare() function,
670 * this method returned the second of these options instead of the first.
674 * @param m 4x4 matrix to be returned
675 * @param centerA return center of rotation; if null, then standard 4x4 matrix is returned
678 public static float getTransformMatrix4(Lst<P3> ptsA, Lst<P3> ptsB, M4 m,
680 P3[] cptsA = getCenterAndPoints(ptsA);
681 P3[] cptsB = getCenterAndPoints(ptsB);
682 float[] retStddev = new float[2];
683 Quat q = calculateQuaternionRotation(new P3[][] { cptsA, cptsB },
685 M3 r = q.getMatrix();
689 centerA.setT(cptsA[0]);
690 V3 t = V3.newVsub(cptsB[0], cptsA[0]);
696 * from a list of points, create an array that includes the center
697 * point as the first point. This array is used as a starting point for
698 * a quaternion analysis of superposition.
701 * @return array of points with first point center
703 public static P3[] getCenterAndPoints(Lst<P3> vPts) {
705 P3[] pts = new P3[n + 1];
708 for (int i = 0; i < n; i++) {
709 pts[0].add(pts[i + 1] = vPts.get(i));
711 pts[0].scale(1f / n);
716 public static float getRmsd(P3[][] centerAndPoints, Quat q) {
718 P3[] ptsA = centerAndPoints[0];
719 P3[] ptsB = centerAndPoints[1];
722 int n = ptsA.length - 1;
723 P3 ptAnew = new P3();
725 for (int i = n + 1; --i >= 1;) {
726 ptAnew.sub2(ptsA[i], cA);
727 q.transform2(ptAnew, ptAnew).add(cB);
728 sum2 += ptAnew.distanceSquared(ptsB[i]);
730 return (float) Math.sqrt(sum2 / n);