Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_UMEYAMA_H
00026 #define EIGEN_UMEYAMA_H
00027
00028
00029
00030
00031
00032
00033
00034 namespace Eigen {
00035
00036 #ifndef EIGEN_PARSED_BY_DOXYGEN
00037
00038
00039
00040
00041 namespace internal {
00042
00043
00044
00045
00046 template<typename MatrixType, typename OtherMatrixType>
00047 struct umeyama_transform_matrix_type
00048 {
00049 enum {
00050 MinRowsAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(MatrixType::RowsAtCompileTime, OtherMatrixType::RowsAtCompileTime),
00051
00052
00053
00054 HomogeneousDimension = int(MinRowsAtCompileTime) == Dynamic ? Dynamic : int(MinRowsAtCompileTime)+1
00055 };
00056
00057 typedef Matrix<typename traits<MatrixType>::Scalar,
00058 HomogeneousDimension,
00059 HomogeneousDimension,
00060 AutoAlign | (traits<MatrixType>::Flags & RowMajorBit ? RowMajor : ColMajor),
00061 HomogeneousDimension,
00062 HomogeneousDimension
00063 > type;
00064 };
00065
00066 }
00067
00068 #endif
00069
00108 template <typename Derived, typename OtherDerived>
00109 typename internal::umeyama_transform_matrix_type<Derived, OtherDerived>::type
00110 umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, bool with_scaling = true)
00111 {
00112 typedef typename internal::umeyama_transform_matrix_type<Derived, OtherDerived>::type TransformationMatrixType;
00113 typedef typename internal::traits<TransformationMatrixType>::Scalar Scalar;
00114 typedef typename NumTraits<Scalar>::Real RealScalar;
00115 typedef typename Derived::Index Index;
00116
00117 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL)
00118 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename internal::traits<OtherDerived>::Scalar>::value),
00119 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00120
00121 enum { Dimension = EIGEN_SIZE_MIN_PREFER_DYNAMIC(Derived::RowsAtCompileTime, OtherDerived::RowsAtCompileTime) };
00122
00123 typedef Matrix<Scalar, Dimension, 1> VectorType;
00124 typedef Matrix<Scalar, Dimension, Dimension> MatrixType;
00125 typedef typename internal::plain_matrix_type_row_major<Derived>::type RowMajorMatrixType;
00126
00127 const Index m = src.rows();
00128 const Index n = src.cols();
00129
00130
00131 const RealScalar one_over_n = 1 / static_cast<RealScalar>(n);
00132
00133
00134 const VectorType src_mean = src.rowwise().sum() * one_over_n;
00135 const VectorType dst_mean = dst.rowwise().sum() * one_over_n;
00136
00137
00138 const RowMajorMatrixType src_demean = src.colwise() - src_mean;
00139 const RowMajorMatrixType dst_demean = dst.colwise() - dst_mean;
00140
00141
00142 const Scalar src_var = src_demean.rowwise().squaredNorm().sum() * one_over_n;
00143
00144
00145 const MatrixType sigma = one_over_n * dst_demean * src_demean.transpose();
00146
00147 JacobiSVD<MatrixType> svd(sigma, ComputeFullU | ComputeFullV);
00148
00149
00150 TransformationMatrixType Rt = TransformationMatrixType::Identity(m+1,m+1);
00151
00152
00153 VectorType S = VectorType::Ones(m);
00154 if (sigma.determinant()<0) S(m-1) = -1;
00155
00156
00157 const VectorType& d = svd.singularValues();
00158 Index rank = 0; for (Index i=0; i<m; ++i) if (!internal::isMuchSmallerThan(d.coeff(i),d.coeff(0))) ++rank;
00159 if (rank == m-1) {
00160 if ( svd.matrixU().determinant() * svd.matrixV().determinant() > 0 ) {
00161 Rt.block(0,0,m,m).noalias() = svd.matrixU()*svd.matrixV().transpose();
00162 } else {
00163 const Scalar s = S(m-1); S(m-1) = -1;
00164 Rt.block(0,0,m,m).noalias() = svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose();
00165 S(m-1) = s;
00166 }
00167 } else {
00168 Rt.block(0,0,m,m).noalias() = svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose();
00169 }
00170
00171
00172 const Scalar c = 1/src_var * svd.singularValues().dot(S);
00173
00174
00175
00176
00177 Rt.col(m).head(m) = dst_mean;
00178 Rt.col(m).head(m).noalias() -= c*Rt.topLeftCorner(m,m)*src_mean;
00179
00180 if (with_scaling) Rt.block(0,0,m,m) *= c;
00181
00182 return Rt;
00183 }
00184
00185 }
00186
00187 #endif // EIGEN_UMEYAMA_H