Merge pull request #4 from josemduarte/superpose

lafita · web-flow · commit df93f78d9e25 · 2016-09-20T11:05:30.000+02:00
TestSuperPosition now tests rmsd of non-centered points.
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/CalcPoint.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/CalcPoint.java
@@ -2,6 +2,7 @@
 
 import javax.vecmath.Matrix4d;
 import javax.vecmath.Point3d;
+import javax.vecmath.Vector3d;
 
 import org.biojava.nbio.structure.jama.Matrix;
 
@@ -28,7 +29,7 @@ private CalcPoint() {
 	public static void center(Point3d[] x) {
 		Point3d center = centroid(x);
 		center.negate();
-		translate(center, x);
+		translate(new Vector3d(center), x);
 	}
 
 	/**
@@ -52,8 +53,8 @@ public static Point3d centroid(Point3d[] x) {
 	 * 
 	 * @param rotTrans
 	 *            4x4 transformation matrix
-	 * @param array
-	 *            of points. Point objects will be modified
+	 * @param x
+	 *            array of points. Point objects will be modified
 	 */
 	public static void transform(Matrix4d rotTrans, Point3d[] x) {
 		for (Point3d p : x) {
@@ -64,12 +65,12 @@ public static void transform(Matrix4d rotTrans, Point3d[] x) {
 	/**
 	 * Translate all points with a translation vector.
 	 * 
-	 * @param rotTrans
-	 *            4x4 transformation matrix
-	 * @param array
-	 *            of points. Point objects will be modified
+	 * @param trans
+	 *            the translation vector to apply
+	 * @param x
+	 *            array of points. Point objects will be modified
 	 */
-	public static void translate(Point3d trans, Point3d[] x) {
+	public static void translate(Vector3d trans, Point3d[] x) {
 		for (Point3d p : x) {
 			p.add(trans);
 		}
@@ -90,7 +91,7 @@ public static Point3d[] clonePoint3dArray(Point3d[] x) {
 		return clone;
 	}
 
-	/**
+	/*
 	 * Peter can you document this method? TODO
 	 * 
 	 * @param moved
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQCP.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQCP.java
@@ -261,13 +261,13 @@ private void calcRmsd(Point3d[] x, Point3d[] y) {
 			xtrans = CalcPoint.centroid(xref);
 			logger.debug("x centroid: " + xtrans);
 			xtrans.negate();
-			CalcPoint.translate(xtrans, xref);
+			CalcPoint.translate(new Vector3d(xtrans), xref);
 
 			yref = CalcPoint.clonePoint3dArray(y);
 			ytrans = CalcPoint.centroid(yref);
 			logger.debug("y centroid: " + ytrans);
 			ytrans.negate();
-			CalcPoint.translate(ytrans, yref);
+			CalcPoint.translate(new Vector3d(ytrans), yref);
 			innerProduct(yref, xref);
 		}
 		calcRmsd(wsum);
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQuat.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionQuat.java
@@ -59,12 +59,12 @@ public Matrix4d superpose(Point3d[] fixed, Point3d[] moved) {
 		Point3d[] xref = CalcPoint.clonePoint3dArray(fixed);
 		Point3d xtrans = CalcPoint.centroid(xref);
 		xtrans.negate();
-		CalcPoint.translate(xtrans, xref);
+		CalcPoint.translate(new Vector3d(xtrans), xref);
 
 		Point3d[] yref = CalcPoint.clonePoint3dArray(moved);
 		Point3d ytrans = CalcPoint.centroid(yref);
 		ytrans.negate();
-		CalcPoint.translate(ytrans, yref);
+		CalcPoint.translate(new Vector3d(ytrans), yref);
 
 		// calculate rotational component (rotation around origin)
 		Quat4d q = UnitQuaternions.relativeOrientation(xref, yref);
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionSVD.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/geometry/SuperPositionSVD.java
@@ -24,6 +24,7 @@
 
 import javax.vecmath.Matrix4d;
 import javax.vecmath.Point3d;
+import javax.vecmath.Vector3d;
 
 import org.biojava.nbio.structure.jama.Matrix;
 import org.biojava.nbio.structure.jama.SingularValueDecomposition;
@@ -73,10 +74,10 @@ public Matrix4d superpose(Point3d[] fixed, Point3d[] moved) {
 		cenb.negate();
 
 		Point3d[] ats1 = CalcPoint.clonePoint3dArray(fixed);
-		CalcPoint.translate(cena, ats1);
+		CalcPoint.translate(new Vector3d(cena), ats1);
 
 		Point3d[] ats2 = CalcPoint.clonePoint3dArray(moved);
-		CalcPoint.translate(cenb, ats2);
+		CalcPoint.translate(new Vector3d(cenb), ats2);
 
 		double[][] coordSet1 = new double[ats1.length][3];
 		double[][] coordSet2 = new double[ats2.length][3];
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/C2RotationSolver.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/symmetry/core/C2RotationSolver.java
@@ -75,9 +75,9 @@ private void solve() {
 //		SuperPosition.center(y);
 
 		Point3d[] x = CalcPoint.clonePoint3dArray(traces.get(0));
-		CalcPoint.translate(new Point3d(trans), x);
+		CalcPoint.translate(trans, x);
 		Point3d[] y = CalcPoint.clonePoint3dArray(traces.get(1));
-		CalcPoint.translate(new Point3d(trans), y);
+		CalcPoint.translate(trans, y);
 
 		// TODO implement this piece of code using at origin superposition
 		Quat4d quat = UnitQuaternions.relativeOrientation(
diff --git a/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestSuperPosition.java b/biojava-structure/src/test/java/org/biojava/nbio/structure/geometry/TestSuperPosition.java
@@ -25,20 +25,24 @@
  * algorithm for RMSD and Superposition calculations.
  * 
  * @author Aleix Lafita
+ * @author Jose Duarte
  * @since 5.0.0
  *
  */
 public class TestSuperPosition {
 
-	private static final Logger logger = LoggerFactory
-			.getLogger(TestSuperPosition.class);
+	private static final Logger LOGGER = LoggerFactory.getLogger(TestSuperPosition.class);
 
-	private List<Point3d[]> cloud1;
-	private List<Point3d[]> cloud2;
+	private static List<Point3d[]> cloud1;
+	private static List<Point3d[]> cloud2;
 
-	private AxisAngle4d rotAxis;
-	private Vector3d translation;
-	private Matrix4d transform;
+	// the transformation to apply to cloud points 1 that needs to be recovered by the superposition code
+	private static final AxisAngle4d rotAxis = new AxisAngle4d(0.440, 0.302, 0.845, 1.570);
+	private static final Vector3d translation = new Vector3d(0.345, 2.453, 5.324);;
+	private static Matrix4d transform;
+		
+	// a translation to apply to cloud point 2 for the rmsd test only
+	private static final Vector3d translation2 = new Vector3d(1.32, -1.03, 6.28);
 
 	/**
 	 * Generate two clouds of random points of different sizes to test
@@ -54,8 +58,6 @@ public void setUp() throws StructureException {
 
 		Random rnd = new Random(0);
 
-		rotAxis = new AxisAngle4d(0.440, 0.302, 0.845, 1.570);
-		translation = new Vector3d(0.345, 2.453, 5.324);
 		transform = new Matrix4d();
 		transform.set(rotAxis);
 		transform.setTranslation(translation);
@@ -82,10 +84,15 @@ public void setUp() throws StructureException {
 			CalcPoint.center(c1);
 			CalcPoint.center(c2);
 
-			CalcPoint.transform(transform, c1);
+			CalcPoint.transform(transform, c1);			
 
 			cloud1.add(c1);
 			cloud2.add(c2);
+					
+			Point3d centroid1 = CalcPoint. centroid(c1);
+			Point3d centroid2 = CalcPoint. centroid(c2);
+			LOGGER.debug("Centroid c1 (size %d): (%.2f, %.2f, %.2f)\n", size, centroid1.x, centroid1.y, centroid1.z);
+			LOGGER.debug("Centroid c2 (size %d): (%.2f, %.2f, %.2f)\n", size, centroid2.x, centroid2.y, centroid2.z);
 		}
 
 	}
@@ -116,7 +123,7 @@ public void testSuperposition() {
 			Matrix4d qcpTransform = qcp.superpose(cloud1.get(c), cloud2.get(c));
 			long qcpTime = (System.nanoTime() - qcpStart) / 1000;
 
-			logger.error(String.format("Transformation Matrix %d points: "
+			LOGGER.info(String.format("Transformation Matrix %d points: "
 					+ "SVD time %d us, SP time: %d us, QCP time: %d us",
 					cloud1.get(c).length, svdTime, quatTime, qcpTime));
 
@@ -133,6 +140,16 @@ public void testSuperposition() {
 	 */
 	@Test
 	public void testRMSD() {
+		
+		// for the rmsd test we first make sure that both cloud points are not centered in origin so that the centering is tested too
+		// first cloud points are already centered, we translate cloud2 only
+		for (int c=0; c<cloud2.size(); c++) {
+			CalcPoint.translate(translation2, cloud2.get(c));
+			
+			Point3d centroid2 = CalcPoint. centroid(cloud2.get(c));			
+			LOGGER.debug("Centroid c2 (index %d): (%.2f, %.2f, %.2f)\n", c, centroid2.x, centroid2.y, centroid2.z);
+		}
+		
 
 		for (int c = 0; c < cloud1.size(); c++) {
 
@@ -154,11 +171,12 @@ public void testRMSD() {
 			double qcprmsd = qcp.getRmsd(cloud1.get(c), cloud2.get(c));
 			long qcpTime = (System.nanoTime() - qcpStart) / 1000;
 
-			logger.error(String.format("RMSD %d points: SVD time %d us, "
+			LOGGER.info(String.format("RMSD %d points: SVD time %d us, "
 					+ "Quat time: %d us, QCP time: %d us", cloud1.get(c).length,
 					svdTime, quatTime, qcpTime));
 
 			// Check that the returned RMSDs are equal
+			assertEquals(quatrmsd, qcprmsd, 0.001);
 			assertEquals(svdrmsd, quatrmsd, 0.001);
 			assertEquals(svdrmsd, qcprmsd, 0.001);
 		}
