Functions
absor Namespace Reference

Functions

int hornAbsOrientation (const Eigen::MatrixXd &refPoints, const Eigen::MatrixXd &targetPoints, std::vector< double > *quat, double *rmsd, std::vector< double > *rmsdList, double *scale)
 Get the absolute orientation using Horn's algorithm (with quaternions) More...
 
Eigen::MatrixXd calcMatrixS (const Eigen::MatrixXd &centeredRefPnts, const Eigen::MatrixXd &centeredTargetPnts, int nop, int dim)
 
Eigen::MatrixXd calcMatrixN (const Eigen::MatrixXd &S)
 
Eigen::MatrixXd centerWRTcentroid (const Eigen::MatrixXd &pointSet)
 Center a point set wrt the centroid. More...
 
double calcScaleFactor (const Eigen::MatrixXd &rightSys, const Eigen::MatrixXd &leftSys, int n)
 Calculate the scale factor from the centered left and right point sets. More...
 
Eigen::MatrixXd quat2RotMatrix (const Eigen::VectorXd &quat)
 Get a rotation matrix from a unit quaternion. More...
 
double getRMSD (const Eigen::MatrixXd &centeredRefPnts, const Eigen::MatrixXd &centeredTargetPnts, const Eigen::VectorXd &quat, std::vector< double > *rmsdList, int nop, double scale)
 Calculate the RMSD. More...
 

Function Documentation

◆ calcMatrixN()

Eigen::MatrixXd absor::calcMatrixN ( const Eigen::MatrixXd &  S)

Compute the matrix N, a 4x4 symmetric matrix, by combining sums (saved as elements in the matrix S)

Compute the matrix \((4 \times 4)\) N, whose largest eigenvector corresponds to the optimal rotation. It is calculated from the elements of the \((3 \times 3)\) matrix S. The matrix N is:

N=[(Sxx+Syy+Szz) (Syz-Szy) (Szx-Sxz) (Sxy-Syx);...

(Syz-Szy) (Sxx-Syy-Szz) (Sxy+Syx) (Szx+Sxz);...

(Szx-Sxz) (Sxy+Syx) (-Sxx+Syy-Szz) (Syz+Szy);...

(Sxy-Syx) (Szx+Sxz) (Syz+Szy) (-Sxx-Syy+Szz)];

Parameters
[in]S(3 \times 3) Eigen matrix whose elements are the sums of products of the coordinates measured in the left and right systems.
Returns
a \( (4 \times 4) \) Eigen matrix, or N.

Definition at line 172 of file absOrientation.cpp.

172  {
173  //
174  Eigen::MatrixXd N(4, 4); // Output matrix N
175  // Components of S
176  double Sxx = S(0, 0);
177  double Sxy = S(0, 1);
178  double Sxz = S(0, 2);
179  double Syx = S(1, 0);
180  double Syy = S(1, 1);
181  double Syz = S(1, 2);
182  double Szx = S(2, 0);
183  double Szy = S(2, 1);
184  double Szz = S(2, 2);
185 
186  // N=[(Sxx+Syy+Szz) (Syz-Szy) (Szx-Sxz) (Sxy-Syx);...
187  // (Syz-Szy) (Sxx-Syy-Szz) (Sxy+Syx) (Szx+Sxz);...
188  // (Szx-Sxz) (Sxy+Syx) (-Sxx+Syy-Szz) (Syz+Szy);...
189  // (Sxy-Syx) (Szx+Sxz) (Syz+Szy) (-Sxx-Syy+Szz)];
190 
191  // ------------------
192  // Calculating N:
193  // Diagonal elements
194  N(0, 0) = Sxx + Syy + Szz;
195  N(1, 1) = Sxx - Syy - Szz;
196  N(2, 2) = -Sxx + Syy - Szz;
197  N(3, 3) = -Sxx - Syy + Szz;
198  // Other elements
199  // First row
200  N(0, 1) = Syz - Szy;
201  N(0, 2) = Szx - Sxz;
202  N(0, 3) = Sxy - Syx;
203  // Second row
204  N(1, 0) = Syz - Szy;
205  N(1, 2) = Sxy + Syx;
206  N(1, 3) = Szx + Sxz;
207  // Third row
208  N(2, 0) = Szx - Sxz;
209  N(2, 1) = Sxy + Syx;
210  N(2, 3) = Syz + Szy;
211  // Fourth row
212  N(3, 0) = Sxy - Syx;
213  N(3, 1) = Szx + Sxz;
214  N(3, 2) = Syz + Szy;
215  // ------------------
216  // Output matrix
217  return N;
218 } // end of function

◆ calcMatrixS()

Eigen::MatrixXd absor::calcMatrixS ( const Eigen::MatrixXd &  centeredRefPnts,
const Eigen::MatrixXd &  centeredTargetPnts,
int  nop,
int  dim 
)

Compute the matrix S, or M, whose elements are the sums of products of coordinates measured in the left and right systems

Compute the matrix S, or M, whose elements are the sums of products of coordinates measured in the left and right systems. The matrix S is : S=[Sxx Sxy Sxz;... Syx Syy Syz;... Szx Szy Szz];

Parameters
[in]centeredRefPntsThe point set of the reference system (or right system), centered with respect to the centroid. This is a \( (n \times 3) \) Eigen matrix. Here, \( n \) is the number of particles.
[in]centeredTargetPnts\( (n \times 3) \) Eigen matrix of the candidate/test system (or left system), centered with respect to the centroid.
[in]nopThe number of particles.
[in]dimThe number of dimensions = 3.
Returns
a \( (3 \times 3) \) Eigen matrix, or S.

Definition at line 131 of file absOrientation.cpp.

133  {
134  Eigen::MatrixXd S(nop, dim); // Output matrix S
135  Eigen::VectorXd targetCoord(nop); // Column of the target point set
136  Eigen::VectorXd refCoord(nop); // Column of the reference point set
137  double Svalue; // Current value being filled (Sxx, Sxy etc.)
138 
139  // Calculate Sxx, Sxy, Sxz etc
140  for (int iCol = 0; iCol < dim; iCol++) {
141  //
142  for (int jCol = 0; jCol < dim; jCol++) {
143  targetCoord =
144  centeredTargetPnts.col(iCol); // iCol^th column of target point set
145  refCoord = centeredRefPnts.col(jCol); // jCol^th of reference point set
146  Svalue = targetCoord.dot(refCoord);
147  S(iCol, jCol) = Svalue;
148  } // end column wise filling
149  } // end of filling
150 
151  // Output matrix
152  return S;
153 } // end of function

◆ calcScaleFactor()

double absor::calcScaleFactor ( const Eigen::MatrixXd &  rightSys,
const Eigen::MatrixXd &  leftSys,
int  n 
)

Calculate the scale factor from the centered left and right point sets.

Calculate the scale factor from the centered right and left point sets.

Parameters
[in]centeredRefPntsThe point set of the reference system (or right system), centered with respect to the centroid. This is a \( (n \times 3) \) Eigen matrix. Here, \( n \) is the number of particles.
[in]centeredTargetPnts\( (n \times 3) \) Eigen matrix of the candidate/test system (or left system), centered with respect to the centroid.
[in]nThe number of points.
Returns
the scale factor.

Definition at line 276 of file absOrientation.cpp.

277  {
278  double scale; // Output scale
279  double v1, v2; // Sum of the length of the vector
280  Eigen::VectorXd rightVec(
281  3); // Vector of the i^th particle in the right system
282  Eigen::VectorXd leftVec(3); // Vector of the i^th particle in the right system
283 
284  // scale = (sigma_to_n ||r_r||^2 / ||r_l||^2)^0.5
285  // ref: http://people.csail.mit.edu/bkph/papers/Absolute_Orientation.pdf
286 
287  // Loop through all the points, and get the dot product of the vector for each
288  // point
289  for (int i = 0; i < n; i++) {
290  //
291  rightVec = rightSys.row(i); // i^th row of the right system
292  leftVec = leftSys.row(i); // i^th row of the left system
293  v1 += rightVec.dot(rightVec);
294  v2 += leftVec.dot(leftVec);
295  } // end of loop through all points
296 
297  // The optimum scale is the ratio of v1 and v2
298  scale = std::sqrt(v1 / v2);
299 
300  return scale;
301 } // end of function
std::sqrt
T sqrt(T... args)

◆ centerWRTcentroid()

Eigen::MatrixXd absor::centerWRTcentroid ( const Eigen::MatrixXd &  pointSet)

Center a point set wrt the centroid.

Centers a point set (which is an Eigen matrix), with respect to the centroid.

Parameters
[in]pointSet\( (n \times 3) \) Eigen matrix for the point set. Here \( n \) is the number of points.
Returns
a \( (n \times 3) \) Eigen matrix of the same size as the input point set.

Definition at line 228 of file absOrientation.cpp.

228  {
229  int nop = pointSet.rows(); // Number of particles
230  int dim = pointSet.cols(); // Number of dimensions
231  Eigen::MatrixXd centeredPointSet(nop, dim); // Output point set
232  Eigen::VectorXd vecOfOnes(nop); // vector of ones
233  std::vector<double> centroid;
234  double coordValue;
235  double centeredVal;
236  //
237  centroid.resize(dim); // Init to zero
238  vecOfOnes = Eigen::VectorXd::Ones(nop);
239  // --------------------------------
240 
241  for (int i = 0; i < nop; i++) {
242  for (int k = 0; k < dim; k++) {
243  coordValue = pointSet(i, k);
244  centroid[k] += coordValue;
245  } // loop through columns
246  } // end of loop through rows
247  // Divide by the total number of particles
248  centroid[0] /= nop; // x
249  centroid[1] /= nop; // y
250  centroid[2] /= nop; // z
251  // --------------------------------
252  // Subtract the centroid from the coordinates to get the centered point set
253  for (int i = 0; i < nop; i++) {
254  for (int k = 0; k < dim; k++) {
255  coordValue = pointSet(i, k);
256  centeredVal = coordValue - centroid[k];
257  centeredPointSet(i, k) = centeredVal;
258  } // end of loop through columns (dimensions)
259  } // end of loop through the rows
260  // --------------------------------
261  return centeredPointSet;
262 } // end of function
std::vector::resize
T resize(T... args)

◆ getRMSD()

double absor::getRMSD ( const Eigen::MatrixXd &  centeredRefPnts,
const Eigen::MatrixXd &  centeredTargetPnts,
const Eigen::VectorXd &  quat,
std::vector< double > *  rmsdList,
int  nop,
double  scale 
)

Calculate the RMSD.

Get the root mean square of the errors from Horn's absolute orientation algorithm.

Parameters
[in]quatThe Eigen vector of length 4, for the input quaternion.
[in]centeredRefPntsThe point set of the reference system (or right system), centered with respect to the centroid. This is a \( (n \times 3) \) Eigen matrix. Here, \( n \) is the number of particles.
[in]centeredTargetPnts\( (n \times 3) \) Eigen matrix of the candidate/test system (or left system), centered with respect to the centroid.
[in]quatThe Eigen vector of length 4, for the quaternion denoting the rotation.
[in]nopThe number of particles.
[in]scaleThe scale factor
Returns
the least RMSD.

Definition at line 359 of file absOrientation.cpp.

362  {
363  Eigen::MatrixXd R(3, 3); // The (3x3) rotation vector
364  // The RMSD per atom is filled in this vector
365  (*rmsdList).resize(nop);
366  //
367  R = absor::quat2RotMatrix(
368  quat); // orthonormal rotation matrix from the quaternion
369  // The total error is:
370  // sum_over_all_n (r'_r - s*rotated_r_l')
371  double rmsd = 0.0; // Error
372  Eigen::VectorXd errorVec(3); // The vector which is the r_r -s*R
373  Eigen::VectorXd rotatedLeft(3);
374  Eigen::VectorXd targetCol(3);
375  Eigen::VectorXd refCol(3);
376  //
377  for (int i = 0; i < nop; i++) {
378  //
379  // Rotate the left system coordinate using the rotation matrix
380  targetCol = (centeredTargetPnts.row(i)).transpose();
381  refCol = (centeredRefPnts.row(i)).transpose();
382  //
383  rotatedLeft = R * targetCol;
384  errorVec = refCol - scale * rotatedLeft;
385  rmsd += errorVec.dot(errorVec); // Total error
386  (*rmsdList)[i] = sqrt(errorVec.dot(errorVec)); // Error per atom
387  } // end of loop through every row
388  //
389  return sqrt(rmsd / nop);
390 } // end of function
absor::quat2RotMatrix
Eigen::MatrixXd quat2RotMatrix(const Eigen::VectorXd &quat)
Get a rotation matrix from a unit quaternion.
Definition: absOrientation.cpp:309

◆ hornAbsOrientation()

int absor::hornAbsOrientation ( const Eigen::MatrixXd &  refPoints,
const Eigen::MatrixXd &  targetPoints,
std::vector< double > *  quat,
double *  rmsd,
std::vector< double > *  rmsdList,
double *  scale 
)

Get the absolute orientation using Horn's algorithm (with quaternions)

Get the absolute orientation using Horn's algorithm (with quaternions). The algorithm is described in the linked paper.

Parameters
[in]refPointsThe point set of the reference system (or right system). This is a \( (n \times 3) \) Eigen matrix. Here, \( n \) is the number of particles.
[in]targetPoints\( (n \times 3) \) Eigen matrix of the candidate/test system (or left system).
[in,out]quatThe quaternion for the optimum rotation of the left system into the right system.
[in,out]scaleThe scale factor obtained from Horn's algorithm.

Definition at line 30 of file absOrientation.cpp.

33  {
34  int nop =
35  refPoints.rows(); // Number of particles (equal to the number of rows)
36  int dim =
37  refPoints.cols(); // Number of dimensions (equal to the number of columns)
38  Eigen::MatrixXd centeredRefPnts(
39  nop, dim); // Reference point set after centering wrt the centroid
40  Eigen::MatrixXd centeredTargetPnts(
41  nop, dim); // Target point set after centering wrt the centroid
42  Eigen::MatrixXd S(dim,
43  dim); // Matrix containing sums of products of coordinates
44  Eigen::MatrixXd N(
45  4, 4); // 4x4 Matrix, whose largest eigenvector must be calculated
46  Eigen::VectorXd calcEigenVec(
47  4); // This should have 4 components (eigen vector calculated from N)
48  // -----
49  // Check that the sizes of the reference point set (right point system) and
50  // the target point set (left point system) are the same
51  if (refPoints.rows() != targetPoints.rows() ||
52  refPoints.cols() != targetPoints.cols()) {
53  // Throw error
54  std::cerr
55  << "The reference and target point sets are not of the same size.\n";
56  return 1;
57  } // unequal size; error!
58  // -----
59  //
60  // ---------------------------------------------------
61  // FINDING THE CENTROIDS AND THE NEW COORDINATES WRT THE CENTROIDS
62  centeredRefPnts = absor::centerWRTcentroid(refPoints);
63  centeredTargetPnts = absor::centerWRTcentroid(targetPoints);
64  // ---------------------------------------------------
65  // FINDING THE ROTATION MATRIX
66  //
67  // Find the 3x3 matrix S
68  S = absor::calcMatrixS(centeredRefPnts, centeredTargetPnts, nop, dim);
69  // Calculate the 4x4 symmetric matrix N, whose largest eigenvector yields the
70  // quaternion in the same direction
71  N = absor::calcMatrixN(S);
72  // --------
73  // Calculate the eigenvector corresponding to the largest eigenvalue
74  //
75  // Construct matrix operation object (op) using the wrapper class
76  // DenseSymMatProd for the matrix N
77  Spectra::DenseSymMatProd<double> op(N);
78  //
79  // Construct eigen solver object, requesting the largest 1 eigenvalue and
80  // eigenvector
81  Spectra::SymEigsSolver<double, Spectra::LARGEST_ALGE,
82  Spectra::DenseSymMatProd<double>>
83  eigs(&op, 1, 4);
84  //
85  // Initialize and compute
86  eigs.init();
87  int nconv = eigs.compute();
88  // Get the eigenvalue and eigenvector
89  if (eigs.info() == Spectra::SUCCESSFUL) {
90  Eigen::VectorXd calcEigenValue = eigs.eigenvalues(); // Eigenvalue
91  calcEigenVec = eigs.eigenvectors();
92  } // end of eigenvector calculation
93  //
94  // --------
95  // Normalize the eigenvector calculated
96  double qNorm = sqrt(calcEigenVec.dot(calcEigenVec));
97  calcEigenVec /= qNorm; // Divide by the square root of the sum
98  // Update the quaternion with the normalized eigenvector
99  (*quat).resize(4); // Output quaternion update
100  for (int i = 0; i < 4; i++) {
101  (*quat)[i] = calcEigenVec(i);
102  } // end of quaternion update
103  // --------
104  // ---------------------------------------------------
105  // COMPUTE THE OPTIMUM SCALE
106  (*scale) = absor::calcScaleFactor(centeredRefPnts, centeredTargetPnts, nop);
107  // ---------------------------------------------------
108  // GETTING THE ERROR
109  (*rmsd) = absor::getRMSD(centeredRefPnts, centeredTargetPnts, calcEigenVec,
110  rmsdList, nop, (*scale));
111  // ---------------------------------------------------
112  return 0;
113 } // end of function
absor::getRMSD
double getRMSD(const Eigen::MatrixXd &centeredRefPnts, const Eigen::MatrixXd &centeredTargetPnts, const Eigen::VectorXd &quat, std::vector< double > *rmsdList, int nop, double scale)
Calculate the RMSD.
Definition: absOrientation.cpp:359
absor::calcMatrixS
Eigen::MatrixXd calcMatrixS(const Eigen::MatrixXd &centeredRefPnts, const Eigen::MatrixXd &centeredTargetPnts, int nop, int dim)
Definition: absOrientation.cpp:131
absor::calcScaleFactor
double calcScaleFactor(const Eigen::MatrixXd &rightSys, const Eigen::MatrixXd &leftSys, int n)
Calculate the scale factor from the centered left and right point sets.
Definition: absOrientation.cpp:276
absor::centerWRTcentroid
Eigen::MatrixXd centerWRTcentroid(const Eigen::MatrixXd &pointSet)
Center a point set wrt the centroid.
Definition: absOrientation.cpp:228
absor::calcMatrixN
Eigen::MatrixXd calcMatrixN(const Eigen::MatrixXd &S)
Definition: absOrientation.cpp:172

◆ quat2RotMatrix()

Eigen::MatrixXd absor::quat2RotMatrix ( const Eigen::VectorXd &  quat)

Get a rotation matrix from a unit quaternion.

Get a \( (3 \times 3) \) rotation matrix from a unit quaternion.

Parameters
[in]quatThe Eigen vector of length 4, for the input quaternion.
Returns
the rotation matrix.

Definition at line 309 of file absOrientation.cpp.

309  {
310  //
311  Eigen::MatrixXd R(3, 3); // Rotation matrix
312  // Components of the quaternion
313  double q0 = quat(0);
314  double qx = quat(1);
315  double qy = quat(2);
316  double qz = quat(3);
317 
318  // Quaternion derived rotation matrix, when q (q=q0+qx*i+qy*j+qz*k)
319  // is a unit quaternion:
320  // R=[(q0^2+qx^2+qy^2+qz^2) 2(qx*qy-q0*qz) 2(qx*qz+q0*qy);...
321  // 2(qy*qx-q0*qz) (q0^2+qx^2+qy^2+qz^2) 2(qy*qz-q0*qx);...
322  // 2(qz*qx-q0*qy) 2(qz*qy+q0*qz) (q0^2+qx^2+qy^2+qz^2)];
323 
324  // Fill up the rotation matrix R according to the above formula
325  //
326  // First row
327  R(0, 0) = q0 * q0 + qx * qx + qy * qy + qz * qz;
328  R(0, 1) = 2 * (qx * qy - q0 * qz);
329  R(0, 2) = 2 * (qx * qz + q0 * qy);
330  // Second row
331  R(1, 0) = 2 * (qy * qx + q0 * qz);
332  R(1, 1) = q0 * q0 + qx * qx + qy * qy + qz * qz;
333  R(1, 2) = 2 * (qy * qz - q0 * qy);
334  // Third row
335  R(2, 0) = 2 * (qz * qx - q0 * qy);
336  R(2, 1) = 2 * (qz * qy + q0 * qz);
337  R(2, 2) = q0 * q0 + qx * qx + qy * qy + qz * qz;
338 
339  // return the rotation matrix
340  return R;
341 } // end of function
Generated by Doxygen 1.9.1 and the doxyYoda theme