polygon laplacian and other polygon functionality

include/igl/adjacency_matrix.cpp

@@ -65,8 +65,57 @@ IGL_INLINE void igl::adjacency_matrix(
   }
 }
 
+template <typename DerivedI, typename DerivedC, typename T>
+IGL_INLINE void igl::adjacency_matrix(
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedC> & C,
+  Eigen::SparseMatrix<T>& A)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  typedef Triplet<T> IJV;
+  vector<IJV > ijv;
+  ijv.reserve(C(C.size()-1)*2);
+  const typename DerivedI::Scalar n = I.maxCoeff()+1;
+  {
+    // loop over polygons
+    for(int p = 0;p<C.size()-1;p++)
+    {
+      // number of edges
+      const int np = C(p+1)-C(p);
+      // loop over edges
+      for(int c = 0;c<np;c++)
+      {
+        const int i = I(C(p)+c);
+        const int j = I(C(p)+((c+1)%np));
+        ijv.emplace_back(i,j,1);
+        ijv.emplace_back(j,i,1);
+      }
+    }
+  }
+
+  A.resize(n,n);
+  A.reserve(6*n);
+  A.setFromTriplets(ijv.begin(),ijv.end());
+
+  // Force all non-zeros to be one
+
+  // Iterate over outside
+  for(int k=0; k<A.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
+    {
+      assert(it.value() != 0);
+      A.coeffRef(it.row(),it.col()) = 1;
+    }
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template void igl::adjacency_matrix<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<int, 0, int>& );
 // generated by autoexplicit.sh
 template void igl::adjacency_matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, bool>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<bool, 0, int>&);
 template void igl::adjacency_matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);

include/igl/adjacency_matrix.h

@@ -21,7 +21,7 @@ namespace igl
   // Inputs:
   //   F  #F by dim list of mesh simplices
   // Outputs: 
-  //   A  max(F) by max(F) adjacency matrix, each row i corresponding to V(i,:)
+  //   A  max(F)+1 by max(F)+1 adjacency matrix, each row i corresponding to V(i,:)
   //
   // Example:
   //   // Mesh in (V,F)
@@ -42,6 +42,21 @@ namespace igl
   IGL_INLINE void adjacency_matrix(
     const Eigen::MatrixBase<DerivedF> & F, 
     Eigen::SparseMatrix<T>& A);
+  // Constructs an vertex adjacency for a polygon mesh.
+  //
+  // Inputs:
+  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+  //   C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) =
+  //     size of the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the
+  //     indices of the ith polygon
+  // Outputs:
+  //   A  max(I)+1 by max(I)+1 adjacency matrix, each row i corresponding to V(i,:)
+  //
+  template <typename DerivedI, typename DerivedC, typename T>
+  IGL_INLINE void adjacency_matrix(
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedC> & C,
+    Eigen::SparseMatrix<T>& A);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/cotmatrix.cpp

@@ -78,9 +78,151 @@ IGL_INLINE void igl::cotmatrix(
   L.setFromTriplets(IJV.begin(),IJV.end());
 }
 
+#include "massmatrix.h"
+#include "pinv.h"
+#include "cotmatrix_entries.h"
+#include "diag.h"
+#include "massmatrix.h"
+#include <Eigen/Geometry>
+
+template <
+  typename DerivedV, 
+  typename DerivedI, 
+  typename DerivedC, 
+  typename Scalar>
+IGL_INLINE void igl::cotmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedI> & I, 
+  const Eigen::MatrixBase<DerivedC> & C, 
+  Eigen::SparseMatrix<Scalar>& L,
+  Eigen::SparseMatrix<Scalar>& M,
+  Eigen::SparseMatrix<Scalar>& P)
+{
+  typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
+  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> MatrixXS;
+  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> VectorXS;
+  // number of vertices
+  const int n = V.rows();
+  // number of polyfaces
+  const int m = C.size()-1;
+  assert(V.cols() == 2 || V.cols() == 3);
+  std::vector<Eigen::Triplet<Scalar> > Lfijv;
+  std::vector<Eigen::Triplet<Scalar> > Mfijv;
+  std::vector<Eigen::Triplet<Scalar> > Pijv;
+  // loop over vertices; set identity for original vertices
+  for(int i = 0;i<V.rows();i++) { Pijv.emplace_back(i,i,1); }
+  // loop over faces
+  for(int p = 0;p<C.size()-1;p++)
+  {
+    // number of faces/vertices in this simple polygon
+    const int np = C(p+1)-C(p);
+    // Working "local" list of vertices; last vertex is new one
+    // this needs to have 3 columns so Eigen doesn't complain about cross
+    // products below.
+    Eigen::Matrix<Scalar,Eigen::Dynamic,3> X = decltype(X)::Zero(np+1,3);
+    for(int i = 0;i<np;i++){ X.row(i).head(V.cols()) = V.row(I(C(p)+i)); };
+    // determine weights definig position of inserted vertex
+    {
+      MatrixXS A = decltype(A)::Zero(np+1,np);
+      // My equation (38) would be A w = b.
+      VectorXS b = decltype(b)::Zero(np+1);
+      for(int k = 0;k<np;k++)
+      { 
+        const RowVector3S Xkp1mk = X.row((k+1)%np)-X.row(k);
+        const RowVector3S Xkp1mkck = Xkp1mk.cross(X.row(k));
+        for(int i = 0;i<np;i++)
+        { 
+          b(i) -= 2.*(X.row(i).cross(Xkp1mk)).dot(Xkp1mkck);
+          for(int j = 0;j<np;j++)
+          { 
+            A(i,j) += 2.*(X.row(j).cross(Xkp1mk)).dot(X.row(i).cross(Xkp1mk));
+          }
+        }
+      }
+      A.row(np).setConstant(1);
+      b(np) = 1;
+      const VectorXS w =
+        Eigen::CompleteOrthogonalDecomposition<Eigen::MatrixXd>(A).solve(b);
+      X.row(np) = w.transpose()*X.topRows(np);
+      // scatter w into new row of P
+      for(int i = 0;i<np;i++) { Pijv.emplace_back(n+p,I(C(p)+i),w(i)); }
+    }
+    // "local" fan of faces. These could be statically cached, but this will
+    // not be the bottleneck.
+    Eigen::MatrixXi F(np,3);
+    for(int i = 0;i<np;i++)
+    { 
+      F(i,0) = i; 
+      F(i,1) = (i+1)%np; 
+      F(i,2) = np; 
+    }
+    // Cotangent contributions
+    MatrixXS K;
+    igl::cotmatrix_entries(X,F,K);
+    // Massmatrix entried
+    VectorXS Mp;
+    {
+      Eigen::SparseMatrix<Scalar> M;
+      igl::massmatrix(X,F,igl::MASSMATRIX_TYPE_DEFAULT,M);
+      Mp = M.diagonal();
+    }
+    // Scatter into fine Laplacian and mass matrices
+    const auto J = [&n,&np,&p,&I,&C](int i)->int{return i==np?n+p:I(C(p)+i);};
+    // Should just build Mf as a vector...
+    for(int i = 0;i<np+1;i++) { Mfijv.emplace_back(J(i),J(i),Mp(i)); }
+    // loop over faces
+    for(int f = 0;f<np;f++)
+    {
+      for(int c = 0;c<3;c++)
+      {
+        const int i = F(f,(c+1)%3);
+        const int j = F(f,(c+2)%3);
+        // symmetric off-diagonal
+        Lfijv.emplace_back(J(i),J(j),K(f,c));
+        Lfijv.emplace_back(J(j),J(i),K(f,c));
+        // diagonal
+        Lfijv.emplace_back(J(i),J(i),-K(f,c));
+        Lfijv.emplace_back(J(j),J(j),-K(f,c));
+      }
+    }
+  }
+  P.resize(n+m,n);
+  P.setFromTriplets(Pijv.begin(),Pijv.end());
+  Eigen::SparseMatrix<Scalar> Lf(n+m,n+m);
+  Lf.setFromTriplets(Lfijv.begin(),Lfijv.end());
+  Eigen::SparseMatrix<Scalar> Mf(n+m,n+m);
+  Mf.setFromTriplets(Mfijv.begin(),Mfijv.end());
+  L = P.transpose() * Lf * P;
+  // "unlumped" M
+  const Eigen::SparseMatrix<Scalar> PTMP = P.transpose() * Mf * P;
+  // Lump M
+  const VectorXS Mdiag = PTMP * VectorXS::Ones(n,1);
+  igl::diag(Mdiag,M);
+
+  MatrixXS Vf = P*V;
+  Eigen::MatrixXi Ff(I.size(),3);
+  {
+    int f = 0;
+    for(int p = 0;p<C.size()-1;p++)
+    {
+      const int np = C(p+1)-C(p);
+      for(int c = 0;c<np;c++)
+      {
+        Ff(f,0) = I(C(p)+c);
+        Ff(f,1) = I(C(p)+(c+1)%np);
+        Ff(f,2) = V.rows()+p;
+        f++;
+      }
+    }
+    assert(f == Ff.rows());
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::SparseMatrix<double, 0, int>&, Eigen::SparseMatrix<double, 0, int>&);
+// generated by autoexplicit.sh
 template void igl::cotmatrix<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
 template void igl::cotmatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::SparseMatrix<double, 0, int>&);
 template void igl::cotmatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<double, 0, int>&);

include/igl/cotmatrix.h

@@ -48,6 +48,31 @@ namespace igl
     const Eigen::MatrixBase<DerivedV> & V, 
     const Eigen::MatrixBase<DerivedF> & F, 
     Eigen::SparseMatrix<Scalar>& L);
+  // Cotangent Laplacian (and mass matrix) for polygon meshes according to
+  // "Polygon Laplacian Made Simple" [Bunge et al. 2020]
+  //
+  // Inputs:
+  //   V  #V by 3 list of mesh vertex positions
+  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+  //   C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) = size of
+  //     the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the indices of
+  //     the ith polygon
+  // Outputs:
+  //   L  #V by #V polygon Laplacian made simple matrix
+  //   M  #V by #V mass matrix
+  //   P  #V+#polygons by #V prolongation operator
+  template <
+    typename DerivedV, 
+    typename DerivedI, 
+    typename DerivedC, 
+    typename Scalar>
+  IGL_INLINE void cotmatrix(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedI> & I, 
+    const Eigen::MatrixBase<DerivedC> & C, 
+    Eigen::SparseMatrix<Scalar>& L,
+    Eigen::SparseMatrix<Scalar>& M,
+    Eigen::SparseMatrix<Scalar>& P);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/cotmatrix_entries.cpp

@@ -136,6 +136,8 @@ IGL_INLINE void igl::cotmatrix_entries(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
 template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);

include/igl/edge_lengths.cpp

@@ -22,6 +22,8 @@ IGL_INLINE void igl::edge_lengths(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+// generated by autoexplicit.sh
 template void igl::edge_lengths<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
 // generated by autoexplicit.sh
 template void igl::edge_lengths<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);

include/igl/edges.cpp

@@ -18,6 +18,26 @@ IGL_INLINE void igl::edges(
   typedef typename DerivedF::Scalar Index;
   Eigen::SparseMatrix<Index> A;
   igl::adjacency_matrix(F,A);
+  igl::edges(A,E);
+}
+
+template <typename DerivedI, typename DerivedC, typename DerivedE>
+IGL_INLINE void igl::edges(
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedE> & E)
+{
+  typedef typename DerivedE::Scalar Index;
+  Eigen::SparseMatrix<Index> A;
+  igl::adjacency_matrix(I,C,A);
+  igl::edges(A,E);
+}
+
+template <typename T, typename DerivedE>
+IGL_INLINE void igl::edges(
+  const Eigen::SparseMatrix<T> & A,
+  Eigen::PlainObjectBase<DerivedE> & E)
+{
   // Number of non zeros should be twice number of edges
   assert(A.nonZeros()%2 == 0);
   // Resize to fit edges
@@ -27,7 +47,7 @@ IGL_INLINE void igl::edges(
   for(int k=0; k<A.outerSize(); ++k)
   {
     // Iterate over inside
-    for(typename Eigen::SparseMatrix<Index>::InnerIterator it (A,k); it; ++it)
+    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
     {
       // only add edge in one direction
       if(it.row()<it.col())
@@ -38,10 +58,13 @@ IGL_INLINE void igl::edges(
       }
     }
   }
+  assert(i == E.rows() && "A should be symmetric");
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::edges<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 template void igl::edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::edges<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);

include/igl/edges.h

@@ -10,6 +10,7 @@
 #include "igl_inline.h"
 
 #include <Eigen/Dense>
+#include <Eigen/Sparse>
 
 namespace igl
 {
@@ -27,6 +28,28 @@ namespace igl
   IGL_INLINE void edges(
     const Eigen::MatrixBase<DerivedF> & F, 
     Eigen::PlainObjectBase<DerivedE> & E);
+  // Constructs a list of unique edges represented in a given polygon mesh.
+  //
+  // Inputs:
+  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+  //   C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) =
+  //     size of the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the
+  //     indices of the ith polygon
+  // Outputs:
+  //   E #E by 2 list of edges in no particular order
+  template <typename DerivedI, typename DerivedC, typename DerivedE>
+  IGL_INLINE void edges(
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedE> & E);
+  // Inputs:
+  //   A  #V by #V symmetric adjacency matrix
+  // Outputs:
+  //   E  #E by 2 list of edges in no particular order
+  template <typename T, typename DerivedE>
+  IGL_INLINE void edges(
+    const Eigen::SparseMatrix<T> & A,
+    Eigen::PlainObjectBase<DerivedE> & E);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/massmatrix.cpp

@@ -84,6 +84,8 @@ IGL_INLINE void igl::massmatrix(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::massmatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
+// generated by autoexplicit.sh
 template void igl::massmatrix<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
 // generated by autoexplicit.sh
 template void igl::massmatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);

include/igl/pinv.cpp

@@ -1,4 +1,5 @@
 #include "pinv.h"
+#include <Eigen/SVD>
 #include <limits>
 #include <cmath>
 
@@ -33,3 +34,9 @@ void igl::pinv(
 {
   return pinv(A,-1,X);
 }
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::pinv<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/pinv.h

@@ -12,6 +12,9 @@ namespace igl
   // Outputs:
   //   X  n by m matrix so that A*X*A = A and X*A*X = X and A*X = (A*X)' and
   //     (X*A) = (X*A)'
+  //
+  // Obsolete: Use Eigen::CompleteOrthogonalDecomposition<Eigen::MatrixXd>
+  // .solve() or .pseudoinverse() instead.
   template <typename DerivedA, typename DerivedX>
   void pinv(
     const Eigen::MatrixBase<DerivedA> & A,

include/igl/polygon_corners.cpp

@@ -0,0 +1,37 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+// 
+// 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 "polygon_corners.h"
+
+template <
+  typename PType, 
+  typename DerivedI,
+  typename DerivedC>
+IGL_INLINE void igl::polygon_corners(
+  const std::vector<std::vector<PType> > & P,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  std::vector<int> vI;vI.reserve(P.size()*4);
+  C.resize(P.size()+1);
+  C(0) = 0;
+  for(int p = 0;p<P.size();p++)
+  {
+    C(p+1) = C(p)+P[p].size();
+    for(int c = 0;c<P[p].size();c++)
+    {
+      vI.push_back(P[p][c]);
+    }
+  }
+  I = Eigen::Map<DerivedI>(vI.data(),vI.size());
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::polygon_corners<int, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/polygon_corners.h

@@ -0,0 +1,41 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+// 
+// 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_POLYGON_CORNERS_H
+#define IGL_POLYGON_CORNERS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  // Convert a list-of-lists polygon mesh faces representation to list of
+  // polygon corners and sizes
+  //
+  // Inputs:
+  //   P  #P list of lists of vertex indices into rows of some matrix V
+  // Outputs:
+  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+  //   C  #P+1 list of cumulative polygon sizes so that C(i+1)-C(i) = size of
+  //     the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the indices of
+  //     the ith polygon
+  //
+  template <
+    typename PType, 
+    typename DerivedI,
+    typename DerivedC>
+  IGL_INLINE void polygon_corners(
+    const std::vector<std::vector<PType> > & P,
+    Eigen::PlainObjectBase<DerivedI> & I,
+    Eigen::PlainObjectBase<DerivedC> & C);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "polygon_corners.cpp"
+#endif
+
+#endif 

include/igl/polygon_mesh_to_triangle_mesh.cpp

@@ -1,77 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2014 Daniele Panozzo <[email protected]>
-//
-// 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 "polygon_mesh_to_triangle_mesh.h"
-#include "matrix_to_list.h"
-
-template <typename Index, typename DerivedF>
-IGL_INLINE void igl::polygon_mesh_to_triangle_mesh(
-  const std::vector<std::vector<Index> > & vF,
-  Eigen::PlainObjectBase<DerivedF>& F)
-{
-  using namespace std;
-  using namespace Eigen;
-  int m = 0;
-  // estimate of size
-  for(typename vector<vector<Index > >::const_iterator fit = vF.begin();
-    fit!=vF.end();
-    fit++)
-  {
-    if(fit->size() >= 3)
-    {
-      m += fit->size() - 2;
-    }
-  }
-  // Resize output
-  F.resize(m,3);
-  {
-    int k = 0;
-    for(typename vector<vector<Index > >::const_iterator fit = vF.begin();
-      fit!=vF.end();
-      fit++)
-    {
-      if(fit->size() >= 3)
-      {
-        typename vector<Index >::const_iterator cit = fit->begin();
-        cit++;
-        typename vector<Index >::const_iterator pit = cit++;
-        for(;
-          cit!=fit->end();
-          cit++,pit++)
-        {
-          F(k,0) = *(fit->begin());
-          F(k,1) = *pit;
-          F(k,2) = *cit;
-          k++;
-        }
-      }
-    }
-    assert(k==m);
-  }
-
-}
-
-template <typename DerivedP, typename DerivedF>
-IGL_INLINE void igl::polygon_mesh_to_triangle_mesh(
-  const Eigen::PlainObjectBase<DerivedP>& P,
-  Eigen::PlainObjectBase<DerivedF>& F)
-{
-  std::vector<std::vector<typename DerivedP::Scalar> > vP;
-  matrix_to_list(P,vP);
-  return polygon_mesh_to_triangle_mesh(vP,F);
-}
-
-#ifdef IGL_STATIC_LIBRARY
-// Explicit template instantiation
-// generated by autoexplicit.sh
-template void igl::polygon_mesh_to_triangle_mesh<int, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
-// generated by autoexplicit.sh
-template void igl::polygon_mesh_to_triangle_mesh<int, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
-template void igl::polygon_mesh_to_triangle_mesh<int, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::polygon_mesh_to_triangle_mesh<int, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
-template void igl::polygon_mesh_to_triangle_mesh<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-#endif

include/igl/polygon_mesh_to_triangle_mesh.h

@@ -1,48 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2013 Daniele Panozzo <[email protected]>
-//
-// 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_POLYGON_MESH_TO_TRIANGLE_MESH_H
-#define IGL_POLYGON_MESH_TO_TRIANGLE_MESH_H
-#include "igl_inline.h"
-
-#ifndef IGL_NO_EIGEN
-#  include <Eigen/Core>
-#endif
-#include <vector>
-
-namespace igl
-{
-  // Triangulate a general polygonal mesh into a triangle mesh.
-  //
-  // Inputs:
-  //   vF  list of polygon index lists
-  // Outputs:
-  //   F  eigen int matrix #F by 3
-  //
-  // Example:
-  //   vector<vector<double > > vV;
-  //   vector<vector<int > > vF;
-  //   read_triangle_mesh("poly.obj",vV,vF);
-  //   MatrixXd V;
-  //   MatrixXi F;
-  //   list_to_matrix(vV,V);
-  //   triangulate(vF,F);
-  template <typename Index, typename DerivedF>
-  IGL_INLINE void polygon_mesh_to_triangle_mesh(
-    const std::vector<std::vector<Index> > & vF,
-    Eigen::PlainObjectBase<DerivedF>& F);
-  template <typename DerivedP, typename DerivedF>
-  IGL_INLINE void polygon_mesh_to_triangle_mesh(
-    const Eigen::PlainObjectBase<DerivedP>& P,
-    Eigen::PlainObjectBase<DerivedF>& F);
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "polygon_mesh_to_triangle_mesh.cpp"
-#endif
-
-#endif

include/igl/polygons_to_triangles.cpp

@@ -0,0 +1,45 @@
+#include "polygons_to_triangles.h"
+
+template <
+  typename DerivedI,
+  typename DerivedC,
+  typename DerivedF,
+  typename DerivedJ>
+IGL_INLINE void igl::polygons_to_triangles(
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  // Each polygon results in #sides-2 triangles. So ∑#sides-2
+  F.resize(C(C.size()-1) - (C.size()-1)*2,3);
+  J.resize(F.rows());
+  {
+    int f = 0;
+    for(int p = 0;p<C.size()-1;p++)
+    {
+      const int np = C(p+1)-C(p);
+      for(int c = 1;c<np-1;c++)
+      {
+        F(f,0) = I(C(p)+0);
+        F(f,1) = I(C(p)+c);
+        F(f,2) = I(C(p)+c+1);
+        J(f) = p;
+        f++;
+      }
+    }
+    assert(f == F.rows());
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::polygons_to_triangles<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::polygons_to_triangles<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::polygons_to_triangles<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::polygons_to_triangles<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/polygons_to_triangles.h

@@ -0,0 +1,46 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+// 
+// 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_POLYGONS_TO_TRIANGLES_H
+#define IGL_POLYGONS_TO_TRIANGLES_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  // Given a polygon mesh, trivially triangulate each polygon with a fan. This
+  // purely combinatorial triangulation will work well for convex/flat polygons
+  // and degrade otherwise.
+  //
+  // Inputs:
+  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+  //   C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) =
+  //     size of the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the
+  //     indices of the ith polygon
+  // Outputs:
+  //   F  #F by 3 list of triangle indices into rows of V
+  //   J  #F list of indices into 0:#P-1 of corresponding polygon
+  //
+  template <
+    typename DerivedI,
+    typename DerivedC,
+    typename DerivedF,
+    typename DerivedJ>
+  IGL_INLINE void polygons_to_triangles(
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedJ> & J);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "polygons_to_triangles.cpp"
+#endif
+
+#endif

include/igl/predicates/polygons_to_triangles.cpp

@@ -0,0 +1,247 @@
+#include "polygons_to_triangles.h"
+#include "ear_clipping.h"
+#include "../sort.h"
+#include "../slice.h"
+#include <Eigen/Eigenvalues>
+
+template <
+  typename DerivedV,
+  typename DerivedI,
+  typename DerivedC,
+  typename DerivedF,
+  typename DerivedJ>
+IGL_INLINE void igl::predicates::polygons_to_triangles(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  // Each polygon results in #sides-2 triangles. So ∑#sides-2
+  F.resize(C(C.size()-1) - (C.size()-1)*2,3);
+  J.resize(F.rows());
+  {
+    int f = 0;
+    for(int p = 0;p<C.size()-1;p++)
+    {
+      const int np = C(p+1)-C(p);
+      Eigen::MatrixXi pF;
+      if(np == 3)
+      {
+        pF = (Eigen::MatrixXi(1,3)<<0,1,2).finished();
+      }else
+      {
+        // Make little copy of this polygon with an initial fan
+        DerivedV pV(np,V.cols());
+        for(int c = 0;c<np;c++)
+        {
+          pV.row(c) = V.row(I(C(p)+c));
+        }
+        // Use PCA to project to 2D
+        Eigen::MatrixXd S;
+        switch(V.cols())
+        {
+          case 2:
+            S = V.template cast<double>();
+            break;
+          case 3:
+          {
+            Eigen::MatrixXd P = (pV.rowwise() - pV.colwise().mean()).template cast<double>();
+            Eigen::Matrix3d O = P.transpose() * P;
+            Eigen::EigenSolver<Eigen::Matrix3d> es(O);
+            Eigen::Matrix3d C = es.eigenvectors().real();
+            {
+              Eigen::Vector3d _1;
+              Eigen::Vector3i I;
+              igl::sort(es.eigenvalues().real().eval(),1,false,_1,I);
+              igl::slice(Eigen::Matrix3d(C),I,2,C);
+            }
+            S = P*C.leftCols(2);
+            break;
+          }
+          default: assert(false && "dim>3 not supported");
+        }
+
+        Eigen::VectorXi RT = Eigen::VectorXi::Zero(S.rows(),1);
+        Eigen::VectorXi _I;
+        Eigen::MatrixXd _nS;
+
+        // compute signed area
+        {
+          double area = 0;
+          for(int c = 0;c<np;c++)
+          {
+            area += S((c+0)%np,0)*S((c+1)%np,1) - S((c+1)%np,0)*S((c+0)%np,1);
+          }
+          //printf("area: %g\n",area);
+          if(area<0)
+          {
+            S.col(0) *= -1;
+          }
+        }
+
+        // This is a really low quality triangulator and will contain nearly
+        // degenerate elements which become degenerate or worse when unprojected
+        // back to 3D.
+        igl::predicates::ear_clipping(S,RT,_I,pF,_nS);
+        // igl::predicates::ear_clipping does not gracefully fail when the input
+        // is not simple. Instead it (tends?) to output too few triangles.
+        if(pF.rows() < np-2)
+        {
+          // Fallback, use a fan
+          //std::cout<<igl::matlab_format(S,"S")<<std::endl;
+          //std::cout<<igl::matlab_format(RT,"RT")<<std::endl;
+          //std::cout<<igl::matlab_format(_I,"I")<<std::endl;
+          //std::cout<<igl::matlab_format(pF,"pF")<<std::endl;
+          //std::cout<<igl::matlab_format(_nS,"nS")<<std::endl;
+          //std::cout<<std::endl;
+
+          pF.resize(np-2,3);
+          for(int c = 0;c<np;c++)
+          {
+            if(c>0 && c<np-1)
+            {
+              pF(c-1,0) = 0;
+              pF(c-1,1) = c;
+              pF(c-1,2) = c+1;
+            }
+          }
+        }
+        assert(pF.rows() == np-2);
+
+        // Could at least flip edges of degenerate edges
+
+        //if(pF.rows()>1)
+        //{
+        //  // Delaunay-ize 
+        //  Eigen::MatrixXd pl;
+        //  igl::edge_lengths(pV,pF,pl);
+
+        //  typedef Eigen::Matrix<int,Eigen::Dynamic,2> MatrixX2I;
+        //  typedef Eigen::Matrix<int,Eigen::Dynamic,1> VectorXI;
+        //  MatrixX2I E,uE;
+        //  VectorXI EMAP;
+        //  std::vector<std::vector<int> > uE2E;
+        //  igl::unique_edge_map(pF, E, uE, EMAP, uE2E);
+        //  typedef int Index;
+        //  typedef double Scalar;
+        //  const Index num_faces = pF.rows();
+        //  std::vector<Index> Q;
+        //  Q.reserve(uE2E.size());
+        //  for (size_t uei=0; uei<uE2E.size(); uei++) 
+        //  {
+        //    Q.push_back(uei);
+        //  }
+        //  while(!Q.empty())
+        //  {
+        //    const Index uei = Q.back();
+        //    Q.pop_back();
+        //    if (uE2E[uei].size() == 2) 
+        //    {
+        //      double w;
+        //      igl::is_intrinsic_delaunay(pl,uE2E,num_faces,uei,w);
+        //      printf("%d : %0.17f\n",uei,w);
+        //      if(w<-1e-7) 
+        //      {
+        //        printf("  flippin'\n");
+        //        //
+        //        //          v1                 v1
+        //        //          /|\                / \
+        //        //        c/ | \b            c/f1 \b
+        //        //     v3 /f2|f1\ v4  =>  v3 /__f__\ v4
+        //        //        \  e  /            \ f2  /
+        //        //        d\ | /a            d\   /a
+        //        //          \|/                \ /
+        //        //          v2                 v2
+        //        //
+        //        // hmm... is the flip actually in the other direction?
+        //        const Index f1 = uE2E[uei][0]%num_faces;
+        //        const Index f2 = uE2E[uei][1]%num_faces;
+        //        const Index c1 = uE2E[uei][0]/num_faces;
+        //        const Index c2 = uE2E[uei][1]/num_faces;
+        //        const size_t e_24 = f1 + ((c1 + 1) % 3) * num_faces;
+        //        const size_t e_41 = f1 + ((c1 + 2) % 3) * num_faces;
+        //        const size_t e_13 = f2 + ((c2 + 1) % 3) * num_faces;
+        //        const size_t e_32 = f2 + ((c2 + 2) % 3) * num_faces;
+        //        const size_t ue_24 = EMAP(e_24);
+        //        const size_t ue_41 = EMAP(e_41);
+        //        const size_t ue_13 = EMAP(e_13);
+        //        const size_t ue_32 = EMAP(e_32);
+        //        // new edge lengths
+        //        const Index v1 = pF(f1, (c1+1)%3);
+        //        const Index v2 = pF(f1, (c1+2)%3);
+        //        const Index v4 = pF(f1, c1);
+        //        const Index v3 = pF(f2, c2);
+        //        {
+        //          const Scalar e = pl(f1,c1);
+        //          const Scalar a = pl(f1,(c1+1)%3);
+        //          const Scalar b = pl(f1,(c1+2)%3);
+        //          const Scalar c = pl(f2,(c2+1)%3);
+        //          const Scalar d = pl(f2,(c2+2)%3);
+        //          const double f = (pV.row(v3)-pV.row(v4)).norm();
+        //          // New order
+        //          pl(f1,0) = f;
+        //          pl(f1,1) = b;
+        //          pl(f1,2) = c;
+        //          pl(f2,0) = f;
+        //          pl(f2,1) = d;
+        //          pl(f2,2) = a;
+        //        }
+        //        printf("%d,%d %d,%d -> %d,%d\n",uE(uei,0),uE(uei,1),v1,v2,v3,v4);
+        //        igl::flip_edge(pF, E, uE, EMAP, uE2E, uei);
+        //        std::cout<<"  "<<pl.row(f1)<<std::endl;
+        //        std::cout<<"  "<<pl.row(f2)<<std::endl;
+        //        //// new edge lengths, slow!
+        //        //igl::edge_lengths(pV,pF,pl);
+        //        // recompute edge lengths of two faces. (extra work on untouched
+        //        // edges)
+        //        for(int f : {f1,f2})
+        //        {
+        //          for(int c=0;c<3;c++)
+        //          {
+        //            pl(f,c) = 
+        //              (pV.row(pF(f,(c+1)%3))-pV.row(pF(f,(c+2)%3))).norm();
+        //          }
+        //        }
+        //        std::cout<<"  "<<pl.row(f1)<<std::endl;
+        //        std::cout<<"  "<<pl.row(f2)<<std::endl;
+        //        std::cout<<std::endl;
+
+        //        Q.push_back(ue_24);
+        //        Q.push_back(ue_41);
+        //        Q.push_back(ue_13);
+        //        Q.push_back(ue_32);
+        //      }
+        //    }
+        //  }
+
+
+        //  // check for self-loops (I claim these cannot happen)
+        //  for(int f = 0;f<pF.rows();f++)
+        //  {
+        //    for(int c =0;c<3;c++)
+        //    {
+        //      assert(pF(f,c) != pF(f,(c+1)%3) && "self loops should not exist");
+        //    }
+        //  }
+        //}
+      }
+      // Copy into global list
+      for(int i = 0;i<pF.rows();i++)
+      {
+        for(int c =0;c<3;c++)
+        {
+          F(f,c) = I(C(p)+pF(i,c));
+        }
+        J(f) = p;
+        f++;
+      }
+
+    }
+    assert(f == F.rows());
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+#endif

include/igl/predicates/polygons_to_triangles.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+// 
+// 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_POLYGONS_TO_TRIANGLES_H
+#define IGL_POLYGONS_TO_TRIANGLES_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  namespace predicates
+  {
+    // Given a polygon mesh, trivially triangulate each polygon with a fan. This
+    // purely combinatorial triangulation will work well for convex/flat polygons
+    // and degrade otherwise.
+    //
+    // Inputs:
+    //   V  #V by dim list of vertex positions
+    //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+    //   C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) =
+    //     size of the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the
+    //     indices of the ith polygon
+    // Outputs:
+    //   F  #F by 3 list of triangle indices into rows of V
+    //   J  #F list of indices into 0:#P-1 of corresponding polygon
+    //
+    template <
+      typename DerivedV,
+      typename DerivedI,
+      typename DerivedC,
+      typename DerivedF,
+      typename DerivedJ>
+    IGL_INLINE void polygons_to_triangles(
+      const Eigen::MatrixBase<DerivedV> & V,
+      const Eigen::MatrixBase<DerivedI> & I,
+      const Eigen::MatrixBase<DerivedC> & C,
+      Eigen::PlainObjectBase<DerivedF> & F,
+      Eigen::PlainObjectBase<DerivedJ> & J);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "polygons_to_triangles.cpp"
+#endif
+
+#endif

include/igl/readOBJ.cpp

@@ -10,6 +10,8 @@
 #include "list_to_matrix.h"
 #include "max_size.h"
 #include "min_size.h"
+#include "polygon_corners.h"
+#include "polygons_to_triangles.h"
 
 #include <iostream>
 #include <cstdio>
@@ -396,9 +398,36 @@ IGL_INLINE bool igl::readOBJ(
   return true;
 }
 
+template <typename DerivedV, typename DerivedI, typename DerivedC>
+IGL_INLINE bool igl::readOBJ(
+  const std::string str,
+  Eigen::PlainObjectBase<DerivedV>& V,
+  Eigen::PlainObjectBase<DerivedI>& I,
+  Eigen::PlainObjectBase<DerivedC>& C)
+{
+  // we should flip this so that the base implementation uses arrays.
+  std::vector<std::vector<double> > vV,vTC,vN;
+  std::vector<std::vector<int> > vF,vFTC,vFN;
+  bool success = igl::readOBJ(str,vV,vTC,vN,vF,vFTC,vFN);
+  if(!success)
+  {
+    // readOBJ(str,vV,vTC,vN,vF,vFTC,vFN) should have already printed an error
+    // message to stderr
+    return false;
+  }
+  if(!igl::list_to_matrix(vV,V))
+  {
+    return false;
+  }
+  igl::polygon_corners(vF,I,C);
+  return true;
+}
+
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::readOBJ<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template bool igl::readOBJ<double, int>(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
 template bool igl::readOBJ<double, int>(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::tuple<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, int, int>, std::allocator<std::tuple<std::basic_string<char, std::char_traits<char>, std::allocator<char> >, int, int> > >&);
 template bool igl::readOBJ<double, int>(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);

include/igl/readOBJ.h

@@ -5,6 +5,7 @@
 // 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_READOBJ_H
 #define IGL_READOBJ_H
 #include "igl_inline.h"
@@ -121,6 +122,17 @@ namespace igl
     const std::string str,
     Eigen::PlainObjectBase<DerivedV>& V,
     Eigen::PlainObjectBase<DerivedF>& F);
+  // Outputs:
+  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+  //   C  #P+1 list of cumulative polygon sizes so that C(i+1)-C(i) = size of
+  //     the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the indices of
+  //     the ith polygon
+  template <typename DerivedV, typename DerivedI, typename DerivedC>
+  IGL_INLINE bool readOBJ(
+    const std::string str,
+    Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedI>& I,
+    Eigen::PlainObjectBase<DerivedC>& C);
 
 }
 

include/igl/read_triangle_mesh.cpp

@@ -17,7 +17,8 @@
 #include "readWRL.h"
 #include "pathinfo.h"
 #include "boundary_facets.h"
-#include "polygon_mesh_to_triangle_mesh.h"
+#include "polygon_corners.h"
+#include "polygons_to_triangles.h"
 
 #include <algorithm>
 #include <iostream>
@@ -192,7 +193,12 @@ IGL_INLINE bool igl::read_triangle_mesh(
     {
       return false;
     }
-    polygon_mesh_to_triangle_mesh(vF,F);
+    {
+      Eigen::VectorXi I,C;
+      igl::polygon_corners(vF,I,C);
+      Eigen::VectorXi J;
+      igl::polygons_to_triangles(I,C,F,J);
+    }
   }
   return true;
 }

include/igl/squared_edge_lengths.cpp

@@ -73,6 +73,8 @@ IGL_INLINE void igl::squared_edge_lengths(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+// generated by autoexplicit.sh
 template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
 // generated by autoexplicit.sh
 template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);

tutorial/207_PolygonLaplacian/CMakeLists.txt

@@ -0,0 +1,5 @@
+get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+project(${PROJECT_NAME})
+
+add_executable(${PROJECT_NAME} main.cpp)
+target_link_libraries(${PROJECT_NAME} igl::core igl::opengl igl::opengl_glfw tutorials)

tutorial/207_PolygonLaplacian/main.cpp

@@ -0,0 +1,128 @@
+#include <igl/readOBJ.h>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/edges.h>
+#include <igl/find.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/polygons_to_triangles.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <iostream>
+#include "tutorial_shared_path.h"
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+
+  // Load a mesh in OBJ format 
+  Eigen::MatrixXd OV,V;
+  Eigen::VectorXi I,C;
+  igl::readOBJ(
+    argc<=1?TUTORIAL_SHARED_PATH "/cylinder.obj" :argv[1],V,I,C);
+  // center
+  V.rowwise() -= V.colwise().mean();
+  OV = V;
+  // Convert polygon representation to triangles
+  Eigen::MatrixXi F;
+  Eigen::VectorXi J;
+  igl::polygons_to_triangles(I,C,F,J);
+  Eigen::SparseMatrix<double> pL,pM,pP;
+  igl::cotmatrix(V,I,C,pL,pM,pP);
+  Eigen::SparseMatrix<double> tL,tM;
+  igl::cotmatrix(V,F,tL);
+  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,tM);
+  const double bbd = (V.colwise().maxCoeff()- V.colwise().minCoeff()).norm();
+  igl::opengl::glfw::Viewer vr;
+  vr.data_list[0].set_mesh(V,F);
+  vr.append_mesh();
+  vr.selected_data_index = 0;
+  vr.data_list[0].set_face_based(true);
+  Eigen::MatrixXi E;
+  igl::edges(I,C,E);
+  bool show_edges = true;
+  bool use_poly = true;
+  Eigen::MatrixXd pHN;
+  Eigen::MatrixXd tHN;
+  const auto update = [&]()
+  {
+    pHN = (pL*V).array().colwise() / pM.diagonal().array();
+    tHN = (tL*V).array().colwise() / tM.diagonal().array();
+    pHN *= 1.0/pHN.rowwise().norm().maxCoeff();
+    tHN *= 1.0/tHN.rowwise().norm().maxCoeff();
+    const auto was_face_based  = vr.data_list[0].face_based;
+    Eigen::MatrixXd QV(V.rows()*2,3);
+    QV.topRows(V.rows()) = V;
+    if(use_poly)
+    {
+      printf("using polygon Laplacian\n");
+      QV.bottomRows(V.rows()) = V-pHN;
+    }else
+    {
+      printf("using triangle Laplacian\n");
+      QV.bottomRows(V.rows()) = V-tHN;
+    }
+    Eigen::MatrixXi QE(V.rows(),2);
+    for(int i = 0;i<V.rows();i++){ QE(i,0)=i;QE(i,1)=i+V.rows();}
+    vr.data_list[1].set_edges(QV,QE,Eigen::RowVector3d(1,1,1));
+
+    if(use_poly)
+    {
+      vr.data_list[0].show_lines = false;
+      if(show_edges)
+      {
+        vr.data_list[0].set_edges(V,E,Eigen::RowVector3d(0,0,0));
+      }else
+      {
+        vr.data_list[0].clear_edges();
+      }
+    }else
+    {
+      vr.data_list[0].clear_edges();
+      vr.data_list[0].show_lines = show_edges;
+    }
+    vr.data_list[0].set_face_based(was_face_based);
+  };
+  const double original_area = pM.diagonal().sum();
+  const auto recompute_M = [&]()
+  {
+    Eigen::SparseMatrix<double> _1,_2;
+    igl::cotmatrix(V,I,C,_1,pM,_2);
+    igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,tM);
+    V *= sqrt(original_area / pM.diagonal().sum());
+  };
+  const auto cmcf_step = [&]()
+  {
+    // take step
+    Eigen::SparseMatrix<double> & L = tL;
+    Eigen::SparseMatrix<double> & M = tM;
+    if(use_poly)
+    {
+      L = pL;
+      M = tM;
+    }
+    const auto & S = (M - 0.05*L);
+    Eigen::SimplicialLLT<Eigen::SparseMatrix<double > > solver(S);
+    assert(solver.info() == Eigen::Success);
+    V = solver.solve(M*V).eval();
+    // recompute just mass matrices
+    recompute_M();
+    // center
+    V.rowwise() -= V.colwise().mean();
+    vr.data_list[0].set_vertices(V);
+    vr.data_list[0].compute_normals();
+    update();
+  };
+  vr.callback_key_pressed = [&](decltype(vr) &,unsigned int key, int mod)
+  {
+    switch(key)
+    {
+      case ' ': cmcf_step(); return true;
+      case 'R': case 'r': V=OV;recompute_M();vr.data_list[0].set_vertices(V);vr.data_list[0].compute_normals(); update();return true;
+      case 'P': case 'p': use_poly=!use_poly; update();return true;
+      case 'L': case 'l': show_edges=!show_edges; update();return true;
+    }
+    return false;
+  };
+  update();
+  vr.launch();
+}

tutorial/CMakeLists.txt

@@ -75,6 +75,7 @@ if(TUTORIALS_CHAPTER2)
   add_subdirectory("204_Gradient")
   add_subdirectory("205_Laplacian")
   add_subdirectory("206_GeodesicDistance")
+  add_subdirectory("207_PolygonLaplacian")
 endif()
 
 # Chapter 3