3 * $Date: 2005-11-10 09:52:44 -0600 (Thu, 10 Nov 2005) $
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.
9 * Copyright (C) 2003-2005 The Jmol Development Team
11 * Contact: jmol-developers@lists.sf.net
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.
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.
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.
29 import javajs.api.EigenInterface;
31 import javajs.api.Interface;
36 //import org.jmol.script.T;
38 final public class Measure {
40 public final static float radiansPerDegree = (float) (2 * Math.PI / 360);
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);
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);
55 public static float computeTorsion(T3 p1, T3 p2, T3 p3, T3 p4, boolean asDegrees) {
57 float ijx = p1.x - p2.x;
58 float ijy = p1.y - p2.y;
59 float ijz = p1.z - p2.z;
61 float kjx = p3.x - p2.x;
62 float kjy = p3.y - p2.y;
63 float kjz = p3.z - p2.z;
65 float klx = p3.x - p4.x;
66 float kly = p3.y - p4.y;
67 float klz = p3.z - p4.z;
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;
76 float ai2 = 1f / (ax * ax + ay * ay + az * az);
77 float ci2 = 1f / (cx * cx + cy * cy + cz * cz);
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;
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);
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.
106 * @return new T3[] { pt_a_prime, n, r, P3.new3(theta, pitch, residuesPerTurn), pt_b_prime };
108 public static T3[] computeHelicalAxis(P3 a, P3 b, Quat dq) {
116 // | / \ Vcb = Vab . n
117 // n | / \d Vda = (Vcb - Vab) / 2
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.
129 * a negative angle implies a left-handed axis (sheets)
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)
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();
143 v_dot_n = PT.FLOAT_MIN_SAFE; // allow for perpendicular axis to vab
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);
152 P3 pt_a_prime = P3.newP(a);
154 // already done this. ??
155 if (v_dot_n != PT.FLOAT_MIN_SAFE)
157 // must calculate directed angle:
158 P3 pt_b_prime = P3.newP(pt_a_prime);
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 };
170 public static P4 getPlaneThroughPoints(T3 pointA,
174 float w = getNormalThroughPoints(pointA, pointB, pointC, vNorm, vAB);
175 plane.set4(vNorm.x, vNorm.y, vNorm.z, w);
179 public static void getPlaneThroughPoint(T3 pt, V3 normal, P4 plane) {
180 plane.set4(normal.x, normal.y, normal.z, -normal.dot(pt));
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)));
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));
194 public static float distanceToPlaneD(P4 plane, float d, P3 pt) {
195 return (plane == null ? Float.NaN : (plane.dot(pt) + plane.w) / d);
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)));
204 * note that if vAB or vAC is dispensible, vNormNorm can be one of them
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();
219 public static float getDirectedNormalThroughPoints(T3 pointA,
220 T3 pointB, T3 pointC, T3 ptRef, V3 vNorm,
222 // for x = plane({atomno=1}, {atomno=2}, {atomno=3}, {atomno=4})
223 float nd = getNormalThroughPoints(pointA, pointB, pointC, vNorm, vAB);
225 P3 pt0 = P3.newP(pointA);
227 float d = pt0.distance(ptRef);
228 pt0.sub2(pointA, vNorm);
229 if (d > pt0.distance(ptRef)) {
238 * if vAC is dispensible vNorm can be vAC
246 public static float getNormalThroughPoints(T3 pointA, T3 pointB,
247 T3 pointC, T3 vNorm, T3 vTemp) {
249 calcNormalizedNormal(pointA, pointB, pointC, vNorm, vTemp);
250 // ax + by + cz + d = 0
251 // so if a point is in the plane, then N dot X = -d
253 return -vTemp.dot(vNorm);
256 public static void getPlaneProjection(P3 pt, P4 plane, P3 ptProj, V3 vNorm) {
257 float dist = distanceToPlane(plane, pt);
258 vNorm.set(plane.x, plane.y, plane.z);
261 ptProj.add2(pt, vNorm);
271 * @param normal set to be opposite to direction of ptCenter from ABC
273 * @return true if winding is CCW; false if CW
275 public static boolean getNormalFromCenter(P3 ptCenter, P3 ptA, P3 ptB, P3 ptC,
276 boolean isOutward, V3 normal, V3 vTemp) {
277 float d = getNormalThroughPoints(ptA, ptB, ptC, normal, vTemp);
278 boolean isReversed = (distanceToPlaneV(normal, d, ptCenter) > 0);
279 if (isReversed == isOutward)
284 public final static V3 axisY = V3.new3(0, 1, 0);
286 public static void getNormalToLine(P3 pointA, P3 pointB,
288 // vector in xy plane perpendicular to a line between two points RMH
289 vNormNorm.sub2(pointA, pointB);
290 vNormNorm.cross(vNormNorm, axisY);
291 vNormNorm.normalize();
292 if (Float.isNaN(vNormNorm.x))
293 vNormNorm.set(1, 0, 0);
296 public static void getBisectingPlane(P3 pointA, V3 vAB,
297 T3 ptTemp, V3 vTemp, P4 plane) {
298 ptTemp.scaleAdd2(0.5f, vAB, pointA);
301 getPlaneThroughPoint(ptTemp, vTemp, plane);
304 public static void projectOntoAxis(P3 point, P3 axisA,
306 V3 vectorProjection) {
307 vectorProjection.sub2(point, axisA);
308 float projectedLength = vectorProjection.dot(axisUnitVector);
309 point.scaleAdd2(projectedLength, axisUnitVector, axisA);
310 vectorProjection.sub2(point, axisA);
313 public static void calcBestAxisThroughPoints(P3[] points, P3 axisA,
317 // just a crude starting point.
319 int nPoints = points.length;
320 axisA.setT(points[0]);
321 axisUnitVector.sub2(points[nPoints - 1], axisA);
322 axisUnitVector.normalize();
325 * We now calculate the least-squares 3D axis
326 * through the helix alpha carbons starting with Vo
327 * as a first approximation.
329 * This uses the simple 0-centered least squares fit:
333 * minimizing R^2 = SUM(|Y - Yi|^2)
335 * where Yi is the vector PERPENDICULAR of the point onto axis Vo
336 * and Xi is the vector PROJECTION of the point onto axis Vo
337 * and M is a vector adjustment
339 * M = SUM_(Xi cross Yi) / sum(|Xi|^2)
341 * from which we arrive at:
343 * V = Vo + (M cross Vo)
345 * Basically, this is just a 3D version of a
346 * standard 2D least squares fit to a line, where we would say:
350 * D = n (sum xi^2) - (sum xi)^2
352 * m = [(n sum xiyi) - (sum xi)(sum yi)] / D
353 * b = [(sum yi) (sum xi^2) - (sum xi)(sum xiyi)] / D
355 * but here we demand that the line go through the center, so we
356 * require (sum xi) = (sum yi) = 0, so b = 0 and
358 * m = (sum xiyi) / (sum xi^2)
360 * In 3D we do the same but
361 * instead of x we have Vo,
362 * instead of multiplication we use cross products
364 * A bit of iteration is necessary.
370 calcAveragePointN(points, nPoints, axisA);
373 while (nTries++ < nTriesMax
374 && findAxis(points, nPoints, axisA, axisUnitVector, vectorProjection) > 0.001) {
378 * Iteration here gets the job done.
379 * We now find the projections of the endpoints onto the axis
383 P3 tempA = P3.newP(points[0]);
384 projectOntoAxis(tempA, axisA, axisUnitVector, vectorProjection);
388 public static float findAxis(P3[] points, int nPoints, P3 axisA,
389 V3 axisUnitVector, V3 vectorProjection) {
390 V3 sumXiYi = new V3();
393 P3 ptProj = new P3();
394 V3 a = V3.newV(axisUnitVector);
397 for (int i = nPoints; --i >= 0;) {
400 projectOntoAxis(ptProj, axisA, axisUnitVector,
402 vTemp.sub2(pt, ptProj);
403 //sum_Yi2 += vTemp.lengthSquared();
404 vTemp.cross(vectorProjection, vTemp);
406 sum_Xi2 += vectorProjection.lengthSquared();
408 V3 m = V3.newV(sumXiYi);
409 m.scale(1 / sum_Xi2);
410 vTemp.cross(m, axisUnitVector);
411 axisUnitVector.add(vTemp);
412 axisUnitVector.normalize();
413 //check for change in direction by measuring vector difference length
414 vTemp.sub2(axisUnitVector, a);
415 return vTemp.length();
419 public static void calcAveragePoint(P3 pointA, P3 pointB,
421 pointC.set((pointA.x + pointB.x) / 2, (pointA.y + pointB.y) / 2,
422 (pointA.z + pointB.z) / 2);
425 public static void calcAveragePointN(P3[] points, int nPoints,
427 averagePoint.setT(points[0]);
428 for (int i = 1; i < nPoints; i++)
429 averagePoint.add(points[i]);
430 averagePoint.scale(1f / nPoints);
433 public static Lst<P3> transformPoints(Lst<P3> vPts, M4 m4, P3 center) {
434 Lst<P3> v = new Lst<P3>();
435 for (int i = 0; i < vPts.size(); i++) {
436 P3 pt = P3.newP(vPts.get(i));
445 public static boolean isInTetrahedron(P3 pt, P3 ptA, P3 ptB,
448 V3 vTemp2, boolean fullyEnclosed) {
449 boolean b = (distanceToPlane(getPlaneThroughPoints(ptC, ptD, ptA, vTemp, vTemp2, plane), pt) >= 0);
450 if (b != (distanceToPlane(getPlaneThroughPoints(ptA, ptD, ptB, vTemp, vTemp2, plane), pt) >= 0))
452 if (b != (distanceToPlane(getPlaneThroughPoints(ptB, ptD, ptC, vTemp, vTemp2, plane), pt) >= 0))
454 float d = distanceToPlane(getPlaneThroughPoints(ptA, ptB, ptC, vTemp, vTemp2, plane), pt);
456 return (b == (d >= 0));
457 float d1 = distanceToPlane(plane, ptD);
458 return d1 * d <= 0 || Math.abs(d1) > Math.abs(d);
466 * @return [ point, vector ] or []
468 public static Lst<Object> getIntersectionPP(P4 plane1, P4 plane2) {
477 V3 norm1 = V3.new3(a1, b1, c1);
478 V3 norm2 = V3.new3(a2, b2, c2);
480 nxn.cross(norm1, norm2);
481 float ax = Math.abs(nxn.x);
482 float ay = Math.abs(nxn.y);
483 float az = Math.abs(nxn.z);
485 int type = (ax > ay ? (ax > az ? 1 : 3) : ay > az ? 2 : 3);
489 diff = (b1 * c2 - b2 * c1);
490 if (Math.abs(diff) < 0.01) return null;
491 y = (c1 * d2 - c2 * d1) / diff;
492 z = (b2 * d1 - d2 * b1) / diff;
495 diff = (a1 * c2 - a2 * c1);
496 if (Math.abs(diff) < 0.01) return null;
497 x = (c1 * d2 - c2 * d1) / diff;
499 z = (a2 * d1 - d2 * a1) / diff;
503 diff = (a1 * b2 - a2 * b1);
504 if (Math.abs(diff) < 0.01) return null;
505 x = (b1 * d2 - b2 * d1) / diff;
506 y = (a2 * d1 - d2 * a1) / diff;
509 Lst<Object>list = new Lst<Object>();
510 list.addLast(P3.new3(x, y, z));
518 * @param pt1 point on line
519 * @param v unit vector of line
521 * @param ptRet point of intersection of line with plane
526 public static P3 getIntersection(P3 pt1, V3 v,
527 P4 plane, P3 ptRet, V3 tempNorm, V3 vTemp) {
528 getPlaneProjection(pt1, plane, ptRet, tempNorm);
529 tempNorm.set(plane.x, plane.y, plane.z);
530 tempNorm.normalize();
532 v = V3.newV(tempNorm);
533 float l_dot_n = v.dot(tempNorm);
534 if (Math.abs(l_dot_n) < 0.01) return null;
535 vTemp.sub2(ptRet, pt1);
536 ptRet.scaleAdd2(vTemp.dot(tempNorm) / l_dot_n, v, pt1);
541 * public static Point3f getTriangleIntersection(Point3f a1, Point3f a2,
542 * Point3f a3, Point4f plane, Point3f b1, Point3f b2, Point3f b3, Vector3f
543 * vNorm, Vector3f vTemp, Point3f ptRet, Point3f ptTemp, Vector3f vTemp2,
544 * Point4f pTemp, Vector3f vTemp3) {
546 * if (getTriangleIntersection(b1, b2, a1, a2, a3, vTemp, plane, vNorm,
547 * vTemp2, vTemp3, ptRet, ptTemp)) return ptRet; if
548 * (getTriangleIntersection(b2, b3, a1, a2, a3, vTemp, plane, vNorm, vTemp2,
549 * vTemp3, ptRet, ptTemp)) return ptRet; if (getTriangleIntersection(b3, b1,
550 * a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp)) return
551 * ptRet; return null; }
554 * public static boolean getTriangleIntersection(Point3f b1, Point3f b2,
555 * Point3f a1, Point3f a2, Point3f a3, Vector3f vTemp, Point4f plane, Vector3f
556 * vNorm, Vector3f vTemp2, Vector3f vTemp3, Point3f ptRet, Point3f ptTemp) {
557 * if (distanceToPlane(plane, b1) * distanceToPlane(plane, b2) >= 0) return
558 * false; vTemp.sub(b2, b1); vTemp.normalize(); if (getIntersection(b1, vTemp,
559 * plane, ptRet, vNorm, vTemp2) != null) { if (isInTriangle(ptRet, a1, a2, a3,
560 * vTemp, vTemp2, vTemp3)) return true; } return false; } private static
561 * boolean isInTriangle(Point3f p, Point3f a, Point3f b, Point3f c, Vector3f
562 * v0, Vector3f v1, Vector3f v2) { // from
563 * http://www.blackpawn.com/texts/pointinpoly/default.html // Compute
564 * barycentric coordinates v0.sub(c, a); v1.sub(b, a); v2.sub(p, a); float
565 * dot00 = v0.dot(v0); float dot01 = v0.dot(v1); float dot02 = v0.dot(v2);
566 * float dot11 = v1.dot(v1); float dot12 = v1.dot(v2); float invDenom = 1 /
567 * (dot00 * dot11 - dot01 * dot01); float u = (dot11 * dot02 - dot01 * dot12)
568 * * invDenom; float v = (dot00 * dot12 - dot01 * dot02) * invDenom; return (u
569 * > 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 * see Berthold K. P. Horn,
587 * "Closed-form solution of absolute orientation using unit quaternions" J.
588 * Opt. Soc. Amer. A, 1987, Vol. 4, pp. 629-642
589 * http://www.opticsinfobase.org/viewmedia.cfm?uri=josaa-4-4-629&seq=0
592 * A similar treatment was developed independently (and later!) by G.
593 * Kramer, in G. R. Kramer,
594 * "Superposition of Molecular Structures Using Quaternions" Molecular
595 * Simulation, 1991, Vol. 7, pp. 113-119.
597 * In that treatment there is a lot of unnecessary calculation along the
598 * trace of matrix M (eqn 20). I'm not sure why the extra x^2 + y^2 + z^2 +
599 * x'^2 + y'^2 + z'^2 is in there, but they are unnecessary and only
600 * contribute to larger numerical averaging errors and additional processing
601 * time, as far as I can tell. Adding aI, where a is a scalar and I is the
602 * 4x4 identity just offsets the eigenvalues but doesn't change the
605 * and Lydia E. Kavraki, "Molecular Distance Measures"
606 * http://cnx.org/content/m11608/latest/
610 retStddev[1] = Float.NaN;
612 P3[] ptsA = centerAndPoints[0];
613 P3[] ptsB = centerAndPoints[1];
614 int nPts = ptsA.length - 1;
615 if (nPts < 2 || ptsA.length != ptsB.length)
617 double Sxx = 0, Sxy = 0, Sxz = 0, Syx = 0, Syy = 0, Syz = 0, Szx = 0, Szy = 0, Szz = 0;
622 for (int i = nPts + 1; --i >= 1;) {
623 ptA.sub2(ptsA[i], ptA0);
624 ptB.sub2(ptsB[i], ptB0);
625 Sxx += (double) ptA.x * (double) ptB.x;
626 Sxy += (double) ptA.x * (double) ptB.y;
627 Sxz += (double) ptA.x * (double) ptB.z;
628 Syx += (double) ptA.y * (double) ptB.x;
629 Syy += (double) ptA.y * (double) ptB.y;
630 Syz += (double) ptA.y * (double) ptB.z;
631 Szx += (double) ptA.z * (double) ptB.x;
632 Szy += (double) ptA.z * (double) ptB.y;
633 Szz += (double) ptA.z * (double) ptB.z;
635 retStddev[0] = getRmsd(centerAndPoints, q);
636 double[][] N = new double[4][4];
637 N[0][0] = Sxx + Syy + Szz;
638 N[0][1] = N[1][0] = Syz - Szy;
639 N[0][2] = N[2][0] = Szx - Sxz;
640 N[0][3] = N[3][0] = Sxy - Syx;
642 N[1][1] = Sxx - Syy - Szz;
643 N[1][2] = N[2][1] = Sxy + Syx;
644 N[1][3] = N[3][1] = Szx + Sxz;
646 N[2][2] = -Sxx + Syy - Szz;
647 N[2][3] = N[3][2] = Syz + Szy;
649 N[3][3] = -Sxx - Syy + Szz;
651 // this construction prevents JavaScript from requiring preloading of Eigen
653 float[] v = ((EigenInterface) Interface.getInterface("javajs.util.Eigen"))
654 .setM(N).getEigenvectorsFloatTransposed()[3];
655 q = Quat.newP4(P4.new4(v[1], v[2], v[3], v[0]));
656 retStddev[1] = getRmsd(centerAndPoints, q);
661 * Fills a 4x4 matrix with rotation-translation of mapped points A to B.
662 * If centerA is null, this is a standard 4x4 rotation-translation matrix;
663 * otherwise, this 4x4 matrix is a rotation around a vector through the center of ptsA,
664 * and centerA is filled with that center;
665 * Prior to Jmol 14.3.12_2014.02.14, when used from the JmolScript compare() function,
666 * this method returned the second of these options instead of the first.
670 * @param m 4x4 matrix to be returned
671 * @param centerA return center of rotation; if null, then standard 4x4 matrix is returned
674 public static float getTransformMatrix4(Lst<P3> ptsA, Lst<P3> ptsB, M4 m,
676 P3[] cptsA = getCenterAndPoints(ptsA);
677 P3[] cptsB = getCenterAndPoints(ptsB);
678 float[] retStddev = new float[2];
679 Quat q = calculateQuaternionRotation(new P3[][] { cptsA, cptsB },
681 M3 r = q.getMatrix();
685 centerA.setT(cptsA[0]);
686 V3 t = V3.newVsub(cptsB[0], cptsA[0]);
692 * from a list of points, create an array that includes the center
693 * point as the first point. This array is used as a starting point for
694 * a quaternion analysis of superposition.
697 * @return array of points with first point center
699 public static P3[] getCenterAndPoints(Lst<P3> vPts) {
701 P3[] pts = new P3[n + 1];
704 for (int i = 0; i < n; i++) {
705 pts[0].add(pts[i + 1] = vPts.get(i));
707 pts[0].scale(1f / n);
712 public static float getRmsd(P3[][] centerAndPoints, Quat q) {
714 P3[] ptsA = centerAndPoints[0];
715 P3[] ptsB = centerAndPoints[1];
718 int n = ptsA.length - 1;
719 P3 ptAnew = new P3();
721 for (int i = n + 1; --i >= 1;) {
722 ptAnew.sub2(ptsA[i], cA);
723 q.transform2(ptAnew, ptAnew).add(cB);
724 sum2 += ptAnew.distanceSquared(ptsB[i]);
726 return (float) Math.sqrt(sum2 / n);