loop, upsample, bug fix + tutorial

Former-commit-id: 8346cf35adb87eaaf6a394eff262245240ac811b

include/igl/loop.cpp

@@ -14,139 +14,156 @@
 
 #include <vector>
 
-namespace igl
+template <
+  typename DerivedF,
+  typename SType,
+  typename DerivedNF>
+IGL_INLINE void igl::loop(
+  const int n_verts,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::SparseMatrix<SType>& S,
+  Eigen::PlainObjectBase<DerivedNF> & NF)
 {
-    
-    IGL_INLINE void loop(const int n_verts,
-                         const Eigen::MatrixXi& F,
-                         Eigen::SparseMatrix<double>& S,
-                         Eigen::MatrixXi& newF)
+  typedef Eigen::SparseMatrix<SType> SparseMat;
+  typedef Eigen::Triplet<SType> Triplet_t;
+  
+  //Ref. https://graphics.stanford.edu/~mdfisher/subdivision.html
+  //Heavily borrowing from igl::upsample
+  
+  DerivedF FF, FFi;
+  triangle_triangle_adjacency(F, FF, FFi);
+  std::vector<std::vector<typename DerivedF::Scalar>> adjacencyList;
+  adjacency_list(F, adjacencyList, true);
+  
+  //Compute the number and positions of the vertices to insert (on edges)
+  Eigen::MatrixXi NI = Eigen::MatrixXi::Constant(FF.rows(), FF.cols(), -1);
+  Eigen::MatrixXi NIdoubles = Eigen::MatrixXi::Zero(FF.rows(), FF.cols());
+  Eigen::VectorXi vertIsOnBdry = Eigen::VectorXi::Zero(n_verts);
+  int counter = 0;
+  for(int i=0; i<FF.rows(); ++i)
+  {
+    for(int j=0; j<3; ++j)
     {
-        
-        typedef Eigen::SparseMatrix<double> SparseMat;
-        typedef Eigen::Triplet<double> Triplet_t;
-        
-        //Ref. https://graphics.stanford.edu/~mdfisher/subdivision.html
-        //Heavily borrowing from igl::upsample
-        
-        Eigen::MatrixXi FF, FFi;
-        triangle_triangle_adjacency(F, FF, FFi);
-        std::vector<std::vector<int>> adjacencyList;
-        adjacency_list(F, adjacencyList, true);
-        
-        //Compute the number and positions of the vertices to insert (on edges)
-        Eigen::MatrixXi NI = Eigen::MatrixXi::Constant(FF.rows(), FF.cols(), -1);
-        Eigen::MatrixXi NIdoubles = Eigen::MatrixXi::Zero(FF.rows(), FF.cols());
-        Eigen::VectorXi vertIsOnBdry = Eigen::VectorXi::Zero(n_verts);
-        int counter = 0;
-        for(int i=0; i<FF.rows(); ++i)
+      if(NI(i,j) == -1)
+      {
+        NI(i,j) = counter;
+        NIdoubles(i,j) = 0;
+        if (FF(i,j) != -1) 
         {
-            for(int j=0; j<3; ++j)
-            {
-                if(NI(i,j) == -1)
-                {
-                    NI(i,j) = counter;
-                    NIdoubles(i,j) = 0;
-                    if (FF(i,j) != -1) {
-                        //If it is not a boundary
-                        NI(FF(i,j), FFi(i,j)) = counter;
-                        NIdoubles(i,j) = 1;
-                    } else {
-                        //Mark boundary vertices for later
-                        vertIsOnBdry(F(i,j)) = 1;
-                        vertIsOnBdry(F(i,(j+1)%3)) = 1;
-                    }
-                    ++counter;
-                }
-            }
-        }
-        
-        const int& n_odd = n_verts;
-        const int& n_even = counter;
-        const int n_newverts = n_odd + n_even;
-        
-        //Construct vertex positions
-        std::vector<Triplet_t> tripletList;
-        for(int i=0; i<n_odd; ++i) {
-            //Old vertices
-            const std::vector<int>& localAdjList = adjacencyList[i];
-            if(vertIsOnBdry(i)==1) {
-                //Boundary vertex
-                tripletList.emplace_back(i, localAdjList.front(), 1./8.);
-                tripletList.emplace_back(i, localAdjList.back(), 1./8.);
-                tripletList.emplace_back(i, i, 3./4.);
-            } else {
-                const int n = localAdjList.size();
-                const double dn = n;
-                double beta;
-                if(n==3)
-                    beta = 3./16.;
-                else
-                    beta = 3./8./dn;
-                for(int j=0; j<n; ++j)
-                    tripletList.emplace_back(i, localAdjList[j], beta);
-                tripletList.emplace_back(i, i, 1.-dn*beta);
-            }
-        }
-        for(int i=0; i<FF.rows(); ++i) {
-            //New vertices
-            for(int j=0; j<3; ++j) {
-                if(NIdoubles(i,j)==0) {
-                    if(FF(i,j)==-1) {
-                        //Boundary vertex
-                        tripletList.emplace_back(NI(i,j) + n_odd, F(i,j), 1./2.);
-                        tripletList.emplace_back(NI(i,j) + n_odd, F(i, (j+1)%3), 1./2.);
-                    } else {
-                        tripletList.emplace_back(NI(i,j) + n_odd, F(i,j), 3./8.);
-                        tripletList.emplace_back(NI(i,j) + n_odd, F(i, (j+1)%3), 3./8.);
-                        tripletList.emplace_back(NI(i,j) + n_odd, F(i, (j+2)%3), 1./8.);
-                        tripletList.emplace_back(NI(i,j) + n_odd, F(FF(i,j), (FFi(i,j)+2)%3), 1./8.);
-                    }
-                }
-            }
-        }
-        S.resize(n_newverts, n_verts);
-        S.setFromTriplets(tripletList.begin(), tripletList.end());
-        
-        // Build the new topology (Every face is replaced by four)
-        newF.resize(F.rows()*4, 3);
-        for(int i=0; i<F.rows();++i)
+          //If it is not a boundary
+          NI(FF(i,j), FFi(i,j)) = counter;
+          NIdoubles(i,j) = 1;
+        } else 
         {
-            Eigen::VectorXi VI(6);
-            VI << F(i,0), F(i,1), F(i,2), NI(i,0) + n_odd, NI(i,1) + n_odd, NI(i,2) + n_odd;
-            
-            Eigen::VectorXi f0(3), f1(3), f2(3), f3(3);
-            f0 << VI(0), VI(3), VI(5);
-            f1 << VI(1), VI(4), VI(3);
-            f2 << VI(3), VI(4), VI(5);
-            f3 << VI(4), VI(2), VI(5);
-            
-            newF.row((i*4)+0) = f0;
-            newF.row((i*4)+1) = f1;
-            newF.row((i*4)+2) = f2;
-            newF.row((i*4)+3) = f3;
+          //Mark boundary vertices for later
+          vertIsOnBdry(F(i,j)) = 1;
+          vertIsOnBdry(F(i,(j+1)%3)) = 1;
         }
-        
+        ++counter;
+      }
     }
-    
-    
-    IGL_INLINE void loop(const Eigen::MatrixXd& V,
-                         const Eigen::MatrixXi& F,
-                         Eigen::MatrixXd& newV,
-                         Eigen::MatrixXi& newF,
-                         const int number_of_subdivs)
+  }
+  
+  const int& n_odd = n_verts;
+  const int& n_even = counter;
+  const int n_newverts = n_odd + n_even;
+  
+  //Construct vertex positions
+  std::vector<Triplet_t> tripletList;
+  for(int i=0; i<n_odd; ++i) 
+  {
+    //Old vertices
+    const std::vector<int>& localAdjList = adjacencyList[i];
+    if(vertIsOnBdry(i)==1) 
     {
-        typedef Eigen::SparseMatrix<double> SparseMat;
-        typedef Eigen::Triplet<double> Triplet_t;
-        
-        newV = V;
-        newF = F;
-        for(int i=0; i<number_of_subdivs; ++i) {
-            Eigen::MatrixXi tempF = newF;
-            SparseMat S;
-            loop(newV.rows(), tempF, S, newF);
-            newV = S*newV;
+      //Boundary vertex
+      tripletList.emplace_back(i, localAdjList.front(), 1./8.);
+      tripletList.emplace_back(i, localAdjList.back(), 1./8.);
+      tripletList.emplace_back(i, i, 3./4.);
+    } else 
+    {
+      const int n = localAdjList.size();
+      const SType dn = n;
+      SType beta;
+      if(n==3)
+      {
+        beta = 3./16.;
+      } else
+      {
+        beta = 3./8./dn;
+      }
+      for(int j=0; j<n; ++j)
+      {
+        tripletList.emplace_back(i, localAdjList[j], beta);
+      }
+      tripletList.emplace_back(i, i, 1.-dn*beta);
+    }
+  }
+  for(int i=0; i<FF.rows(); ++i) 
+  {
+    //New vertices
+    for(int j=0; j<3; ++j) 
+    {
+      if(NIdoubles(i,j)==0) 
+      {
+        if(FF(i,j)==-1) 
+        {
+          //Boundary vertex
+          tripletList.emplace_back(NI(i,j) + n_odd, F(i,j), 1./2.);
+          tripletList.emplace_back(NI(i,j) + n_odd, F(i, (j+1)%3), 1./2.);
+        } else 
+        {
+          tripletList.emplace_back(NI(i,j) + n_odd, F(i,j), 3./8.);
+          tripletList.emplace_back(NI(i,j) + n_odd, F(i, (j+1)%3), 3./8.);
+          tripletList.emplace_back(NI(i,j) + n_odd, F(i, (j+2)%3), 1./8.);
+          tripletList.emplace_back(NI(i,j) + n_odd, F(FF(i,j), (FFi(i,j)+2)%3), 1./8.);
         }
+      }
     }
+  }
+  S.resize(n_newverts, n_verts);
+  S.setFromTriplets(tripletList.begin(), tripletList.end());
+  
+  // Build the new topology (Every face is replaced by four)
+  NF.resize(F.rows()*4, 3);
+  for(int i=0; i<F.rows();++i)
+  {
+    Eigen::VectorXi VI(6);
+    VI << F(i,0), F(i,1), F(i,2), NI(i,0) + n_odd, NI(i,1) + n_odd, NI(i,2) + n_odd;
+    
+    Eigen::VectorXi f0(3), f1(3), f2(3), f3(3);
+    f0 << VI(0), VI(3), VI(5);
+    f1 << VI(1), VI(4), VI(3);
+    f2 << VI(3), VI(4), VI(5);
+    f3 << VI(4), VI(2), VI(5);
     
+    NF.row((i*4)+0) = f0;
+    NF.row((i*4)+1) = f1;
+    NF.row((i*4)+2) = f2;
+    NF.row((i*4)+3) = f3;
+  }
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF,
+  typename DerivedNV,
+  typename DerivedNF>
+IGL_INLINE void igl::loop(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedNV>& NV,
+  Eigen::PlainObjectBase<DerivedNF>& NF,
+  const int number_of_subdivs)
+{
+  NV = V;
+  NF = F;
+  for(int i=0; i<number_of_subdivs; ++i) 
+  {
+    DerivedNF tempF = NF;
+    Eigen::SparseMatrix<typename DerivedV::Scalar> S;
+    loop(NV.rows(), tempF, S, NF);
+    // This .eval is super important
+    NV = (S*NV).eval();
+  }
 }

include/igl/loop.h

@@ -15,38 +15,49 @@
 
 namespace igl
 {
-    // LOOP Given the triangle mesh [V, F], where n_verts = V.rows(), computes newV and a sparse matrix S s.t. [newV, newF] is the subdivided mesh where newV = S*V.
-    //
-    // Inputs:
-    //  n_verts an integer (number of mesh vertices)
-    //  F an m by 3 matrix of integers of triangle faces
-    // Outputs:
-    //  S a sparse matrix (will become the subdivision matrix)
-    //  newF a matrix containing the new faces
-    IGL_INLINE void loop(const int n_verts,
-                         const Eigen::MatrixXi& F,
-                         Eigen::SparseMatrix<double>& S,
-                         Eigen::MatrixXi& newF);
-    
-    // LOOP Given the triangle mesh [V, F], computes number_of_subdivs steps of loop subdivision and outputs the new mesh [newV, newF]
-    //
-    // Inputs:
-    //  V an n by 3 matrix of vertices
-    //  F an m by 3 matrix of integers of triangle faces
-    //  number_of_subdivs an integer that specifies how many subdivision steps to do
-    // Outputs:
-    //  newV a matrix containing the new vertices
-    //  newF a matrix containing the new faces
-    IGL_INLINE void loop(const Eigen::MatrixXd& V,
-                         const Eigen::MatrixXi& F,
-                         Eigen::MatrixXd& newV,
-                         Eigen::MatrixXi& newF,
-                         const int number_of_subdivs = 1);
-    
+  // LOOP Given the triangle mesh [V, F], where n_verts = V.rows(), computes
+  // newV and a sparse matrix S s.t. [newV, newF] is the subdivided mesh where
+  // newV = S*V.
+  //
+  // Inputs:
+  //   n_verts  an integer (number of mesh vertices)
+  //   F  an m by 3 matrix of integers of triangle faces
+  // Outputs:
+  //   S  a sparse matrix (will become the subdivision matrix)
+  //   newF  a matrix containing the new faces
+  template <
+    typename DerivedF,
+    typename SType,
+    typename DerivedNF>
+  IGL_INLINE void loop(
+    const int n_verts,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::SparseMatrix<SType>& S,
+    Eigen::PlainObjectBase<DerivedNF> & NF);
+  // LOOP Given the triangle mesh [V, F], computes number_of_subdivs steps of loop subdivision and outputs the new mesh [newV, newF]
+  //
+  // Inputs:
+  //  V an n by 3 matrix of vertices
+  //  F an m by 3 matrix of integers of triangle faces
+  //  number_of_subdivs an integer that specifies how many subdivision steps to do
+  // Outputs:
+  //  NV a matrix containing the new vertices
+  //  NF a matrix containing the new faces
+  template <
+    typename DerivedV, 
+    typename DerivedF,
+    typename DerivedNV,
+    typename DerivedNF>
+  IGL_INLINE void loop(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedNV>& NV,
+    Eigen::PlainObjectBase<DerivedNF>& NF,
+    const int number_of_subdivs = 1);
 }
 
 #ifndef IGL_STATIC_LIBRARY
-#include "loop.cpp"
+#  include "loop.cpp"
 #endif
 
 #endif

include/igl/upsample.cpp

@@ -12,21 +12,21 @@
 
 template <
   typename DerivedF,
-  typename DerivedS,
+  typename SType,
   typename DerivedNF>
-IGL_INLINE void igl::upsample(const int n_verts,
-                              const Eigen::PlainObjectBase<DerivedF>& F,
-                              Eigen::SparseMatrix<DerivedS>& S,
-                              Eigen::PlainObjectBase<DerivedNF>& NF)
+IGL_INLINE void igl::upsample(
+  const int n_verts,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::SparseMatrix<SType>& S,
+  Eigen::PlainObjectBase<DerivedNF>& NF)
 {
   using namespace std;
   using namespace Eigen;
 
-  typedef Eigen::Triplet<DerivedS> Triplet_t;
+  typedef Eigen::Triplet<SType> Triplet_t;
 
-  Eigen::Matrix<
-  typename DerivedF::Scalar,Eigen::Dynamic,Eigen::Dynamic>
-  FF,FFi;
+  Eigen::Matrix< typename DerivedF::Scalar,Eigen::Dynamic,Eigen::Dynamic>
+    FF,FFi;
   triangle_triangle_adjacency(F,FF,FFi);
 
   // TODO: Cache optimization missing from here, it is a mess
@@ -106,31 +106,31 @@ template <
   typename DerivedNV,
   typename DerivedNF>
 IGL_INLINE void igl::upsample(
-                              const Eigen::PlainObjectBase<DerivedV>& V,
-                              const Eigen::PlainObjectBase<DerivedF>& F,
-                              Eigen::PlainObjectBase<DerivedNV>& NV,
-                              Eigen::PlainObjectBase<DerivedNF>& NF,
-                              const int number_of_subdivs)
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedNV>& NV,
+  Eigen::PlainObjectBase<DerivedNF>& NF,
+  const int number_of_subdivs)
 {
-  typedef Eigen::SparseMatrix<double> SparseMat;
-  typedef Eigen::Triplet<double> Triplet_t;
-
   NV = V;
   NF = F;
-  for(int i=0; i<number_of_subdivs; ++i) {
+  for(int i=0; i<number_of_subdivs; ++i) 
+  {
     DerivedNF tempF = NF;
-    SparseMat S;
+    Eigen::SparseMatrix<typename DerivedV::Scalar >S;
     upsample(NV.rows(), tempF, S, NF);
-    NV = S*NV;
+    // This .eval is super important
+    NV = (S*NV).eval();
   }
 }
 
 template <
   typename MatV,
   typename MatF>
-IGL_INLINE void igl::upsample(MatV& V,
-                              MatF& F,
-                              const int number_of_subdivs)
+IGL_INLINE void igl::upsample(
+  MatV& V,
+  MatF& F,
+  const int number_of_subdivs)
 {
   const MatV V_copy = V;
   const MatF F_copy = F;

include/igl/upsample.h

@@ -21,20 +21,20 @@ namespace igl
   // [newV, newF] is the subdivided mesh where newV = S*V.
   //
   // Inputs:
-  //  n_verts an integer (number of mesh vertices)
-  //  F an m by 3 matrix of integers of triangle faces
+  //   n_verts  an integer (number of mesh vertices)
+  //   F  an m by 3 matrix of integers of triangle faces
   // Outputs:
-  //  S a sparse matrix (will become the subdivision matrix)
-  //  newF a matrix containing the new faces
+  //   S  a sparse matrix (will become the subdivision matrix)
+  //   newF  a matrix containing the new faces
   template <
     typename DerivedF,
-    typename DerivedS,
+    typename SType,
     typename DerivedNF>
-  IGL_INLINE void upsample(const int n_verts,
-                           const Eigen::PlainObjectBase<DerivedF>& F,
-                           Eigen::SparseMatrix<DerivedS>& S,
-                           Eigen::PlainObjectBase<DerivedNF>& NF);
-
+  IGL_INLINE void upsample(
+    const int n_verts,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::SparseMatrix<SType>& S,
+    Eigen::PlainObjectBase<DerivedNF>& NF);
   // Subdivide a mesh without moving vertices: loop subdivision but odd
   // vertices stay put and even vertices are just edge midpoints
   // 
@@ -58,19 +58,21 @@ namespace igl
     typename DerivedF,
     typename DerivedNV,
     typename DerivedNF>
-  IGL_INLINE void upsample(const Eigen::PlainObjectBase<DerivedV>& V,
-                           const Eigen::PlainObjectBase<DerivedF>& F,
-                           Eigen::PlainObjectBase<DerivedNV>& NV,
-                           Eigen::PlainObjectBase<DerivedNF>& NF,
-                           const int number_of_subdivs = 1);
+  IGL_INLINE void upsample(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedNV>& NV,
+    Eigen::PlainObjectBase<DerivedNF>& NF,
+    const int number_of_subdivs = 1);
 
   // Virtually in place wrapper
   template <
     typename MatV, 
     typename MatF>
-  IGL_INLINE void upsample(MatV& V,
-                           MatF& F,
-                           const int number_of_subdivs = 1);
+  IGL_INLINE void upsample(
+    MatV& V,
+    MatF& F,
+    const int number_of_subdivs = 1);
 }
 
 #ifndef IGL_STATIC_LIBRARY

tutorial/710_SLIM/CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 2.8.12)
-project(611_SLIM)
+project(710_SLIM)
 
 add_executable(${PROJECT_NAME}_bin
   main.cpp)

tutorial/711_Subdivision/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(711_Subdivision)
+
+add_executable(${PROJECT_NAME}_bin
+  main.cpp)
+target_include_directories(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_INCLUDE_DIRS})
+target_compile_definitions(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_DEFINITIONS})
+target_link_libraries(${PROJECT_NAME}_bin ${LIBIGL_LIBRARIES} ${LIBIGL_VIEWER_EXTRA_LIBRARIES})

tutorial/711_Subdivision/main.cpp

@@ -0,0 +1,69 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/loop.h>
+#include <igl/upsample.h>
+#include <igl/false_barycentric_subdivision.h>
+#include <igl/viewer/Viewer.h>
+#include <Eigen/Core>
+#include <iostream>
+
+#include "tutorial_shared_path.h"
+
+int main(int argc, char * argv[])
+{
+  using namespace std;
+  using namespace igl;
+  Eigen::MatrixXi F,UF,LF,FF;
+  Eigen::MatrixXd V,UV,LV,FV;
+  bool show_swept_volume = false;
+  read_triangle_mesh(
+      TUTORIAL_SHARED_PATH "/decimated-knight.off",V,F);
+  cout<<R"(Usage:
+1  Original mesh
+2  In-plane upsampled mesh
+3  Loop subdivided mesh
+4  False barycentric subdivision
+)";
+  igl::viewer::Viewer viewer;
+  viewer.data.set_mesh(V,F);
+  viewer.data.set_face_based(true);
+  igl::upsample(V,F,UV,UF);
+  igl::loop(V,F,LV,LF);
+  igl::false_barycentric_subdivision(V,F,FV,FF);
+
+  viewer.callback_key_down =
+    [&](igl::viewer::Viewer & viewer, unsigned char key, int mod)->bool
+    {
+      switch(key)
+      {
+        default:
+          return false;
+        case '1':
+        {
+          viewer.data.clear();
+          viewer.data.set_mesh(V,F);
+          break;
+        }
+        case '2':
+        {
+          viewer.data.clear();
+          viewer.data.set_mesh(UV,UF);
+          break;
+        }
+        case '3':
+        {
+          viewer.data.clear();
+          viewer.data.set_mesh(LV,LF);
+          break;
+        }
+        case '4':
+        {
+          viewer.data.clear();
+          viewer.data.set_mesh(FV,FF);
+          break;
+        }
+      }
+      viewer.data.set_face_based(true);
+      return true;
+    };
+  viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -164,7 +164,6 @@ if(TUTORIALS_CHAPTER6)
     add_subdirectory("609_Boolean")
     add_subdirectory("610_CSGTree")
   endif()
-  add_subdirectory("611_SLIM")
 endif()
 
 # Chapter 7
@@ -180,5 +179,6 @@ if(TUTORIALS_CHAPTER7)
   add_subdirectory("707_SweptVolume")
   add_subdirectory("708_Picking")
   add_subdirectory("709_VectorFieldVisualizer")
-
+  add_subdirectory("710_SLIM")
+  add_subdirectory("711_Subdivision")
 endif()