Added loop subdivision functionality.

Former-commit-id: 9d542d7ce3626ac3b2a1581465e4e14c35216ef6

include/igl/loop.cpp

@@ -0,0 +1,152 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Oded Stein <oded.stein@columbia.edu>
+//
+// 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 "loop.h"
+
+#include <igl/adjacency_list.h>
+#include <igl/triangle_triangle_adjacency.h>
+#include <igl/unique.h>
+
+#include <vector>
+
+namespace igl
+{
+    
+    IGL_INLINE void loop(const int n_verts,
+                         const Eigen::MatrixXi& F,
+                         Eigen::SparseMatrix<double>& S,
+                         Eigen::MatrixXi& newF)
+    {
+        
+        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)
+        {
+            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)
+        {
+            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;
+        }
+        
+    }
+    
+    
+    IGL_INLINE void loop(const Eigen::MatrixXd& V,
+                         const Eigen::MatrixXi& F,
+                         Eigen::MatrixXd& newV,
+                         Eigen::MatrixXi& newF,
+                         const int number_of_subdivs)
+    {
+        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;
+        }
+    }
+    
+}

include/igl/loop.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Oded Stein <oded.stein@columbia.edu>
+//
+// 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_LOOP_H
+#define IGL_LOOP_H
+
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+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);
+    
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "loop.cpp"
+#endif
+
+#endif