Added Umeyama implementation.

Eigen/Geometry

@@ -44,6 +44,7 @@ namespace Eigen {
 #include "src/Geometry/Hyperplane.h"
 #include "src/Geometry/ParametrizedLine.h"
 #include "src/Geometry/AlignedBox.h"
+#include "src/Geometry/Umeyama.h"
 
 #if defined EIGEN_VECTORIZE_SSE
   #include "src/Geometry/arch/Geometry_SSE.h"

Eigen/src/Geometry/Umeyama.h

@@ -0,0 +1,205 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Hauke Heibel <hauke.heibel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_UMEYAMA_H
+#define EIGEN_UMEYAMA_H
+
+// This file requires the user to include 
+// * Eigen/Core
+// * Eigen/LU 
+// * Eigen/SVD
+// * Eigen/Array
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+
+// These helpers are required since it allows to use mixed types as parameters
+// for the Umeyama. The problem with mixed parameters is that the return type
+// cannot trivially be deduced when float and double types are mixed.
+namespace
+{
+  // Compile time return type deduction for different MatrixBase types.
+  // Different means here different alignment and parameters but the same underlying
+  // real scalar type.
+  template<typename MatrixType, typename OtherMatrixType>
+  struct ei_umeyama_transform_matrix_type
+  {    
+    enum {
+      MinRowsAtCompileTime = EIGEN_ENUM_MIN(MatrixType::RowsAtCompileTime, OtherMatrixType::RowsAtCompileTime),
+      MinMaxRowsAtCompileTime = EIGEN_ENUM_MIN(MatrixType::MaxRowsAtCompileTime, OtherMatrixType::MaxRowsAtCompileTime),
+
+      // When possible we want to choose some small fixed size value since the result
+      // is likely to fit on the stack.
+      HomogeneousDimension = EIGEN_ENUM_MIN(MinRowsAtCompileTime+1, Dynamic),
+      MaxRowsAtCompileTime = EIGEN_ENUM_MIN(MinMaxRowsAtCompileTime+1, Dynamic),
+      MaxColsAtCompileTime = EIGEN_ENUM_MIN(MatrixType::MaxColsAtCompileTime, OtherMatrixType::MaxColsAtCompileTime)
+    };
+
+    typedef Matrix<typename ei_traits<MatrixType>::Scalar,
+      HomogeneousDimension,
+      HomogeneousDimension,
+      AutoAlign | (ei_traits<MatrixType>::Flags & RowMajorBit ? RowMajor : ColMajor),
+      MaxRowsAtCompileTime, 
+      MaxColsAtCompileTime
+    > type;
+  };
+}
+
+#endif
+
+/**
+* \geometry_module \ingroup Geometry_Module
+*
+* \brief Returns the transformation between two point sets.
+*
+* The algorithm is based on:
+* "Least-squares estimation of transformation parameters between two point patterns",
+* Shinji Umeyama, PAMI 1991, DOI: 10.1109/34.88573
+*
+* It estimates parameters \f$ c, \mathbf{R}, \f$ and \f$ \mathbf{t} \f$ such that
+* \f{align*}
+*   \frac{1}{n} \sum_{i=1}^n \vert\vert y_i - (c\mathbf{R}x_i + \mathbf{t}) \vert\vert_2^2
+* \f}
+* is minimized.
+*
+* The algorithm is based on the analysis of the covariance matrix
+* \f$ \Sigma_{\mathbf{x}\mathbf{y}} \in \mathbb{R}^{d \times d} \f$
+* of the input point sets \f$ \mathbf{x} \f$ and \f$ \mathbf{y} \f$ where 
+* \f$d\f$ is corresponding to the dimension (which is typically small).
+* The analysis is involving the SVD having a complexity of \f$O(d^3)\f$
+* though the actual bottleneck usually lies in the computation of the covariance
+* matrix which has an asymptotic lower bound of \f$O(dm)\f$ when the input point 
+* sets have dimension \f$d \times m\f$.
+*
+* Currently the method is working only for floating point matrices.
+*
+* \todo Should the return type of umeyama() become a Transform?
+*
+* \param src Source points \f$ \mathbf{x} = \left( x_1, \hdots, x_n \right) \f$.
+* \param dst Destination points \f$ \mathbf{y} = \left( y_1, \hdots, y_n \right) \f$.
+* \param with_scaling Sets \f$ c=1 \f$ when <code>false</code> is passed.
+* \return The homogeneous transformation 
+* \f{align*}
+*   T = \begin{bmatrix} c\mathbf{R} & \mathbf{t} \\ \mathbf{0} & 1 \end{bmatrix}
+* \f}
+* minimizing the resudiual above. This transformation is always returned as an 
+* Eigen::Matrix.
+*/
+template <typename Derived, typename OtherDerived>
+typename ei_umeyama_transform_matrix_type<Derived, OtherDerived>::type
+umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, bool with_scaling = true)
+{
+  typedef typename ei_umeyama_transform_matrix_type<Derived, OtherDerived>::type TransformationMatrixType;
+  typedef typename ei_traits<TransformationMatrixType>::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+
+  EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
+  EIGEN_STATIC_ASSERT((ei_is_same_type<Scalar, typename ei_traits<OtherDerived>::Scalar>::ret),
+    YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
+
+  enum { Dimension = EIGEN_ENUM_MIN(Derived::RowsAtCompileTime, OtherDerived::RowsAtCompileTime) };
+
+  typedef Matrix<Scalar, Dimension, 1> VectorType;
+  typedef Matrix<Scalar, Dimension, Dimension> MatrixType;
+
+  const int m = src.rows(); // dimension
+  const int n = src.cols(); // number of measurements
+
+  // required for demeaning ...
+  const RealScalar one_over_n = 1 / static_cast<RealScalar>(n);
+
+  // computation of mean
+  const VectorType src_mean = src.rowwise().sum() * one_over_n;
+  const VectorType dst_mean = dst.rowwise().sum() * one_over_n;
+
+  // demeaning of src and dst points
+  MatrixType src_demean(m,n);
+  MatrixType dst_demean(m,n);
+  for (int i=0; i<n; ++i)
+  {
+    src_demean.col(i) = src.col(i) - src_mean;
+    dst_demean.col(i) = dst.col(i) - dst_mean;
+  }
+
+  // Eq. (36)-(37)
+  const Scalar src_var = src_demean.rowwise().squaredNorm().sum() * one_over_n;
+  const Scalar dst_var = dst_demean.rowwise().squaredNorm().sum() * one_over_n;
+
+  // Eq. (38)
+  const MatrixType sigma = (dst_demean*src_demean.transpose()).lazy() * one_over_n;
+
+  SVD<MatrixType> svd(sigma);
+
+  // Initialize the resulting transformation with an identity matrix...
+  TransformationMatrixType Rt = TransformationMatrixType::Identity(m+1,m+1);
+
+  // Eq. (39)
+  VectorType S = VectorType::Ones(m);
+  if (sigma.determinant()<0) S(m-1) = -1;
+
+  // Eq. (40) and (43)
+  const VectorType& d = svd.singularValues();
+  int rank = 0; for (int i=0; i<m; ++i) if (!ei_isMuchSmallerThan(d.coeff(i),d.coeff(0))) ++rank;
+  if (rank == m-1) {
+    if ( svd.matrixU().determinant() * svd.matrixV().determinant() > 0 ) {
+      Rt.block(0,0,m,m) = (svd.matrixU()*svd.matrixV().transpose()).lazy();
+    } else {
+      const Scalar s = S(m-1); S(m-1) = -1;
+      Rt.block(0,0,m,m) = (svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose()).lazy();
+      S(m-1) = s;
+    }
+  } else {
+    Rt.block(0,0,m,m) = (svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose()).lazy();
+  }
+
+  // Eq. (42)
+  const Scalar c = 1/src_var * svd.singularValues().dot(S);
+
+  // Eq. (41)
+  // TODO: lazyness does not make much sense over here, right?
+  Rt.col(m).segment(0,m) = dst_mean - c*Rt.block(0,0,m,m)*src_mean;
+
+  if (with_scaling) Rt.block(0,0,m,m) *= c;
+
+  return Rt;
+}
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+
+/**
+* This is simply here to prevent the creation of dozens compile time errors for
+* std::complex types...
+*/
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols,
+         typename _OtherScalar, int _OtherRows, int _OtherCols, int _OtherOptions, int _OtherMaxRows, int _OtherMaxCols>
+         typename ei_umeyama_transform_matrix_type<Matrix<std::complex<_Scalar>,_Rows,_Cols,_Options,_MaxRows,_MaxCols>, 
+                                                   Matrix<std::complex<_OtherScalar>,_OtherRows,_OtherCols,_OtherOptions,_OtherMaxRows,_OtherMaxCols> >::type
+umeyama(const MatrixBase<Matrix<std::complex<_Scalar>,_Rows,_Cols,_Options,_MaxRows,_MaxCols> >& src, 
+        const MatrixBase<Matrix<std::complex<_OtherScalar>,_OtherRows,_OtherCols,_OtherOptions,_OtherMaxRows,_OtherMaxCols> >& dst, bool with_scaling = true)
+{
+  EIGEN_STATIC_ASSERT(false, NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
+}
+
+#endif
+
+#endif // EIGEN_UMEYAMA_H

doc/Doxyfile.in

@@ -962,7 +962,8 @@ PAPER_TYPE             = a4wide
 # The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX
 # packages that should be included in the LaTeX output.
 
-EXTRA_PACKAGES         = amssymb
+EXTRA_PACKAGES         = amssymb \
+                         amsmath
 
 # The LATEX_HEADER tag can be used to specify a personal LaTeX header for
 # the generated latex document. The header should contain everything until

test/CMakeLists.txt

@@ -139,6 +139,7 @@ ei_add_test(sparse_vector)
 ei_add_test(sparse_basic)
 ei_add_test(sparse_product)
 ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
+ei_add_test(umeyama)
 
 
 ei_add_property(EIGEN_TESTING_SUMMARY "CXX:               ${CMAKE_CXX_COMPILER}\n")

test/umeyama.cpp

@@ -0,0 +1,203 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2009 Hauke Heibel <hauke.heibel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or1 FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+#include <Eigen/Core>
+#include <Eigen/Array>
+#include <Eigen/Geometry>
+
+#include <Eigen/LU> // required for MatrixBase::determinant
+#include <Eigen/SVD> // required for SVD
+
+using namespace Eigen;
+
+#define VAR_CALL_SUBTEST(...) do { \
+  g_test_stack.push_back(EI_PP_MAKE_STRING(__VA_ARGS__)); \
+  __VA_ARGS__; \
+  g_test_stack.pop_back(); \
+} while (0)
+
+//  Constructs a random matrix from the unitary group U(size).
+template <typename T>
+Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> randMatrixUnitary(int size)
+{
+  typedef typename T Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> MatrixType;
+
+  MatrixType Q;
+
+  int max_tries = 40;
+  double is_unitary = false;
+
+  while (!is_unitary && max_tries > 0)
+  {
+	  // initialize random matrix
+    Q = MatrixType::Random(size, size);
+
+	  // orthogonalize columns using the Gram-Schmidt algorithm
+	  for (int col = 0; col < size; ++col)
+	  {
+      MatrixType::ColXpr colVec = Q.col(col);
+		  for (int prevCol = 0; prevCol < col; ++prevCol)
+		  {
+			  MatrixType::ColXpr prevColVec = Q.col(prevCol);
+        colVec -= colVec.dot(prevColVec)*prevColVec;
+		  }
+		  Q.col(col) = colVec.normalized();
+	  }
+
+	  // this additional orthogonalization is not necessary in theory but should enhance 
+    // the numerical orthogonality of the matrix
+	  for (int row = 0; row < size; ++row)
+	  {
+      MatrixType::RowXpr rowVec = Q.row(row);
+		  for (int prevRow = 0; prevRow < row; ++prevRow)
+		  {
+			  MatrixType::RowXpr prevRowVec = Q.row(prevRow);
+			  rowVec -= rowVec.dot(prevRowVec)*prevRowVec;
+		  }
+		  Q.row(row) = rowVec.normalized();
+	  }
+
+	  // final check
+    is_unitary = Q.isUnitary();
+    --max_tries;
+  }
+
+  if (max_tries == 0)
+    throw std::runtime_error("randMatrixUnitary: Could not construct unitary matrix!");
+
+	return Q;
+}
+
+//  Constructs a random matrix from the special unitary group SU(size).
+template <typename T>
+Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> randMatrixSpecialUnitary(int size)
+{
+  typedef typename T Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> MatrixType;
+
+	// initialize unitary matrix
+	MatrixType Q = randMatrixUnitary<Scalar>(size);
+
+	// tweak the first column to make the determinant be 1
+  Q.col(0) *= ei_conj(Q.determinant());
+
+  return Q;
+}
+
+template <typename MatrixType>
+void run_test(int dim, int num_elements)
+{
+  typedef typename ei_traits<MatrixType>::Scalar Scalar;
+  typedef Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> MatrixX;
+  typedef Matrix<Scalar, Eigen::Dynamic, 1> VectorX;
+
+  // MUST be positive because in any other case det(cR_t) may become negative for
+  // odd dimensions!
+  const Scalar c = ei_abs(ei_random<Scalar>());
+
+  MatrixX R = randMatrixSpecialUnitary<Scalar>(dim);
+  VectorX t = Scalar(50)*VectorX::Random(dim,1);
+
+  MatrixX cR_t = MatrixX::Identity(dim+1,dim+1);
+  cR_t.block(0,0,dim,dim) = c*R;
+  cR_t.block(0,dim,dim,1) = t;
+
+  MatrixX src = MatrixX::Random(dim+1, num_elements);
+  src.row(dim) = Matrix<Scalar, 1, Dynamic>::Constant(num_elements, Scalar(1));
+
+  MatrixX dst = (cR_t*src).lazy();
+
+  MatrixX cR_t_umeyama = umeyama(src.block(0,0,dim,num_elements), dst.block(0,0,dim,num_elements));
+
+  const Scalar error = ( cR_t_umeyama*src - dst ).cwise().square().sum();
+
+  VERIFY(error < Scalar(10)*std::numeric_limits<Scalar>::epsilon());
+}
+
+template<typename Scalar, int Dimension>
+void run_fixed_size_test(int num_elements)
+{
+  typedef Matrix<Scalar, Dimension+1, Dynamic> MatrixX;
+  typedef Matrix<Scalar, Dimension+1, Dimension+1> HomMatrix;
+  typedef Matrix<Scalar, Dimension, Dimension> FixedMatrix;
+  typedef Matrix<Scalar, Dimension, 1> FixedVector;
+
+  const int dim = Dimension;
+
+  // MUST be positive because in any other case det(cR_t) may become negative for
+  // odd dimensions!
+  const Scalar c = ei_abs(ei_random<Scalar>());
+
+  FixedMatrix R = randMatrixSpecialUnitary<Scalar>(dim);
+  FixedVector t = Scalar(50)*FixedVector::Random(dim,1);
+
+  HomMatrix cR_t = HomMatrix::Identity(dim+1,dim+1);
+  cR_t.block(0,0,dim,dim) = c*R;
+  cR_t.block(0,dim,dim,1) = t;
+
+  MatrixX src = MatrixX::Random(dim+1, num_elements);
+  src.row(dim) = Matrix<Scalar, 1, Dynamic>::Constant(num_elements, Scalar(1));
+
+  MatrixX dst = (cR_t*src).lazy();
+
+  HomMatrix cR_t_umeyama = umeyama(src.block(0,0,dim,num_elements), dst.block(0,0,dim,num_elements));
+
+  const Scalar error = ( cR_t_umeyama*src - dst ).cwise().square().sum();
+
+  VERIFY(error < Scalar(10)*std::numeric_limits<Scalar>::epsilon());
+}
+
+void test_umeyama()
+{
+  for (int i=0; i<g_repeat; ++i)
+  {
+    const int num_elements = ei_random<int>(40,500);
+
+    // works also for dimensions bigger than 3...
+    for (int dim=2; dim<8; ++dim)
+    {
+      CALL_SUBTEST(run_test<MatrixXd>(dim, num_elements));
+      CALL_SUBTEST(run_test<MatrixXf>(dim, num_elements));
+    }
+
+    VAR_CALL_SUBTEST(run_fixed_size_test<float, 2>(num_elements));
+    VAR_CALL_SUBTEST(run_fixed_size_test<float, 3>(num_elements));
+    VAR_CALL_SUBTEST(run_fixed_size_test<float, 4>(num_elements));
+
+    VAR_CALL_SUBTEST(run_fixed_size_test<double, 2>(num_elements));
+    VAR_CALL_SUBTEST(run_fixed_size_test<double, 3>(num_elements));
+    VAR_CALL_SUBTEST(run_fixed_size_test<double, 4>(num_elements));
+  }
+
+  // Those two calls don't compile and result in meaningful error messages!
+  // umeyama(MatrixXcf(),MatrixXcf());
+  // umeyama(MatrixXcd(),MatrixXcd());
+}