GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-29 15:24:51 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source 
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.
This commit is contained in:
Antonio Sanchez 2021-01-11 09:27:25 -08:00
parent 587fd6ab70
commit 3daf92c7a5
1 changed files with 6 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Eigen/src/Geometry/Transform.h
View File

@ -1097,16 +1097,17 @@ template<typename Scalar, int Dim, int Mode, int Options>
template<typename RotationMatrixType, typename ScalingMatrixType> template<typename RotationMatrixType, typename ScalingMatrixType>
EIGEN_DEVICE_FUNC void Transform<Scalar,Dim,Mode,Options>::computeRotationScaling(RotationMatrixType *rotation, ScalingMatrixType *scaling) const 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); 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()); VectorType sv(svd.singularValues());
sv.coeffRef(Dim-1) *= x; sv.coeffRef(Dim-1) *= x;
if(scaling) *scaling = svd.matrixV() * sv.asDiagonal() * svd.matrixV().adjoint(); if(scaling) *scaling = svd.matrixV() * sv.asDiagonal() * svd.matrixV().adjoint();
if(rotation) if(rotation)
{ {
LinearMatrixType m(svd.matrixU()); LinearMatrixType m(svd.matrixU());
m.col(Dim-1) /= x; m.col(Dim-1) *= x;
*rotation = m * svd.matrixV().adjoint(); *rotation = m * svd.matrixV().adjoint();
} }
} }
@ -1126,16 +1127,17 @@ template<typename Scalar, int Dim, int Mode, int Options>
template<typename ScalingMatrixType, typename RotationMatrixType> template<typename ScalingMatrixType, typename RotationMatrixType>
EIGEN_DEVICE_FUNC void Transform<Scalar,Dim,Mode,Options>::computeScalingRotation(ScalingMatrixType *scaling, RotationMatrixType *rotation) const 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); 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()); VectorType sv(svd.singularValues());
sv.coeffRef(Dim-1) *= x; sv.coeffRef(Dim-1) *= x;
if(scaling) *scaling = svd.matrixU() * sv.asDiagonal() * svd.matrixU().adjoint(); if(scaling) *scaling = svd.matrixU() * sv.asDiagonal() * svd.matrixU().adjoint();
if(rotation) if(rotation)
{ {
LinearMatrixType m(svd.matrixU()); LinearMatrixType m(svd.matrixU());
m.col(Dim-1) /= x; m.col(Dim-1) *= x;
*rotation = m * svd.matrixV().adjoint(); *rotation = m * svd.matrixV().adjoint();
} }
} }