|
Chemical Data Processing Library C++ API - Version 1.1.1
|
Go to the documentation of this file.
27 #ifndef CDPL_MATH_KABSCHALGORITHM_HPP
28 #define CDPL_MATH_KABSCHALGORITHM_HPP
85 template <
typename M1,
typename M2,
typename V>
87 bool do_center =
true, std::size_t max_svd_iter = 0)
91 typename V::SizeType>::Type SizeType;
93 SizeType dim = points().getSize1();
94 SizeType num_pts = points().getSize2();
96 CDPL_MATH_CHECK(dim == SizeType(ref_points().getSize1()) && num_pts == SizeType(ref_points().getSize2()),
100 "KabschAlgorithm: Number of points != number of weights",
Base::SizeError);
104 for (SizeType i = 0; i < num_pts; i++) {
106 w_sum += weights()(i);
112 prod(points, weights, centroid1);
113 prod(ref_points, weights, centroid2);
118 tmpPoints.
resize(dim, num_pts,
false);
121 tmpRefPoints.
resize(dim, num_pts,
false);
122 tmpRefPoints.
assign(ref_points);
124 for (SizeType i = 0; i < num_pts; i++) {
125 column(tmpPoints, i).minusAssign(centroid1) *= weights()(i) / w_sum;
126 column(tmpRefPoints, i).minusAssign(centroid2);
130 tmpPoints.
resize(dim, num_pts,
false);
133 for (SizeType i = 0; i < num_pts; i++)
134 column(tmpPoints, i) *= weights()(i) / w_sum;
137 covarMatrix.
resize(dim, dim,
false);
140 prod(tmpPoints,
trans(tmpRefPoints), covarMatrix);
142 prod(tmpPoints,
trans(ref_points), covarMatrix);
144 return align(dim, do_center, max_svd_iter);
160 template <
typename M1,
typename M2>
162 bool do_center =
true, std::size_t max_svd_iter = 0)
167 SizeType dim = points().getSize1();
168 SizeType num_pts = points().getSize2();
170 CDPL_MATH_CHECK(dim == SizeType(ref_points().getSize1()) && num_pts == SizeType(ref_points().getSize2()),
177 tmpPoints.
resize(dim, num_pts,
false);
180 tmpRefPoints.
resize(dim, num_pts,
false);
181 tmpRefPoints.
assign(ref_points);
183 for (SizeType i = 0; i < num_pts; i++) {
184 column(tmpPoints, i).minusAssign(centroid1);
185 column(tmpRefPoints, i).minusAssign(centroid2);
189 covarMatrix.
resize(dim, dim,
false);
192 prod(tmpPoints,
trans(tmpRefPoints), covarMatrix);
194 prod(points,
trans(ref_points), covarMatrix);
196 return align(dim, do_center, max_svd_iter);
205 template <
typename SizeType>
206 bool align(SizeType dim,
bool do_center, std::size_t max_svd_iter)
209 svdV.
resize(dim, dim,
false);
211 if (!
svDecompose(covarMatrix, svdW, svdV, max_svd_iter))
217 SizeType xform_dim = dim + 1;
219 transform.
resize(xform_dim, xform_dim,
false);
221 range(transform, 0, dim, 0, dim).assign(
prod(svdV,
trans(covarMatrix)));
229 range(last_col, 0, dim).assign(centroid2 -
prod(
range(transform, 0, dim, 0, dim), centroid1));
250 #endif // CDPL_MATH_KABSCHALGORITHM_HPP
T ValueType
Definition: KabschAlgorithm.hpp:65
std::common_type< T1, T2 >::type Type
Definition: CommonType.hpp:43
MatrixTranspose< E > trans(MatrixExpression< E > &e)
Definition: MatrixExpression.hpp:941
Definition: Vector.hpp:1470
MatrixColumn< M > column(MatrixExpression< M > &e, typename MatrixColumn< M >::SizeType j)
Definition: MatrixProxy.hpp:730
Matrix< T > MatrixType
Definition: KabschAlgorithm.hpp:66
Definition: Expression.hpp:54
Matrix & assign(const MatrixExpression< E > &e)
Definition: Matrix.hpp:459
Definition: Expression.hpp:76
#define CDPL_MATH_CHECK(expr, msg, e)
Definition: Check.hpp:36
Definition: CommonType.hpp:41
Definition of vector proxy types.
E::ValueType det(const MatrixExpression< E > &e)
Definition: Matrix.hpp:1721
Thrown to indicate that the size of a (multidimensional) array is not correct.
Definition: Base/Exceptions.hpp:133
Definition: MatrixProxy.hpp:49
Definition: MatrixProxy.hpp:196
Definition of type traits.
Matrix1VectorBinaryTraits< E1, E2, MatrixVectorProduct< E1, E2 > >::ResultType prod(const MatrixExpression< E1 > &e1, const VectorExpression< E2 > &e2)
Definition: MatrixExpression.hpp:833
void resize(SizeType m, SizeType n, bool preserve=true, const ValueType &v=ValueType())
Definition: Matrix.hpp:519
Definition: Vector.hpp:1292
Thrown to indicate errors caused by some invalid value.
Definition: Base/Exceptions.hpp:76
Implementation of matrix singular value decomposition and associated operations.
bool align(const MatrixExpression< M1 > &points, const MatrixExpression< M2 > &ref_points, bool do_center=true, std::size_t max_svd_iter=0)
Computes the rigid body transformation that aligns a set of -dimensional points points with a corres...
Definition: KabschAlgorithm.hpp:161
const unsigned int T
Specifies Hydrogen (Tritium).
Definition: AtomType.hpp:67
bool align(const MatrixExpression< M1 > &points, const MatrixExpression< M2 > &ref_points, const VectorExpression< V > &weights, bool do_center=true, std::size_t max_svd_iter=0)
Computes the rigid body transformation that aligns a set of -dimensional points points with a corres...
Definition: KabschAlgorithm.hpp:86
Implementation of the Kabsch algorithm [KABA].
Definition: KabschAlgorithm.hpp:62
Definition of exception classes.
Vector< T > VectorType
Definition: KabschAlgorithm.hpp:67
The namespace of the Chemical Data Processing Library.
void resize(SizeType n, const ValueType &v=ValueType())
Definition: Vector.hpp:491
Definition of matrix proxy types.
Definition of matrix data types.
const MatrixType & getTransform() const
Definition: KabschAlgorithm.hpp:199
Definition of various preprocessor macros for error checking.
MatrixRange< E > range(MatrixExpression< E > &e, const typename MatrixRange< E >::RangeType &r1, const typename MatrixRange< E >::RangeType &r2)
Definition: MatrixProxy.hpp:744
bool svDecompose(MatrixExpression< A > &a, VectorExpression< W > &w, MatrixExpression< V > &v, std::size_t max_iter=0)
Computes the Singular Value Decomposition [WSVD] of a -dimensional matrix a.
Definition: SVDecomposition.hpp:70
Definition of vector data types.