arap works in 2D (bug fix in fit_rotations_planar)

Former-commit-id: d63950c2d212e604cd80567a58e1ce84afe44438

examples/arap/Makefile

@@ -13,6 +13,7 @@ LIBIGL_INC=-I$(LIBIGL)/include
 LIBIGL_LIB=-L$(LIBIGL)/lib -ligl -liglpng -liglsvd3x3
 
 EIGEN3_INC=-I/opt/local/include/eigen3 -I/opt/local/include/eigen3/unsupported
+CFLAGS+=-DIGL_HEADER_ONLY
 
 
 #EMBREE=$(LIBIGL)/external/embree
@@ -37,7 +38,11 @@ ifeq ($(UNAME), Darwin)
 	ANTTWEAKBAR_LIB+=-framework AppKit
 endif
 
-INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(YIMG_INC) 
+# SVD 
+SINGULAR_VALUE_DECOMPOSITION_INC=\
+	-I$(LIBIGL)/external/Singular_Value_Decomposition/
+
+INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(YIMG_INC) $(SINGULAR_VALUE_DECOMPOSITION_INC)
 LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(YIMG_LIB) 
 
 example: example.o

include/igl/matlab_format.cpp

@@ -110,4 +110,12 @@ template Eigen::WithFormat<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const igl::m
 template Eigen::WithFormat<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const igl::matlab_format<Eigen::Matrix<double, 2, 1, 0, 2, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 template Eigen::WithFormat<Eigen::Matrix<int, 1, -1, 1, 1, -1> > const igl::matlab_format<Eigen::Matrix<int, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, 1, -1, 1, 1, -1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 template Eigen::WithFormat<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const igl::matlab_format<Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<double, 3, 3, 0, 3, 3> > const igl::matlab_format<Eigen::Matrix<double, 3, 3, 0, 3, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 3, 0, 3, 3> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<double, 2, 2, 0, 2, 2> > const igl::matlab_format<Eigen::Matrix<double, 2, 2, 0, 2, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 2, 0, 2, 2> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<float, 2, 2, 0, 2, 2> > const igl::matlab_format<Eigen::Matrix<float, 2, 2, 0, 2, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<float, 2, 2, 0, 2, 2> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<double, 2, 2, 1, 2, 2> > const igl::matlab_format<Eigen::Matrix<double, 2, 2, 1, 2, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 2, 1, 2, 2> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const igl::matlab_format<Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<float, 2, 1, 0, 2, 1> > const igl::matlab_format<Eigen::Matrix<float, 2, 1, 0, 2, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, 2, 1, 0, 2, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<float, 2, 2, 1, 2, 2> > const igl::matlab_format<Eigen::Matrix<float, 2, 2, 1, 2, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<float, 2, 2, 1, 2, 2> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const igl::matlab_format<Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 #endif

include/igl/min_quad_with_fixed.h

@@ -13,7 +13,6 @@
 #include <Eigen/Core>
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
-#include <Eigen/SparseExtra>
 // Bug in unsupported/Eigen/SparseExtra needs iostream first
 #include <iostream>
 #include <unsupported/Eigen/SparseExtra>

include/igl/polar_svd.cpp

@@ -39,6 +39,7 @@ IGL_INLINE void igl::polar_svd(
   Eigen::PlainObjectBase<DerivedS> & S,
   Eigen::PlainObjectBase<DerivedV> & V)
 {
+  using namespace std;
   Eigen::JacobiSVD<DerivedA> svd;
   svd.compute(A, Eigen::ComputeFullU | Eigen::ComputeFullV );
   U = svd.matrixU();

include/igl/svd3x3/Makefile

@@ -20,7 +20,6 @@ INC+=-I../../../include/
 
 # EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
 
-
 # SVD 
 SINGULAR_VALUE_DECOMPOSITION_INC=\
 	-I../../../external/Singular_Value_Decomposition/

include/igl/svd3x3/arap.cpp

@@ -21,6 +21,7 @@
 #include <cassert>
 #include <iostream>
 
+
 template <
   typename DerivedV,
   typename DerivedF,
@@ -32,6 +33,7 @@ IGL_INLINE bool igl::arap_precomputation(
   ARAPData & data)
 {
   using namespace igl;
+  using namespace std;
   using namespace Eigen;
   typedef typename DerivedV::Scalar Scalar;
   // number of vertices
@@ -44,7 +46,7 @@ IGL_INLINE bool igl::arap_precomputation(
   //assert(F.cols() == 3 && "For now only triangles");
   // dimension
   const int dim = V.cols();
-  assert(dim == 3 && "Only 3d supported");
+  //assert(dim == 3 && "Only 3d supported");
   // Defaults
   data.f_ext = Eigen::MatrixXd::Zero(n,dim);
 
@@ -154,7 +156,10 @@ IGL_INLINE bool igl::arap_solve(
   if(U.size() == 0)
   {
     // terrible initial guess.. should at least copy input mesh
+#ifndef NDEBUG
+    cerr<<"arap_solve: Using terrible initial guess for U. Try U = V."<<endl;
     U = MatrixXd::Zero(data.n,dim);
+#endif
   }
   // changes each arap iteration
   MatrixXd U_prev = U;
@@ -173,13 +178,21 @@ IGL_INLINE bool igl::arap_solve(
       U.row(data.b(bi)) = bc.row(bi);
     }
 
-    MatrixXd S = data.CSM * U.replicate(dim,1);
+    const auto & Udim = U.replicate(dim,1);
+    MatrixXd S = data.CSM * Udim;
+
     MatrixXd R(dim,data.CSM.rows());
+    if(dim == 2)
+    {
+      fit_rotations_planar(S,R);
+    }else
+    {
 #ifdef __SSE__ // fit_rotations_SSE will convert to float if necessary
-    fit_rotations_SSE(S,R);
+      fit_rotations_SSE(S,R);
 #else
-    fit_rotations(S,R);
+      fit_rotations(S,R);
 #endif
+    }
     //for(int k = 0;k<(data.CSM.rows()/dim);k++)
     //{
     //  R.block(0,dim*k,dim,dim) = MatrixXd::Identity(dim,dim);
@@ -207,7 +220,7 @@ IGL_INLINE bool igl::arap_solve(
     MatrixXd Dl;
     if(data.with_dynamics)
     {
-      assert(M.rows() == n && 
+      assert(data.M.rows() == n && 
         "No mass matrix. Call arap_precomputation if changing with_dynamics");
       const double h = data.h;
       assert(h != 0);

include/igl/svd3x3/arap.h

@@ -75,7 +75,10 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedF> & F,
     const Eigen::PlainObjectBase<Derivedb> & b,
     ARAPData & data);
-
+  // Inputs:
+  //   bc  #b by dim list of boundary conditions
+  //   data  struct containing necessary precomputation and parameters
+  //   U  #V by dim initial guess
   template <
     typename Derivedbc,
     typename DerivedU>

include/igl/svd3x3/fit_rotations.cpp

@@ -11,6 +11,7 @@
 #include <igl/verbose.h>
 #include <igl/polar_dec.h>
 #include <igl/polar_svd.h>
+#include <igl/matlab_format.h>
 #include <iostream>
 
 template <typename DerivedS, typename DerivedD>
@@ -42,21 +43,8 @@ IGL_INLINE void igl::fit_rotations(
     }
     Eigen::Matrix<typename DerivedD::Scalar,3,3> ri;
     Eigen::Matrix<typename DerivedD::Scalar,3,3> ti;
-    //polar_dec(si,ri,ti);
-    //polar_svd(si,ri,ti);
     polar_svd3x3(si, ri);
     assert(ri.determinant() >= 0);
-#ifndef FIT_ROTATIONS_ALLOW_FLIPS
-    // Check for reflection
-    if(ri.determinant() < 0)
-    {
-      cerr<<"Error: Warning: flipping is wrong..."<<endl;
-      assert(false && "This is wrong. Need to flip column in U and recompute R = U*V'");
-      // flip sign of last row
-      ri.row(2) *= -1;
-    }
-    assert(ri.determinant() >= 0);
-#endif  
     // Not sure why polar_dec computes transpose...
     R.block(0,r*dim,dim,dim) = ri.block(0,0,dim,dim).transpose();
   }
@@ -70,6 +58,7 @@ IGL_INLINE void igl::fit_rotations_planar(
   using namespace std;
   const int dim = S.cols();
   const int nr = S.rows()/dim;
+  assert(dim == 2 && "_planar input should be 2D");
   assert(nr * dim == S.rows());
 
   // resize output
@@ -87,28 +76,24 @@ IGL_INLINE void igl::fit_rotations_planar(
         si(i,j) = S(i*nr+r,j);
       }
     }
-    Eigen::Matrix<typename DerivedD::Scalar,2,2> ri;
-    Eigen::Matrix<typename DerivedD::Scalar,2,2> ti;
-    igl::polar_svd(si,ri,ti);
+    typedef Eigen::Matrix<typename DerivedD::Scalar,2,2> Mat2;
+    typedef Eigen::Matrix<typename DerivedD::Scalar,2,1> Vec2;
+    Mat2 ri,ti,ui,vi;
+    Vec2 _;
+    igl::polar_svd(si,ri,ti,ui,_,vi);
+
 #ifndef FIT_ROTATIONS_ALLOW_FLIPS
     // Check for reflection
     if(ri.determinant() < 0)
     {
-      cerr<<"Error: Warning: flipping is wrong..."<<endl;
-      assert(false && "This is wrong. Need to flip column in U and recompute R = U*V'");
-      // flip sign of last row
-      ri.row(1) *= -1;
+      vi.col(1) *= -1.;
+      ri = ui * vi.transpose();
     }
     assert(ri.determinant() >= 0);
 #endif  
+
     // Not sure why polar_dec computes transpose...
-    R.block(0,r*dim,2,2) = ri.block(0,0,2,2).transpose();
-    // Set remaining part to identity
-    R(0,r*3+2) = 0;
-    R(1,r*3+2) = 0;
-    R(2,r*3+0) = 0;
-    R(2,r*3+1) = 0;
-    R(2,r*3+2) = 1;
+    R.block(0,r*dim,2,2) = ri.transpose();
   }
 }
 

include/igl/svd3x3/polar_svd3x3.h

@@ -22,7 +22,8 @@ namespace igl
   //
   //  This means that det(R) = 1. Technically it's not polar decomposition
   //  which guarantees positive semidefinite stretch factor (at the cost of
-  //  having det(R) = -1).
+  //  having det(R) = -1). "• The orthogonal factors U and V will be true
+  //  rotation matrices..." [McAdams, Selle, Tamstorf, Teran, Sefakis 2011]
   //
   template<typename Mat>
   IGL_INLINE void polar_svd3x3(const Mat& A, Mat& R);

include/igl/svd3x3/svd3x3.cpp

@@ -23,7 +23,7 @@ IGL_INLINE void igl::svd3x3(const Eigen::Matrix<T, 3, 3>& A, Eigen::Matrix<T, 3,
 {
   // this code only supports the scalar version (otherwise we'd need to pass arrays of matrices)  
 
-#include "Singular_Value_Decomposition_Kernel_Declarations.hpp"
+#include <Singular_Value_Decomposition_Kernel_Declarations.hpp>
 
   ENABLE_SCALAR_IMPLEMENTATION(Sa11.f=A(0,0);)                                      ENABLE_SSE_IMPLEMENTATION(Va11=_mm_loadu_ps(a11);)                                  ENABLE_AVX_IMPLEMENTATION(Va11=_mm256_loadu_ps(a11);)
     ENABLE_SCALAR_IMPLEMENTATION(Sa21.f=A(1,0);)                                      ENABLE_SSE_IMPLEMENTATION(Va21=_mm_loadu_ps(a21);)                                  ENABLE_AVX_IMPLEMENTATION(Va21=_mm256_loadu_ps(a21);)
@@ -35,7 +35,7 @@ IGL_INLINE void igl::svd3x3(const Eigen::Matrix<T, 3, 3>& A, Eigen::Matrix<T, 3,
     ENABLE_SCALAR_IMPLEMENTATION(Sa23.f=A(1,2);)                                      ENABLE_SSE_IMPLEMENTATION(Va23=_mm_loadu_ps(a23);)                                  ENABLE_AVX_IMPLEMENTATION(Va23=_mm256_loadu_ps(a23);)
     ENABLE_SCALAR_IMPLEMENTATION(Sa33.f=A(2,2);)                                      ENABLE_SSE_IMPLEMENTATION(Va33=_mm_loadu_ps(a33);)                                  ENABLE_AVX_IMPLEMENTATION(Va33=_mm256_loadu_ps(a33);)
 
-#include "Singular_Value_Decomposition_Main_Kernel_Body.hpp"
+#include <Singular_Value_Decomposition_Main_Kernel_Body.hpp>
 
     ENABLE_SCALAR_IMPLEMENTATION(U(0,0)=Su11.f;)                                      ENABLE_SSE_IMPLEMENTATION(_mm_storeu_ps(u11,Vu11);)                                 ENABLE_AVX_IMPLEMENTATION(_mm256_storeu_ps(u11,Vu11);)
     ENABLE_SCALAR_IMPLEMENTATION(U(1,0)=Su21.f;)                                      ENABLE_SSE_IMPLEMENTATION(_mm_storeu_ps(u21,Vu21);)                                 ENABLE_AVX_IMPLEMENTATION(_mm256_storeu_ps(u21,Vu21);)