From 8eb0fc1e72687ad32e1bbc8a635f3834e3592ed3 Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Tue, 12 Oct 2010 10:19:59 -0400
Subject: remove SVD class (was bad code taked from elsewhere) Use JacobiSVD
 for now. We do plan to reintroduce a bidiagonalizing SVD asap.

---
 Eigen/src/Geometry/Transform.h | 26 ++++++++++++++++++++++++--
 Eigen/src/Geometry/Umeyama.h   |  2 +-
 2 files changed, 25 insertions(+), 3 deletions(-)

(limited to 'Eigen/src/Geometry')

diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h
index 58bb2d0c0..e8099495d 100644
--- a/Eigen/src/Geometry/Transform.h
+++ b/Eigen/src/Geometry/Transform.h
@@ -928,7 +928,18 @@ template<typename Scalar, int Dim, int Mode>
 template<typename RotationMatrixType, typename ScalingMatrixType>
 void Transform<Scalar,Dim,Mode>::computeRotationScaling(RotationMatrixType *rotation, ScalingMatrixType *scaling) const
 {
-  linear().svd().computeRotationScaling(rotation, scaling);
+  JacobiSVD<LinearMatrixType> svd(linear(), ComputeFullU | ComputeFullV);
+
+  Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant(); // so x has absolute value 1
+  VectorType sv(svd.singularValues());
+  sv.coeffRef(0) *= x;
+  if(scaling) scaling->lazyAssign(svd.matrixV() * sv.asDiagonal() * svd.matrixV().adjoint());
+  if(rotation)
+  {
+    LinearMatrixType m(svd.matrixU());
+    m.col(0) /= x;
+    rotation->lazyAssign(m * svd.matrixV().adjoint());
+  }
 }
 
 /** decomposes the linear part of the transformation as a product rotation x scaling, the scaling being
@@ -946,7 +957,18 @@ template<typename Scalar, int Dim, int Mode>
 template<typename ScalingMatrixType, typename RotationMatrixType>
 void Transform<Scalar,Dim,Mode>::computeScalingRotation(ScalingMatrixType *scaling, RotationMatrixType *rotation) const
 {
-  linear().svd().computeScalingRotation(scaling, rotation);
+  JacobiSVD<LinearMatrixType> svd(linear(), ComputeFullU | ComputeFullV);
+
+  Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant(); // so x has absolute value 1
+  VectorType sv(svd.singularValues());
+  sv.coeffRef(0) *= x;
+  if(scaling) scaling->lazyAssign(svd.matrixU() * sv.asDiagonal() * svd.matrixU().adjoint());
+  if(rotation)
+  {
+    LinearMatrixType m(svd.matrixU());
+    m.col(0) /= x;
+    rotation->lazyAssign(m * svd.matrixV().adjoint());
+  }
 }
 
 /** Convenient method to set \c *this from a position, orientation and scale
diff --git a/Eigen/src/Geometry/Umeyama.h b/Eigen/src/Geometry/Umeyama.h
index fa6bfcbee..04449f623 100644
--- a/Eigen/src/Geometry/Umeyama.h
+++ b/Eigen/src/Geometry/Umeyama.h
@@ -141,7 +141,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo
   // Eq. (38)
   const MatrixType sigma = one_over_n * dst_demean * src_demean.transpose();
 
-  SVD<MatrixType> svd(sigma);
+  JacobiSVD<MatrixType> svd(sigma, ComputeFullU | ComputeFullV);
 
   // Initialize the resulting transformation with an identity matrix...
   TransformationMatrixType Rt = TransformationMatrixType::Identity(m+1,m+1);
-- 
cgit v1.2.3