Merge branch 'master' of https://github.com/libigl/libigl

Former-commit-id: 3a307021edecbc8cf919cf01a58b856965f56a39

documentation/implemented-papers.tex

@@ -40,6 +40,9 @@ choosing skinning transformations \cite{Jacobson:FAST:2012}.
 ``skeletal subspace deformation'', or ``enveloping''. This technique is often
 attributed to \cite{Magnenat-Thalmann:1988:JLD}.
 
+\paragraph{\texttt{winding\_number}} implements ``Generalized Winding Numbers''
+\cite{Jacobson:WN:2013}
+
 \bibliographystyle{acmsiggraph}
 \bibliography{references} 
 

documentation/references.bib.REMOVED.git-id

@@ -1 +1 @@
-92b6b20756fd71dd62d7b52066a7bcf07e4acd28
+346cd2c8206ee2c83e7aef50adb832e483320d09

examples/ambient-occlusion-mex/compile.m

@@ -4,6 +4,6 @@ mex -v -o ambient_occlusion -DMEX -largeArrayDims ...
   CXXFLAGS="\$CXXFLAGS -m64 -msse4.2 -fopenmp" ...
   -I/usr/local/igl/libigl/external/embree/embree/ ...
   -I/usr/local/igl/libigl/external/embree/ ...
-  -L/usr/local/igl/libigl/external/embree/bin -lembree -lsys ... 
+  -L/usr/local/igl/libigl/external/embree/build -lembree -lsys ... 
   mexFunction.cpp parse_rhs.cpp;
 

examples/arap/Makefile

@@ -14,6 +14,7 @@ LIBIGL_LIB=-L$(LIBIGL)/lib -ligl -liglpng -liglsvd3x3
 
 EIGEN3_INC=-I/opt/local/include/eigen3 -I/opt/local/include/eigen3/unsupported
 
+
 #EMBREE=$(LIBIGL)/external/embree
 #EMBREE_INC=-I$(EMBREE)/rtcore -I$(EMBREE)/common
 #EMBREE_LIB=-L$(EMBREE)/build -lrtcore -lsys
@@ -36,8 +37,8 @@ ifeq ($(UNAME), Darwin)
 	ANTTWEAKBAR_LIB+=-framework AppKit
 endif
 
-INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(YIMG_INC)
-LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(YIMG_LIB)
+INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(YIMG_INC) 
+LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(YIMG_LIB) 
 
 example: example.o
 	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) $(LIB) -o example $^

examples/arap/example.cpp

@@ -28,7 +28,6 @@
 #include <igl/material_colors.h>
 #include <igl/barycenter.h>
 #include <igl/matlab_format.h>
-#include <igl/material_colors.h>
 #include <igl/ReAntTweakBar.h>
 #include <igl/pathinfo.h>
 #include <igl/Camera.h>
@@ -70,6 +69,9 @@ Eigen::Quaterniond animation_to_quat;
 // Use vector for range-based `for`
 std::vector<State> undo_stack;
 std::vector<State> redo_stack;
+bool paused = false;
+double t = 0;
+double pause_time = 0.0;
 
 void push_undo()
 {
@@ -133,7 +135,7 @@ float light_pos[4] = {0.1,0.1,-0.9,0};
 Eigen::MatrixXd V,U,N,C,mid;
 Eigen::VectorXi S;
 igl::ARAPData arap_data;
-Eigen::MatrixXi F;
+Eigen::MatrixXi F,T;
 int selected_col = 0;
 // Faces
 // Bounding box diagonal length
@@ -164,6 +166,14 @@ bool init_arap()
   assert(S.rows() == V.rows());
   C.resize(S.rows(),3);
   MatrixXd bc = MatrixXd::Zero(b.size(),S.maxCoeff()+1);
+  MatrixXi * Ele;
+  if(T.rows()>0)
+  {
+    Ele = &T;
+  }else
+  {
+    Ele = &F;
+  }
   // get b from S
   {
     int bi = 0;
@@ -200,12 +210,15 @@ bool init_arap()
   VectorXi _S;
   VectorXd _D;
   MatrixXd W;
-  if(!harmonic(V,F,b,bc,1,W))
+  if(!harmonic(V,*Ele,b,bc,1,W))
   {
     return false;
   }
-  partition(W,100,arap_data.G,_S,_D);
-  return arap_precomputation(V,F,b,arap_data);
+  //arap_data.with_dynamics = true;
+  //arap_data.h = 0.5;
+  //arap_data.max_iter = 100;
+  //partition(W,100,arap_data.G,_S,_D);
+  return arap_precomputation(V,*Ele,b,arap_data);
 }
 
 bool update_arap()
@@ -216,6 +229,10 @@ bool update_arap()
   MatrixXd bc(num_in_selection(S),V.cols());
   // get b from S
   {
+    if(!paused)
+    {
+      t = get_seconds()-pause_time;
+    }
     int bi = 0;
     for(int v = 0;v<S.rows(); v++)
     {
@@ -227,15 +244,15 @@ bool update_arap()
           case 0:
           {
             const double r = mid(0)*0.25;
-            bc(bi,0) += r*cos(0.5*get_seconds()*2.*PI);
-            bc(bi,1) -= r+r*sin(0.5*get_seconds()*2.*PI);
+            bc(bi,0) += r*cos(0.5*t*2.*PI);
+            bc(bi,1) -= r+r*sin(0.5*t*2.*PI);
             break;
           }
           case 1:
           {
             const double r = mid(1)*0.15;
-            bc(bi,1) += r+r*cos(0.15*get_seconds()*2.*PI);
-            bc(bi,2) -= r*sin(0.15*get_seconds()*2.*PI);
+            bc(bi,1) += r+r*cos(0.15*t*2.*PI);
+            bc(bi,2) -= r*sin(0.15*t*2.*PI);
 
             //// Pull-up
             //bc(bi,0) += 0.42;//mid(0)*0.5;
@@ -579,6 +596,7 @@ void mouse_drag(int mouse_x, int mouse_y)
 void key(unsigned char key, int mouse_x, int mouse_y)
 {
   using namespace std;
+  using namespace igl;
   switch(key)
   {
     // ESC
@@ -587,6 +605,20 @@ void key(unsigned char key, int mouse_x, int mouse_y)
     // ^C
     case char(3):
       exit(0);
+    case ' ':
+      {
+        static double pause_start,pause_stop;
+        paused = !paused;
+        if(paused)
+        {
+          pause_start = get_seconds();
+        }else
+        {
+          pause_stop = get_seconds();
+          pause_time += (pause_stop-pause_start);
+        }
+        break;
+      }
     default:
       if(!TwEventKeyboardGLUT(key,mouse_x,mouse_y))
       {
@@ -606,6 +638,8 @@ int main(int argc, char * argv[])
   // init mesh
   string filename = "../shared/decimated-knight.obj";
   string sfilename = "../shared/decimated-knight-selection.dmat";
+  //string filename = "../shared/decimated-knight.mesh";
+  //string sfilename = "../shared/decimated-knight-1-selection.dmat";
   if(argc < 3)
   {
     cerr<<"Usage:"<<endl<<"    ./example input.obj selection.dmat"<<endl;
@@ -616,11 +650,27 @@ int main(int argc, char * argv[])
     filename = argv[1];
     sfilename = argv[2];
   }
+  string d,b,ext,f;
+  pathinfo(filename,d,b,ext,f);
+  // Convert extension to lower case
+  transform(ext.begin(), ext.end(), ext.begin(), ::tolower);
 
   vector<vector<double > > vV,vN,vTC;
   vector<vector<int > > vF,vTF,vFN;
   // Convert extension to lower case
-  if(!igl::readOBJ(filename,vV,vTC,vN,vF,vTF,vFN))
+  if(ext == "obj")
+  {
+    if(!igl::readOBJ(filename,vV,vTC,vN,vF,vTF,vFN))
+    {
+      return 1;
+    }
+  }else if(ext == "mesh")
+  {
+    if(!igl::readMESH(filename,V,T,F))
+    {
+      return 1;
+    }
+  }else
   {
     return 1;
   }

include/igl/WindingNumberTree.h

@@ -322,6 +322,7 @@ inline double igl::WindingNumberTree<Point>::winding_number_boundary(const Point
   using namespace std;
 
   double w = 0;
+  // `cap` is already flipped inside out, so we don't need to flip sign of w
   igl::winding_number_3(
     V.data(),
     V.rows(),

include/igl/cgal/SelfIntersectMesh.h

@@ -85,6 +85,8 @@ namespace igl
     public:
       // Constructs (VV,FF) a new mesh with self-intersections of (V,F)
       // subdivided
+      //
+      // See also: selfintersect.h
       inline SelfIntersectMesh(
         const Eigen::MatrixXd & V,
         const Eigen::MatrixXi & F,
@@ -92,7 +94,9 @@ namespace igl
         Eigen::MatrixXd & VV,
         Eigen::MatrixXi & FF,
         Eigen::MatrixXi & IF,
-        Eigen::VectorXi & J);
+        Eigen::VectorXi & J,
+        Eigen::VectorXi & IM
+        );
     private:
       // Helper function to mark a face as offensive
       //
@@ -252,7 +256,8 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
   Eigen::MatrixXd & VV,
   Eigen::MatrixXi & FF,
   Eigen::MatrixXi & IF,
-  Eigen::VectorXi & J):
+  Eigen::VectorXi & J,
+  Eigen::VectorXi & IM):
   V(V),
   F(F),
   count(0),
@@ -357,8 +362,8 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
           //  P[o].to_3d(vit->point())<<endl;
 #ifndef NDEBUG
           // I want to be sure that the original corners really show up as the
-          // original corners of the CDT. I.e. I don't trust CGAL to maintain the
-          // order
+          // original corners of the CDT. I.e. I don't trust CGAL to maintain
+          // the order
           assert(T[f].vertex(i) == P[o].to_3d(vit->point()));
 #endif
           // For first three, use original index in F
@@ -496,6 +501,38 @@ inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
       i++;
     }
   }
+  IM.resize(VV.rows(),1);
+  map<Point_3,int> vv2i;
+  // Safe to check for duplicates using double for original vertices: if
+  // incoming reps are different then the points are unique.
+  for(int v = 0;v<V.rows();v++)
+  {
+    const Point_3 p(V(v,0),V(v,1),V(v,2));
+    if(vv2i.count(p)==0)
+    {
+      vv2i[p] = v;
+    }
+    assert(vv2i.count(p) == 1);
+    IM(v) = vv2i[p];
+  }
+  // Must check for duplicates of new vertices using exact.
+  {
+    int v = V.rows();
+    for(
+      typename list<Point_3>::const_iterator nvit = NV.begin();
+      nvit != NV.end();
+      nvit++)
+    {
+      const Point_3 & p = *nvit;
+      if(vv2i.count(p)==0)
+      {
+        vv2i[p] = v;
+      }
+      assert(vv2i.count(p) == 1);
+      IM(v) = vv2i[p];
+      v++;
+    }
+  }
 
   // Q: Does this give the same result as TETGEN?
   // A: For the cow and beast, yes.

include/igl/cgal/intersect_other.cpp

@@ -128,3 +128,7 @@ IGL_INLINE void igl::intersect_other(
   }
 
 }
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

include/igl/cgal/mesh_to_cgal_triangle_list.h

@@ -14,6 +14,9 @@ namespace igl
 {
   // Convert a mesh (V,F) to a list of CGAL triangles
   //
+  // Templates:
+  //   Kernal  CGAL computation and construction kernel (e.g.
+  //     CGAL::Exact_predicates_exact_constructions_kernel)
   // Inputs:
   //   V  #V by 3 list of vertex positions
   //   F  #F by 3 list of triangle indices

include/igl/cgal/point_mesh_squared_distance.cpp

@@ -0,0 +1,63 @@
+// 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/.
+#include "point_mesh_squared_distance.h"
+#include "mesh_to_cgal_triangle_list.h"
+
+template <typename Kernel>
+IGL_INLINE void igl::point_mesh_squared_distance(
+  const Eigen::MatrixXd & P,
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  Eigen::VectorXd & sqrD,
+  Eigen::VectorXi & I,
+  Eigen::MatrixXd & C)
+{
+  using namespace std;
+  typedef CGAL::Point_3<Kernel>    Point_3;
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+
+  // Must be 3D
+  assert(V.cols() == 3);
+  assert(P.cols() == 3);
+  // Must be triangles
+  assert(F.cols() == 3);
+  // Make list of cgal triangles
+  Tree tree;
+  vector<Triangle_3> T;
+  mesh_to_cgal_triangle_list(V,F,T);
+  tree.insert(T.begin(),T.end());
+
+  tree.accelerate_distance_queries();
+  const int n = P.rows();
+  sqrD.resize(n,1);
+  I.resize(n,1);
+  C.resize(n,P.cols());
+  for(int p = 0;p < n;p++)
+  {
+    Point_3 query(P(p,0),P(p,1),P(p,2));
+    // Find closest point and primitive id
+    Point_and_primitive_id pp = tree.closest_point_and_primitive(query);
+    Point_3 closest_point = pp.first;
+    C(p,0) = CGAL::to_double(closest_point[0]);
+    C(p,1) = CGAL::to_double(closest_point[1]);
+    C(p,2) = CGAL::to_double(closest_point[2]);
+    sqrD(p) = CGAL::to_double((closest_point-query).squared_length());
+    I(p) = pp.second - T.begin();
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
+
+#endif

include/igl/cgal/point_mesh_squared_distance.h

@@ -0,0 +1,44 @@
+// 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_POINT_MESH_SQUARED_DISTANCE_H
+#define IGL_POINT_MESH_SQUARED_DISTANCE_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include "CGAL_includes.hpp"
+namespace igl
+{
+  // Compute distances from a set of points P to a triangle mesh (V,F)
+  //
+  // Templates:
+  //   Kernal  CGAL computation and construction kernel (e.g.
+  //     CGAL::Simple_cartesian<double>)
+  // Inputs:
+  //   P  #P by 3 list of query point positions
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   sqrD  #P list of smallest squared distances
+  //   I  #P list of facet indices corresponding to smallest distances
+  //   C  #P by 3 list of closest points
+  //
+  // Known bugs: This only computes distances to triangles. So unreferenced
+  // vertices are ignored.
+  template <typename Kernel>
+  IGL_INLINE void point_mesh_squared_distance(
+    const Eigen::MatrixXd & P,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    Eigen::VectorXd & sqrD,
+    Eigen::VectorXi & I,
+    Eigen::MatrixXd & C);
+}
+#ifdef IGL_HEADER_ONLY
+#  include "point_mesh_squared_distance.cpp"
+#endif
+
+#endif

include/igl/cgal/selfintersect.cpp

@@ -17,7 +17,8 @@ IGL_INLINE void igl::selfintersect(
   Eigen::MatrixXd & VV,
   Eigen::MatrixXi & FF,
   Eigen::MatrixXi & IF,
-  Eigen::VectorXi & J)
+  Eigen::VectorXi & J,
+  Eigen::VectorXi & IM)
 {
   using namespace std;
   if(params.detect_only)
@@ -35,7 +36,7 @@ IGL_INLINE void igl::selfintersect(
 //#endif
 
     typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF,J);
+    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF,J,IM);
 
 //#ifdef __APPLE__
 //    signal(SIGFPE,SIG_DFL);
@@ -44,6 +45,6 @@ IGL_INLINE void igl::selfintersect(
   }else
   {
     typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
-    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF,J);
+    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF,J,IM);
   }
 }

include/igl/cgal/selfintersect.h

@@ -47,6 +47,12 @@ namespace igl
   //   IF  #intersecting face pairs by 2  list of intersecting face pairs,
   //     indexing F
   //   J  #FF list of indices into F denoting birth triangle
+  //   IM  #VV list of indices into VV of unique vertices.
+  //
+  // Known bugs: If an existing edge in (V,F) lies exactly on another face then
+  // any resulting additional vertices along that edge may not get properly
+  // connected so that the output mesh has the same global topology. This is
+  // because 
   IGL_INLINE void selfintersect(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
@@ -54,7 +60,8 @@ namespace igl
     Eigen::MatrixXd & VV,
     Eigen::MatrixXi & FF,
     Eigen::MatrixXi & IF,
-    Eigen::VectorXi & J);
+    Eigen::VectorXi & J,
+    Eigen::VectorXi & IM);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/harmonic.cpp

@@ -24,7 +24,16 @@ IGL_INLINE bool igl::harmonic(
   using namespace Eigen;
   SparseMatrix<double> L,M,Mi;
   cotmatrix(V,F,L);
-  massmatrix(V,F,MASSMATRIX_VORONOI,M);
+  switch(F.cols())
+  {
+    case 3:
+      massmatrix(V,F,MASSMATRIX_VORONOI,M);
+      break;
+    case 4:
+    default:
+      massmatrix(V,F,MASSMATRIX_BARYCENTRIC,M);
+      break;
+  }
   invert_diag(M,Mi);
   SparseMatrix<double> Q = -L;
   for(int p = 1;p<k;p++)

include/igl/is_symmetric.cpp

@@ -15,6 +15,7 @@ IGL_INLINE bool igl::is_symmetric(const Eigen::SparseMatrix<T>& A)
   {
     return false;
   }
+  assert(A.size() != 0);
   Eigen::SparseMatrix<T> AT = A.transpose();
   Eigen::SparseMatrix<T> AmAT = A-AT;
   //// Eigen screws up something with LLT if you try to do
@@ -32,6 +33,7 @@ IGL_INLINE bool igl::is_symmetric(
   {
     return false;
   }
+  assert(A.size() != 0);
   const typename Eigen::PlainObjectBase<DerivedA>& AT = A.transpose();
   const typename Eigen::PlainObjectBase<DerivedA>& AmAT = A-AT;
   //// Eigen screws up something with LLT if you try to do
@@ -52,6 +54,7 @@ IGL_INLINE bool igl::is_symmetric(
   {
     return false;
   }
+  assert(A.size() != 0);
   SparseMatrix<AType> AT = A.transpose();
   SparseMatrix<AType> AmAT = A-AT;
   VectorXi AmATI,AmATJ;

include/igl/massmatrix.cpp

@@ -21,12 +21,21 @@ IGL_INLINE void igl::massmatrix(
 {
   using namespace Eigen;
   using namespace std;
-  assert(type!=MASSMATRIX_FULL);
 
   const int n = V.rows();
   const int m = F.rows();
   const int simplex_size = F.cols();
 
+  MassMatrixType eff_type = type;
+  // Use voronoi of for triangles by default, otherwise barycentric
+  if(type == MASSMATRIX_DEFAULT)
+  {
+    eff_type = (simplex_size == 3?MASSMATRIX_VORONOI:MASSMATRIX_BARYCENTRIC);
+  }
+
+  // Not yet supported
+  assert(type!=MASSMATRIX_FULL);
+
   Matrix<int,Dynamic,1> MI;
   Matrix<int,Dynamic,1> MJ;
   Matrix<Scalar,Dynamic,1> MV;
@@ -52,7 +61,7 @@ IGL_INLINE void igl::massmatrix(
       dblA(i) = 2.0*sqrt(s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2)));
     }
 
-    switch(type)
+    switch(eff_type)
     {
       case MASSMATRIX_BARYCENTRIC:
         // diagonal entries for each face corner
@@ -115,13 +124,13 @@ IGL_INLINE void igl::massmatrix(
         assert(false && "Implementation incomplete");
         break;
       default:
-        assert(false && "Unknown Mass matrix type");
+        assert(false && "Unknown Mass matrix eff_type");
     }
 
   }else if(simplex_size == 4)
   {
     assert(V.cols() == 3);
-    assert(type == MASSMATRIX_BARYCENTRIC);
+    assert(eff_type == MASSMATRIX_BARYCENTRIC);
     MI.resize(m*4,1); MJ.resize(m*4,1); MV.resize(m*4,1);
     MI.block(0*m,0,m,1) = F.col(0);
     MI.block(1*m,0,m,1) = F.col(1);

include/igl/massmatrix.h

@@ -18,11 +18,12 @@ namespace igl
 
   enum MassMatrixType
   {
-    MASSMATRIX_BARYCENTRIC,
-    MASSMATRIX_VORONOI,
-    MASSMATRIX_FULL
+    MASSMATRIX_BARYCENTRIC = 0,
+    MASSMATRIX_VORONOI = 1,
+    MASSMATRIX_FULL = 2,
+    MASSMATRIX_DEFAULT = 3,
+    NUM_MASSMATRIX = 4
   };
-#define NUM_MASSMATRIXTYPE 3
 
   // Constructs the mass (area) matrix for a given mesh (V,F).
   //

include/igl/project_to_line_segment.h

@@ -43,7 +43,7 @@ namespace igl
 }
 
 #ifdef IGL_HEADER_ONLY
-#  include "project_to_line.cpp"
+#  include "project_to_line_segment.cpp"
 #endif
 
 #endif

include/igl/readOFF.cpp

@@ -29,7 +29,7 @@ IGL_INLINE bool igl::readOFF(
   char header[1000];
   const std::string OFF("OFF");
   const std::string NOFF("NOFF");
-  if(!fscanf(off_file,"%s\n",header)==1
+  if(fscanf(off_file,"%s\n",header)!=1
      || !(OFF == header || NOFF == header))
   {
     printf("Error: %s's first line should be OFF or NOFF not %s...",off_file_name.c_str(),header);

include/igl/svd3x3/arap.cpp

@@ -8,6 +8,7 @@
 #include "arap.h"
 #include <igl/colon.h>
 #include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
 #include <igl/group_sum_matrix.h>
 #include <igl/covariance_scatter_matrix.h>
 #include <igl/speye.h>
@@ -40,7 +41,7 @@ IGL_INLINE bool igl::arap_precomputation(
   assert((b.size() == 0 || b.minCoeff() >=0) && "b out of bounds");
   // remember b
   data.b = b;
-  assert(F.cols() == 3 && "For now only triangles");
+  //assert(F.cols() == 3 && "For now only triangles");
   // dimension
   const int dim = V.cols();
   assert(dim == 3 && "Only 3d supported");
@@ -66,14 +67,15 @@ IGL_INLINE bool igl::arap_precomputation(
         break;
       case 4:
         eff_energy = ARAP_ENERGY_TYPE_ELEMENTS;
+        break;
       default:
         assert(false);
     }
   }
 
 
-  // Get covariance scatter matrix, when applied collects the covariance matrices
-  // used to fit rotations to during optimization
+  // 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);
 
   // Get group sum scatter matrix, when applied sums all entries of the same
@@ -81,13 +83,19 @@ IGL_INLINE bool igl::arap_precomputation(
   SparseMatrix<double> G_sum;
   if(data.G.size() == 0)
   {
-    speye(n,G_sum);
+    if(eff_energy == ARAP_ENERGY_TYPE_ELEMENTS)
+    {
+      speye(F.rows(),G_sum);
+    }else
+    {
+      speye(n,G_sum);
+    }
   }else
   {
     // groups are defined per vertex, convert to per face using mode
-    Eigen::Matrix<int,Eigen::Dynamic,1> GG;
     if(eff_energy == ARAP_ENERGY_TYPE_ELEMENTS)
     {
+      Eigen::Matrix<int,Eigen::Dynamic,1> GG;
       MatrixXi GF(F.rows(),F.cols());
       for(int j = 0;j<F.cols();j++)
       {
@@ -96,20 +104,34 @@ IGL_INLINE bool igl::arap_precomputation(
         GF.col(j) = GFj;
       }
       mode<int>(GF,2,GG);
-    }else
-    {
-      GG=data.G;
+      data.G=GG;
     }
     //printf("group_sum_matrix()\n");
-    group_sum_matrix(GG,G_sum);
+    group_sum_matrix(data.G,G_sum);
   }
   SparseMatrix<double> G_sum_dim;
   repdiag(G_sum,dim,G_sum_dim);
+  assert(G_sum_dim.cols() == data.CSM.rows());
   data.CSM = (G_sum_dim * data.CSM).eval();
 
+
   arap_rhs(V,F,eff_energy,data.K);
 
   SparseMatrix<double> Q = (-0.5*L).eval();
+
+  if(data.with_dynamics)
+  {
+    const double h = data.h;
+    assert(h != 0);
+    SparseMatrix<double> M;
+    massmatrix(V,F,MASSMATRIX_DEFAULT,data.M);
+    SparseMatrix<double> DQ = 0.5/(h*h)*data.M;
+    Q += DQ;
+    // Dummy external forces
+    data.f_ext = MatrixXd::Zero(n,dim);
+    data.vel = MatrixXd::Zero(n,dim);
+  }
+
   return min_quad_with_fixed_precompute(
     Q,b,SparseMatrix<double>(),true,data.solver_data);
 }
@@ -134,7 +156,14 @@ IGL_INLINE bool igl::arap_solve(
     // terrible initial guess.. should at least copy input mesh
     U = MatrixXd::Zero(data.n,dim);
   }
+  // changes each arap iteration
   MatrixXd U_prev = U;
+  // doesn't change for fixed with_dynamics timestep
+  MatrixXd U0;
+  if(data.with_dynamics)
+  {
+    U0 = U_prev;
+  }
   while(iter < data.max_iter)
   {
     U_prev = U;
@@ -157,6 +186,8 @@ IGL_INLINE bool igl::arap_solve(
     //}
 
 
+    // Number of rotations: #vertices or #elements
+    int num_rots = data.K.cols()/dim/dim;
     // distribute group rotations to vertices in each group
     MatrixXd eff_R;
     if(data.G.size() == 0)
@@ -165,21 +196,40 @@ IGL_INLINE bool igl::arap_solve(
       eff_R = R;
     }else
     {
-      eff_R.resize(dim,dim*n);
-      for(int v = 0;v<n;v++)
+      eff_R.resize(dim,num_rots*dim);
+      for(int r = 0;r<num_rots;r++)
       {
-        eff_R.block(0,dim*v,dim,dim) = 
-          R.block(0,dim*data.G(v),dim,dim);
+        eff_R.block(0,dim*r,dim,dim) = 
+          R.block(0,dim*data.G(r),dim,dim);
       }
     }
 
+    MatrixXd Dl;
+    if(data.with_dynamics)
+    {
+      assert(M.rows() == n && 
+        "No mass matrix. Call arap_precomputation if changing with_dynamics");
+      const double h = data.h;
+      assert(h != 0);
+      //Dl = 1./(h*h*h)*M*(-2.*V0 + Vm1) - fext;
+      // data.vel = (V0-Vm1)/h
+      // h*data.vel = (V0-Vm1)
+      // -h*data.vel = -V0+Vm1)
+      // -V0-h*data.vel = -2V0+Vm1
+      Dl = 1./(h*h)*data.M*(-U0 - h*data.vel) - data.f_ext;
+    }
+
     VectorXd Rcol;
-    columnize(eff_R,n,2,Rcol);
+    columnize(eff_R,num_rots,2,Rcol);
     VectorXd Bcol = -data.K * Rcol;
     for(int c = 0;c<dim;c++)
     {
       VectorXd Uc,Bc,bcc,Beq;
       Bc = Bcol.block(c*n,0,n,1);
+      if(data.with_dynamics)
+      {
+        Bc += Dl.col(c);
+      }
       bcc = bc.col(c);
       min_quad_with_fixed_solve(
         data.solver_data,
@@ -187,10 +237,14 @@ IGL_INLINE bool igl::arap_solve(
         Uc);
       U.col(c) = Uc;
     }
-    
 
     iter++;
   }
+  if(data.with_dynamics)
+  {
+    // Keep track of velocity for next time
+    data.vel = (U-U0)/data.h;
+  }
   return true;
 }
 

include/igl/svd3x3/arap.h

@@ -21,21 +21,24 @@ namespace igl
     // G  #V list of group indices (1 to k) for each vertex, such that vertex i
     //    is assigned to group G(i)
     // energy  type of energy to use
-    // with_dynamics  whether using dynamics 
+    // with_dynamics  whether using dynamics (need to call arap_precomputation
+    //   after changing)
     // f_ext  #V by dim list of external forces
+    // vel  #V by dim list of velocities
     // h  dynamics time step
     // max_iter  maximum inner iterations
     // K  rhs pre-multiplier
+    // M  mass matrix
     // solver_data  quadratic solver data
     // b  list of boundary indices into V
     int n;
     Eigen::VectorXi G;
     ARAPEnergyType energy;
     bool with_dynamics;
-    Eigen::MatrixXd f_ext;
+    Eigen::MatrixXd f_ext,vel;
     double h;
     int max_iter;
-    Eigen::SparseMatrix<double> K;
+    Eigen::SparseMatrix<double> K,M;
     Eigen::SparseMatrix<double> CSM;
     min_quad_with_fixed_data<double> solver_data;
     Eigen::VectorXi b;

include/igl/xml/XMLSerializer.h.REMOVED.git-id

@@ -1 +1 @@
-1dc638aa63c3131bfb705a8549e5b7f47c780223
+8732db96cb04680209315d0fd383a4e392159b2e