Skip to content

Commit 586acf0

Browse files
BobHansonBobHanson
committed
Measure upgrade for getPlaneProjection
1 parent 4fd65bd commit 586acf0

File tree

6 files changed

+140
-138
lines changed

6 files changed

+140
-138
lines changed

‎sources/net.sf.j2s.core/dist/swingjs/SwingJS-site.zip‎

242 Bytes
Binary file not shown.

‎sources/net.sf.j2s.core/dist/swingjs/timestamp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
20220131135154
1+
20220223044927

‎sources/net.sf.j2s.core/dist/swingjs/ver/3.3.1/SwingJS-site.zip‎

242 Bytes
Binary file not shown.

‎sources/net.sf.j2s.core/dist/swingjs/ver/3.3.1/timestamp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
20220131135154
1+
20220223044927

‎sources/net.sf.j2s.java.core/dist/SwingJS-site.zip‎

3 Bytes
Binary file not shown.

‎sources/net.sf.j2s.java.core/src/javajs/util/Measure.java‎

Lines changed: 138 additions & 136 deletions
Original file line numberDiff line numberDiff line change
@@ -251,19 +251,22 @@ public static float getNormalThroughPoints(T3 pointA, T3 pointB,
251251
}
252252

253253
/**
254-
* Project a point onto a plane, also returning the normal vector.
254+
* Project a point onto a plane, also returning the normal vector and the directed distance to the plane.
255255
*
256256
* @param pt
257257
* @param plane
258258
* @param retPtProj returned pt (can be pt)
259259
* @param retNorm returned normal vector
260+
* @return directed distance to plane
260261
*/
261-
public static void getPlaneProjection(T3 pt, P4 plane, T3 retPtProj, V3 retNorm) {
262+
public static float getPlaneProjection(T3 pt, P4 plane, T3 retPtProj, V3 retNorm) {
262263
float dist = distanceToPlane(plane, pt);
263264
retNorm.set(plane.x, plane.y, plane.z);
264265
retNorm.normalize();
265-
retNorm.scale(-dist);
266-
retPtProj.add2(pt, retNorm);
266+
if (dist > 0)
267+
retNorm.scale(-1);
268+
retPtProj.scaleAdd2(Math.abs(dist), retNorm, pt);
269+
return dist;
267270
}
268271

269272
/**
@@ -561,125 +564,125 @@ public static P3 getIntersection(P3 pt1, V3 v,
561564
return ptRet;
562565
}
563566

564-
/*
565-
* public static Point3f getTriangleIntersection(Point3f a1, Point3f a2,
566-
* Point3f a3, Point4f plane, Point3f b1, Point3f b2, Point3f b3, Vector3f
567-
* vNorm, Vector3f vTemp, Point3f ptRet, Point3f ptTemp, Vector3f vTemp2,
568-
* Point4f pTemp, Vector3f vTemp3) {
569-
*
570-
* if (getTriangleIntersection(b1, b2, a1, a2, a3, vTemp, plane, vNorm,
571-
* vTemp2, vTemp3, ptRet, ptTemp)) return ptRet; if
572-
* (getTriangleIntersection(b2, b3, a1, a2, a3, vTemp, plane, vNorm, vTemp2,
573-
* vTemp3, ptRet, ptTemp)) return ptRet; if (getTriangleIntersection(b3, b1,
574-
* a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp)) return
575-
* ptRet; return null; }
576-
*/
577-
/*
578-
* public static boolean getTriangleIntersection(Point3f b1, Point3f b2,
579-
* Point3f a1, Point3f a2, Point3f a3, Vector3f vTemp, Point4f plane, Vector3f
580-
* vNorm, Vector3f vTemp2, Vector3f vTemp3, Point3f ptRet, Point3f ptTemp) {
581-
* if (distanceToPlane(plane, b1) * distanceToPlane(plane, b2) >= 0) return
582-
* false; vTemp.sub(b2, b1); vTemp.normalize(); if (getIntersection(b1, vTemp,
583-
* plane, ptRet, vNorm, vTemp2) != null) { if (isInTriangle(ptRet, a1, a2, a3,
584-
* vTemp, vTemp2, vTemp3)) return true; } return false; } private static
585-
* boolean isInTriangle(Point3f p, Point3f a, Point3f b, Point3f c, Vector3f
586-
* v0, Vector3f v1, Vector3f v2) { // from
587-
* http://www.blackpawn.com/texts/pointinpoly/default.html // Compute
588-
* barycentric coordinates v0.sub(c, a); v1.sub(b, a); v2.sub(p, a); float
589-
* dot00 = v0.dot(v0); float dot01 = v0.dot(v1); float dot02 = v0.dot(v2);
590-
* float dot11 = v1.dot(v1); float dot12 = v1.dot(v2); float invDenom = 1 /
591-
* (dot00 * dot11 - dot01 * dot01); float u = (dot11 * dot02 - dot01 * dot12)
592-
* * invDenom; float v = (dot00 * dot12 - dot01 * dot02) * invDenom; return (u
593-
* > 0 && v > 0 && u + v < 1); }
594-
*/
595-
596-
/**
597-
* Closed-form solution of absolute orientation requiring 1:1 mapping of
598-
* positions.
599-
*
600-
* @param centerAndPoints
601-
* @param retStddev
602-
* @return unit quaternion representation rotation
603-
*
604-
* @author hansonr Bob Hanson
605-
*
606-
*/
607-
public static Quat calculateQuaternionRotation(P3[][] centerAndPoints,
608-
float[] retStddev) {
609-
/*
610-
* see Berthold K. P. Horn,
611-
* "Closed-form solution of absolute orientation using unit quaternions" J.
612-
* Opt. Soc. Amer. A, 1987, Vol. 4, pp. 629-642
613-
* http://www.opticsinfobase.org/viewmedia.cfm?uri=josaa-4-4-629&seq=0
614-
*
615-
*
616-
* A similar treatment was developed independently (and later!) by G.
617-
* Kramer, in G. R. Kramer,
618-
* "Superposition of Molecular Structures Using Quaternions" Molecular
619-
* Simulation, 1991, Vol. 7, pp. 113-119.
620-
*
621-
* In that treatment there is a lot of unnecessary calculation along the
622-
* trace of matrix M (eqn 20). I'm not sure why the extra x^2 + y^2 + z^2 +
623-
* x'^2 + y'^2 + z'^2 is in there, but they are unnecessary and only
624-
* contribute to larger numerical averaging errors and additional processing
625-
* time, as far as I can tell. Adding aI, where a is a scalar and I is the
626-
* 4x4 identity just offsets the eigenvalues but doesn't change the
627-
* eigenvectors.
628-
*
629-
* and Lydia E. Kavraki, "Molecular Distance Measures"
630-
* http://cnx.org/content/m11608/latest/
631-
*/
632-
633-
634-
retStddev[1] = Float.NaN;
635-
Quat q = new Quat();
636-
P3[] ptsA = centerAndPoints[0];
637-
P3[] ptsB = centerAndPoints[1];
638-
int nPts = ptsA.length - 1;
639-
if (nPts < 2 || ptsA.length != ptsB.length)
640-
return q;
641-
double Sxx = 0, Sxy = 0, Sxz = 0, Syx = 0, Syy = 0, Syz = 0, Szx = 0, Szy = 0, Szz = 0;
642-
P3 ptA = new P3();
643-
P3 ptB = new P3();
644-
P3 ptA0 = ptsA[0];
645-
P3 ptB0 = ptsB[0];
646-
for (int i = nPts + 1; --i >= 1;) {
647-
ptA.sub2(ptsA[i], ptA0);
648-
ptB.sub2(ptsB[i], ptB0);
649-
Sxx += (double) ptA.x * (double) ptB.x;
650-
Sxy += (double) ptA.x * (double) ptB.y;
651-
Sxz += (double) ptA.x * (double) ptB.z;
652-
Syx += (double) ptA.y * (double) ptB.x;
653-
Syy += (double) ptA.y * (double) ptB.y;
654-
Syz += (double) ptA.y * (double) ptB.z;
655-
Szx += (double) ptA.z * (double) ptB.x;
656-
Szy += (double) ptA.z * (double) ptB.y;
657-
Szz += (double) ptA.z * (double) ptB.z;
658-
}
659-
retStddev[0] = getRmsd(centerAndPoints, q);
660-
double[][] N = new double[4][4];
661-
N[0][0] = Sxx + Syy + Szz;
662-
N[0][1] = N[1][0] = Syz - Szy;
663-
N[0][2] = N[2][0] = Szx - Sxz;
664-
N[0][3] = N[3][0] = Sxy - Syx;
665-
666-
N[1][1] = Sxx - Syy - Szz;
667-
N[1][2] = N[2][1] = Sxy + Syx;
668-
N[1][3] = N[3][1] = Szx + Sxz;
669-
670-
N[2][2] = -Sxx + Syy - Szz;
671-
N[2][3] = N[3][2] = Syz + Szy;
672-
673-
N[3][3] = -Sxx - Syy + Szz;
674-
675-
// this construction prevents JavaScript from requiring preloading of Eigen
676-
677-
float[] v = ((EigenInterface) Interface.getInterface("javajs.util.Eigen"))
678-
.setM(N).getEigenvectorsFloatTransposed()[3];
679-
q = Quat.newP4(P4.new4(v[1], v[2], v[3], v[0]));
680-
retStddev[1] = getRmsd(centerAndPoints, q);
681-
return q;
682-
}
567+
/*
568+
* public static Point3f getTriangleIntersection(Point3f a1, Point3f a2,
569+
* Point3f a3, Point4f plane, Point3f b1, Point3f b2, Point3f b3, Vector3f
570+
* vNorm, Vector3f vTemp, Point3f ptRet, Point3f ptTemp, Vector3f vTemp2,
571+
* Point4f pTemp, Vector3f vTemp3) {
572+
*
573+
* if (getTriangleIntersection(b1, b2, a1, a2, a3, vTemp, plane, vNorm,
574+
* vTemp2, vTemp3, ptRet, ptTemp)) return ptRet; if
575+
* (getTriangleIntersection(b2, b3, a1, a2, a3, vTemp, plane, vNorm, vTemp2,
576+
* vTemp3, ptRet, ptTemp)) return ptRet; if (getTriangleIntersection(b3, b1,
577+
* a1, a2, a3, vTemp, plane, vNorm, vTemp2, vTemp3, ptRet, ptTemp)) return
578+
* ptRet; return null; }
579+
*/
580+
/*
581+
* public static boolean getTriangleIntersection(Point3f b1, Point3f b2,
582+
* Point3f a1, Point3f a2, Point3f a3, Vector3f vTemp, Point4f plane, Vector3f
583+
* vNorm, Vector3f vTemp2, Vector3f vTemp3, Point3f ptRet, Point3f ptTemp) {
584+
* if (distanceToPlane(plane, b1) * distanceToPlane(plane, b2) >= 0) return
585+
* false; vTemp.sub(b2, b1); vTemp.normalize(); if (getIntersection(b1, vTemp,
586+
* plane, ptRet, vNorm, vTemp2) != null) { if (isInTriangle(ptRet, a1, a2, a3,
587+
* vTemp, vTemp2, vTemp3)) return true; } return false; } private static
588+
* boolean isInTriangle(Point3f p, Point3f a, Point3f b, Point3f c, Vector3f
589+
* v0, Vector3f v1, Vector3f v2) { // from
590+
* http://www.blackpawn.com/texts/pointinpoly/default.html // Compute
591+
* barycentric coordinates v0.sub(c, a); v1.sub(b, a); v2.sub(p, a); float
592+
* dot00 = v0.dot(v0); float dot01 = v0.dot(v1); float dot02 = v0.dot(v2);
593+
* float dot11 = v1.dot(v1); float dot12 = v1.dot(v2); float invDenom = 1 /
594+
* (dot00 * dot11 - dot01 * dot01); float u = (dot11 * dot02 - dot01 * dot12)
595+
* * invDenom; float v = (dot00 * dot12 - dot01 * dot02) * invDenom; return (u
596+
* > 0 && v > 0 && u + v < 1); }
597+
*/
598+
599+
/**
600+
* Closed-form solution of absolute orientation requiring 1:1 mapping of
601+
* positions.
602+
*
603+
* @param centerAndPoints
604+
* @param retStddev
605+
* @return unit quaternion representation rotation
606+
*
607+
* @author hansonr Bob Hanson
608+
*
609+
*/
610+
public static Quat calculateQuaternionRotation(P3[][] centerAndPoints,
611+
float[] retStddev) {
612+
/*
613+
* see Berthold K. P. Horn,
614+
* "Closed-form solution of absolute orientation using unit quaternions" J.
615+
* Opt. Soc. Amer. A, 1987, Vol. 4, pp. 629-642
616+
* http://www.opticsinfobase.org/viewmedia.cfm?uri=josaa-4-4-629&seq=0
617+
*
618+
*
619+
* A similar treatment was developed independently (and later!) by G.
620+
* Kramer, in G. R. Kramer,
621+
* "Superposition of Molecular Structures Using Quaternions" Molecular
622+
* Simulation, 1991, Vol. 7, pp. 113-119.
623+
*
624+
* In that treatment there is a lot of unnecessary calculation along the
625+
* trace of matrix M (eqn 20). I'm not sure why the extra x^2 + y^2 + z^2 +
626+
* x'^2 + y'^2 + z'^2 is in there, but they are unnecessary and only
627+
* contribute to larger numerical averaging errors and additional processing
628+
* time, as far as I can tell. Adding aI, where a is a scalar and I is the
629+
* 4x4 identity just offsets the eigenvalues but doesn't change the
630+
* eigenvectors.
631+
*
632+
* and Lydia E. Kavraki, "Molecular Distance Measures"
633+
* http://cnx.org/content/m11608/latest/
634+
*/
635+
636+
637+
retStddev[1] = Float.NaN;
638+
Quat q = new Quat();
639+
P3[] ptsA = centerAndPoints[0];
640+
P3[] ptsB = centerAndPoints[1];
641+
int nPts = ptsA.length - 1;
642+
if (nPts < 2 || ptsA.length != ptsB.length)
643+
return q;
644+
double Sxx = 0, Sxy = 0, Sxz = 0, Syx = 0, Syy = 0, Syz = 0, Szx = 0, Szy = 0, Szz = 0;
645+
P3 ptA = new P3();
646+
P3 ptB = new P3();
647+
P3 ptA0 = ptsA[0];
648+
P3 ptB0 = ptsB[0];
649+
for (int i = nPts + 1; --i >= 1;) {
650+
ptA.sub2(ptsA[i], ptA0);
651+
ptB.sub2(ptsB[i], ptB0);
652+
Sxx += (double) ptA.x * (double) ptB.x;
653+
Sxy += (double) ptA.x * (double) ptB.y;
654+
Sxz += (double) ptA.x * (double) ptB.z;
655+
Syx += (double) ptA.y * (double) ptB.x;
656+
Syy += (double) ptA.y * (double) ptB.y;
657+
Syz += (double) ptA.y * (double) ptB.z;
658+
Szx += (double) ptA.z * (double) ptB.x;
659+
Szy += (double) ptA.z * (double) ptB.y;
660+
Szz += (double) ptA.z * (double) ptB.z;
661+
}
662+
retStddev[0] = getRmsd(centerAndPoints, q);
663+
double[][] N = new double[4][4];
664+
N[0][0] = Sxx + Syy + Szz;
665+
N[0][1] = N[1][0] = Syz - Szy;
666+
N[0][2] = N[2][0] = Szx - Sxz;
667+
N[0][3] = N[3][0] = Sxy - Syx;
668+
669+
N[1][1] = Sxx - Syy - Szz;
670+
N[1][2] = N[2][1] = Sxy + Syx;
671+
N[1][3] = N[3][1] = Szx + Sxz;
672+
673+
N[2][2] = -Sxx + Syy - Szz;
674+
N[2][3] = N[3][2] = Syz + Szy;
675+
676+
N[3][3] = -Sxx - Syy + Szz;
677+
678+
// this construction prevents JavaScript from requiring preloading of Eigen
679+
680+
float[] v = ((EigenInterface) Interface.getInterface("javajs.util.Eigen"))
681+
.setM(N).getEigenvectorsFloatTransposed()[3];
682+
q = Quat.newP4(P4.new4(v[1], v[2], v[3], v[0]));
683+
retStddev[1] = getRmsd(centerAndPoints, q);
684+
return q;
685+
}
683686

684687
/**
685688
* Fills a 4x4 matrix with rotation-translation of mapped points A to B.
@@ -720,18 +723,18 @@ public static float getTransformMatrix4(Lst<P3> ptsA, Lst<P3> ptsB, M4 m,
720723
* @param vPts
721724
* @return array of points with first point center
722725
*/
723-
public static P3[] getCenterAndPoints(Lst<P3> vPts) {
724-
int n = vPts.size();
725-
P3[] pts = new P3[n + 1];
726-
pts[0] = new P3();
727-
if (n > 0) {
728-
for (int i = 0; i < n; i++) {
729-
pts[0].add(pts[i + 1] = vPts.get(i));
730-
}
731-
pts[0].scale(1f / n);
732-
}
733-
return pts;
734-
}
726+
public static P3[] getCenterAndPoints(Lst<P3> vPts) {
727+
int n = vPts.size();
728+
P3[] pts = new P3[n + 1];
729+
pts[0] = new P3();
730+
if (n > 0) {
731+
for (int i = 0; i < n; i++) {
732+
pts[0].add(pts[i + 1] = vPts.get(i));
733+
}
734+
pts[0].scale(1f / n);
735+
}
736+
return pts;
737+
}
735738

736739
public static float getRmsd(P3[][] centerAndPoints, Quat q) {
737740
double sum2 = 0;
@@ -978,7 +981,6 @@ static void testRnd() {
978981
pts[i] = new P3();
979982
P3 p = rndPt();
980983
getPlaneProjection(p, plane, ptProj, vNorm );
981-
vNorm.scale(1/vNorm.length());
982984
pts[i].setT(ptProj);
983985
float d = (float)Math.random()*0.1f;
984986
pts[i].scaleAdd2(d, vNorm, ptProj);

0 commit comments

Comments
 (0)