fast exmaple

Former-commit-id: e92defde6dc453444d68ec6c7989a1a2c08c5767

.gitignore

@@ -59,6 +59,7 @@ external/tinyxml2/tinyxml2.pc
 external/yimg/showpng
 README.html
 tutorial/readme.html
+tutorial/tutorial.html
 tutorial/*/build/*
 tutorial/*/Makefile
 external/glew/build

build/Makefile

@@ -125,7 +125,7 @@ obj/%.o: ../include/igl/%.cpp ../include/igl/%.h
 
 clean:
 	rm -rf ../lib/igl.framework/
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libigl.a
 	make -C ../examples clean
 	for p in  $(EXTRA_DIRS); \

build/Makefile_bbw

@@ -54,5 +54,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	$(GG) $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libiglbbw.a

build/Makefile_boost

@@ -40,5 +40,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	g++ $(AFLAGS) $(OPENMP) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libiglboost.a

build/Makefile_cgal

@@ -41,5 +41,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libiglcgal.a

build/Makefile_embree

@@ -43,5 +43,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	g++ $(AFLAGS) $(OPENMP) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libiglembree.a

build/Makefile_matlab

@@ -42,5 +42,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libiglmatlab.a

build/Makefile_mosek

@@ -49,5 +49,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libiglmosek.a

build/Makefile_png

@@ -38,5 +38,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libiglpng.a

build/Makefile_svd3x3

@@ -48,5 +48,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	$(GG) $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libiglsvd3x3.a

build/Makefile_tetgen

@@ -40,5 +40,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libigltetgen.a

build/Makefile_triangle

@@ -40,5 +40,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libigltriangle.a

build/Makefile_viewer

@@ -46,5 +46,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	$(GG) $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libiglmatlab.a

build/Makefile_xml

@@ -53,5 +53,5 @@ obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
 	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
 
 clean:
-	rm -f obj/*.o
+	rm -f $(OBJ_FILES)
 	rm -f ../lib/libiglxml.a

include/igl/arap_rhs.cpp

@@ -10,6 +10,7 @@
 #include "verbose.h"
 #include "repdiag.h"
 #include "cat.h"
+#include <iostream>
 
 IGL_INLINE void igl::arap_rhs(
   const Eigen::MatrixXd & V, 
@@ -19,6 +20,7 @@ IGL_INLINE void igl::arap_rhs(
   Eigen::SparseMatrix<double>& K)
 {
   using namespace igl;
+  using namespace std;
   using namespace Eigen;
   // Number of dimensions
   int Vdim = V.cols();

include/igl/repdiag.cpp

@@ -6,11 +6,7 @@
 // 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 "repdiag.h"
-#ifndef IGL_NO_OPENGL
-
-// Bug in unsupported/Eigen/SparseExtra needs iostream first
-#include <iostream>
-#include <unsupported/Eigen/SparseExtra>
+#include <vector>
 
 template <typename T>
 IGL_INLINE void igl::repdiag(
@@ -18,27 +14,52 @@ IGL_INLINE void igl::repdiag(
   const int d,
   Eigen::SparseMatrix<T>& B)
 {
+  using namespace std;
+  using namespace Eigen;
   int m = A.rows();
   int n = A.cols();
 
-  B.resize(m*d,n*d);
-  // Reserve enough space for new non zeros
-  B.reserve(d*A.nonZeros());
-
-  // loop over reps
-  for(int i=0;i<d;i++)
+  vector<Triplet<T> > IJV;
+  IJV.reserve(A.nonZeros()*d);
+  // Loop outer level
+  for (int k=0; k<A.outerSize(); ++k)
   {
-    // Loop outer level
-    for (int k=0; k<A.outerSize(); ++k)
+    // loop inner level
+    for (typename Eigen::SparseMatrix<T>::InnerIterator it(A,k); it; ++it)
     {
-      // loop inner level
-      for (typename Eigen::SparseMatrix<T>::InnerIterator it(A,k); it; ++it)
+      for(int i = 0;i<d;i++)
       {
-        B.insert(i*m+it.row(),i*n+it.col()) = it.value();
+        IJV.push_back(Triplet<T>(i*m+it.row(),i*n+it.col(),it.value()));
       }
     }
   }
-  B.makeCompressed();
+  B.resize(m*d,n*d);
+  B.setFromTriplets(IJV.begin(),IJV.end());
+  
+
+  // Q: Why is this **Very** slow?
+
+  //int m = A.rows();
+  //int n = A.cols();
+
+  //B.resize(m*d,n*d);
+  //// Reserve enough space for new non zeros
+  //B.reserve(d*A.nonZeros());
+
+  //// loop over reps
+  //for(int i=0;i<d;i++)
+  //{
+  //  // Loop outer level
+  //  for (int k=0; k<A.outerSize(); ++k)
+  //  {
+  //    // loop inner level
+  //    for (typename Eigen::SparseMatrix<T>::InnerIterator it(A,k); it; ++it)
+  //    {
+  //      B.insert(i*m+it.row(),i*n+it.col()) = it.value();
+  //    }
+  //  }
+  //}
+  //B.makeCompressed();
 }
 
 template <typename T>
@@ -73,4 +94,3 @@ template void igl::repdiag<double>(Eigen::SparseMatrix<double, 0, int> const&, i
 // generated by autoexplicit.sh
 template Eigen::SparseMatrix<double, 0, int> igl::repdiag<Eigen::SparseMatrix<double, 0, int> >(Eigen::SparseMatrix<double, 0, int> const&, int);
 #endif
-#endif

include/igl/svd3x3/arap.cpp

@@ -18,6 +18,7 @@
 #include <igl/arap_rhs.h>
 #include <igl/repdiag.h>
 #include <igl/columnize.h>
+#include <igl/matlab/MatlabWorkspace.h>
 #include "fit_rotations.h"
 #include <cassert>
 #include <iostream>
@@ -95,13 +96,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(ref_V,ref_F,eff_energy,data.CSM);
-  assert(data.CSM.rows() == ref_F.rows()*data.dim);
   if(flat)
   {
     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
@@ -216,6 +215,9 @@ IGL_INLINE bool igl::arap_solve(
     assert(U.cols() == data.dim);
     // As if U.col(2) was 0
     MatrixXd S = data.CSM * Udim;
+    // THIS NORMALIZATION IS IMPORTANT TO GET SINGLE PRECISION SVD CODE TO WORK
+    // CORRECTLY.
+    S /= S.array().abs().maxCoeff();
 
     const int Rdim = data.dim;
     MatrixXd R(Rdim,data.CSM.rows());
@@ -224,11 +226,12 @@ IGL_INLINE bool igl::arap_solve(
       fit_rotations_planar(S,R);
     }else
     {
-#ifdef __SSE__ // fit_rotations_SSE will convert to float if necessary
-    fit_rotations_SSE(S,R);
-#else
-    fit_rotations(S,true,R);
-#endif
+      fit_rotations(S,true,R);
+//#ifdef __SSE__ // fit_rotations_SSE will convert to float if necessary
+//      fit_rotations_SSE(S,R);
+//#else
+//      fit_rotations(S,true,R);
+//#endif
     }
     //for(int k = 0;k<(data.CSM.rows()/dim);k++)
     //{
@@ -299,6 +302,7 @@ IGL_INLINE bool igl::arap_solve(
     // Keep track of velocity for next time
     data.vel = (U-U0)/data.h;
   }
+
   return true;
 }
 

include/igl/viewer/Viewer.cpp

@@ -940,6 +940,13 @@ namespace igl
     data.set_texture(R,G,B);
   }
 
+  void Viewer::set_points(
+    const Eigen::MatrixXd& P,
+    const Eigen::MatrixXd& C)
+  {
+    data.set_points(P,C);
+  }
+
   void Viewer::add_points(const Eigen::MatrixXd& P,  const Eigen::MatrixXd& C)
   {
     data.add_points(P,C);

include/igl/viewer/Viewer.h

@@ -93,8 +93,28 @@ namespace igl
                       const Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& G,
                       const Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& B);
 
+    // Sets points given a list of point vertices. In constrast to `set_points`
+    // this will (purposefully) clober existing points.
+    //
+    // Inputs:
+    //   P  #P by 3 list of vertex positions
+    //   C  #P|1 by 3 color(s)
+    void set_points(
+      const Eigen::MatrixXd& P,  
+      const Eigen::MatrixXd& C);
     void add_points(const Eigen::MatrixXd& P,  const Eigen::MatrixXd& C);
-    void set_edges (const Eigen::MatrixXd& P, const Eigen::MatrixXi& E, const Eigen::MatrixXd& C);
+    // Sets edges given a list of edge vertices and edge indices. In
+    // constrast
+    // to `add_edges` this will (purposefully) clober existing edges.
+    //
+    // Inputs:
+    //   P  #P by 3 list of vertex positions
+    //   E  #E by 2 list of edge indices into P
+    //   C  #E|1 by 3 color(s)
+    void set_edges(
+      const Eigen::MatrixXd& P, 
+      const Eigen::MatrixXi& E, 
+      const Eigen::MatrixXd& C);
     void add_edges (const Eigen::MatrixXd& P1, const Eigen::MatrixXd& P2, const Eigen::MatrixXd& C);
     void add_label (const Eigen::VectorXd& P,  const std::string& str);
 

include/igl/viewer/ViewerData.cpp

@@ -213,6 +213,15 @@ IGL_INLINE void igl::ViewerData::set_texture(
   dirty |= DIRTY_TEXTURE;
 }
 
+IGL_INLINE void igl::ViewerData::set_points(
+  const Eigen::MatrixXd& P,  
+  const Eigen::MatrixXd& C)
+{
+  // clear existing points
+  points.resize(0,0);
+  add_points(P,C);
+}
+
 IGL_INLINE void igl::ViewerData::add_points(const Eigen::MatrixXd& P,  const Eigen::MatrixXd& C)
 {
   Eigen::MatrixXd P_temp;

include/igl/viewer/ViewerData.h

@@ -66,6 +66,15 @@ public:
                     const Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& G,
                     const Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& B);
 
+  // Sets points given a list of point vertices. In constrast to `set_points`
+  // this will (purposefully) clober existing points.
+  //
+  // Inputs:
+  //   P  #P by 3 list of vertex positions
+  //   C  #P|1 by 3 color(s)
+  IGL_INLINE void set_points(
+    const Eigen::MatrixXd& P,  
+    const Eigen::MatrixXd& C);
   IGL_INLINE void add_points(const Eigen::MatrixXd& P,  const Eigen::MatrixXd& C);
   // Sets edges given a list of edge vertices and edge indices. In constrast
   // to `add_edges` this will (purposefully) clober existing edges.

tutorial/406_FastAutomaticSkinningTransformations/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(406_FastAutomaticSkinningTransformations)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${PROJECT_NAME}_bin ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES})

tutorial/406_FastAutomaticSkinningTransformations/main.cpp

@@ -0,0 +1,216 @@
+#include <igl/colon.h>
+#include <igl/directed_edge_orientations.h>
+#include <igl/directed_edge_parents.h>
+#include <igl/forward_kinematics.h>
+#include <igl/PI.h>
+#include <igl/partition.h>
+#include <igl/mat_max.h>
+#include <igl/lbs_matrix.h>
+#include <igl/slice.h>
+#include <igl/deform_skeleton.h>
+#include <igl/dqs.h>
+#include <igl/lbs_matrix.h>
+#include <igl/columnize.h>
+#include <igl/readDMAT.h>
+#include <igl/readOBJ.h>
+#include <igl/svd3x3/arap.h>
+#include <igl/svd3x3/arap_dof.h>
+#include <igl/viewer/Viewer.h>
+
+#include <Eigen/Geometry>
+#include <Eigen/StdVector>
+#include <vector>
+#include <algorithm>
+#include <iostream>
+
+typedef 
+  std::vector<Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> >
+  RotationList;
+
+const Eigen::RowVector3d sea_green(70./255.,252./255.,167./255.);
+Eigen::MatrixXd V,U,M;
+Eigen::MatrixXi F;
+Eigen::VectorXi S,b;
+Eigen::MatrixXd L;
+Eigen::RowVector3d mid;
+double anim_t = 0.0;
+double anim_t_dir = 0.03;
+double bbd = 1.0;
+bool resolve = true;
+igl::ARAPData arap_data,arap_grouped_data;
+igl::ArapDOFData<Eigen::MatrixXd,double> arap_dof_data;
+Eigen::SparseMatrix<double> Aeq;
+
+enum ModeType
+{
+  MODE_TYPE_ARAP = 0,
+  MODE_TYPE_ARAP_GROUPED = 1,
+  MODE_TYPE_ARAP_DOF = 2,
+  NUM_MODE_TYPES = 4
+} mode = MODE_TYPE_ARAP;
+
+bool pre_draw(igl::Viewer & viewer)
+{
+  using namespace Eigen;
+  using namespace std;
+  if(resolve)
+  {
+    MatrixXd bc(b.size(),V.cols());
+    VectorXd Beq(3*b.size());
+    for(int i = 0;i<b.size();i++)
+    {
+      bc.row(i) = V.row(b(i));
+      switch(i%4)
+      {
+        case 2:
+          bc(i,0) += 0.15*bbd*sin(0.5*anim_t);
+          bc(i,1) += 0.15*bbd*(1.-cos(0.5*anim_t));
+          break;
+        case 1:
+          bc(i,1) += 0.10*bbd*sin(1.*anim_t*(i+1));
+          bc(i,2) += 0.10*bbd*(1.-cos(1.*anim_t*(i+1)));
+          break;
+        case 0:
+          bc(i,0) += 0.20*bbd*sin(2.*anim_t*(i+1));
+          break;
+      }
+      Beq(3*i+0) = bc(i,0);
+      Beq(3*i+1) = bc(i,1);
+      Beq(3*i+2) = bc(i,2);
+    }
+    switch(mode)
+    {
+      case MODE_TYPE_ARAP:
+        igl::arap_solve(bc,arap_data,U);
+        break;
+      case MODE_TYPE_ARAP_GROUPED:
+        igl::arap_solve(bc,arap_grouped_data,U);
+        break;
+      case MODE_TYPE_ARAP_DOF:
+      {
+        VectorXd L0 = L;
+        arap_dof_update(arap_dof_data,Beq,L0,30,0,L);
+        const auto & Ucol = M*L;
+        U.col(0) = Ucol.block(0*U.rows(),0,U.rows(),1);
+        U.col(1) = Ucol.block(1*U.rows(),0,U.rows(),1);
+        U.col(2) = Ucol.block(2*U.rows(),0,U.rows(),1);
+        break;
+      }
+    }
+    viewer.set_vertices(U);
+    viewer.set_points(bc,sea_green);
+    viewer.compute_normals();
+    if(viewer.core.is_animating)
+    {
+      anim_t += anim_t_dir;
+    }else
+    {
+      resolve = false;
+    }
+  }
+  return false;
+}
+
+bool key_down(igl::Viewer &viewer, unsigned char key, int mods)
+{
+  switch(key)
+  {
+    case '0':
+      anim_t = 0;
+      resolve = true;
+      return true;
+    case '.':
+      mode = (ModeType)(((int)mode+1)%((int)NUM_MODE_TYPES-1));
+      resolve = true;
+      return true;
+    case ',':
+      mode = (ModeType)(((int)mode-1)%((int)NUM_MODE_TYPES-1));
+      resolve = true;
+      return true;
+    case ' ':
+      viewer.core.is_animating = !viewer.core.is_animating;
+      if(viewer.core.is_animating)
+      {
+        resolve = true;
+      }
+      return true;
+  }
+  return false;
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  igl::readOBJ("../shared/armadillo.obj",V,F);
+  U=V;
+  MatrixXd W;
+  igl::readDMAT("../shared/armadillo-weights.dmat",W);
+  igl::lbs_matrix_column(V,W,M);
+  VectorXi G;
+  {
+    VectorXi S;
+    VectorXd D;
+    igl::partition(W,50,G,S,D);
+  }
+
+  // vertices in selection (in regions with weight = delta)
+  //Matrix<bool,Dynamic,1> S = W.array().rowwise().maxCoeff()==1;
+  //igl::colon<int>(0,V.rows()-1,b);
+  //b.conservativeResize(stable_partition( b.data(), b.data()+b.size(), 
+  // [&](int i)->bool{return S(i);})-b.data());
+
+  {
+    VectorXd maxW;
+    igl::mat_max(W,1,maxW,b);
+  }
+
+  // Precomputation for FAST
+  // number of weights
+  const int m = W.cols();
+  Aeq.resize(m*3,m*3*(3+1));
+  vector<Triplet<double> > ijv;
+  for(int i = 0;i<m;i++)
+  {
+    RowVector4d homo;
+    homo << V.row(b(i)),1.;
+    for(int d = 0;d<3;d++)
+    {
+      for(int c = 0;c<(3+1);c++)
+      {
+        ijv.push_back(Triplet<double>(3*i + d,i + c*m*3 + d*m, homo(c)));
+      }
+    }
+  }
+  Aeq.setFromTriplets(ijv.begin(),ijv.end());
+  igl::arap_dof_precomputation(V,F,M,G,arap_dof_data);
+  igl::arap_dof_recomputation(VectorXi(),Aeq,arap_dof_data);
+  // Initialize
+  MatrixXd Istack = MatrixXd::Identity(3,3+1).replicate(1,m);
+  igl::columnize(Istack,m,2,L);
+
+  // Precomputation for ARAP
+  arap_data.max_iter = 1;
+  igl::arap_precomputation(V,F,V.cols(),b,arap_data);
+  // Grouped arap
+  arap_grouped_data.max_iter = 2;
+  arap_grouped_data.G = G;
+  igl::arap_precomputation(V,F,V.cols(),b,arap_grouped_data);
+
+
+  // bounding box diagonal
+  bbd = (V.colwise().maxCoeff()- V.colwise().minCoeff()).norm();
+
+  // Plot the mesh with pseudocolors
+  igl::Viewer viewer;
+  viewer.set_mesh(U, F);
+  viewer.add_points(igl::slice(V,b,1),sea_green);
+  viewer.core.show_lines = false;
+  viewer.callback_pre_draw = &pre_draw;
+  viewer.callback_key_down = &key_down;
+  viewer.core.is_animating = false;
+  viewer.core.animation_max_fps = 30.;
+  cout<<
+    "Press [space] to toggle animation"<<endl;
+  viewer.launch();
+}

tutorial/images/armadillo-fast.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+15c64de24bc8c6f038afc3975355364907f1a8b6

tutorial/shared/armadillo-weights.dmat.REMOVED.git-id

@@ -0,0 +1 @@
+29ec0f5e3f445b96d480d3e53d712946b459caf1

tutorial/shared/armadillo.obj.REMOVED.git-id

@@ -0,0 +1 @@
+b7df6d6c21f2d8002037f48d4156d57d1dcf0148

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-d373d835336654fc66f128ce06a64f6d043f1e87
+d87566b894b52e6aa2ad61a4d18afe555186f7da