arap parameterization no longer allowing reflections, using ref triangles

Former-commit-id: 9492e2df91da2ee1d178fcdca3c632c004d90299

include/igl/cotmatrix.h

@@ -9,7 +9,6 @@
 #define IGL_COTMATRIX_H
 #include "igl_inline.h"
 
-#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 

include/igl/embree/project_mesh.cpp

@@ -0,0 +1,56 @@
+#include "plane_project.h"
+#include "edge_lengths.h"
+
+template <
+  typename DerivedV, 
+  typename DerivedF,
+  typename DerivedU,
+  typename DerivedUF,
+  typename Scalar>
+IGL_INLINE void igl::plane_project(
+  const Eigen::PlainObjectBase<DerivedV> & V, 
+  const Eigen::PlainObjectBase<DerivedF> & F, 
+  Eigen::PlainObjectBase<DerivedU> & U,
+  Eigen::PlainObjectBase<DerivedUF> & UF, 
+  Eigen::SparseMatrix<Scalar>& I)
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(F.rows() == 3 && "F should contain triangles");
+  typedef Eigen::PlainObjectBase<DerivedV> MatrixX;
+  MatrixX l;
+  edge_lengths(V,F,l);
+  // Number of faces
+  const int m = F.rows();
+
+  // First corner at origin
+  U = Eigen::PlainObjectBase<DerivedU>::Zero(m*3,2);
+  // Second corner along x-axis
+  U.block(m,0,m,1) = l.col(2);
+  // Third corner rotated onto plane
+  U.block(m*2,0,m,1) = 
+    (-l.col(0).array().square() + 
+     l.col(1).array().square() + 
+     l.col(2).array().square())/(2.*l.col(2).array());
+  U.block(m*2,1,m,1) =
+    (l.col(1).array().square()-U.block(m*2,0,m,1).array().square()).sqrt();
+
+  typedef Triplet<Scalar> IJV;
+  vector<IJV > ijv;
+  ijv.reserve(3*m);
+  UF.resize(m,3);
+  for(int f = 0;f<m;f++)
+  {
+    for(int c = 0;c<3;c++)
+    {
+      UF(f,c) = c*m+f;
+      ijv.push_back(IJV(F(f,c),c*m+f,1));
+    }
+  }
+  I.resize(V.rows(),m*3);
+  I.setFromTriplets(ijv.begin(),ijv.end());
+}
+
+#ifndef IGL_HEADER_ONLY
+template void igl::plane_project<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/plane_project.h

@@ -0,0 +1,48 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_PLANE_PROJECT_H
+#define IGL_PLANE_PROJECT_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl 
+{
+  // Project each triangle to the plane
+  //
+  // [U,UF,I] = plane_project(V,F)
+  //
+  // Inputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of mesh indices
+  // Outputs:
+  //   U  #F*3 by 2 list of triangle positions
+  //   UF  #F by 3 list of mesh indices into U
+  //   I  #V by #F such that I(i,j) = 1 implies U(j,:) corresponds to V(i,:)
+  //
+  template <
+    typename DerivedV, 
+    typename DerivedF,
+    typename DerivedU,
+    typename DerivedUF,
+    typename Scalar>
+  IGL_INLINE void plane_project(
+    const Eigen::PlainObjectBase<DerivedV> & V, 
+    const Eigen::PlainObjectBase<DerivedF> & F, 
+    Eigen::PlainObjectBase<DerivedU> & U,
+    Eigen::PlainObjectBase<DerivedUF> & UF, 
+    Eigen::SparseMatrix<Scalar>& I);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "plane_project.cpp"
+#endif
+
+#endif
+

include/igl/svd3x3/arap.cpp

@@ -13,6 +13,7 @@
 #include <igl/covariance_scatter_matrix.h>
 #include <igl/speye.h>
 #include <igl/mode.h>
+#include <igl/plane_project.h>
 #include <igl/slice.h>
 #include <igl/arap_rhs.h>
 #include <igl/repdiag.h>
@@ -32,7 +33,6 @@ IGL_INLINE bool igl::arap_precomputation(
   const Eigen::PlainObjectBase<Derivedb> & b,
   ARAPData & data)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   typedef typename DerivedV::Scalar Scalar;
@@ -48,12 +48,25 @@ IGL_INLINE bool igl::arap_precomputation(
   //const int dim = V.cols();
   assert((dim == 3 || dim ==2) && "dim should be 2 or 3");
   data.dim = dim;
-  data.Vdim = V.cols();
   //assert(dim == 3 && "Only 3d supported");
   // Defaults
-  data.f_ext = Eigen::MatrixXd::Zero(n,data.dim);
+  data.f_ext = MatrixXd::Zero(n,data.dim);
+
+  assert(data.dim <= V.cols() && "solve dim should be <= embedding");
+  bool flat = (V.cols() - data.dim)==1;
 
-  SparseMatrix<Scalar> L;
+  PlainObjectBase<DerivedV> plane_V;
+  PlainObjectBase<DerivedF> plane_F;
+  typedef SparseMatrix<Scalar> SparseMatrixS;
+  SparseMatrixS ref_map,ref_map_dim;
+  if(flat)
+  {
+    plane_project(V,F,plane_V,plane_F,ref_map);
+    repdiag(ref_map,dim,ref_map_dim);
+  }
+  const PlainObjectBase<DerivedV>& ref_V = (flat?plane_V:V);
+  const PlainObjectBase<DerivedF>& ref_F = (flat?plane_F:F);
+  SparseMatrixS L;
   cotmatrix(V,F,L);
 
   ARAPEnergyType eff_energy = data.energy;
@@ -81,7 +94,11 @@ IGL_INLINE bool igl::arap_precomputation(
 
   // Get covariance scatter matrix, when applied collects the covariance
   // matrices used to fit rotations to during optimization
-  covariance_scatter_matrix(V,F,eff_energy,data.CSM);
+  covariance_scatter_matrix(ref_V,ref_F,eff_energy,data.CSM);
+  assert(data.CSM.rows() == ref_F.rows()*data.dim);
+  data.CSM = (data.CSM * ref_map_dim.transpose()).eval();
+  assert(data.CSM.cols() == V.rows()*data.dim);
+  assert(data.CSM.rows() == ref_F.rows()*data.dim);
 
   // Get group sum scatter matrix, when applied sums all entries of the same
   // group according to G
@@ -115,12 +132,13 @@ IGL_INLINE bool igl::arap_precomputation(
     group_sum_matrix(data.G,G_sum);
   }
   SparseMatrix<double> G_sum_dim;
-  repdiag(G_sum,data.Vdim,G_sum_dim);
+  repdiag(G_sum,data.dim,G_sum_dim);
   assert(G_sum_dim.cols() == data.CSM.rows());
   data.CSM = (G_sum_dim * data.CSM).eval();
 
 
-  arap_rhs(V,F,data.dim,eff_energy,data.K);
+  arap_rhs(ref_V,ref_F,data.dim,eff_energy,data.K);
+  data.K = (ref_map_dim * data.K).eval();
   assert(data.K.rows() == data.n*data.dim);
 
   SparseMatrix<double> Q = (-0.5*L).eval();
@@ -188,19 +206,16 @@ IGL_INLINE bool igl::arap_solve(
       U.row(data.b(bi)) = bc.row(bi);
     }
 
-    const auto & Udim = U.replicate(data.Vdim,1);
+    const auto & Udim = U.replicate(data.dim,1);
     assert(U.cols() == data.dim);
-    assert(data.Vdim >= data.dim);
     // As if U.col(2) was 0
-    MatrixXd S = MatrixXd::Zero(data.CSM.rows(),data.Vdim);
-    S.leftCols(data.dim) = data.CSM * Udim;
+    MatrixXd S = data.CSM * Udim;
 
-    const int Rdim = data.Vdim;
+    const int Rdim = data.dim;
     MatrixXd R(Rdim,data.CSM.rows());
     if(R.rows() == 2)
     {
       fit_rotations_planar(S,R);
-      cerr<<"PLANAR"<<endl;
     }else
     {
 #ifdef __SSE__ // fit_rotations_SSE will convert to float if necessary

include/igl/svd3x3/arap.h

@@ -32,7 +32,6 @@ namespace igl
     // solver_data  quadratic solver data
     // b  list of boundary indices into V
     // dim  dimension being used for solving
-    // Vdim  dimension of original V (determines dimension of rotation fitting)
     int n;
     Eigen::VectorXi G;
     ARAPEnergyType energy;
@@ -45,7 +44,6 @@ namespace igl
     min_quad_with_fixed_data<double> solver_data;
     Eigen::VectorXi b;
     int dim;
-    int Vdim;
       ARAPData():
         n(0),
         G(),
@@ -58,8 +56,7 @@ namespace igl
         CSM(),
         solver_data(),
         b(),
-        dim(-1), // force this to be set by _precomputation
-        Vdim(-1) // force this to be set by _precomputation
+        dim(-1) // force this to be set by _precomputation
     {
     };
   };

tutorial/readme.md

@@ -1,5 +1,8 @@
-xhtml header:   <script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
 css: style.css
+html header:   <script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
+<link rel="stylesheet" href="http://yandex.st/highlightjs/7.3/styles/default.min.css">
+<script src="http://yandex.st/highlightjs/7.3/highlight.min.js"></script>
+<script>hljs.initHighlightingOnLoad();</script>
 
 # Introduction
 
@@ -18,7 +21,7 @@ All examples depends on glfw, glew and anttweakbar. A copy
 of the sourcecode of each library is provided together with libigl
 and they can be precompiled using:
 
-** Alec: Is this just compiling the dependencies? Then perhaps rename `compile_dependencies_*`**
+**Alec: Is this just compiling the dependencies? Then perhaps rename `compile_dependencies_*`**
 
     sh compile_macosx.sh (MACOSX)
     sh compile_linux.sh (LINUX)
@@ -29,7 +32,13 @@ On Linux and MacOSX, you can use the provided bash script:
 
     sh ../compile_example.sh
 
-## Optional compilation with libigl as static library
+## (Optional: compilation with libigl as static library)
+
+By default, libigl is a _headers only_ library, thus it does not require
+compilation. However, one can precompile libigl as a statically linked library.
+See `../README.md` in the main directory for compilations instructions to
+produce `libigl.a` and other libraries. Once compiled, these examples can be
+compiled using the `CMAKE` flag `-DLIBIGL_USE_STATIC_LIBRARY=ON`:
 
     ../compile_example.sh -DLIBIGL_USE_STATIC_LIBRARY=ON
 
@@ -66,3 +75,13 @@ elliptic, hyperbolic and parabolic vertices on the domain.
 
 ![The `GaussianCurvature` example computes discrete Gaussian curvature and visualizes it in
 pseudocolor.](images/bumpy-gaussian-curvature.jpg)
+
+
+This is an example of syntax highlighted code:
+```cpp
+#include <foo.html>
+int main(int argc, char * argv[])
+{
+  return 0;
+}
+```