commit 3daf92c7a5e8288d47839e47a461b8d249d206f1
author Antonio Sanchez <cantonios@google.com> Mon Jan 11 09:27:25 2021 -0800
committer Antonio Sanchez <cantonios@google.com> Mon Jan 11 10:13:38 2021 -0800
tree 9e4d495f68eac302ef26c5781d6a3cb9aa829d0c
parent 587fd6ab707b0743f48593e401eb546d709eca4d

Transform::computeScalingRotation flush determinant to +/- 1.

In the previous code, in attempting to correct for a negative
determinant, we end up multiplying and dividing by a number that
is often very near, but not exactly +/-1.  By flushing to +/-1,
we can replace a division with a multiplication, and results
are more numerically consistent.

Eigen/src/Geometry/Transform.h

diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h
index 547b366..7497b31 100644
--- a/Eigen/src/Geometry/Transform.h
+++ b/Eigen/src/Geometry/Transform.h
@@ -1097,16 +1097,17 @@
 template<typename RotationMatrixType, typename ScalingMatrixType>
 EIGEN_DEVICE_FUNC void Transform<Scalar,Dim,Mode,Options>::computeRotationScaling(RotationMatrixType *rotation, ScalingMatrixType *scaling) const
 {
+  // TODO: investigate BDCSVD implementation.
   JacobiSVD<LinearMatrixType> svd(linear(), ComputeFullU | ComputeFullV);
 
-  Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant(); // so x has absolute value 1
+  Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant() < Scalar(0) ? Scalar(-1) : Scalar(1); // so x has absolute value 1
   VectorType sv(svd.singularValues());
   sv.coeffRef(Dim-1) *= x;
   if(scaling) *scaling = svd.matrixV() * sv.asDiagonal() * svd.matrixV().adjoint();
   if(rotation)
   {
     LinearMatrixType m(svd.matrixU());
-    m.col(Dim-1) /= x;
+    m.col(Dim-1) *= x;
     *rotation = m * svd.matrixV().adjoint();
   }
 }
@@ -1126,16 +1127,17 @@
 template<typename ScalingMatrixType, typename RotationMatrixType>
 EIGEN_DEVICE_FUNC void Transform<Scalar,Dim,Mode,Options>::computeScalingRotation(ScalingMatrixType *scaling, RotationMatrixType *rotation) const
 {
+  // TODO: investigate BDCSVD implementation.
   JacobiSVD<LinearMatrixType> svd(linear(), ComputeFullU | ComputeFullV);
 
-  Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant(); // so x has absolute value 1
+  Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant() < Scalar(0) ? Scalar(-1) : Scalar(1); // so x has absolute value 1
   VectorType sv(svd.singularValues());
   sv.coeffRef(Dim-1) *= x;
   if(scaling) *scaling = svd.matrixU() * sv.asDiagonal() * svd.matrixU().adjoint();
   if(rotation)
   {
     LinearMatrixType m(svd.matrixU());
-    m.col(Dim-1) /= x;
+    m.col(Dim-1) *= x;
     *rotation = m * svd.matrixV().adjoint();
   }
 }