polar_dec produces rotations

Former-commit-id: 89cb1f36b59cfc9d24c5c945676a80f585e638bf

include/igl/polar_dec.cpp

@@ -32,34 +32,46 @@ IGL_INLINE void igl::polar_dec(
   Eigen::PlainObjectBase<DerivedV> & V)
 {
   using namespace std;
- typedef typename DerivedA::Scalar Scalar;
+  using namespace Eigen;
+  typedef typename DerivedA::Scalar Scalar;
 
- const Scalar th = std::sqrt(Eigen::NumTraits<Scalar>::dummy_precision());
+  const Scalar th = std::sqrt(Eigen::NumTraits<Scalar>::dummy_precision());
 
- Eigen::SelfAdjointEigenSolver<DerivedA> eig;
- feclearexcept(FE_UNDERFLOW);
- eig.computeDirect(A.transpose()*A);
- if(fetestexcept(FE_UNDERFLOW) || eig.eigenvalues()(0)/eig.eigenvalues()(2)<th)
- {
-   cout<<"resorting to svd 1..."<<endl;
-   return polar_svd(A,R,T,U,S,V);
- }
+  Eigen::SelfAdjointEigenSolver<DerivedA> eig;
+  feclearexcept(FE_UNDERFLOW);
+  eig.computeDirect(A.transpose()*A);
+  if(fetestexcept(FE_UNDERFLOW) || eig.eigenvalues()(0)/eig.eigenvalues()(2)<th)
+  {
+    cout<<"resorting to svd 1..."<<endl;
+    return polar_svd(A,R,T,U,S,V);
+  }
 
- S = eig.eigenvalues().cwiseSqrt();
+  S = eig.eigenvalues().cwiseSqrt();
 
- T = eig.eigenvectors() * S.asDiagonal() * eig.eigenvectors().transpose();
- U = A * eig.eigenvectors();
- V = eig.eigenvectors();
- R = U * S.asDiagonal().inverse() * V.transpose();
- S = S.reverse().eval();
- V = V.rowwise().reverse().eval();
- U = U.rowwise().reverse().eval() * S.asDiagonal().inverse();
+  V = eig.eigenvectors();
+  U = A * V;
+  R = U * S.asDiagonal().inverse() * V.transpose();
+  T = V * S.asDiagonal() * V.transpose();
 
- if(std::fabs(R.squaredNorm()-3.) > th)
- {
-   cout<<"resorting to svd 2..."<<endl;
-   return polar_svd(A,R,T,U,S,V);
- }
+  S = S.reverse().eval();
+  V = V.rowwise().reverse().eval();
+  U = U.rowwise().reverse().eval() * S.asDiagonal().inverse();
+
+  if(R.determinant() < 0)
+  {
+    // Annoyingly the .eval() is necessary
+    auto W = V.eval();
+    const auto & SVT = S.asDiagonal() * V.adjoint();
+    W.col(V.cols()-1) *= -1.;
+    R = U*W.transpose();
+    T = W*SVT;
+  }
+
+  if(std::fabs(R.squaredNorm()-3.) > th)
+  {
+    cout<<"resorting to svd 2..."<<endl;
+    return polar_svd(A,R,T,U,S,V);
+  }
 }
 
 template <

include/igl/polar_svd.cpp

@@ -8,6 +8,7 @@
 #include "polar_svd.h"
 #include <Eigen/SVD>
 #include <Eigen/Geometry>
+#include <iostream>
 
 // Adapted from Olga's CGAL mentee's ARAP code
 template <
@@ -46,16 +47,20 @@ IGL_INLINE void igl::polar_svd(
   U = svd.matrixU();
   V = svd.matrixV();
   S = svd.singularValues();
-  // TODO: Check for sign and switch here.
   R = U*V.transpose();
-  const auto & SV = S.asDiagonal() * V.adjoint();
+  const auto & SVT = S.asDiagonal() * V.adjoint();
   // Check for reflection
   if(R.determinant() < 0)
   {
-    V.col(V.cols()-1) *= -1.;
-    R = U*V.transpose();
+    // Annoyingly the .eval() is necessary
+    auto W = V.eval();
+    W.col(V.cols()-1) *= -1.;
+    R = U*W.transpose();
+    T = W*SVT;
+  }else
+  {
+    T = V*SVT;
   }
-  T = V*SV;
 }
 
 #ifndef IGL_HEADER_ONLY