Merge branch 'master' of https://source.jalview.org/git/jalviewjs.git
[jalviewjs.git] / site / j2s / JU / Measure.js
1 Clazz.declarePackage ("JU");
2 Clazz.load (["JU.V3"], "JU.Measure", ["java.lang.Float", "javajs.api.Interface", "JU.Lst", "$.P3", "$.P4", "$.Quat"], function () {
3 c$ = Clazz.declareType (JU, "Measure");
4 c$.computeAngle = Clazz.defineMethod (c$, "computeAngle", 
5 function (pointA, pointB, pointC, vectorBA, vectorBC, asDegrees) {
6 vectorBA.sub2 (pointA, pointB);
7 vectorBC.sub2 (pointC, pointB);
8 var angle = vectorBA.angle (vectorBC);
9 return (asDegrees ? angle / 0.017453292 : angle);
10 }, "JU.T3,JU.T3,JU.T3,JU.V3,JU.V3,~B");
11 c$.computeAngleABC = Clazz.defineMethod (c$, "computeAngleABC", 
12 function (pointA, pointB, pointC, asDegrees) {
13 var vectorBA =  new JU.V3 ();
14 var vectorBC =  new JU.V3 ();
15 return JU.Measure.computeAngle (pointA, pointB, pointC, vectorBA, vectorBC, asDegrees);
16 }, "JU.T3,JU.T3,JU.T3,~B");
17 c$.computeTorsion = Clazz.defineMethod (c$, "computeTorsion", 
18 function (p1, p2, p3, p4, asDegrees) {
19 var ijx = p1.x - p2.x;
20 var ijy = p1.y - p2.y;
21 var ijz = p1.z - p2.z;
22 var kjx = p3.x - p2.x;
23 var kjy = p3.y - p2.y;
24 var kjz = p3.z - p2.z;
25 var klx = p3.x - p4.x;
26 var kly = p3.y - p4.y;
27 var klz = p3.z - p4.z;
28 var ax = ijy * kjz - ijz * kjy;
29 var ay = ijz * kjx - ijx * kjz;
30 var az = ijx * kjy - ijy * kjx;
31 var cx = kjy * klz - kjz * kly;
32 var cy = kjz * klx - kjx * klz;
33 var cz = kjx * kly - kjy * klx;
34 var ai2 = 1 / (ax * ax + ay * ay + az * az);
35 var ci2 = 1 / (cx * cx + cy * cy + cz * cz);
36 var ai = Math.sqrt (ai2);
37 var ci = Math.sqrt (ci2);
38 var denom = ai * ci;
39 var cross = ax * cx + ay * cy + az * cz;
40 var cosang = cross * denom;
41 if (cosang > 1) {
42 cosang = 1;
43 }if (cosang < -1) {
44 cosang = -1;
45 }var torsion = Math.acos (cosang);
46 var dot = ijx * cx + ijy * cy + ijz * cz;
47 var absDot = Math.abs (dot);
48 torsion = (dot / absDot > 0) ? torsion : -torsion;
49 return (asDegrees ? torsion / 0.017453292 : torsion);
50 }, "JU.T3,JU.T3,JU.T3,JU.T3,~B");
51 c$.computeHelicalAxis = Clazz.defineMethod (c$, "computeHelicalAxis", 
52 function (a, b, dq) {
53 var vab =  new JU.V3 ();
54 vab.sub2 (b, a);
55 var theta = dq.getTheta ();
56 var n = dq.getNormal ();
57 var v_dot_n = vab.dot (n);
58 if (Math.abs (v_dot_n) < 0.0001) v_dot_n = 0;
59 var va_prime_d =  new JU.V3 ();
60 va_prime_d.cross (vab, n);
61 if (va_prime_d.dot (va_prime_d) != 0) va_prime_d.normalize ();
62 var vda =  new JU.V3 ();
63 var vcb = JU.V3.newV (n);
64 if (v_dot_n == 0) v_dot_n = 1.4E-45;
65 vcb.scale (v_dot_n);
66 vda.sub2 (vcb, vab);
67 vda.scale (0.5);
68 va_prime_d.scale (theta == 0 ? 0 : (vda.length () / Math.tan (theta / 2 / 180 * 3.141592653589793)));
69 var r = JU.V3.newV (va_prime_d);
70 if (theta != 0) r.add (vda);
71 var pt_a_prime = JU.P3.newP (a);
72 pt_a_prime.sub (r);
73 if (v_dot_n != 1.4E-45) n.scale (v_dot_n);
74 var pt_b_prime = JU.P3.newP (pt_a_prime);
75 pt_b_prime.add (n);
76 theta = JU.Measure.computeTorsion (a, pt_a_prime, pt_b_prime, b, true);
77 if (Float.isNaN (theta) || r.length () < 0.0001) theta = dq.getThetaDirectedV (n);
78 var residuesPerTurn = Math.abs (theta == 0 ? 0 : 360 / theta);
79 var pitch = Math.abs (v_dot_n == 1.4E-45 ? 0 : n.length () * (theta == 0 ? 1 : 360 / theta));
80 return  Clazz.newArray (-1, [pt_a_prime, n, r, JU.P3.new3 (theta, pitch, residuesPerTurn), pt_b_prime]);
81 }, "JU.P3,JU.P3,JU.Quat");
82 c$.getPlaneThroughPoints = Clazz.defineMethod (c$, "getPlaneThroughPoints", 
83 function (pointA, pointB, pointC, vNorm, vAB, plane) {
84 var w = JU.Measure.getNormalThroughPoints (pointA, pointB, pointC, vNorm, vAB);
85 plane.set4 (vNorm.x, vNorm.y, vNorm.z, w);
86 return plane;
87 }, "JU.T3,JU.T3,JU.T3,JU.V3,JU.V3,JU.P4");
88 c$.getPlaneThroughPoint = Clazz.defineMethod (c$, "getPlaneThroughPoint", 
89 function (pt, normal, plane) {
90 plane.set4 (normal.x, normal.y, normal.z, -normal.dot (pt));
91 }, "JU.T3,JU.V3,JU.P4");
92 c$.distanceToPlane = Clazz.defineMethod (c$, "distanceToPlane", 
93 function (plane, pt) {
94 return (plane == null ? NaN : (plane.x * pt.x + plane.y * pt.y + plane.z * pt.z + plane.w) / Math.sqrt (plane.x * plane.x + plane.y * plane.y + plane.z * plane.z));
95 }, "JU.P4,JU.T3");
96 c$.distanceToPlaneD = Clazz.defineMethod (c$, "distanceToPlaneD", 
97 function (plane, d, pt) {
98 return (plane == null ? NaN : (plane.x * pt.x + plane.y * pt.y + plane.z * pt.z + plane.w) / d);
99 }, "JU.P4,~N,JU.P3");
100 c$.distanceToPlaneV = Clazz.defineMethod (c$, "distanceToPlaneV", 
101 function (norm, w, pt) {
102 return (norm == null ? NaN : (norm.x * pt.x + norm.y * pt.y + norm.z * pt.z + w) / Math.sqrt (norm.x * norm.x + norm.y * norm.y + norm.z * norm.z));
103 }, "JU.V3,~N,JU.P3");
104 c$.calcNormalizedNormal = Clazz.defineMethod (c$, "calcNormalizedNormal", 
105 function (pointA, pointB, pointC, vNormNorm, vAB) {
106 vAB.sub2 (pointB, pointA);
107 vNormNorm.sub2 (pointC, pointA);
108 vNormNorm.cross (vAB, vNormNorm);
109 vNormNorm.normalize ();
110 }, "JU.T3,JU.T3,JU.T3,JU.V3,JU.V3");
111 c$.getDirectedNormalThroughPoints = Clazz.defineMethod (c$, "getDirectedNormalThroughPoints", 
112 function (pointA, pointB, pointC, ptRef, vNorm, vAB) {
113 var nd = JU.Measure.getNormalThroughPoints (pointA, pointB, pointC, vNorm, vAB);
114 if (ptRef != null) {
115 var pt0 = JU.P3.newP (pointA);
116 pt0.add (vNorm);
117 var d = pt0.distance (ptRef);
118 pt0.sub2 (pointA, vNorm);
119 if (d > pt0.distance (ptRef)) {
120 vNorm.scale (-1);
121 nd = -nd;
122 }}return nd;
123 }, "JU.T3,JU.T3,JU.T3,JU.T3,JU.V3,JU.V3");
124 c$.getNormalThroughPoints = Clazz.defineMethod (c$, "getNormalThroughPoints", 
125 function (pointA, pointB, pointC, vNorm, vTemp) {
126 JU.Measure.calcNormalizedNormal (pointA, pointB, pointC, vNorm, vTemp);
127 vTemp.setT (pointA);
128 return -vTemp.dot (vNorm);
129 }, "JU.T3,JU.T3,JU.T3,JU.V3,JU.V3");
130 c$.getPlaneProjection = Clazz.defineMethod (c$, "getPlaneProjection", 
131 function (pt, plane, ptProj, vNorm) {
132 var dist = JU.Measure.distanceToPlane (plane, pt);
133 vNorm.set (plane.x, plane.y, plane.z);
134 vNorm.normalize ();
135 vNorm.scale (-dist);
136 ptProj.add2 (pt, vNorm);
137 }, "JU.P3,JU.P4,JU.P3,JU.V3");
138 c$.getNormalFromCenter = Clazz.defineMethod (c$, "getNormalFromCenter", 
139 function (ptCenter, ptA, ptB, ptC, isOutward, normal) {
140 var vAB =  new JU.V3 ();
141 var d = JU.Measure.getNormalThroughPoints (ptA, ptB, ptC, normal, vAB);
142 var isReversed = (JU.Measure.distanceToPlaneV (normal, d, ptCenter) > 0);
143 if (isReversed == isOutward) normal.scale (-1.0);
144 return !isReversed;
145 }, "JU.P3,JU.P3,JU.P3,JU.P3,~B,JU.V3");
146 c$.getNormalToLine = Clazz.defineMethod (c$, "getNormalToLine", 
147 function (pointA, pointB, vNormNorm) {
148 vNormNorm.sub2 (pointA, pointB);
149 vNormNorm.cross (vNormNorm, JU.Measure.axisY);
150 vNormNorm.normalize ();
151 if (Float.isNaN (vNormNorm.x)) vNormNorm.set (1, 0, 0);
152 }, "JU.P3,JU.P3,JU.V3");
153 c$.getBisectingPlane = Clazz.defineMethod (c$, "getBisectingPlane", 
154 function (pointA, vAB, ptTemp, vTemp, plane) {
155 ptTemp.scaleAdd2 (0.5, vAB, pointA);
156 vTemp.setT (vAB);
157 vTemp.normalize ();
158 JU.Measure.getPlaneThroughPoint (ptTemp, vTemp, plane);
159 }, "JU.P3,JU.V3,JU.T3,JU.V3,JU.P4");
160 c$.projectOntoAxis = Clazz.defineMethod (c$, "projectOntoAxis", 
161 function (point, axisA, axisUnitVector, vectorProjection) {
162 vectorProjection.sub2 (point, axisA);
163 var projectedLength = vectorProjection.dot (axisUnitVector);
164 point.scaleAdd2 (projectedLength, axisUnitVector, axisA);
165 vectorProjection.sub2 (point, axisA);
166 }, "JU.P3,JU.P3,JU.V3,JU.V3");
167 c$.calcBestAxisThroughPoints = Clazz.defineMethod (c$, "calcBestAxisThroughPoints", 
168 function (points, axisA, axisUnitVector, vectorProjection, nTriesMax) {
169 var nPoints = points.length;
170 axisA.setT (points[0]);
171 axisUnitVector.sub2 (points[nPoints - 1], axisA);
172 axisUnitVector.normalize ();
173 JU.Measure.calcAveragePointN (points, nPoints, axisA);
174 var nTries = 0;
175 while (nTries++ < nTriesMax && JU.Measure.findAxis (points, nPoints, axisA, axisUnitVector, vectorProjection) > 0.001) {
176 }
177 var tempA = JU.P3.newP (points[0]);
178 JU.Measure.projectOntoAxis (tempA, axisA, axisUnitVector, vectorProjection);
179 axisA.setT (tempA);
180 }, "~A,JU.P3,JU.V3,JU.V3,~N");
181 c$.findAxis = Clazz.defineMethod (c$, "findAxis", 
182 function (points, nPoints, axisA, axisUnitVector, vectorProjection) {
183 var sumXiYi =  new JU.V3 ();
184 var vTemp =  new JU.V3 ();
185 var pt =  new JU.P3 ();
186 var ptProj =  new JU.P3 ();
187 var a = JU.V3.newV (axisUnitVector);
188 var sum_Xi2 = 0;
189 for (var i = nPoints; --i >= 0; ) {
190 pt.setT (points[i]);
191 ptProj.setT (pt);
192 JU.Measure.projectOntoAxis (ptProj, axisA, axisUnitVector, vectorProjection);
193 vTemp.sub2 (pt, ptProj);
194 vTemp.cross (vectorProjection, vTemp);
195 sumXiYi.add (vTemp);
196 sum_Xi2 += vectorProjection.lengthSquared ();
197 }
198 var m = JU.V3.newV (sumXiYi);
199 m.scale (1 / sum_Xi2);
200 vTemp.cross (m, axisUnitVector);
201 axisUnitVector.add (vTemp);
202 axisUnitVector.normalize ();
203 vTemp.sub2 (axisUnitVector, a);
204 return vTemp.length ();
205 }, "~A,~N,JU.P3,JU.V3,JU.V3");
206 c$.calcAveragePoint = Clazz.defineMethod (c$, "calcAveragePoint", 
207 function (pointA, pointB, pointC) {
208 pointC.set ((pointA.x + pointB.x) / 2, (pointA.y + pointB.y) / 2, (pointA.z + pointB.z) / 2);
209 }, "JU.P3,JU.P3,JU.P3");
210 c$.calcAveragePointN = Clazz.defineMethod (c$, "calcAveragePointN", 
211 function (points, nPoints, averagePoint) {
212 averagePoint.setT (points[0]);
213 for (var i = 1; i < nPoints; i++) averagePoint.add (points[i]);
214
215 averagePoint.scale (1 / nPoints);
216 }, "~A,~N,JU.P3");
217 c$.transformPoints = Clazz.defineMethod (c$, "transformPoints", 
218 function (vPts, m4, center) {
219 var v =  new JU.Lst ();
220 for (var i = 0; i < vPts.size (); i++) {
221 var pt = JU.P3.newP (vPts.get (i));
222 pt.sub (center);
223 m4.rotTrans (pt);
224 pt.add (center);
225 v.addLast (pt);
226 }
227 return v;
228 }, "JU.Lst,JU.M4,JU.P3");
229 c$.isInTetrahedron = Clazz.defineMethod (c$, "isInTetrahedron", 
230 function (pt, ptA, ptB, ptC, ptD, plane, vTemp, vTemp2, fullyEnclosed) {
231 var b = (JU.Measure.distanceToPlane (JU.Measure.getPlaneThroughPoints (ptC, ptD, ptA, vTemp, vTemp2, plane), pt) >= 0);
232 if (b != (JU.Measure.distanceToPlane (JU.Measure.getPlaneThroughPoints (ptA, ptD, ptB, vTemp, vTemp2, plane), pt) >= 0)) return false;
233 if (b != (JU.Measure.distanceToPlane (JU.Measure.getPlaneThroughPoints (ptB, ptD, ptC, vTemp, vTemp2, plane), pt) >= 0)) return false;
234 var d = JU.Measure.distanceToPlane (JU.Measure.getPlaneThroughPoints (ptA, ptB, ptC, vTemp, vTemp2, plane), pt);
235 if (fullyEnclosed) return (b == (d >= 0));
236 var d1 = JU.Measure.distanceToPlane (plane, ptD);
237 return d1 * d <= 0 || Math.abs (d1) > Math.abs (d);
238 }, "JU.P3,JU.P3,JU.P3,JU.P3,JU.P3,JU.P4,JU.V3,JU.V3,~B");
239 c$.getIntersectionPP = Clazz.defineMethod (c$, "getIntersectionPP", 
240 function (plane1, plane2) {
241 var a1 = plane1.x;
242 var b1 = plane1.y;
243 var c1 = plane1.z;
244 var d1 = plane1.w;
245 var a2 = plane2.x;
246 var b2 = plane2.y;
247 var c2 = plane2.z;
248 var d2 = plane2.w;
249 var norm1 = JU.V3.new3 (a1, b1, c1);
250 var norm2 = JU.V3.new3 (a2, b2, c2);
251 var nxn =  new JU.V3 ();
252 nxn.cross (norm1, norm2);
253 var ax = Math.abs (nxn.x);
254 var ay = Math.abs (nxn.y);
255 var az = Math.abs (nxn.z);
256 var x;
257 var y;
258 var z;
259 var diff;
260 var type = (ax > ay ? (ax > az ? 1 : 3) : ay > az ? 2 : 3);
261 switch (type) {
262 case 1:
263 x = 0;
264 diff = (b1 * c2 - b2 * c1);
265 if (Math.abs (diff) < 0.01) return null;
266 y = (c1 * d2 - c2 * d1) / diff;
267 z = (b2 * d1 - d2 * b1) / diff;
268 break;
269 case 2:
270 diff = (a1 * c2 - a2 * c1);
271 if (Math.abs (diff) < 0.01) return null;
272 x = (c1 * d2 - c2 * d1) / diff;
273 y = 0;
274 z = (a2 * d1 - d2 * a1) / diff;
275 break;
276 case 3:
277 default:
278 diff = (a1 * b2 - a2 * b1);
279 if (Math.abs (diff) < 0.01) return null;
280 x = (b1 * d2 - b2 * d1) / diff;
281 y = (a2 * d1 - d2 * a1) / diff;
282 z = 0;
283 }
284 var list =  new JU.Lst ();
285 list.addLast (JU.P3.new3 (x, y, z));
286 nxn.normalize ();
287 list.addLast (nxn);
288 return list;
289 }, "JU.P4,JU.P4");
290 c$.getIntersection = Clazz.defineMethod (c$, "getIntersection", 
291 function (pt1, v, plane, ptRet, tempNorm, vTemp) {
292 JU.Measure.getPlaneProjection (pt1, plane, ptRet, tempNorm);
293 tempNorm.set (plane.x, plane.y, plane.z);
294 tempNorm.normalize ();
295 if (v == null) v = JU.V3.newV (tempNorm);
296 var l_dot_n = v.dot (tempNorm);
297 if (Math.abs (l_dot_n) < 0.01) return null;
298 vTemp.sub2 (ptRet, pt1);
299 ptRet.scaleAdd2 (vTemp.dot (tempNorm) / l_dot_n, v, pt1);
300 return ptRet;
301 }, "JU.P3,JU.V3,JU.P4,JU.P3,JU.V3,JU.V3");
302 c$.calculateQuaternionRotation = Clazz.defineMethod (c$, "calculateQuaternionRotation", 
303 function (centerAndPoints, retStddev) {
304 retStddev[1] = NaN;
305 var q =  new JU.Quat ();
306 if (centerAndPoints[0].length == 1 || centerAndPoints[0].length != centerAndPoints[1].length) return q;
307 var n = centerAndPoints[0].length - 1;
308 if (n < 2) return q;
309 var Sxx = 0;
310 var Sxy = 0;
311 var Sxz = 0;
312 var Syx = 0;
313 var Syy = 0;
314 var Syz = 0;
315 var Szx = 0;
316 var Szy = 0;
317 var Szz = 0;
318 var ptA =  new JU.P3 ();
319 var ptB =  new JU.P3 ();
320 for (var i = n + 1; --i >= 1; ) {
321 var aij = centerAndPoints[0][i];
322 var bij = centerAndPoints[1][i];
323 ptA.sub2 (aij, centerAndPoints[0][0]);
324 ptB.sub2 (bij, centerAndPoints[0][1]);
325 Sxx += ptA.x * ptB.x;
326 Sxy += ptA.x * ptB.y;
327 Sxz += ptA.x * ptB.z;
328 Syx += ptA.y * ptB.x;
329 Syy += ptA.y * ptB.y;
330 Syz += ptA.y * ptB.z;
331 Szx += ptA.z * ptB.x;
332 Szy += ptA.z * ptB.y;
333 Szz += ptA.z * ptB.z;
334 }
335 retStddev[0] = JU.Measure.getRmsd (centerAndPoints, q);
336 var N =  Clazz.newDoubleArray (4, 4, 0);
337 N[0][0] = Sxx + Syy + Szz;
338 N[0][1] = N[1][0] = Syz - Szy;
339 N[0][2] = N[2][0] = Szx - Sxz;
340 N[0][3] = N[3][0] = Sxy - Syx;
341 N[1][1] = Sxx - Syy - Szz;
342 N[1][2] = N[2][1] = Sxy + Syx;
343 N[1][3] = N[3][1] = Szx + Sxz;
344 N[2][2] = -Sxx + Syy - Szz;
345 N[2][3] = N[3][2] = Syz + Szy;
346 N[3][3] = -Sxx - Syy + Szz;
347 var v = (javajs.api.Interface.getInterface ("JU.Eigen")).setM (N).getEigenvectorsFloatTransposed ()[3];
348 q = JU.Quat.newP4 (JU.P4.new4 (v[1], v[2], v[3], v[0]));
349 retStddev[1] = JU.Measure.getRmsd (centerAndPoints, q);
350 return q;
351 }, "~A,~A");
352 c$.getTransformMatrix4 = Clazz.defineMethod (c$, "getTransformMatrix4", 
353 function (ptsA, ptsB, m, centerA) {
354 var cptsA = JU.Measure.getCenterAndPoints (ptsA);
355 var cptsB = JU.Measure.getCenterAndPoints (ptsB);
356 var retStddev =  Clazz.newFloatArray (2, 0);
357 var q = JU.Measure.calculateQuaternionRotation ( Clazz.newArray (-1, [cptsA, cptsB]), retStddev);
358 var r = q.getMatrix ();
359 if (centerA == null) r.rotate (cptsA[0]);
360  else centerA.setT (cptsA[0]);
361 var t = JU.V3.newVsub (cptsB[0], cptsA[0]);
362 m.setMV (r, t);
363 return retStddev[1];
364 }, "JU.Lst,JU.Lst,JU.M4,JU.P3");
365 c$.getCenterAndPoints = Clazz.defineMethod (c$, "getCenterAndPoints", 
366 function (vPts) {
367 var n = vPts.size ();
368 var pts =  new Array (n + 1);
369 pts[0] =  new JU.P3 ();
370 if (n > 0) {
371 for (var i = 0; i < n; i++) {
372 pts[0].add (pts[i + 1] = vPts.get (i));
373 }
374 pts[0].scale (1 / n);
375 }return pts;
376 }, "JU.Lst");
377 c$.getRmsd = Clazz.defineMethod (c$, "getRmsd", 
378 function (centerAndPoints, q) {
379 var sum2 = 0;
380 var ptsA = centerAndPoints[0];
381 var ptsB = centerAndPoints[1];
382 var cA = ptsA[0];
383 var cB = ptsB[0];
384 var n = ptsA.length - 1;
385 var ptAnew =  new JU.P3 ();
386 for (var i = n + 1; --i >= 1; ) {
387 ptAnew.sub2 (ptsA[i], cA);
388 q.transform2 (ptAnew, ptAnew).add (cB);
389 sum2 += ptAnew.distanceSquared (ptsB[i]);
390 }
391 return Math.sqrt (sum2 / n);
392 }, "~A,JU.Quat");
393 Clazz.defineStatics (c$,
394 "radiansPerDegree", (0.017453292519943295));
395 c$.axisY = c$.prototype.axisY = JU.V3.new3 (0, 1, 0);
396 });