Make usample work like loop subdivision: returning sparse matrix

Former-commit-id: ad4de0e50cbe4500383c956e6c4d0588384eb53d

include/igl/upsample.cpp

@@ -11,15 +11,13 @@
 
 
 template <
-  typename DerivedV,
   typename DerivedF,
-  typename DerivedNV,
+  typename DerivedS,
   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)
+IGL_INLINE void igl::upsample(const int n_verts,
+                              const Eigen::PlainObjectBase<DerivedF>& F,
+                              Eigen::SparseMatrix<DerivedS>& S,
+                              Eigen::PlainObjectBase<DerivedNF>& NF)
 {
   // Use "in place" wrapper instead
   assert(&V != &NV);
@@ -27,15 +25,18 @@ IGL_INLINE void igl::upsample(
   using namespace std;
   using namespace Eigen;
 
+  typedef Eigen::Triplet<DerivedS> Triplet_t;
+
   Eigen::Matrix<
-    typename DerivedF::Scalar,Eigen::Dynamic,Eigen::Dynamic>
-    FF,FFi;
+  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
 
   // 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());
   int counter = 0;
 
   for(int i=0;i<FF.rows();++i)
@@ -45,33 +46,45 @@ IGL_INLINE void igl::upsample(
       if(NI(i,j) == -1)
       {
         NI(i,j) = counter;
-        if (FF(i,j) != -1) // If it is not a border
-          NI(FF(i,j),FFi(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;
+        }
         ++counter;
       }
     }
   }
 
-  int n_odd = V.rows();
-  int n_even = counter;
+  const int& n_odd = n_verts;
+  const int& n_even = counter;
+  const int n_newverts = n_odd + n_even;
 
-  // Preallocate NV and NF
-  NV.resize(V.rows()+n_even,V.cols());
-  NF.resize(F.rows()*4,3);
+  //Construct vertex positions
+  std::vector<Triplet_t> tripletList;
 
   // Fill the odd vertices position
-  NV.block(0,0,V.rows(),V.cols()) = V;
+  for (int i=0; i<n_odd; ++i)
+  {
+    tripletList.emplace_back(i, i, 1.);
+  }
 
-  // Fill the even vertices position
   for(int i=0;i<FF.rows();++i)
   {
     for(int j=0;j<3;++j)
     {
-      NV.row(NI(i,j) + n_odd) = 0.5 * V.row(F(i,j)) + 0.5 * V.row(F(i,(j+1)%3));
+      if(NIdoubles(i,j)==0) {
+        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.);
+      }
     }
   }
+  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)
   {
     VectorXi VI(6);
@@ -88,19 +101,43 @@ IGL_INLINE void igl::upsample(
     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::upsample(
+                              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) {
+    DerivedNF tempF = NF;
+    SparseMat S;
+    upsample(NV.rows(), tempF, S, NF);
+    NV = S*NV;
+  }
 }
 
 template <
   typename MatV,
   typename MatF>
-IGL_INLINE void igl::upsample(
-  MatV& V,
-  MatF& F)
+IGL_INLINE void igl::upsample(MatV& V,
+                              MatF& F,
+                              const int number_of_subdivs)
 {
   const MatV V_copy = V;
   const MatF F_copy = F;
-  return upsample(V_copy,F_copy,V,F);
+  return upsample(V_copy,F_copy,V,F,number_of_subdivs);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/upsample.h

@@ -10,11 +10,31 @@
 #include "igl_inline.h"
 
 #include <Eigen/Core>
+#include <Eigen/Sparse>
 
 // History:
 //  changed templates from generic matrices to PlainObjectBase Alec May 7, 2011
 namespace igl
 {
+  // Subdivide without moving vertices: 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 DerivedS,
+    typename DerivedNF>
+  IGL_INLINE void upsample(const int n_verts,
+                           const Eigen::PlainObjectBase<DerivedF>& F,
+                           Eigen::SparseMatrix<DerivedS>& 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
   // 
@@ -38,18 +58,19 @@ 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);
+  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);
+  IGL_INLINE void upsample(MatV& V,
+                           MatF& F,
+                           const int number_of_subdivs = 1);
 }
 
 #ifndef IGL_STATIC_LIBRARY