diff options
author | Antonio Sanchez <cantonios@google.com> | 2021-01-11 09:27:25 -0800 |
---|---|---|
committer | Antonio Sanchez <cantonios@google.com> | 2021-01-11 10:13:38 -0800 |
commit | 3daf92c7a5e8288d47839e47a461b8d249d206f1 (patch) | |
tree | 9e4d495f68eac302ef26c5781d6a3cb9aa829d0c /Eigen/src | |
parent | 587fd6ab707b0743f48593e401eb546d709eca4d (diff) |
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.
Diffstat (limited to 'Eigen/src')
-rw-r--r-- | Eigen/src/Geometry/Transform.h | 10 |
1 files changed, 6 insertions, 4 deletions
diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 547b3661f..7497b31ac 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -1097,16 +1097,17 @@ template<typename Scalar, int Dim, int Mode, int Options> 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 Scalar, int Dim, int Mode, int Options> 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(); } } |