Update class version with newest Code

Newest version of ASOPA code. Changed from row to column vectors.

Common/Transforms/vtkAnisotropicLandmarkTransform.cxx

@@ -114,8 +114,8 @@ void vtkAnisotropicLandmarkTransform::InternalUpdate()
   }
 
   //put the points into eigen matrices
-  Eigen::MatrixXd X(N_PTS, 3);
-  Eigen::MatrixXd Y(N_PTS, 3);
+  Eigen::MatrixXd X(3, N_PTS);
+  Eigen::MatrixXd Y(3, N_PTS);
   double p1[3];
   double p2[3];
   for(int i=0; i < N_PTS; i++)
@@ -124,41 +124,42 @@ void vtkAnisotropicLandmarkTransform::InternalUpdate()
     this->TargetLandmarks->GetPoint(i, p2);
     for(int j=0; j < 3; j++)
     {
-      X(i, j) = p1[j];
-      Y(i, j) = p2[j];
+      X(j, i) = p1[j];
+      Y(j, i) = p2[j];
     }
   }
 
   //find centroids
-  Eigen::Vector3d X_centroid = X.colwise().mean();
-  Eigen::Vector3d Y_centroid = Y.colwise().mean();
+  Eigen::Vector3d X_centroid = X.rowwise().mean();
+  Eigen::Vector3d Y_centroid = Y.rowwise().mean();
 
   //translate input by centroids
-  Eigen::MatrixXd X_trans = X.rowwise() - X_centroid.transpose();
-  Eigen::MatrixXd Y_trans = Y.rowwise() - Y_centroid.transpose();
+  Eigen::MatrixXd X_trans = X.colwise() - X_centroid;
+  Eigen::MatrixXd Y_trans = Y.colwise() - Y_centroid;
 
   //normalise translated source input
-  Eigen::MatrixXd X_norm = X_trans.colwise().normalized();
+  Eigen::MatrixXd X_norm = X_trans.rowwise().normalized();
 
   //Find the cross covariance matrix
-  Eigen::Matrix3d B = Y_trans.transpose() * X_norm;
+  Eigen::Matrix3d B = Y_trans * X_norm.transpose();
 
   //use SVD to decompose B such that B=USV^T
   Eigen::JacobiSVD<Eigen::Matrix3d> svd(B, Eigen::ComputeFullU | Eigen::ComputeFullV);
   Eigen::Matrix3d U = svd.matrixU();
   Eigen::Matrix3d V = svd.matrixV().transpose();
 
-  //create D where the last element is the determinant of UV^T
+  //create D where the last element is the determinant of UV
   Eigen::Matrix3d D(3, 3);
   D = Eigen::Matrix3d::Identity();
   D(2, 2) = (U*V).determinant();
 
   //Find the rotation
-  Eigen::Matrix3d Q = U * D * V;
+  Eigen::Matrix3d Q(3, 3);
+  Q = U * D * V;
 
   //Calculate the FRE for the given rotation
-  Eigen::MatrixXd FRE_vect(N_PTS, 3);
-  FRE_vect = X_norm * Q.transpose() - Y_trans;
+  Eigen::MatrixXd FRE_vect(3, N_PTS);
+  FRE_vect = Y_trans - Q.transpose() * X_norm;
   double FRE = sqrt(FRE_vect.squaredNorm()/N_PTS);
 
   //starting FRE value to compare to
@@ -181,29 +182,30 @@ void vtkAnisotropicLandmarkTransform::InternalUpdate()
     Q = U * D * V;
 
     //recompute FRE value
-    FRE_vect = X_norm * Q.transpose() - Y_trans;
+    FRE_vect = Y_trans - Q.transpose() * X_norm;
     FRE_orig = FRE;
     FRE = 0.0f;
     FRE = sqrt(FRE_vect.squaredNorm()/N_PTS);
   }
 
   //calculate final scaling
-  B = Y_trans.transpose() * X_trans;
+  B = Y_trans * X_trans.transpose();
   U = B.transpose() * Q;
-  V = X_trans.transpose() * X_trans;
-  Eigen::Matrix3d A = Eigen::Matrix3d::Zero();
+  V = X_trans * X_trans.transpose();
+  Eigen::Matrix3d A(3, 3);
+  A = Eigen::Matrix3d::Zero();
   A(0, 0) = U(0, 0)/V(0, 0);
   A(1, 1) = U(1, 1)/V(1, 1);
   A(2, 2) = U(2, 2)/V(2, 2);
 
+  //calculate final translation
+  Eigen::Vector3d t(3);
+  t =  Y_centroid - (Q * (A * X_centroid));
+
   //calculate final FRE values
-  FRE_vect = Y_trans - (X_trans * A * Q.transpose());
-  Eigen::MatrixXd FRE_mag = FRE_vect.rowwise().squaredNorm();
+  FRE_vect = Y - ((Q * (A * X)).colwise() + t);
   FRE = sqrt(FRE_vect.squaredNorm()/N_PTS);
 
-  //calculate final translation
-  Eigen::Vector3d t =  Y_centroid - (Q * A * X_centroid);
-
   //set scaling and rotation matrices sperate to regular
   //transformation matrix output
   this->RotationMatrix->Identity();