Merge pull request #988 from libigl/dev-intrinsic-delaunay

Heat Geodesics and Intrinsic Delaunay Triangulation

include/igl/boundary_facets.cpp

@@ -103,7 +103,7 @@ IGL_INLINE void igl::boundary_facets(
 
 template <typename DerivedT, typename DerivedF>
 IGL_INLINE void igl::boundary_facets(
-  const Eigen::PlainObjectBase<DerivedT>& T,
+  const Eigen::MatrixBase<DerivedT>& T,
   Eigen::PlainObjectBase<DerivedF>& F)
 {
   assert(T.cols() == 0 || T.cols() == 4 || T.cols() == 3);
@@ -119,7 +119,7 @@ IGL_INLINE void igl::boundary_facets(
 
 template <typename DerivedT, typename Ret>
 Ret igl::boundary_facets(
-  const Eigen::PlainObjectBase<DerivedT>& T)
+  const Eigen::MatrixBase<DerivedT>& T)
 {
   Ret F;
   igl::boundary_facets(T,F);
@@ -130,12 +130,12 @@ Ret igl::boundary_facets(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
+template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
 // generated by autoexplicit.sh
-template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
-template void igl::boundary_facets<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> >&);
+template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::boundary_facets<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::boundary_facets<int, int>(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
-//template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > igl::boundary_facets(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
-template Eigen::Matrix<int, -1, -1, 0, -1, -1> igl::boundary_facets<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&);
-template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
+//template Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > igl::boundary_facets(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
+template Eigen::Matrix<int, -1, -1, 0, -1, -1> igl::boundary_facets<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&);
+template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
 #endif

include/igl/boundary_facets.h

@@ -37,12 +37,12 @@ namespace igl
   //   DerivedF  integer-value: i.e. from MatrixXi
   template <typename DerivedT, typename DerivedF>
   IGL_INLINE void boundary_facets(
-    const Eigen::PlainObjectBase<DerivedT>& T,
+    const Eigen::MatrixBase<DerivedT>& T,
     Eigen::PlainObjectBase<DerivedF>& F);
   // Same as above but returns F
   template <typename DerivedT, typename Ret>
   Ret boundary_facets(
-    const Eigen::PlainObjectBase<DerivedT>& T);
+    const Eigen::MatrixBase<DerivedT>& T);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/circulation.cpp

@@ -7,12 +7,11 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "circulation.h"
 #include "list_to_matrix.h"
+#include <cassert>
 
 IGL_INLINE std::vector<int> igl::circulation(
   const int e,
   const bool ccw,
-  const Eigen::MatrixXi & F,
-  const Eigen::MatrixXi & E,
   const Eigen::VectorXi & EMAP,
   const Eigen::MatrixXi & EF,
   const Eigen::MatrixXi & EI)
@@ -20,7 +19,8 @@ IGL_INLINE std::vector<int> igl::circulation(
   // prepare output
   std::vector<int> N;
   N.reserve(6);
-  const int m = F.rows();
+  const int m = EMAP.size()/3;
+  assert(m*3 == EMAP.size());
   const auto & step = [&](
     const int e, 
     const int ff,
@@ -59,13 +59,11 @@ IGL_INLINE std::vector<int> igl::circulation(
 IGL_INLINE void igl::circulation(
   const int e,
   const bool ccw,
-  const Eigen::MatrixXi & F,
-  const Eigen::MatrixXi & E,
   const Eigen::VectorXi & EMAP,
   const Eigen::MatrixXi & EF,
   const Eigen::MatrixXi & EI,
   Eigen::VectorXi & vN)
 {
-  std::vector<int> N = circulation(e,ccw,F,E,EMAP,EF,EI);
+  std::vector<int> N = circulation(e,ccw,EMAP,EF,EI);
   igl::list_to_matrix(N,vN);
 }

include/igl/circulation.h

@@ -20,8 +20,6 @@ namespace igl
   //   e  index into E of edge to circulate
   //   ccw  whether to _continue_ in ccw direction of edge (circulate around
   //     E(e,1))
-  //   F  #F by 3 list of face indices
-  //   E  #E by 2 list of edge indices
   //   EMAP #F*3 list of indices into E, mapping each directed edge to unique
   //     unique edge in E
   //   EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
@@ -30,11 +28,10 @@ namespace igl
   //   EI  #E by 2 list of edge flap corners (see above).
   // Returns list of faces touched by circulation (in cyclically order).
   //   
+  // See also: edge_flaps
   IGL_INLINE std::vector<int> circulation(
     const int e,
     const bool ccw,
-    const Eigen::MatrixXi & F,
-    const Eigen::MatrixXi & E,
     const Eigen::VectorXi & EMAP,
     const Eigen::MatrixXi & EF,
     const Eigen::MatrixXi & EI);
@@ -42,8 +39,6 @@ namespace igl
   IGL_INLINE void circulation(
     const int e,
     const bool ccw,
-    const Eigen::MatrixXi & F,
-    const Eigen::MatrixXi & E,
     const Eigen::VectorXi & EMAP,
     const Eigen::MatrixXi & EF,
     const Eigen::MatrixXi & EI,

include/igl/collapse_edge.cpp

@@ -40,7 +40,7 @@ IGL_INLINE bool igl::collapse_edge(
   }
 
   // Important to grab neighbors of d before monkeying with edges
-  const std::vector<int> nV2Fd = circulation(e,!eflip,F,E,EMAP,EF,EI);
+  const std::vector<int> nV2Fd = circulation(e,!eflip,EMAP,EF,EI);
 
   // The following implementation strongly relies on s<d
   assert(s<d && "s should be less than d");
@@ -339,8 +339,8 @@ IGL_INLINE bool igl::collapse_edge(
   Q.erase(Q.begin());
   e = p.second;
   Qit[e] = Q.end();
-  std::vector<int> N  = circulation(e, true,F,E,EMAP,EF,EI);
-  std::vector<int> Nd = circulation(e,false,F,E,EMAP,EF,EI);
+  std::vector<int> N  = circulation(e, true,EMAP,EF,EI);
+  std::vector<int> Nd = circulation(e,false,EMAP,EF,EI);
   N.insert(N.begin(),Nd.begin(),Nd.end());
   bool collapsed = true;
   if(pre_collapse(V,F,E,EMAP,EF,EI,Q,Qit,C,e))

include/igl/copyleft/progressive_hulls_cost_and_placement.cpp

@@ -32,8 +32,8 @@ IGL_INLINE void igl::copyleft::progressive_hulls_cost_and_placement(
 
   assert(V.cols() == 3 && "V.cols() should be 3");
   // Gather list of unique face neighbors
-  vector<int> Nall =  circulation(e, true,F,E,EMAP,EF,EI);
-  vector<int> Nother= circulation(e,false,F,E,EMAP,EF,EI);
+  vector<int> Nall =  circulation(e, true,EMAP,EF,EI);
+  vector<int> Nother= circulation(e,false,EMAP,EF,EI);
   Nall.insert(Nall.end(),Nother.begin(),Nother.end());
   vector<int> N;
   igl::unique(Nall,N);

include/igl/cotmatrix_entries.cpp

@@ -1,6 +1,7 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
 // 
 // 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 
@@ -50,6 +51,7 @@ IGL_INLINE void igl::cotmatrix_entries(
       C.resize(m,3);
       for(int i = 0;i<m;i++)
       {
+        // Alec: I'm doubtful that using l2 here is actually improving numerics.
         C(i,0) = (l2(i,1) + l2(i,2) - l2(i,0))/dblA(i)/4.0;
         C(i,1) = (l2(i,2) + l2(i,0) - l2(i,1))/dblA(i)/4.0;
         C(i,2) = (l2(i,0) + l2(i,1) - l2(i,2))/dblA(i)/4.0;
@@ -97,14 +99,49 @@ IGL_INLINE void igl::cotmatrix_entries(
   }
 }
 
+template <typename Derivedl, typename DerivedC>
+IGL_INLINE void igl::cotmatrix_entries(
+  const Eigen::MatrixBase<Derivedl>& l,
+  Eigen::PlainObjectBase<DerivedC>& C)
+{
+  using namespace Eigen;
+  const int m = l.rows();
+  assert(l.cols() == 3 && "Only triangles accepted");
+  //Compute squared Edge lengths 
+  Matrix<typename DerivedC::Scalar,Dynamic,3> l2;
+  l2 = l.array().square();
+  // Alec: It's a little annoying that there's duplicate code here. The
+  // "extrinic" version above is first computing squared edge lengths, taking
+  // the square root and calling this. We can't have a cotmatrix_entries(l,l2,C)
+  // overload because it will confuse Eigen with the cotmatrix_entries(V,F,C)
+  // overload. In the end, I'd like to be convinced that using l2 directly above
+  // is actually better numerically (or significantly faster) than just calling
+  // edge_lengths and this cotmatrix_entries(l,C);
+  //
+  // double area
+  Matrix<typename DerivedC::Scalar,Dynamic,1> dblA;
+  doublearea(l,0.,dblA);
+  // cotangents and diagonal entries for element matrices
+  // correctly divided by 4 (alec 2010)
+  C.resize(m,3);
+  for(int i = 0;i<m;i++)
+  {
+    // Alec: I'm doubtful that using l2 here is actually improving numerics.
+    C(i,0) = (l2(i,1) + l2(i,2) - l2(i,0))/dblA(i)/4.0;
+    C(i,1) = (l2(i,2) + l2(i,0) - l2(i,1))/dblA(i)/4.0;
+    C(i,2) = (l2(i,0) + l2(i,1) - l2(i,2))/dblA(i)/4.0;
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+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> >&);
+template void igl::cotmatrix_entries<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> >&);
+template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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::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, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, 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, 4, 0, -1, 4> > 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<int, -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::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/cotmatrix_entries.h

@@ -1,6 +1,7 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
 // 
 // 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 
@@ -28,6 +29,17 @@ namespace igl
     const Eigen::MatrixBase<DerivedV>& V,
     const Eigen::MatrixBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedC>& C);
+  // Intrinsic version.
+  //
+  // Inputs:
+  //   l  #F by 3 list of triangle edge lengths (see edge_lengths)
+  // Outputs:
+  //   C  #F by 3 list of 1/2*cotangents corresponding angles
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  template <typename Derivedl, typename DerivedC>
+  IGL_INLINE void cotmatrix_entries(
+    const Eigen::MatrixBase<Derivedl>& l,
+    Eigen::PlainObjectBase<DerivedC>& C);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/cotmatrix_intrinsic.cpp

@@ -0,0 +1,67 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "cotmatrix_intrinsic.h"
+#include "cotmatrix_entries.h"
+#include <iostream>
+
+template <typename Derivedl, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::cotmatrix_intrinsic(
+  const Eigen::MatrixBase<Derivedl> & l, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L)
+{
+  using namespace Eigen;
+  using namespace std;
+  // Cribbed from cotmatrix
+
+  const int nverts = F.maxCoeff()+1;
+  L.resize(nverts,nverts);
+  Matrix<int,Dynamic,2> edges;
+  int simplex_size = F.cols();
+  // 3 for triangles, 4 for tets
+  assert(simplex_size == 3);
+  // This is important! it could decrease the comptuation time by a factor of 2
+  // Laplacian for a closed 2d manifold mesh will have on average 7 entries per
+  // row
+  L.reserve(10*nverts);
+  edges.resize(3,2);
+  edges << 
+    1,2,
+    2,0,
+    0,1;
+  // Gather cotangents
+  Matrix<Scalar,Dynamic,Dynamic> C;
+  cotmatrix_entries(l,C);
+  
+  vector<Triplet<Scalar> > IJV;
+  IJV.reserve(F.rows()*edges.rows()*4);
+  // Loop over triangles
+  for(int i = 0; i < F.rows(); i++)
+  {
+    // loop over edges of element
+    for(int e = 0;e<edges.rows();e++)
+    {
+      int source = F(i,edges(e,0));
+      int dest = F(i,edges(e,1));
+      IJV.push_back(Triplet<Scalar>(source,dest,C(i,e)));
+      IJV.push_back(Triplet<Scalar>(dest,source,C(i,e)));
+      IJV.push_back(Triplet<Scalar>(source,source,-C(i,e)));
+      IJV.push_back(Triplet<Scalar>(dest,dest,-C(i,e)));
+    }
+  }
+  L.setFromTriplets(IJV.begin(),IJV.end());
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::cotmatrix_intrinsic<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&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cotmatrix_intrinsic<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_intrinsic<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_intrinsic<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>&);
+template void igl::cotmatrix_intrinsic<Eigen::Matrix<double, -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::SparseMatrix<double, 0, int>&);
+#endif

include/igl/cotmatrix_intrinsic.h

@@ -0,0 +1,40 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_COTMATRIX_INTRINSIC_H
+#define IGL_COTMATRIX_INTRINSIC_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl 
+{
+  // Constructs the cotangent stiffness matrix (discrete laplacian) for a given
+  // mesh with faces F and edge lengths l.
+  //
+  // Inputs:
+  //   l  #F by 3 list of (half-)edge lengths
+  //   F  #F by 3 list of face indices into some (not necessarily
+  //     determined/embedable) list of vertex positions V. It is assumed #V ==
+  //     F.maxCoeff()+1
+  // Outputs:
+  //   L  #V by #V sparse Laplacian matrix
+  //
+  // See also: cotmatrix, intrinsic_delaunay_cotmatrix
+  template <typename Derivedl, typename DerivedF, typename Scalar>
+  IGL_INLINE void cotmatrix_intrinsic(
+    const Eigen::MatrixBase<Derivedl> & l, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    Eigen::SparseMatrix<Scalar>& L);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cotmatrix.cpp"
+#endif
+
+#endif

include/igl/cumprod.cpp

@@ -0,0 +1,76 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "cumprod.h"
+#include <numeric>
+#include <iostream>
+
+template <typename DerivedX, typename DerivedY>
+IGL_INLINE void igl::cumprod(
+  const Eigen::MatrixBase<DerivedX > & X,
+  const int dim,
+  Eigen::PlainObjectBase<DerivedY > & Y)
+{
+  using namespace Eigen;
+  using namespace std;
+  Y.resizeLike(X);
+  // get number of columns (or rows)
+  int num_outer = (dim == 1 ? X.cols() : X.rows() );
+  // get number of rows (or columns)
+  int num_inner = (dim == 1 ? X.rows() : X.cols() );
+  // This has been optimized so that dim = 1 or 2 is roughly the same cost.
+  // (Optimizations assume ColMajor order)
+  if(dim == 1)
+  {
+#pragma omp parallel for
+    for(int o = 0;o<num_outer;o++)
+    {
+      typename DerivedX::Scalar prod = 1;
+      for(int i = 0;i<num_inner;i++)
+      {
+        prod *= X(i,o);
+        Y(i,o) = prod;
+      }
+    }
+  }else
+  {
+    for(int i = 0;i<num_inner;i++)
+    {
+      // Notice that it is *not* OK to put this above the inner loop
+      // Though here it doesn't seem to pay off...
+//#pragma omp parallel for
+      for(int o = 0;o<num_outer;o++)
+      {
+        if(i == 0)
+        {
+          Y(o,i) = X(o,i);
+        }else
+        {
+          Y(o,i) = Y(o,i-1) * X(o,i);
+        }
+      }
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cumprod<Eigen::Matrix<double, 4, 1, 0, 4, 1>, Eigen::Matrix<double, 4, 1, 0, 4, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> >&);
+// generated by autoexplicit.sh
+template void igl::cumprod<Eigen::Matrix<double, 1, 4, 1, 1, 4>, Eigen::Matrix<double, 1, 4, 1, 1, 4> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 4, 1, 1, 4> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 4, 1, 1, 4> >&);
+template void igl::cumprod<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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::cumprod<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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::cumprod<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);
+template void igl::cumprod<Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1>, Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1> >(Eigen::MatrixBase<Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1> >&);
+template void igl::cumprod<Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1>, Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1> >&);
+template void igl::cumprod<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&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#ifdef WIN32
+template void igl::cumprod<class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>, class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>>(class Eigen::MatrixBase<class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>> &);
+template void igl::cumprod<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>, class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>>(class Eigen::MatrixBase<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>> &);
+#endif
+#endif

include/igl/cumprod.h

@@ -0,0 +1,44 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_CUMPROD_H
+#define IGL_CUMPROD_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Computes a cumulative product of the columns of X, like matlab's `cumprod`.
+  //
+  // Templates:
+  //   DerivedX  Type of matrix X
+  //   DerivedY  Type of matrix Y
+  // Inputs:
+  //   X  m by n Matrix to be cumulatively multiplied.
+  //   dim  dimension to take cumulative product (1 or 2)
+  // Output:
+  //   Y  m by n Matrix containing cumulative product.
+  //
+  template <typename DerivedX, typename DerivedY>
+  IGL_INLINE void cumprod(
+    const Eigen::MatrixBase<DerivedX > & X,
+    const int dim,
+    Eigen::PlainObjectBase<DerivedY > & Y);
+  //template <typename DerivedX, typename DerivedY>
+  //IGL_INLINE void cumprod(
+  //  const Eigen::MatrixBase<DerivedX > & X,
+  //  Eigen::PlainObjectBase<DerivedY > & Y);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cumprod.cpp"
+#endif
+
+#endif
+
+

include/igl/cumsum.cpp

@@ -59,6 +59,10 @@ IGL_INLINE void igl::cumsum(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cumsum<Eigen::Matrix<double, 4, 1, 0, 4, 1>, Eigen::Matrix<double, 4, 1, 0, 4, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> >&);
+// generated by autoexplicit.sh
+template void igl::cumsum<Eigen::Matrix<double, 1, 4, 1, 1, 4>, Eigen::Matrix<double, 1, 4, 1, 1, 4> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 4, 1, 1, 4> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 4, 1, 1, 4> >&);
 template void igl::cumsum<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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::cumsum<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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::cumsum<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);

include/igl/doublearea.cpp

@@ -231,13 +231,12 @@ IGL_INLINE void igl::doublearea_quad(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
+template void igl::doublearea<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::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::doublearea<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(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, 1, 0, -1, 1> >&);
+template void igl::doublearea<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::Matrix<double, -1, -1, 1, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 1, 0, 1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 1, 0, 1, 1> >&);

include/igl/edge_collapse_is_valid.cpp

@@ -49,7 +49,7 @@ IGL_INLINE bool igl::edge_collapse_is_valid(
       const Eigen::MatrixXi & EI) 
     {
       vector<int> N,uN;
-      vector<int> V2Fe = circulation(e, ccw,F,E,EMAP,EF,EI);
+      vector<int> V2Fe = circulation(e, ccw,EMAP,EF,EI);
       for(auto f : V2Fe)
       {
         N.push_back(F(f,0));

include/igl/edge_exists_near.cpp

@@ -0,0 +1,96 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "edge_exists_near.h"
+
+template <
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename uE2EType,
+  typename Index>
+IGL_INLINE bool igl::edge_exists_near(
+  const Eigen::MatrixBase<DeriveduE> & uE,
+  const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+  const std::vector<std::vector< uE2EType> > & uE2E,
+  const Index & a,
+  const Index & b,
+  const Index & uei)
+{
+  std::vector<Index> face_queue;
+  face_queue.reserve(32);
+  std::vector<Index> pushed;
+  // 32 is faster than 8
+  pushed.reserve(32);
+  assert(a!=b);
+  // starting with the (2) faces incident on e, consider all faces
+  // incident on edges containing either a or b.
+  //
+  // face_queue  Queue containing faces incident on exactly one of a/b
+  // Using a vector seems mildly faster
+  assert(EMAP.size()%3 == 0);
+  const Index num_faces = EMAP.size()/3;
+  const Index f1 = uE2E[uei][0]%num_faces;
+  const Index f2 = uE2E[uei][1]%num_faces;
+  // map is faster than unordered_map here, and vector + brute force
+  // is_member check is even faster
+  face_queue.push_back(f1);
+  pushed.push_back(f1);
+  face_queue.push_back(f2);
+  pushed.push_back(f2);
+  while(!face_queue.empty())
+  {
+    const Index f = face_queue.back();
+    face_queue.pop_back();
+    // consider each edge of this face
+    for(int c = 0;c<3;c++)
+    {
+      // Unique edge id
+      const Index uec = EMAP(c*num_faces+f);
+      const Index s = uE(uec,0);
+      const Index d = uE(uec,1);
+      const bool ona = s == a || d == a;
+      const bool onb = s == b || d == b;
+      // Is this the edge we're looking for?
+      if(ona && onb)
+      {
+        return true;
+      }
+      // not incident on either?
+      if(!ona && !onb)
+      {
+        continue;
+      }
+      // loop over all incident half-edges
+      for(const auto & he : uE2E[uec])
+      {
+        // face of this he
+        const Index fhe = he%num_faces;
+        bool already_pushed = false;
+        for(const auto & fp : pushed)
+        {
+          if(fp == fhe)
+          {
+            already_pushed = true;
+            break;
+          }
+        }
+        if(!already_pushed)
+        {
+          pushed.push_back(fhe);
+          face_queue.push_back(fhe);
+        }
+      }
+    }
+  }
+  return false;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::edge_exists_near<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, int const&, int const&, int const&);
+#endif 

include/igl/edge_exists_near.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_EDGE_EXISTS_NEAR_H
+#define IGL_EDGE_EXISTS_NEAR_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // Does edge (a,b) exist in the edges of all faces incident on
+  // existing unique edge uei.
+  //
+  // Inputs:
+  //   uE  #uE by 2 list of unique undirected edges
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //     undirected edge
+  //   uE2E  #uE list of lists of indices into E of coexisting edges
+  //   E  #F*3 by 2 list of half-edges
+  //   a  1st end-point of query edge
+  //   b  2nd end-point of query edge
+  //   uei  index into uE/uE2E of unique edge
+  // Returns true if edge exists near uei.
+  //
+  // See also: unique_edge_map
+  template <
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType,
+    typename Index>
+  IGL_INLINE bool edge_exists_near(
+    const Eigen::MatrixBase<DeriveduE> & uE,
+    const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+    const std::vector<std::vector< uE2EType> > & uE2E,
+    const Index & a,
+    const Index & b,
+    const Index & uei);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "edge_exists_near.cpp"
+#endif
+#endif
+

include/igl/edge_flaps.cpp

@@ -12,14 +12,14 @@
 
 IGL_INLINE void igl::edge_flaps(
   const Eigen::MatrixXi & F,
-  const Eigen::MatrixXi & E,
+  const Eigen::MatrixXi & uE,
   const Eigen::VectorXi & EMAP,
   Eigen::MatrixXi & EF,
   Eigen::MatrixXi & EI)
 {
   // Initialize to boundary value
-  EF.setConstant(E.rows(),2,-1);
-  EI.setConstant(E.rows(),2,-1);
+  EF.setConstant(uE.rows(),2,-1);
+  EI.setConstant(uE.rows(),2,-1);
   // loop over all faces
   for(int f = 0;f<F.rows();f++)
   {
@@ -29,13 +29,13 @@ IGL_INLINE void igl::edge_flaps(
       // get edge id
       const int e = EMAP(v*F.rows()+f);
       // See if this is left or right flap w.r.t. edge orientation
-      if( F(f,(v+1)%3) == E(e,0) && F(f,(v+2)%3) == E(e,1))
+      if( F(f,(v+1)%3) == uE(e,0) && F(f,(v+2)%3) == uE(e,1))
       {
         EF(e,0) = f;
         EI(e,0) = v;
       }else
       {
-        assert(F(f,(v+1)%3) == E(e,1) && F(f,(v+2)%3) == E(e,0));
+        assert(F(f,(v+1)%3) == uE(e,1) && F(f,(v+2)%3) == uE(e,0));
         EF(e,1) = f;
         EI(e,1) = v;
       }
@@ -45,16 +45,16 @@ IGL_INLINE void igl::edge_flaps(
 
 IGL_INLINE void igl::edge_flaps(
   const Eigen::MatrixXi & F,
-  Eigen::MatrixXi & E,
+  Eigen::MatrixXi & uE,
   Eigen::VectorXi & EMAP,
   Eigen::MatrixXi & EF,
   Eigen::MatrixXi & EI)
 {
   Eigen::MatrixXi allE;
   std::vector<std::vector<int> > uE2E;
-  igl::unique_edge_map(F,allE,E,EMAP,uE2E);
+  igl::unique_edge_map(F,allE,uE,EMAP,uE2E);
   // Const-ify to call overload
-  const auto & cE = E;
+  const auto & cuE = uE;
   const auto & cEMAP = EMAP;
-  return edge_flaps(F,cE,cEMAP,EF,EI);
+  return edge_flaps(F,cuE,cEMAP,EF,EI);
 }

include/igl/edge_flaps.h

@@ -16,15 +16,17 @@ namespace igl
   //
   // Inputs:
   //   F  #F by 3 list of face indices
-  //   E  #E by 2 list of edge indices into V.
-  //   EMAP #F*3 list of indices into E, mapping each directed edge to unique
-  //     unique edge in E
+  //   uE  #uE by 2 list of edge indices into V.
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //     unique edge in uE
   // Outputs:
   //   EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
   //     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
   //     e=(j->i)
   //   EI  #E by 2 list of edge flap corners (see above).
   //
+  // See also: unique_edge_map
+  //
   // TODO: This seems to be a duplicate of edge_topology.h
   // igl::edge_topology(V,F,etEV,etFE,etEF);
   // igl::edge_flaps(F,efE,efEMAP,efEF,efEI);
@@ -34,14 +36,14 @@ namespace igl
   // all(efEMAP(sub2ind(size(F),repmat(1:size(F,1),3,1)',repmat([1 2 3],size(F,1),1))) == etFE(:,[2 3 1]))
   IGL_INLINE void edge_flaps(
     const Eigen::MatrixXi & F,
-    const Eigen::MatrixXi & E,
+    const Eigen::MatrixXi & uE,
     const Eigen::VectorXi & EMAP,
     Eigen::MatrixXi & EF,
     Eigen::MatrixXi & EI);
   // Only faces as input
   IGL_INLINE void edge_flaps(
     const Eigen::MatrixXi & F,
-    Eigen::MatrixXi & E,
+    Eigen::MatrixXi & uE,
     Eigen::VectorXi & EMAP,
     Eigen::MatrixXi & EF,
     Eigen::MatrixXi & EI);

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, 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> >&);
+// generated by autoexplicit.sh
 template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > 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<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);

include/igl/flip_edge.cpp

@@ -25,6 +25,14 @@ IGL_INLINE void igl::flip_edge(
   typedef typename DerivedF::Scalar Index;
   const size_t num_faces = F.rows();
   assert(F.cols() == 3);
+  // Edge to flip [v1,v2] --> [v3,v4]
+  // Before:
+  // F(f1,:) = [v1,v2,v4] // in some cyclic order
+  // F(f2,:) = [v1,v3,v2] // in some cyclic order
+  // After: 
+  // F(f1,:) = [v1,v3,v4] // in *this* order 
+  // F(f2,:) = [v2,v4,v3] // in *this* order
+  //
   //          v1                 v1
   //          /|\                / \
   //         / | \              /f1 \
@@ -155,5 +163,6 @@ template void igl::flip_edge<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matri
 template void igl::flip_edge<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>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, unsigned long);
 #ifdef WIN32
 template void igl::flip_edge<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>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, unsigned __int64);
+template void igl::flip_edge<class Eigen::Matrix<int,-1,-1,0,-1,-1>,class Eigen::Matrix<int,-1,2,0,-1,2>,class Eigen::Matrix<int,-1,2,0,-1,2>,class Eigen::Matrix<int,-1,1,0,-1,1>,int>(class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,-1,0,-1,-1> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,2,0,-1,2> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,2,0,-1,2> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,1,0,-1,1> > &,class std::vector<class std::vector<int,class std::allocator<int> >,class std::allocator<class std::vector<int,class std::allocator<int> > > > &,unsigned __int64);
 #endif
 #endif

include/igl/grad.cpp

@@ -15,10 +15,12 @@
 #include "doublearea.h"
 
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE void grad_tet(const Eigen::PlainObjectBase<DerivedV>&V,
-                     const Eigen::PlainObjectBase<DerivedF>&T,
-                            Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
-                            bool uniform) {
+IGL_INLINE void grad_tet(
+  const Eigen::MatrixBase<DerivedV>&V,
+  const Eigen::MatrixBase<DerivedF>&T,
+  Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+  bool uniform) 
+{
   using namespace Eigen;
   assert(T.cols() == 4);
   const int n = V.rows(); int m = T.rows();
@@ -116,15 +118,22 @@ IGL_INLINE void grad_tet(const Eigen::PlainObjectBase<DerivedV>&V,
 }
 
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE void grad_tri(const Eigen::PlainObjectBase<DerivedV>&V,
-                     const Eigen::PlainObjectBase<DerivedF>&F,
-                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
-                    bool uniform)
+IGL_INLINE void grad_tri(
+  const Eigen::MatrixBase<DerivedV>&V,
+  const Eigen::MatrixBase<DerivedF>&F,
+  Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+  bool uniform)
 {
+  // Number of faces
+  const int m = F.rows();
+  // Number of vertices
+  const int nv = V.rows();
+  // Number of dimensions
+  const int dims = V.cols();
   Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3>
-    eperp21(F.rows(),3), eperp13(F.rows(),3);
+    eperp21(m,3), eperp13(m,3);
 
-  for (int i=0;i<F.rows();++i)
+  for (int i=0;i<m;++i)
   {
     // renaming indices of vertices of triangles for convenience
     int i1 = F(i,0);
@@ -132,10 +141,14 @@ IGL_INLINE void grad_tri(const Eigen::PlainObjectBase<DerivedV>&V,
     int i3 = F(i,2);
 
     // #F x 3 matrices of triangle edge vectors, named after opposite vertices
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v32 = V.row(i3) - V.row(i2);
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v13 = V.row(i1) - V.row(i3);
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v21 = V.row(i2) - V.row(i1);
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> n = v32.cross(v13);
+    typedef Eigen::Matrix<typename DerivedV::Scalar, 1, 3> RowVector3S;
+    RowVector3S v32 = RowVector3S::Zero(1,3);
+    RowVector3S v13 = RowVector3S::Zero(1,3);
+    RowVector3S v21 = RowVector3S::Zero(1,3);
+    v32.head(V.cols()) = V.row(i3) - V.row(i2);
+    v13.head(V.cols()) = V.row(i1) - V.row(i3);
+    v21.head(V.cols()) = V.row(i2) - V.row(i1);
+    RowVector3S n = v32.cross(v13);
     // area of parallelogram is twice area of triangle
     // area of parallelogram is || v1 x v2 ||
     // This does correct l2 norm of rows, so that it contains #F list of twice
@@ -174,70 +187,46 @@ IGL_INLINE void grad_tri(const Eigen::PlainObjectBase<DerivedV>&V,
     eperp13.row(i) *= norm13 / dblA;
   }
 
-  std::vector<int> rs;
-  rs.reserve(F.rows()*4*3);
-  std::vector<int> cs;
-  cs.reserve(F.rows()*4*3);
-  std::vector<double> vs;
-  vs.reserve(F.rows()*4*3);
-
-  // row indices
-  for(int r=0;r<3;r++)
+  // create sparse gradient operator matrix
+  G.resize(dims*m,nv);
+  std::vector<Eigen::Triplet<typename DerivedV::Scalar> > Gijv;
+  Gijv.reserve(4*dims*m);
+  for(int f = 0;f<F.rows();f++)
   {
-    for(int j=0;j<4;j++)
+    for(int d = 0;d<dims;d++)
     {
-      for(int i=r*F.rows();i<(r+1)*F.rows();i++) rs.push_back(i);
+      Gijv.emplace_back(f+d*m,F(f,1), eperp13(f,d));
+      Gijv.emplace_back(f+d*m,F(f,0),-eperp13(f,d));
+      Gijv.emplace_back(f+d*m,F(f,2), eperp21(f,d));
+      Gijv.emplace_back(f+d*m,F(f,0),-eperp21(f,d));
     }
   }
-
-  // column indices
-  for(int r=0;r<3;r++)
-  {
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,1));
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,2));
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));
-  }
-
-  // values
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,2));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,2));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,2));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,2));
-
-  // create sparse gradient operator matrix
-  G.resize(3*F.rows(),V.rows());
-  std::vector<Eigen::Triplet<typename DerivedV::Scalar> > triplets;
-  for (int i=0;i<(int)vs.size();++i)
-  {
-    triplets.push_back(Eigen::Triplet<typename DerivedV::Scalar>(rs[i],cs[i],vs[i]));
-  }
-  G.setFromTriplets(triplets.begin(), triplets.end());
+  G.setFromTriplets(Gijv.begin(), Gijv.end());
 }
 
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
-                     const Eigen::PlainObjectBase<DerivedF>&F,
-                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
-                    bool uniform)
+IGL_INLINE void igl::grad(
+  const Eigen::MatrixBase<DerivedV>&V,
+  const Eigen::MatrixBase<DerivedF>&F,
+  Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+  bool uniform)
 {
   assert(F.cols() == 3 || F.cols() == 4);
-  if (F.cols() == 3)
-    return grad_tri(V,F,G,uniform);
-  if (F.cols() == 4)
-    return grad_tet(V,F,G,uniform);
+  switch(F.cols())
+  {
+    case 3:
+      return grad_tri(V,F,G,uniform);
+    case 4:
+      return grad_tet(V,F,G,uniform);
+    default:
+      assert(false);
+  }
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template void igl::grad<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 0, int>&, bool);
-template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&, bool);
+// generated by autoexplicit.sh
+template void igl::grad<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 2, 0, -1, 2>::Scalar, 0, int>&, bool);
+template void igl::grad<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 0, int>&, bool);
+template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&, bool);
 #endif

include/igl/grad.h

@@ -5,8 +5,8 @@
 // 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_GRAD_MAT_H
-#define IGL_GRAD_MAT_H
+#ifndef IGL_GRAD_H
+#define IGL_GRAD_H
 #include "igl_inline.h"
 
 #include <Eigen/Core>
@@ -33,11 +33,12 @@ namespace igl {
   // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
   // 90 degrees
   //
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void grad(const Eigen::PlainObjectBase<DerivedV>&V,
-                     const Eigen::PlainObjectBase<DerivedF>&F,
-                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
-                    bool uniform = false);
+  template <typename DerivedV, typename DerivedF>
+  IGL_INLINE void grad(
+    const Eigen::MatrixBase<DerivedV>&V,
+    const Eigen::MatrixBase<DerivedF>&F,
+    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+    bool uniform = false);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "grad.cpp"

include/igl/grad_intrinsic.cpp

@@ -0,0 +1,69 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+//
+// 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 "grad_intrinsic.h"
+#include "grad.h"
+
+template <typename Derivedl, typename DerivedF, typename Gtype>
+IGL_INLINE void igl::grad_intrinsic(
+  const Eigen::MatrixBase<Derivedl>&l,
+  const Eigen::MatrixBase<DerivedF>&F,
+  Eigen::SparseMatrix<Gtype> &G)
+{
+  assert(F.cols() ==3 && "Only triangles supported");
+  // number of vertices
+  const int n = F.maxCoeff()+1;
+  // number of faces
+  const int m = F.rows();
+  typedef Eigen::Matrix<Gtype,Eigen::Dynamic,2> MatrixX2S;
+  MatrixX2S V2 = MatrixX2S::Zero(3*m,2);
+  //     1=[x,y]
+  //     /\
+  // l3 /   \ l2
+  //   /      \
+  //  /         \
+  // 2-----------3
+  //       l1
+  //
+  // x = (l2²-l1²-l3²)/(-2*l1)
+  // y = sqrt(l3² - x²)
+  //
+  //
+  // Place 3rd vertex at [l(:,1) 0]
+  V2.block(2*m,0,m,1) = l.col(0);
+  // Place second vertex at [0 0]
+  // Place third vertex at [x y]
+  V2.block(0,0,m,1) = 
+    (l.col(1).cwiseAbs2()-l.col(0).cwiseAbs2()-l.col(2).cwiseAbs2()).array()/
+    (-2.*l.col(0)).array();
+  V2.block(0,1,m,1) = 
+    (l.col(2).cwiseAbs2() - V2.block(0,0,m,1).cwiseAbs2()).array().sqrt();
+  DerivedF F2(F.rows(),F.cols());
+  std::vector<Eigen::Triplet<Gtype> > Pijv;
+  Pijv.reserve(F.size());
+  for(int f = 0;f<m;f++)
+  {
+    for(int c = 0;c<F.cols();c++)
+    {
+      F2(f,c) = f+c*m;
+      Pijv.emplace_back(F2(f,c),F(f,c),1);
+    }
+  }
+  Eigen::SparseMatrix<Gtype> P(m*3,n);
+  P.setFromTriplets(Pijv.begin(),Pijv.end());
+  Eigen::SparseMatrix<Gtype> G2;
+  grad(V2,F2,G2);
+  G = G2*P;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::grad_intrinsic<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&, Eigen::SparseMatrix<double, 0, int>&);
+// generated by autoexplicit.sh
+template void igl::grad_intrinsic<Eigen::Matrix<double, -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::SparseMatrix<double, 0, int>&);
+#endif

include/igl/grad_intrinsic.h

@@ -0,0 +1,35 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+//
+// 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_GRAD_INTRINSIC_H
+#define IGL_GRAD_INTRINSIC_H
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+namespace igl {
+  // GRAD_INTRINSIC Construct an intrinsic gradient operator.
+  //
+  // Inputs:
+  //  l  #F by 3 list of edge lengths
+  //  F  #F by 3 list of triangle indices into some vertex list V
+  // Outputs:
+  //  G  #F*2 by #V gradient matrix: G=[Gx;Gy] where x runs along the 23 edge and
+  //    y runs in the counter-clockwise 90° rotation.
+  template <typename Derivedl, typename DerivedF, typename Gtype>
+  IGL_INLINE void grad_intrinsic(
+    const Eigen::MatrixBase<Derivedl>&l,
+    const Eigen::MatrixBase<DerivedF>&F,
+    Eigen::SparseMatrix<Gtype> &G);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "grad_intrinsic.cpp"
+#endif
+
+#endif
+

include/igl/grid.cpp

@@ -6,6 +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 "grid.h"
+#include <cassert>
 
 template <
   typename Derivedres,
@@ -16,23 +17,30 @@ IGL_INLINE void igl::grid(
 {
   using namespace Eigen;
   typedef typename DerivedGV::Scalar Scalar;
-  GV.resize(res(0)*res(1)*res(2),3);
-  for(int zi = 0;zi<res(2);zi++)
+  GV.resize(res.array().prod(),res.size());
+  const auto lerp = 
+    [&res](const Scalar di, const int d)->Scalar{return di/(Scalar)(res(d)-1);};
+  int gi = 0;
+  Derivedres sub;
+  sub.resizeLike(res);
+  sub.setConstant(0);
+  for(int gi = 0;gi<GV.rows();gi++)
   {
-    const auto lerp = 
-      [&](const Scalar di, const int d)->Scalar{return di/(Scalar)(res(d)-1);};
-    const Scalar z = lerp(zi,2);
-    for(int yi = 0;yi<res(1);yi++)
+    // omg, I'm implementing addition...
+    for(int c = 0;c<res.size()-1;c++)
     {
-      const Scalar y = lerp(yi,1);
-      for(int xi = 0;xi<res(0);xi++)
+      if(sub(c)>=res(c))
       {
-        const Scalar x = lerp(xi,0);
-        const int gi = xi+res(0)*(yi + res(1)*zi);
-        GV.row(gi) = 
-          Eigen::Matrix<Scalar,1,3>(x,y,z);
+        sub(c) = 0;
+        // roll over
+        sub(c+1)++;
       }
     }
+    for(int c = 0;c<res.size();c++)
+    {
+      GV(gi,c) = lerp(sub(c),c);
+    }
+    sub(0)++;
   }
 }
 
@@ -40,6 +48,8 @@ IGL_INLINE void igl::grid(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::grid<Eigen::Matrix<int, 2, 1, 0, 2, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::grid<Eigen::Matrix<int, 1, 3, 1, 1, 3>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
 template void igl::grid<Eigen::Matrix<int, 3, 1, 0, 3, 1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, 3, 1, 0, 3, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);

include/igl/grid.h

@@ -1,6 +1,7 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
 // 
 // 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 
@@ -15,9 +16,9 @@ namespace igl
   // `igl::marching_cubes`
   //
   // Inputs:
-  //   res  3-long list of number of vertices along the x y and z dimensions
+  //   res  #res list of number of vertices along each dimension
   // Outputs:
-  //   GV  res(0)*res(1)*res(2) by 3 list of mesh vertex positions.
+  //   GV  res.array().prod() by #res list of mesh vertex positions.
   //   
   template <
     typename Derivedres,

include/igl/heat_geodesics.cpp

@@ -0,0 +1,159 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+//
+// 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 "heat_geodesics.h"
+#include "grad.h"
+#include "doublearea.h"
+#include "cotmatrix.h"
+#include "intrinsic_delaunay_cotmatrix.h"
+#include "massmatrix.h"
+#include "massmatrix_intrinsic.h"
+#include "grad_intrinsic.h"
+#include "boundary_facets.h"
+#include "unique.h"
+#include "slice.h"
+#include "avg_edge_length.h"
+
+
+template < typename DerivedV, typename DerivedF, typename Scalar >
+IGL_INLINE bool igl::heat_geodesics_precompute(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  HeatGeodesicsData<Scalar> & data)
+{
+  // default t value
+  const Scalar h = avg_edge_length(V,F);
+  const Scalar t = h*h;
+  return heat_geodesics_precompute(V,F,t,data);
+}
+
+template < typename DerivedV, typename DerivedF, typename Scalar >
+IGL_INLINE bool igl::heat_geodesics_precompute(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Scalar t,
+  HeatGeodesicsData<Scalar> & data)
+{
+  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> VectorXS;
+  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> MatrixXS;
+  Eigen::SparseMatrix<Scalar> L,M;
+  Eigen::Matrix<Scalar,Eigen::Dynamic,3> l_intrinsic;
+  DerivedF F_intrinsic;
+  VectorXS dblA;
+  if(data.use_intrinsic_delaunay)
+  {
+    igl::intrinsic_delaunay_cotmatrix(V,F,L,l_intrinsic,F_intrinsic);
+    igl::massmatrix_intrinsic(l_intrinsic,F_intrinsic,MASSMATRIX_TYPE_DEFAULT,M);
+    igl::doublearea(l_intrinsic,0,dblA);
+    igl::grad_intrinsic(l_intrinsic,F_intrinsic,data.Grad);
+  }else
+  {
+    igl::cotmatrix(V,F,L);
+    igl::massmatrix(V,F,MASSMATRIX_TYPE_DEFAULT,M);
+    igl::doublearea(V,F,dblA);
+    igl::grad(V,F,data.Grad);
+  }
+  // div
+  assert(F.cols() == 3 && "Only triangles are supported");
+  // number of gradient components
+  data.ng = data.Grad.rows() / F.rows();
+  assert(data.ng == 3 || data.ng == 2);
+  data.Div = -0.25*data.Grad.transpose()*dblA.replicate(data.ng,1).asDiagonal();
+
+  Eigen::SparseMatrix<Scalar> Q = M - t*L;
+  Eigen::MatrixXi O;
+  igl::boundary_facets(F,O);
+  igl::unique(O,data.b);
+  {
+    Eigen::SparseMatrix<Scalar> _;
+    if(!igl::min_quad_with_fixed_precompute(
+      Q,Eigen::VectorXi(),_,true,data.Neumann))
+    {
+      return false;
+    }
+    // Only need if there's a boundary
+    if(data.b.size()>0)
+    {
+      if(!igl::min_quad_with_fixed_precompute(Q,data.b,_,true,data.Dirichlet))
+      {
+        return false;
+      }
+    }
+    const Eigen::SparseMatrix<double> Aeq = M.diagonal().transpose().sparseView();
+    if(!igl::min_quad_with_fixed_precompute(
+      (-L*0.5).eval(),Eigen::VectorXi(),Aeq,true,data.Poisson))
+    {
+      return false;
+    }
+  }
+  return true;
+}
+
+template < typename Scalar, typename Derivedgamma, typename DerivedD>
+IGL_INLINE void igl::heat_geodesics_solve(
+  const HeatGeodesicsData<Scalar> & data,
+  const Eigen::MatrixBase<Derivedgamma> & gamma,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  // number of mesh vertices
+  const int n = data.Grad.cols();
+  // Set up delta at gamma
+  DerivedD u0 = DerivedD::Zero(n,1);
+  for(int g = 0;g<gamma.size();g++)
+  {
+    u0(gamma(g)) = 1;
+  }
+  // Neumann solution
+  DerivedD u;
+  igl::min_quad_with_fixed_solve(
+    data.Neumann,u0,DerivedD(),DerivedD(),u);
+  if(data.b.size()>0)
+  {
+    // Average Dirichelt and Neumann solutions
+    DerivedD uD;
+    igl::min_quad_with_fixed_solve(
+      data.Dirichlet,u0,DerivedD::Zero(data.b.size()),DerivedD(),uD);
+    u += uD;
+    u *= 0.5;
+  }
+  DerivedD grad_u = data.Grad*u;
+  const int m = data.Grad.rows()/data.ng;
+  for(int i = 0;i<m;i++)
+  {
+    Scalar norm = 0;
+    for(int d = 0;d<data.ng;d++)
+    {
+      norm += grad_u(d*m+i)*grad_u(d*m+i);
+    }
+    norm = sqrt(norm);
+    if(norm == 0)
+    {
+      for(int d = 0;d<data.ng;d++) { grad_u(d*m+i) = 0; }
+    }else
+    {
+      for(int d = 0;d<data.ng;d++) { grad_u(d*m+i) /= norm; }
+    }
+  }
+  const DerivedD div_X = -data.Div*grad_u;
+  const DerivedD Beq = (DerivedD(1,1)<<0).finished();
+  igl::min_quad_with_fixed_solve(data.Poisson,-2.0*div_X,DerivedD(),Beq,D);
+  DerivedD Dgamma;
+  igl::slice(D,gamma,Dgamma);
+  D.array() -= Dgamma.mean();
+  if(D.mean() < 0)
+  {
+    D = -D;
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::heat_geodesics_solve<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::HeatGeodesicsData<double> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template bool igl::heat_geodesics_precompute<Eigen::Matrix<double, -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&, double, igl::HeatGeodesicsData<double>&);
+template bool igl::heat_geodesics_precompute<Eigen::Matrix<double, -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&, igl::HeatGeodesicsData<double>&);
+#endif

include/igl/heat_geodesics.h

@@ -0,0 +1,69 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+//
+// 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_HEAT_GEODESICS_H
+#define IGL_HEAT_GEODESICS_H
+#include "igl_inline.h"
+#include "min_quad_with_fixed.h"
+#include <Eigen/Sparse>
+#include <Eigen/Sparse>
+namespace igl
+{
+  template <typename Scalar>
+  struct HeatGeodesicsData
+  {
+    // Gradient and Divergence operators
+    Eigen::SparseMatrix<Scalar> Grad,Div;
+    // Number of gradient components
+    int ng;
+    // List of boundary vertex indices
+    Eigen::VectorXi b;
+    // Solvers for Dirichet, Neumann problems
+    min_quad_with_fixed_data<Scalar> Dirichlet,Neumann,Poisson;
+    bool use_intrinsic_delaunay = false;
+  };
+  // Precompute factorized solvers for computing a fast approximation of
+  // geodesic distances on a mesh (V,F). [Crane et al. 2013]
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by 3 list of mesh face indices into V
+  // Outputs:
+  //   data  precomputation data (see heat_geodesics_solve)
+  template < typename DerivedV, typename DerivedF, typename Scalar >
+  IGL_INLINE bool heat_geodesics_precompute(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    HeatGeodesicsData<Scalar> & data);
+  // Inputs:
+  //   t  "heat" parameter (smaller --> more accurate, less stable)
+  template < typename DerivedV, typename DerivedF, typename Scalar >
+  IGL_INLINE bool heat_geodesics_precompute(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Scalar t,
+    HeatGeodesicsData<Scalar> & data);
+  // Compute fast approximate geodesic distances using precomputed data from a
+  // set of selected source vertices (gamma)
+  //
+  // Inputs: 
+  //   data  precomputation data (see heat_geodesics_precompute)
+  //   gamma  #gamma list of indices into V of source vertices
+  // Outputs:
+  //   D  #V list of distances to gamma 
+  template < typename Scalar, typename Derivedgamma, typename DerivedD>
+  IGL_INLINE void heat_geodesics_solve(
+    const HeatGeodesicsData<Scalar> & data,
+    const Eigen::MatrixBase<Derivedgamma> & gamma,
+    Eigen::PlainObjectBase<DerivedD> & D);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "heat_geodesics.cpp"
+#endif
+
+#endif

include/igl/intrinsic_delaunay_cotmatrix.cpp

@@ -0,0 +1,54 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "intrinsic_delaunay_cotmatrix.h"
+#include "edge_lengths.h"
+#include "intrinsic_delaunay_triangulation.h"
+#include "cotmatrix_intrinsic.h"
+#include <cassert>
+
+template <typename DerivedV, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::intrinsic_delaunay_cotmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L)
+{
+  Eigen::Matrix<Scalar, Eigen::Dynamic, 3> l_intrinsic;
+  DerivedF F_intrinsic;
+  return igl::intrinsic_delaunay_cotmatrix(V,F,L,l_intrinsic,F_intrinsic);
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename Scalar,
+  typename Derivedl_intrinsic,
+  typename DerivedF_intrinsic>
+IGL_INLINE void igl::intrinsic_delaunay_cotmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L,
+  Eigen::PlainObjectBase<Derivedl_intrinsic> & l_intrinsic,
+  Eigen::PlainObjectBase<DerivedF_intrinsic> & F_intrinsic)
+{
+  assert(F.cols() == 3 && "Only triangles are supported");
+  Eigen::Matrix<Scalar, Eigen::Dynamic, 3> l;
+  igl::edge_lengths(V,F,l);
+  igl::intrinsic_delaunay_triangulation(l,F,l_intrinsic,F_intrinsic);
+  igl::cotmatrix_intrinsic(l_intrinsic,F_intrinsic,L);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_cotmatrix<Eigen::Matrix<double, -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::SparseMatrix<double, 0, int>&);
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+#endif

include/igl/intrinsic_delaunay_cotmatrix.h

@@ -0,0 +1,51 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_INTRINSIC_DELAUNAY_COTMATRIX_H
+#define IGL_INTRINSIC_DELAUNAY_COTMATRIX_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // INTRINSIC_DELAUNAY_COTMATRIX Computes the discrete cotangent Laplacian of a
+  // mesh after converting it into its intrinsic Delaunay triangulation (see,
+  // e.g., [Fisher et al. 2007].
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by 3 list of mesh elements (triangles or tetrahedra)
+  // Outputs: 
+  //   L  #V by #V cotangent matrix, each row i corresponding to V(i,:)
+  //   l_intrinsic  #F by 3 list of intrinsic edge-lengths used to compute L
+  //   F_intrinsic  #F by 3 list of intrinsic face indices used to compute L
+  //
+  // See also: intrinsic_delaunay_triangulation, cotmatrix, cotmatrix_intrinsic
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename Scalar,
+    typename Derivedl_intrinsic,
+    typename DerivedF_intrinsic>
+  IGL_INLINE void intrinsic_delaunay_cotmatrix(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    Eigen::SparseMatrix<Scalar>& L,
+    Eigen::PlainObjectBase<Derivedl_intrinsic> & l_intrinsic,
+    Eigen::PlainObjectBase<DerivedF_intrinsic> & F_intrinsic);
+  template <typename DerivedV, typename DerivedF, typename Scalar>
+  IGL_INLINE void intrinsic_delaunay_cotmatrix(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    Eigen::SparseMatrix<Scalar>& L);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "intrinsic_delaunay_cotmatrix.cpp"
+#endif
+
+#endif

include/igl/intrinsic_delaunay_triangulation.cpp

@@ -0,0 +1,199 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "intrinsic_delaunay_triangulation.h"
+#include "is_intrinsic_delaunay.h"
+#include "tan_half_angle.h"
+#include "unique_edge_map.h"
+#include "flip_edge.h"
+#include "EPS.h"
+#include "matlab_format.h"
+#include <iostream>
+#include <queue>
+#include <map>
+
+template <
+  typename Derivedl_in,
+  typename DerivedF_in,
+  typename Derivedl,
+  typename DerivedF>
+IGL_INLINE void igl::intrinsic_delaunay_triangulation(
+  const Eigen::MatrixBase<Derivedl_in> & l_in,
+  const Eigen::MatrixBase<DerivedF_in> & F_in,
+  Eigen::PlainObjectBase<Derivedl> & l,
+  Eigen::PlainObjectBase<DerivedF> & F)
+{
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> VectorXI;
+  MatrixX2I E,uE;
+  VectorXI EMAP;
+  std::vector<std::vector<typename DerivedF::Scalar> > uE2E;
+  return intrinsic_delaunay_triangulation(l_in,F_in,l,F,E,uE,EMAP,uE2E);
+}
+
+template <
+  typename Derivedl_in,
+  typename DerivedF_in,
+  typename Derivedl,
+  typename DerivedF,
+  typename DerivedE,
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename uE2EType>
+IGL_INLINE void igl::intrinsic_delaunay_triangulation(
+  const Eigen::MatrixBase<Derivedl_in> & l_in,
+  const Eigen::MatrixBase<DerivedF_in> & F_in,
+  Eigen::PlainObjectBase<Derivedl> & l,
+  Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DeriveduE> & uE,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+  std::vector<std::vector<uE2EType> > & uE2E)
+{
+  igl::unique_edge_map(F_in, E, uE, EMAP, uE2E);
+  // We're going to work in place
+  l = l_in;
+  F = F_in;
+  typedef typename DerivedF::Scalar Index;
+  typedef typename Derivedl::Scalar Scalar;
+  const Index num_faces = F.rows();
+
+  // Vector is faster than queue...
+  std::vector<Index> Q;
+  Q.reserve(uE2E.size());
+  for (size_t uei=0; uei<uE2E.size(); uei++) 
+  {
+    Q.push_back(uei);
+  }
+  // I tried using a "delaunay_since = iter" flag to avoid duplicates, but there
+  // was no speed up.
+
+  while(!Q.empty())
+  {
+#ifdef IGL_INTRINSIC_DELAUNAY_TRIANGULATION_DEBUG
+    // Expensive sanity check
+    {
+      Eigen::Matrix<bool,Eigen::Dynamic,1> inQ(uE2E.size(),1);
+      inQ.setConstant(false);
+      for(const auto uei : Q)
+      {
+        inQ(uei) = true;
+      }
+      for (Index uei=0; uei<uE2E.size(); uei++) 
+      {
+        if(!inQ(uei) && !is_intrinsic_delaunay(l,uE2E,num_faces,uei))
+        {
+          std::cout<<"  "<<uei<<" might never be fixed!"<<std::endl;
+        }
+      }
+    }
+#endif
+
+    const Index uei = Q.back();
+    Q.pop_back();
+    if (uE2E[uei].size() == 2) 
+    {
+      if(!is_intrinsic_delaunay(l,uE2E,num_faces,uei)) 
+      {
+        // update l just before flipping edge
+        //      .        //
+        //     /|\       //
+        //   a/ | \d     //
+        //   /  e  \     //
+        //  /   |   \    //
+        // .----|-f--.   //
+        //  \   |   /    //
+        //   \  |  /     //
+        //   b\α|δ/c     //
+        //     \|/       //
+        //      .        //
+        // Annotated from flip_edge:
+        // Edge to flip [v1,v2] --> [v3,v4]
+        // Before:
+        // F(f1,:) = [v1,v2,v4] // in some cyclic order
+        // F(f2,:) = [v1,v3,v2] // in some cyclic order
+        // After: 
+        // F(f1,:) = [v1,v3,v4] // in *this* order 
+        // F(f2,:) = [v2,v4,v3] // in *this* order
+        //
+        //          v1                 v1
+        //          /|\                / \
+        //        c/ | \b            c/f1 \b
+        //     v3 /f2|f1\ v4  =>  v3 /__f__\ v4
+        //        \  e  /            \ f2  /
+        //        d\ | /a            d\   /a
+        //          \|/                \ /
+        //          v2                 v2
+        //
+        // Compute intrinsic length of oppposite edge
+        assert(uE2E[uei].size() == 2 && "edge should have 2 incident faces");
+        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;
+        assert(c1 < 3);
+        assert(c2 < 3);
+        assert(f1 != f2);
+        const Index v1 = F(f1, (c1+1)%3);
+        const Index v2 = F(f1, (c1+2)%3);
+        const Index v4 = F(f1, c1);
+        const Index v3 = F(f2, c2);
+        assert(F(f2, (c2+2)%3) == v1);
+        assert(F(f2, (c2+1)%3) == v2);
+        assert( std::abs(l(f1,c1)-l(f2,c2)) < igl::EPS<Scalar>() );
+        const Scalar e = l(f1,c1);
+        const Scalar a = l(f1,(c1+1)%3);
+        const Scalar b = l(f1,(c1+2)%3);
+        const Scalar c = l(f2,(c2+1)%3);
+        const Scalar d = l(f2,(c2+2)%3);
+        // tan(α/2)
+        const Scalar tan_a_2= tan_half_angle(a,b,e);
+        // tan(δ/2)
+        const Scalar tan_d_2 = tan_half_angle(d,e,c);
+        // tan((α+δ)/2)
+        const Scalar tan_a_d_2 = (tan_a_2 + tan_d_2)/(1.0-tan_a_2*tan_d_2);
+        // cos(α+δ)
+        const Scalar cos_a_d = 
+          (1.0 - tan_a_d_2*tan_a_d_2)/(1.0+tan_a_d_2*tan_a_d_2);
+        const Scalar f = sqrt(b*b + c*c - 2.0*b*c*cos_a_d);
+        l(f1,0) = f;
+        l(f1,1) = b;
+        l(f1,2) = c;
+        l(f2,0) = f;
+        l(f2,1) = d;
+        l(f2,2) = a;
+        // Important to grab these indices _before_ calling flip_edges (they
+        // will be correct after)
+        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);
+        flip_edge(F, E, uE, EMAP, uE2E, uei);
+        Q.push_back(ue_24);
+        Q.push_back(ue_41);
+        Q.push_back(ue_13);
+        Q.push_back(ue_32);
+      }
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_triangulation<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::Matrix<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_triangulation<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -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, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_triangulation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -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::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, 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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_triangulation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/intrinsic_delaunay_triangulation.h

@@ -0,0 +1,79 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_INTRINSIC_DELAUNAY_TRIANGULATION_H
+#define IGL_INTRINSIC_DELAUNAY_TRIANGULATION_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // INTRINSIC_DELAUNAY_TRIANGULATION Flip edges _intrinsically_ until all are
+  // "intrinsic Delaunay". See "An algorithm for the construction of intrinsic
+  // delaunay triangulations with applications to digital geometry processing"
+  // [Fisher et al. 2007].
+  //
+  // Inputs:
+  //   l_in  #F_in by 3 list of edge lengths (see edge_lengths)
+  //   F_in  #F_in by 3 list of face indices into some unspecified vertex list V
+  // Outputs:
+  //   l  #F by 3 list of edge lengths
+  //   F  #F by 3 list of new face indices. Note: Combinatorially F may contain
+  //     non-manifold edges, duplicate faces and self-loops (e.g., an edge [1,1]
+  //     or a face [1,1,1]). However, the *intrinsic geometry* is still
+  //     well-defined and correct. See [Fisher et al. 2007] Figure 3 and 2nd to
+  //     last paragraph of 1st page. Since F may be "non-eddge-manifold" in the
+  //     usual combinatorial sense, it may be useful to call the more verbose
+  //     overload below if disentangling edges will be necessary later on.
+  //     Calling unique_edge_map on this F will give a _different_ result than
+  //     those outputs.
+  //
+  // See also: is_intrinsic_delaunay
+  template <
+    typename Derivedl_in,
+    typename DerivedF_in,
+    typename Derivedl,
+    typename DerivedF>
+  IGL_INLINE void intrinsic_delaunay_triangulation(
+    const Eigen::MatrixBase<Derivedl_in> & l_in,
+    const Eigen::MatrixBase<DerivedF_in> & F_in,
+    Eigen::PlainObjectBase<Derivedl> & l,
+    Eigen::PlainObjectBase<DerivedF> & F);
+  // Outputs:
+  //   E  #F*3 by 2 list of all directed edges, such that E.row(f+#F*c) is the
+  //     edge opposite F(f,c)
+  //   uE  #uE by 2 list of unique undirected edges
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //     undirected edge
+  //   uE2E  #uE list of lists of indices into E of coexisting edges
+  //
+  // See also: unique_edge_map
+  template <
+    typename Derivedl_in,
+    typename DerivedF_in,
+    typename Derivedl,
+    typename DerivedF,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType>
+  IGL_INLINE void intrinsic_delaunay_triangulation(
+    const Eigen::MatrixBase<Derivedl_in> & l_in,
+    const Eigen::MatrixBase<DerivedF_in> & F_in,
+    Eigen::PlainObjectBase<Derivedl> & l,
+    Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DeriveduE> & uE,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    std::vector<std::vector<uE2EType> > & uE2E);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "intrinsic_delaunay_triangulation.cpp"
+#endif
+
+#endif

include/igl/is_border_vertex.cpp

@@ -10,14 +10,13 @@
 
 #include "triangle_triangle_adjacency.h"
 
-template <typename DerivedV, typename DerivedF>
+template <typename DerivedF>
 IGL_INLINE std::vector<bool> igl::is_border_vertex(
-    const Eigen::PlainObjectBase<DerivedV> &V,
-    const Eigen::PlainObjectBase<DerivedF> &F)
+  const Eigen::MatrixBase<DerivedF> &F)
 {
   DerivedF FF;
   igl::triangle_triangle_adjacency(F,FF);
-  std::vector<bool> ret(V.rows());
+  std::vector<bool> ret(F.maxCoeff()+1);
   for(unsigned i=0; i<ret.size();++i)
     ret[i] = false;
 
@@ -31,9 +30,17 @@ IGL_INLINE std::vector<bool> igl::is_border_vertex(
   return ret;
 }
 
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE std::vector<bool> igl::is_border_vertex(
+  const Eigen::MatrixBase<DerivedV> &/*V*/,
+  const Eigen::MatrixBase<DerivedF> &F)
+{
+  return igl::is_border_vertex(F);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
-template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
-template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
+template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
+template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
+template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 #endif

include/igl/is_border_vertex.h

@@ -22,13 +22,16 @@ namespace igl
   //   F  #F by 3 list of triangle indices
   // Returns #V vector of bools revealing whether vertices are on boundary
   //
-  // Known Bugs: - does not depend on V
-  //             - assumes mesh is edge manifold
+  // Known Bugs: - assumes mesh is edge manifold
   // 
+  template <typename DerivedF>
+  IGL_INLINE std::vector<bool> is_border_vertex(
+   const Eigen::MatrixBase<DerivedF> &F);
+  // Deprecated:
   template <typename DerivedV, typename DerivedF>
   IGL_INLINE std::vector<bool> is_border_vertex(
-   const Eigen::PlainObjectBase<DerivedV> &V,
-   const Eigen::PlainObjectBase<DerivedF> &F);
+   const Eigen::MatrixBase<DerivedV> &V,
+   const Eigen::MatrixBase<DerivedF> &F);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/is_delaunay.cpp

@@ -1,3 +1,10 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "is_delaunay.h"
 #include "unique_edge_map.h"
 #include <cassert>
@@ -11,20 +18,6 @@ IGL_INLINE void igl::is_delaunay(
   const Eigen::MatrixBase<DerivedF> & F,
   Eigen::PlainObjectBase<DerivedD> & D)
 {
-  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
-  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> VectorXI;
-  MatrixX2I E,uE;
-  VectorXI EMAP;
-  std::vector<std::vector<typename DerivedF::Scalar> > uE2E;
-  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
-  const int num_faces = F.rows();
-  D.setConstant(F.rows(),F.cols(),false);
-  const auto D_at = [&D,&num_faces](const int he)->typename DerivedD::Scalar&
-  {
-    const int f = he%num_faces;
-    const int c = he/num_faces;
-    return D(f,c);
-  };
   typedef typename DerivedV::Scalar Scalar;
   // Should use Shewchuk's predicates instead.
   const auto float_incircle = [](
@@ -33,6 +26,7 @@ IGL_INLINE void igl::is_delaunay(
     const Scalar pc[2],
     const Scalar pd[2])->short
   {
+    // I acknowledge that I am cating to double
     const Eigen::Matrix3d A = (Eigen::Matrix3d(3,3)<<
       pa[0]-pd[0], pa[1]-pd[1],(pa[0]-pd[0])*(pa[0]-pd[0])+(pa[1]-pd[1])*(pa[1]-pd[1]),
       pb[0]-pd[0], pb[1]-pd[1],(pb[0]-pd[0])*(pb[0]-pd[0])+(pb[1]-pd[1])*(pb[1]-pd[1]),
@@ -41,28 +35,23 @@ IGL_INLINE void igl::is_delaunay(
     const Scalar detA = A.determinant();
     return (Scalar(0) < detA) - (detA < Scalar(0));
   };
+
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> VectorXI;
+  MatrixX2I E,uE;
+  VectorXI EMAP;
+  std::vector<std::vector<typename DerivedF::Scalar> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  const int num_faces = F.rows();
+  D.setConstant(F.rows(),F.cols(),false);
   // loop over all unique edges
   for(int ue = 0;ue < uE2E.size(); ue++)
   {
-    bool ue_is_d = false;
-    // Is boundary?
-    switch(uE2E[ue].size())
-    {
-      case 1:
-        ue_is_d = true;
-        break;
-      case 2:
-      {
-        ue_is_d = is_delaunay(V,F,uE2E,float_incircle,ue);
-        break;
-      }
-      default:
-        ue_is_d = false;
-        break;
-    }
+    const bool ue_is_d = is_delaunay(V,F,uE2E,float_incircle,ue);
+    // Set for all instances
     for(int e = 0;e<uE2E[ue].size();e++)
     {
-      D_at(uE2E[ue][e]) = ue_is_d;
+      D( uE2E[ue][e]%num_faces, uE2E[ue][e]/num_faces) = ue_is_d;
     }
   }
 }
@@ -80,6 +69,8 @@ IGL_INLINE bool igl::is_delaunay(
   const InCircle incircle,
   const ueiType uei)
 {
+  if(uE2E[uei].size() == 1) return true;
+  if(uE2E[uei].size() > 2) return false;
   const int num_faces = F.rows();
   typedef typename DerivedV::Scalar Scalar;
   const auto& half_edges = uE2E[uei];
@@ -106,6 +97,8 @@ IGL_INLINE bool igl::is_delaunay(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::is_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
 template void igl::is_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
 template bool igl::is_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, short (*)(double const*, double const*, double const*, double const*), unsigned long>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, short (*)(double const*, double const*, double const*, double const*), unsigned long);

include/igl/is_intrinsic_delaunay.cpp

@@ -0,0 +1,147 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "is_intrinsic_delaunay.h"
+#include "unique_edge_map.h"
+#include "tan_half_angle.h"
+#include "EPS.h"
+#include "matlab_format.h"
+#include "cotmatrix_entries.h"
+#include <iostream>
+#include <cassert>
+template <
+  typename Derivedl,
+  typename DerivedF,
+  typename DerivedD>
+IGL_INLINE void igl::is_intrinsic_delaunay(
+  const Eigen::MatrixBase<Derivedl> & l,
+  const Eigen::MatrixBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> VectorXI;
+  MatrixX2I E,uE;
+  VectorXI EMAP;
+  std::vector<std::vector<typename DerivedF::Scalar> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  return is_intrinsic_delaunay(l,F,uE2E,D);
+}
+
+template <
+  typename Derivedl,
+  typename DerivedF,
+  typename uE2EType,
+  typename DerivedD>
+IGL_INLINE void igl::is_intrinsic_delaunay(
+  const Eigen::MatrixBase<Derivedl> & l,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const std::vector<std::vector<uE2EType> > & uE2E,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  const int num_faces = F.rows();
+  D.setConstant(F.rows(),F.cols(),false);
+  // loop over all unique edges
+  for(int ue = 0;ue < uE2E.size(); ue++)
+  {
+    const bool ue_is_d = is_intrinsic_delaunay(l,uE2E,num_faces,ue);
+    // Set for all instances
+    for(int e = 0;e<uE2E[ue].size();e++)
+    {
+      D( uE2E[ue][e]%num_faces, uE2E[ue][e]/num_faces) = ue_is_d;
+    }
+  }
+};
+
+template <
+  typename Derivedl,
+  typename uE2EType,
+  typename Index>
+IGL_INLINE bool igl::is_intrinsic_delaunay(
+  const Eigen::MatrixBase<Derivedl> & l,
+  const std::vector<std::vector<uE2EType> > & uE2E,
+  const Index num_faces,
+  const Index uei)
+{
+  if(uE2E[uei].size() == 1) return true;
+  if(uE2E[uei].size() > 2) return false;
+  typedef typename Derivedl::Scalar Scalar;
+
+  const auto cot_alpha = [](
+    const Scalar & a,
+    const Scalar & b,
+    const Scalar & c)->Scalar
+  {
+    // Fisher 2007
+    const Scalar t = tan_half_angle(a,b,c);
+    return (1.0-t*t)/(2*t);
+  };
+
+  //      2        //
+  //     /|\       //
+  //   a/ | \d     //
+  //   /  e  \     //
+  //  /   |   \    //
+  //0.α---|-f-β.3  //
+  //  \   |   /    //
+  //   \  |  /     //
+  //   b\ | /c     //
+  //     \|/       //
+  //      .        //
+  //      1
+
+  // Fisher 2007
+  assert(uE2E[uei].size() == 2 && "edge should have 2 incident faces");
+  const Index he1 = uE2E[uei][0];
+  const Index he2 = uE2E[uei][1];
+  const Index f1 = he1%num_faces;
+  const Index c1 = he1/num_faces;
+  const Index f2 = he2%num_faces;
+  const Index c2 = he2/num_faces;
+  assert( std::abs(l(f1,c1)-l(f2,c2)) < igl::EPS<Scalar>());
+  const Scalar e = l(f1,c1);
+  const Scalar a = l(f1,(c1+1)%3);
+  const Scalar b = l(f1,(c1+2)%3);
+  const Scalar c = l(f2,(c2+1)%3);
+  const Scalar d = l(f2,(c2+2)%3);
+  const Scalar w = cot_alpha(e,a,b) + cot_alpha(e,c,d);
+
+  //// Test
+  //{
+  //  Eigen::MatrixXd l(2,3);
+  //  l<<e,a,b,
+  //     d,e,c;
+  //  //Eigen::MatrixXi F(2,3);
+  //  //F<<0,1,2,
+  //  //   1,3,2;
+  //  Eigen::MatrixXd C;
+  //  cotmatrix_entries(l,C);
+  //  if(std::abs(w-(2.0*(C(0,0)+C(1,1))))>1e-10)
+  //  {
+  //    std::cout<<matlab_format(l,"l")<<std::endl;
+  //    std::cout<<w<<std::endl;
+  //    std::cout<<(2.0*(C(0,0)+C(1,1)))<<std::endl;
+  //    std::cout<<w-(2.0*(C(0,0)+C(1,1)))<<std::endl;
+  //    std::cout<<std::endl;
+  //  }
+  //}
+
+  return w >= 0;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, 3, 0, -1, 3>, int, int>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, int, int);
+// generated by autoexplicit.sh
+template bool igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, int>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, int, int);
+// generated by autoexplicit.sh
+template void igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, Eigen::Matrix<bool, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 3, 0, -1, 3> >&);
+#endif

include/igl/is_intrinsic_delaunay.h

@@ -0,0 +1,69 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_IS_INTRINSIC_DELAUNAY_H
+#define IGL_IS_INTRINSIC_DELAUNAY_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // IS_INTRINSIC_DELAUNAY Determine if each edge in the mesh (V,F) is Delaunay.
+  //
+  // Inputs:
+  //   l  #l by dim list of edge lengths
+  //   F  #F by 3 list of triangles indices
+  // Outputs:
+  //   D  #F by 3 list of bools revealing whether edges corresponding 23 31 12
+  //     are locally Delaunay. Boundary edges are by definition Delaunay.
+  //     Non-Manifold edges are by definition not Delaunay.
+  template <
+    typename Derivedl,
+    typename DerivedF,
+    typename DerivedD>
+  IGL_INLINE void is_intrinsic_delaunay(
+    const Eigen::MatrixBase<Derivedl> & l,
+    const Eigen::MatrixBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedD> & D);
+  // Inputs:
+  //   uE2E  #uE list of lists mapping unique edges to (half-)edges
+  template <
+    typename Derivedl,
+    typename DerivedF,
+    typename uE2EType,
+    typename DerivedD>
+  IGL_INLINE void is_intrinsic_delaunay(
+    const Eigen::MatrixBase<Derivedl> & l,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const std::vector<std::vector<uE2EType> > & uE2E,
+    Eigen::PlainObjectBase<DerivedD> & D);
+  // Determine whether a single edge is Delaunay using a provided (extrinsic) incirle
+  // test.
+  //
+  // Inputs:
+  //   l  #l by dim list of edge lengths
+  //   uE2E  #uE list of lists of indices into E of coexisting edges (see
+  //     unique_edge_map)
+  //   num_faces  number of faces (==#F)
+  //   uei  index into uE2E of edge to check
+  // Returns true iff edge is Delaunay
+  template <
+    typename Derivedl,
+    typename uE2EType,
+    typename Index>
+  IGL_INLINE bool is_intrinsic_delaunay(
+    const Eigen::MatrixBase<Derivedl> & l,
+    const std::vector<std::vector<uE2EType> > & uE2E,
+    const Index num_faces,
+    const Index uei);
+
+}
+#ifndef IGL_STATIC_LIBRARY
+#include "is_intrinsic_delaunay.cpp"
+#endif
+#endif
+

include/igl/massmatrix.cpp

@@ -6,6 +6,8 @@
 // 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 "massmatrix.h"
+#include "massmatrix_intrinsic.h"
+#include "edge_lengths.h"
 #include "normalize_row_sums.h"
 #include "sparse.h"
 #include "doublearea.h"
@@ -37,92 +39,18 @@ IGL_INLINE void igl::massmatrix(
   // Not yet supported
   assert(type!=MASSMATRIX_TYPE_FULL);
 
-  Matrix<int,Dynamic,1> MI;
-  Matrix<int,Dynamic,1> MJ;
-  Matrix<Scalar,Dynamic,1> MV;
   if(simplex_size == 3)
   {
     // Triangles
     // edge lengths numbered same as opposite vertices
-    Matrix<Scalar,Dynamic,3> l(m,3);
-    // loop over faces
-    for(int i = 0;i<m;i++)
-    {
-      l(i,0) = (V.row(F(i,1))-V.row(F(i,2))).norm();
-      l(i,1) = (V.row(F(i,2))-V.row(F(i,0))).norm();
-      l(i,2) = (V.row(F(i,0))-V.row(F(i,1))).norm();
-    }
-    Matrix<Scalar,Dynamic,1> dblA;
-    doublearea(l,0.,dblA);
-
-    switch(eff_type)
-    {
-      case MASSMATRIX_TYPE_BARYCENTRIC:
-        // diagonal entries for each face corner
-        MI.resize(m*3,1); MJ.resize(m*3,1); MV.resize(m*3,1);
-        MI.block(0*m,0,m,1) = F.col(0);
-        MI.block(1*m,0,m,1) = F.col(1);
-        MI.block(2*m,0,m,1) = F.col(2);
-        MJ = MI;
-        repmat(dblA,3,1,MV);
-        MV.array() /= 6.0;
-        break;
-      case MASSMATRIX_TYPE_VORONOI:
-        {
-          // diagonal entries for each face corner
-          // http://www.alecjacobson.com/weblog/?p=874
-          MI.resize(m*3,1); MJ.resize(m*3,1); MV.resize(m*3,1);
-          MI.block(0*m,0,m,1) = F.col(0);
-          MI.block(1*m,0,m,1) = F.col(1);
-          MI.block(2*m,0,m,1) = F.col(2);
-          MJ = MI;
-
-          // Holy shit this needs to be cleaned up and optimized
-          Matrix<Scalar,Dynamic,3> cosines(m,3);
-          cosines.col(0) = 
-            (l.col(2).array().pow(2)+l.col(1).array().pow(2)-l.col(0).array().pow(2))/(l.col(1).array()*l.col(2).array()*2.0);
-          cosines.col(1) = 
-            (l.col(0).array().pow(2)+l.col(2).array().pow(2)-l.col(1).array().pow(2))/(l.col(2).array()*l.col(0).array()*2.0);
-          cosines.col(2) = 
-            (l.col(1).array().pow(2)+l.col(0).array().pow(2)-l.col(2).array().pow(2))/(l.col(0).array()*l.col(1).array()*2.0);
-          Matrix<Scalar,Dynamic,3> barycentric = cosines.array() * l.array();
-          normalize_row_sums(barycentric,barycentric);
-          Matrix<Scalar,Dynamic,3> partial = barycentric;
-          partial.col(0).array() *= dblA.array() * 0.5;
-          partial.col(1).array() *= dblA.array() * 0.5;
-          partial.col(2).array() *= dblA.array() * 0.5;
-          Matrix<Scalar,Dynamic,3> quads(partial.rows(),partial.cols());
-          quads.col(0) = (partial.col(1)+partial.col(2))*0.5;
-          quads.col(1) = (partial.col(2)+partial.col(0))*0.5;
-          quads.col(2) = (partial.col(0)+partial.col(1))*0.5;
-
-          quads.col(0) = (cosines.col(0).array()<0).select( 0.25*dblA,quads.col(0));
-          quads.col(1) = (cosines.col(0).array()<0).select(0.125*dblA,quads.col(1));
-          quads.col(2) = (cosines.col(0).array()<0).select(0.125*dblA,quads.col(2));
-
-          quads.col(0) = (cosines.col(1).array()<0).select(0.125*dblA,quads.col(0));
-          quads.col(1) = (cosines.col(1).array()<0).select(0.25*dblA,quads.col(1));
-          quads.col(2) = (cosines.col(1).array()<0).select(0.125*dblA,quads.col(2));
-
-          quads.col(0) = (cosines.col(2).array()<0).select(0.125*dblA,quads.col(0));
-          quads.col(1) = (cosines.col(2).array()<0).select(0.125*dblA,quads.col(1));
-          quads.col(2) = (cosines.col(2).array()<0).select( 0.25*dblA,quads.col(2));
-
-          MV.block(0*m,0,m,1) = quads.col(0);
-          MV.block(1*m,0,m,1) = quads.col(1);
-          MV.block(2*m,0,m,1) = quads.col(2);
-          
-          break;
-        }
-      case MASSMATRIX_TYPE_FULL:
-        assert(false && "Implementation incomplete");
-        break;
-      default:
-        assert(false && "Unknown Mass matrix eff_type");
-    }
-
+    Matrix<Scalar,Dynamic,3> l;
+    igl::edge_lengths(V,F,l);
+    return massmatrix_intrinsic(l,F,type,M);
   }else if(simplex_size == 4)
   {
+    Matrix<int,Dynamic,1> MI;
+    Matrix<int,Dynamic,1> MJ;
+    Matrix<Scalar,Dynamic,1> MV;
     assert(V.cols() == 3);
     assert(eff_type == MASSMATRIX_TYPE_BARYCENTRIC);
     MI.resize(m*4,1); MJ.resize(m*4,1); MV.resize(m*4,1);
@@ -145,12 +73,12 @@ IGL_INLINE void igl::massmatrix(
       MV(i+2*m) = v/4.0;
       MV(i+3*m) = v/4.0;
     }
+    sparse(MI,MJ,MV,n,n,M);
   }else
   {
     // Unsupported simplex size
     assert(false && "Unsupported simplex size");
   }
-  sparse(MI,MJ,MV,n,n,M);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/massmatrix.h

@@ -5,11 +5,10 @@
 // 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_MASSMATRIX_TYPE_H
-#define IGL_MASSMATRIX_TYPE_H
+#ifndef IGL_MASSMATRIX_H
+#define IGL_MASSMATRIX_H
 #include "igl_inline.h"
 
-#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 

include/igl/massmatrix_intrinsic.cpp

@@ -0,0 +1,128 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "massmatrix_intrinsic.h"
+#include "edge_lengths.h"
+#include "normalize_row_sums.h"
+#include "sparse.h"
+#include "doublearea.h"
+#include "repmat.h"
+#include <Eigen/Geometry>
+#include <iostream>
+#include <cassert>
+
+template <typename Derivedl, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::massmatrix_intrinsic(
+  const Eigen::MatrixBase<Derivedl> & l, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  const MassMatrixType type,
+  Eigen::SparseMatrix<Scalar>& M)
+{
+  const int n = F.maxCoeff()+1;
+  return massmatrix_intrinsic(l,F,type,n,M);
+}
+
+template <typename Derivedl, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::massmatrix_intrinsic(
+  const Eigen::MatrixBase<Derivedl> & l, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  const MassMatrixType type,
+  const int n,
+  Eigen::SparseMatrix<Scalar>& M)
+{
+  using namespace Eigen;
+  using namespace std;
+  MassMatrixType eff_type = type;
+  const int m = F.rows();
+  const int simplex_size = F.cols();
+  // Use voronoi of for triangles by default, otherwise barycentric
+  if(type == MASSMATRIX_TYPE_DEFAULT)
+  {
+    eff_type = (simplex_size == 3?MASSMATRIX_TYPE_VORONOI:MASSMATRIX_TYPE_BARYCENTRIC);
+  }
+  assert(F.cols() == 3 && "only triangles supported");
+  Matrix<Scalar,Dynamic,1> dblA;
+  doublearea(l,0.,dblA);
+  Matrix<int,Dynamic,1> MI;
+  Matrix<int,Dynamic,1> MJ;
+  Matrix<Scalar,Dynamic,1> MV;
+
+  switch(eff_type)
+  {
+    case MASSMATRIX_TYPE_BARYCENTRIC:
+      // diagonal entries for each face corner
+      MI.resize(m*3,1); MJ.resize(m*3,1); MV.resize(m*3,1);
+      MI.block(0*m,0,m,1) = F.col(0);
+      MI.block(1*m,0,m,1) = F.col(1);
+      MI.block(2*m,0,m,1) = F.col(2);
+      MJ = MI;
+      repmat(dblA,3,1,MV);
+      MV.array() /= 6.0;
+      break;
+    case MASSMATRIX_TYPE_VORONOI:
+      {
+        // diagonal entries for each face corner
+        // http://www.alecjacobson.com/weblog/?p=874
+        MI.resize(m*3,1); MJ.resize(m*3,1); MV.resize(m*3,1);
+        MI.block(0*m,0,m,1) = F.col(0);
+        MI.block(1*m,0,m,1) = F.col(1);
+        MI.block(2*m,0,m,1) = F.col(2);
+        MJ = MI;
+
+        // Holy shit this needs to be cleaned up and optimized
+        Matrix<Scalar,Dynamic,3> cosines(m,3);
+        cosines.col(0) = 
+          (l.col(2).array().pow(2)+l.col(1).array().pow(2)-l.col(0).array().pow(2))/(l.col(1).array()*l.col(2).array()*2.0);
+        cosines.col(1) = 
+          (l.col(0).array().pow(2)+l.col(2).array().pow(2)-l.col(1).array().pow(2))/(l.col(2).array()*l.col(0).array()*2.0);
+        cosines.col(2) = 
+          (l.col(1).array().pow(2)+l.col(0).array().pow(2)-l.col(2).array().pow(2))/(l.col(0).array()*l.col(1).array()*2.0);
+        Matrix<Scalar,Dynamic,3> barycentric = cosines.array() * l.array();
+        normalize_row_sums(barycentric,barycentric);
+        Matrix<Scalar,Dynamic,3> partial = barycentric;
+        partial.col(0).array() *= dblA.array() * 0.5;
+        partial.col(1).array() *= dblA.array() * 0.5;
+        partial.col(2).array() *= dblA.array() * 0.5;
+        Matrix<Scalar,Dynamic,3> quads(partial.rows(),partial.cols());
+        quads.col(0) = (partial.col(1)+partial.col(2))*0.5;
+        quads.col(1) = (partial.col(2)+partial.col(0))*0.5;
+        quads.col(2) = (partial.col(0)+partial.col(1))*0.5;
+
+        quads.col(0) = (cosines.col(0).array()<0).select( 0.25*dblA,quads.col(0));
+        quads.col(1) = (cosines.col(0).array()<0).select(0.125*dblA,quads.col(1));
+        quads.col(2) = (cosines.col(0).array()<0).select(0.125*dblA,quads.col(2));
+
+        quads.col(0) = (cosines.col(1).array()<0).select(0.125*dblA,quads.col(0));
+        quads.col(1) = (cosines.col(1).array()<0).select(0.25*dblA,quads.col(1));
+        quads.col(2) = (cosines.col(1).array()<0).select(0.125*dblA,quads.col(2));
+
+        quads.col(0) = (cosines.col(2).array()<0).select(0.125*dblA,quads.col(0));
+        quads.col(1) = (cosines.col(2).array()<0).select(0.125*dblA,quads.col(1));
+        quads.col(2) = (cosines.col(2).array()<0).select( 0.25*dblA,quads.col(2));
+
+        MV.block(0*m,0,m,1) = quads.col(0);
+        MV.block(1*m,0,m,1) = quads.col(1);
+        MV.block(2*m,0,m,1) = quads.col(2);
+        
+        break;
+      }
+    case MASSMATRIX_TYPE_FULL:
+      assert(false && "Implementation incomplete");
+      break;
+    default:
+      assert(false && "Unknown Mass matrix eff_type");
+  }
+  sparse(MI,MJ,MV,n,n,M);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::massmatrix_intrinsic<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>&);
+template void igl::massmatrix_intrinsic<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>&);
+template void igl::massmatrix_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::massmatrix_intrinsic<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&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/massmatrix_intrinsic.h

@@ -0,0 +1,56 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_MASSMATRIX_INTRINSIC_H
+#define IGL_MASSMATRIX_INTRINSIC_H
+#include "igl_inline.h"
+#include "massmatrix.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl 
+{
+
+  // Constructs the mass (area) matrix for a given mesh (V,F).
+  //
+  // Inputs:
+  //   l  #l by simplex_size list of mesh edge lengths
+  //   F  #F by simplex_size list of mesh elements (triangles or tetrahedra)
+  //   type  one of the following ints:
+  //     MASSMATRIX_TYPE_BARYCENTRIC  barycentric
+  //     MASSMATRIX_TYPE_VORONOI voronoi-hybrid {default}
+  //     MASSMATRIX_TYPE_FULL full {not implemented}
+  // Outputs: 
+  //   M  #V by #V mass matrix
+  //
+  // See also: adjacency_matrix
+  //
+  template <typename Derivedl, typename DerivedF, typename Scalar>
+  IGL_INLINE void massmatrix_intrinsic(
+    const Eigen::MatrixBase<Derivedl> & l, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    const MassMatrixType type,
+    Eigen::SparseMatrix<Scalar>& M);
+  // Inputs:
+  //   n  number of vertices (>= F.maxCoeff()+1)
+  template <typename Derivedl, typename DerivedF, typename Scalar>
+  IGL_INLINE void massmatrix_intrinsic(
+    const Eigen::MatrixBase<Derivedl> & l, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    const MassMatrixType type,
+    const int n,
+    Eigen::SparseMatrix<Scalar>& M);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "massmatrix_intrinsic.cpp"
+#endif
+
+#endif
+
+

include/igl/matlab_format.cpp

@@ -114,6 +114,10 @@ IGL_INLINE const std::string igl::matlab_format(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<bool, int>, Eigen::Matrix<bool, -1, 3, 0, -1, 3> const> > const igl::matlab_format<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<bool, int>, Eigen::Matrix<bool, -1, 3, 0, -1, 3> const> >(Eigen::DenseBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<bool, int>, Eigen::Matrix<bool, -1, 3, 0, -1, 3> const> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> const> const> > const igl::matlab_format<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> const> const> >(Eigen::DenseBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> const> const> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
 template Eigen::WithFormat<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const> > const igl::matlab_format<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const> >(Eigen::DenseBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 // generated by autoexplicit.sh
 template Eigen::WithFormat<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const igl::matlab_format<Eigen::Matrix<float, 1, 3, 1, 1, 3> >(Eigen::DenseBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);

include/igl/min_quad_with_fixed.cpp

@@ -582,7 +582,8 @@ IGL_INLINE bool igl::min_quad_with_fixed(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
+template bool igl::min_quad_with_fixed_solve<double, Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> > > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template bool igl::min_quad_with_fixed<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, 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<double, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template bool igl::min_quad_with_fixed_precompute<double, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int> const&, bool, igl::min_quad_with_fixed_data<double>&);

include/igl/next_filename.h

@@ -17,7 +17,7 @@ namespace igl
   // Inputs:
   //   prefix  path to containing dir and filename prefix
   //   zeros number of leading zeros as if digit printed with printf
-  //   suffix  suffix of filename and extension (should included dot)
+  //   suffix  suffix of filename and extension (should include dot)
   // Outputs:
   //   next  path to next file
   // Returns true if found, false if exceeding range in zeros

include/igl/per_face_normals.cpp

@@ -105,6 +105,8 @@ IGL_INLINE void igl::per_face_normals_stable(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// Nonsense template
+namespace igl{template<> void per_face_normals<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(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, -1, 0, -1, -1> >&){} }
 // generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
 // generated by autoexplicit.sh

include/igl/simplify_polyhedron.cpp

@@ -55,7 +55,7 @@ IGL_INLINE void igl::simplify_polyhedron(
       const auto vi = E(e,positive);
       const auto vj = E(e,!positive);
       p = V.row(vj);
-      std::vector<int> faces = igl::circulation(e,positive,F,E,EMAP,EF,EI);
+      std::vector<int> faces = igl::circulation(e,positive,EMAP,EF,EI);
       cost = 0;
       for(auto f : faces)
       {

include/igl/slice.cpp

@@ -248,6 +248,8 @@ IGL_INLINE DerivedX igl::slice(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::slice<Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+// generated by autoexplicit.sh
 template void igl::slice<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
 // generated by autoexplicit.sh
 template void igl::slice<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);

include/igl/sort.cpp

@@ -317,6 +317,8 @@ if(!ascending)
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, int, bool, 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::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

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, 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> >&);
+// generated by autoexplicit.sh
 template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > 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<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);

include/igl/tan_half_angle.cpp

@@ -0,0 +1,37 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "tan_half_angle.h"
+#include <cmath>
+template < typename Scalar>
+IGL_INLINE Scalar igl::tan_half_angle(
+  const Scalar & a,
+  const Scalar & b,
+  const Scalar & c)
+{
+  //      .
+  //     /|
+  //   c/ |
+  //   /  |
+  //  /   |
+  // .α   | a
+  //  \   |
+  //   \  |
+  //   b\ |  
+  //     \| 
+  //
+  // tan(α/2)
+  // Fisher 2007
+  return sqrt(((a-b+c)*(a+b-c))/((a+b+c)*(-a+b+c)));
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template double igl::tan_half_angle<double>(double const&, double const&, double const&);
+#endif

include/igl/tan_half_angle.h

@@ -0,0 +1,35 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_TAN_HALF_ANGLE_H
+#define IGL_TAN_HALF_ANGLE_H
+#include "igl_inline.h"
+namespace igl
+{
+  // TAN_HALF_ANGLE Compute the tangent of half of the angle opposite the side
+  // with length a, in a triangle with side lengths (a,b,c).
+  //
+  // Inputs:
+  //   a  scalar edge length of first side of triangle
+  //   b  scalar edge length of second side of triangle
+  //   c  scalar edge length of third side of triangle
+  // Returns tangent of half of the angle opposite side with length a
+  //
+  // See also: is_intrinsic_delaunay
+  template < typename Scalar>
+  IGL_INLINE Scalar tan_half_angle(
+    const Scalar & a,
+    const Scalar & b,
+    const Scalar & c);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "tan_half_angle.cpp"
+#endif
+
+#endif
+

include/igl/triangulated_grid.cpp

@@ -0,0 +1,53 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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 "triangulated_grid.h"
+#include "grid.h"
+#include <cassert>
+
+template <
+  typename XType,
+  typename YType,
+  typename DerivedGV,
+  typename DerivedGF>
+IGL_INLINE void igl::triangulated_grid(
+  const XType & nx,
+  const YType & ny,
+  Eigen::PlainObjectBase<DerivedGV> & GV,
+  Eigen::PlainObjectBase<DerivedGF> & GF)
+{
+  using namespace Eigen;
+  Eigen::Matrix<XType,2,1> res(nx,ny);
+  igl::grid(res,GV);
+  GF.resize((nx-1)*(ny-1)*2,3);
+  for(int y = 0;y<ny-1;y++)
+  {
+    for(int x = 0;x<nx-1;x++)
+    {
+      // index of southwest corner
+      const XType sw = (x  +nx*(y+0));
+      const XType se = (x+1+nx*(y+0));
+      const XType ne = (x+1+nx*(y+1));
+      const XType nw = (x  +nx*(y+1));
+      // Index of first triangle in this square
+      const XType gf = 2*(x+(nx-1)*y);
+      GF(gf+0,0) = sw;
+      GF(gf+0,1) = se;
+      GF(gf+0,2) = nw;
+      GF(gf+1,0) = se;
+      GF(gf+1,1) = ne;
+      GF(gf+1,2) = nw;
+    }
+  }
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::triangulated_grid<int, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(int const&, int const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/triangulated_grid.h

@@ -0,0 +1,38 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// 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_TRIANGULATEGRID_H
+#define IGL_TRIANGULATEGRID_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Create a regular grid of elements (only 2D supported, currently)
+  // Vertex position order is compatible with `igl::grid`
+  //
+  // Inputs:
+  //   nx  number of vertices in the x direction
+  //   ny  number of vertices in the y direction
+  // Outputs:
+  //   GV  nx*ny by 2 list of mesh vertex positions.
+  //   GF  2*(nx-1)*(ny-1) by 3  list of triangle indices
+  template <
+    typename XType,
+    typename YType,
+    typename DerivedGV,
+    typename DerivedGF>
+  IGL_INLINE void triangulated_grid(
+    const XType & nx,
+    const YType & ny,
+    Eigen::PlainObjectBase<DerivedGV> & GV,
+    Eigen::PlainObjectBase<DerivedGF> & GF);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "triangulated_grid.cpp"
+#endif
+#endif 
+

include/igl/volume.cpp

@@ -105,8 +105,12 @@ IGL_INLINE void igl::volume(
        (192.*u*v*w);
   }
 }
+
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// Nonsense template
+namespace igl{ template<> void volume<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(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, 1, 0, -1, 1> >&){} }
 // generated by autoexplicit.sh
 template void igl::volume<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh

tests/include/igl/circulation.cpp

@@ -0,0 +1,69 @@
+#include <test_common.h>
+#include <igl/circulation.h>
+#include <igl/edge_flaps.h>
+#include <igl/unique_edge_map.h>
+#include <igl/matlab_format.h>
+
+TEST(circulation,single_edge)
+{
+  //       7    
+  //     /₆|₇\
+  //   4 - 5 - 6
+  //   |₂/₃|₄\₅|
+  //   1 - 2 - 3
+  //     \₀|₁/  
+  //       0     
+  const Eigen::MatrixXi F = (Eigen::MatrixXi(8,3)<<
+    0,2,1,
+    0,3,2,
+    1,5,4,
+    1,2,5,
+    2,3,5,
+    3,6,5,
+    4,5,7,
+    5,6,7).finished();
+  Eigen::MatrixXi E,uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<int> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  Eigen::MatrixXi EI,EF;
+  {
+    const auto & cuE = uE;
+    const auto & cEMAP = EMAP;
+    igl::edge_flaps(F,cuE,cEMAP,EF,EI);
+  }
+  // Find (2,5) in uE
+  int ei = 0;
+  bool flip = false;
+  for(;ei<E.rows();ei++)
+  {
+    if(uE(ei,0) == 2 && uE(ei,1) == 5){flip=false;break;}
+    if(uE(ei,1) == 2 && uE(ei,0) == 5){flip=true;break;}
+  }
+  Eigen::VectorXi Nccw;
+  igl::circulation(ei,!flip,EMAP,EF,EI,Nccw);
+  Eigen::VectorXi Nccwgt = 
+    (Eigen::VectorXi(6)<<
+     4,
+     5,
+     7,
+     6,
+     2,
+     3).finished();
+  Eigen::VectorXi Ncwgt = 
+    (Eigen::VectorXi(4)<<
+     4,
+     1,
+     0,
+     3).finished();
+  if(flip)
+  {
+    Nccwgt = Nccwgt.reverse().eval();
+    Ncwgt = Ncwgt.reverse().eval();
+  }
+  test_common::assert_eq(Nccw,Nccwgt);
+  Eigen::VectorXi Ncw;
+  igl::circulation(ei, flip,EMAP,EF,EI,Ncw);
+  test_common::assert_eq(Ncw,Ncwgt);
+}
+

tests/include/igl/cotmatrix_entries.cpp

@@ -1,5 +1,7 @@
 #include <test_common.h>
 #include <igl/cotmatrix_entries.h>
+#include <igl/edge_lengths.h>
+#include <igl/EPS.h>
 
 TEST(cotmatrix_entries, simple)
 {
@@ -30,7 +32,7 @@ TEST(cotmatrix_entries, simple)
   for(int f = 0;f<C1.rows();f++)
   {
 #ifdef IGL_EDGE_LENGTHS_SQUARED_H
-      //Hard assert if we have edge_lenght_squared
+      //Hard assert if we have edge_length_squared
       for(int v = 0;v<3;v++)
         if (C1(f,v) > 0.1)
             ASSERT_EQ(0.5, C1(f,v));
@@ -39,7 +41,7 @@ TEST(cotmatrix_entries, simple)
        //All cotangents sum 1.0 for those triangles
        ASSERT_EQ(1.0, C1.row(f).sum());
 #else
-      //Soft assert if we have not edge_lenght_squared
+      //Soft assert if we have not edge_length_squared
       for(int v = 0;v<3;v++)
         if (C1(f,v) > 0.1)
             ASSERT_NEAR(0.5, C1(f,v), epsilon);
@@ -72,7 +74,7 @@ TEST(cotmatrix_entries, simple)
   for(int f = 0;f<C1.rows();f++)
   {    
 #ifdef IGL_EDGE_LENGTHS_SQUARED_H
-      //Hard assert if we have edge_lenght_squared
+      //Hard assert if we have edge_length_squared
       for(int v = 0;v<3;v++)
         if (C1(f,v) > 0.1)
             ASSERT_EQ(0.5, C1(f,v));
@@ -81,7 +83,7 @@ TEST(cotmatrix_entries, simple)
        //All cotangents sum 1.0 for those triangles
        ASSERT_EQ(1.0, C1.row(f).sum());
 #else
-      //Soft assert if we have not edge_lenght_squared
+      //Soft assert if we have not edge_length_squared
       for(int v = 0;v<3;v++)
         if (C1(f,v) > 0.1)
             ASSERT_NEAR(0.5, C1(f,v), epsilon);
@@ -134,3 +136,19 @@ TEST(cotmatrix_entries, simple)
   }
 
 }//TEST(cotmatrix_entries, simple)
+
+TEST(cotmatrix_entries, intrinsic)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //This is a cube of dimensions 1.0x1.0x1.0
+  test_common::load_mesh("cube.obj", V, F);
+  Eigen::MatrixXd Cext,Cint;
+  // compute C extrinsically
+  igl::cotmatrix_entries(V,F,Cext);
+  // compute C intrinsically
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V,F,l);
+  igl::cotmatrix_entries(l,Cint);
+  test_common::assert_near(Cext,Cint,igl::EPS<double>());
+}

tests/include/igl/cotmatrix_intrinsic.cpp

@@ -0,0 +1,90 @@
+#include <test_common.h>
+#include <igl/cotmatrix_intrinsic.h>
+#include <igl/cotmatrix.h>
+#include <igl/edge_lengths.h>
+#include <igl/matlab_format.h>
+#include <igl/EPS.h>
+
+class cotmatrix_intrinsic : public ::testing::TestWithParam<std::string> {};
+
+TEST(cotmatrix_intrinsic, periodic)
+{
+  const Eigen::MatrixXi F = (Eigen::MatrixXi(18,3)<<
+    0,3,1,
+    3,4,1,
+    1,4,2,
+    4,5,2,
+    2,5,0,
+    5,3,0,
+    3,6,4,
+    6,7,4,
+    4,7,5,
+    7,8,5,
+    5,8,3,
+    8,6,3,
+    6,0,7,
+    0,1,7,
+    7,1,8,
+    1,2,8,
+    8,2,6,
+    2,0,6).finished();
+  const Eigen::MatrixXd l = (Eigen::MatrixXd(18,3)<<
+    0.47140452079103168,0.33333333333333331,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333331,
+    0.47140452079103168,0.33333333333333331,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333331,
+    0.47140452079103168,0.33333333333333337,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333337,
+    0.47140452079103168,0.33333333333333331,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333331,
+    0.47140452079103168,0.33333333333333331,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333331,
+    0.47140452079103168,0.33333333333333337,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333337,
+    0.47140452079103168,0.33333333333333331,0.33333333333333337,
+    0.33333333333333337,0.47140452079103168,0.33333333333333331,
+    0.47140452079103168,0.33333333333333331,0.33333333333333337,
+    0.33333333333333337,0.47140452079103168,0.33333333333333331,
+    0.47140452079103173,0.33333333333333337,0.33333333333333337,
+    0.33333333333333337,0.47140452079103173,0.33333333333333337).finished();
+  Eigen::SparseMatrix<double> L;
+  igl::cotmatrix_intrinsic(l,F,L);
+  const Eigen::MatrixXd L_d = L;
+  const Eigen::MatrixXd L_gt = (Eigen::MatrixXd(9,9)<<
+    -4,1,1,1,0,0,1,0,0,
+    1,-4,1,0,1,0,0,1,0,
+    1,1,-4,0,0,1,0,0,1,
+    1,0,0,-4,1,1,1,0,0,
+    0,1,0,1,-4,1,0,1,0,
+    0,0,1,1,1,-4,0,0,1,
+    1,0,0,1,0,0,-4,1,1,
+    0,1,0,0,1,0,1,-4,1,
+    0,0,1,0,0,1,1,1,-4).finished();
+  test_common::assert_near(L_d,L_gt,igl::EPS<double>());
+}
+
+TEST_P(cotmatrix_intrinsic , manifold_meshes)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  test_common::load_mesh(GetParam(), V, F);
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V,F,l);
+  Eigen::SparseMatrix<double> L,Li;
+  igl::cotmatrix(V,F,L);
+  igl::cotmatrix_intrinsic(l,F,Li);
+  // Augh, we don't have assert_near for sparse matrices...
+  // Instead test that bilinear form is near equal for 2 random vectors
+  const Eigen::VectorXd u = Eigen::VectorXd::Random(V.rows(),1);
+  const Eigen::VectorXd v = Eigen::VectorXd::Random(V.rows(),1);
+  const double uv = u.norm()*v.norm();
+  ASSERT_NEAR( u.dot(L*v)/uv, u.dot(Li*v)/uv, igl::EPS<double>());
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+ manifold_meshes,
+ cotmatrix_intrinsic,
+ ::testing::ValuesIn(test_common::manifold_meshes()),
+ test_common::string_test_name
+);

tests/include/igl/cumprod.cpp

@@ -0,0 +1,20 @@
+#include <test_common.h>
+#include <igl/cumprod.h>
+
+TEST(cumprod,col_factorial)
+{
+  Eigen::Vector4d X(1,2,3,4);
+  Eigen::Vector4d Y;
+  igl::cumprod(X,1,Y);
+  Eigen::Vector4d Ygt(1,2,6,24);
+  test_common::assert_eq(Y,Ygt);
+}
+
+TEST(cumprod,row_factorial)
+{
+  Eigen::RowVector4d X(1,2,3,4);
+  Eigen::RowVector4d Y;
+  igl::cumprod(X,2,Y);
+  Eigen::RowVector4d Ygt(1,2,6,24);
+  test_common::assert_eq(Y,Ygt);
+}

tests/include/igl/cumsum.cpp

@@ -0,0 +1,21 @@
+#include <test_common.h>
+#include <igl/cumsum.h>
+
+TEST(cumsum,col)
+{
+  Eigen::Vector4d X(1,2,3,4);
+  Eigen::Vector4d Y;
+  igl::cumsum(X,1,Y);
+  Eigen::Vector4d Ygt(1,3,6,10);
+  test_common::assert_eq(Y,Ygt);
+}
+
+TEST(cumsum,row)
+{
+  Eigen::RowVector4d X(1,2,3,4);
+  Eigen::RowVector4d Y;
+  igl::cumsum(X,2,Y);
+  Eigen::RowVector4d Ygt(1,3,6,10);
+  test_common::assert_eq(Y,Ygt);
+}
+

tests/include/igl/doublearea.cpp

@@ -32,7 +32,7 @@ TEST_P(doublearea, VF_vs_ABC )
 
 INSTANTIATE_TEST_CASE_P
 (
- all_meshes,
+ VF_vs_ABC,
  doublearea,
  ::testing::ValuesIn(test_common::all_meshes()),
  test_common::string_test_name

tests/include/igl/edge_exists_near.cpp

@@ -0,0 +1,29 @@
+#include <test_common.h>
+#include <igl/edge_exists_near.h>
+#include <igl/unique_edge_map.h>
+
+TEST(edge_exists_near,tet)
+{
+  const Eigen::MatrixXi F = (Eigen::MatrixXi(4,3)<<
+     0,1,2,
+     0,2,3,
+     0,3,1,
+     1,3,2).finished();
+  Eigen::MatrixXi E,uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<int> > uE2E;
+  igl::unique_edge_map(F,E,uE,EMAP,uE2E);
+  for(int uei = 0;uei<uE2E.size();uei++)
+  {
+    for(int i = 0;i<4;i++)
+    {
+      for(int j = 0;j<4;j++)
+      {
+        if(i != j)
+        {
+          ASSERT_TRUE(igl::edge_exists_near(uE,EMAP,uE2E,i,j,uei));
+        }
+      }
+    }
+  }
+}

tests/include/igl/grad.cpp

@@ -0,0 +1,28 @@
+#include <test_common.h>
+#include <igl/grad.h>
+#include <igl/triangulated_grid.h>
+#include <igl/edge_lengths.h>
+#include <igl/cotmatrix.h>
+#include <igl/doublearea.h>
+#include <igl/EPS.h>
+
+TEST(grad,laplace_grid)
+{
+  Eigen::MatrixXd V2;
+  Eigen::MatrixXi F;
+  igl::triangulated_grid(3,3,V2,F);
+  Eigen::MatrixXd V = Eigen::MatrixXd::Zero(V2.rows(),3);
+  V.topLeftCorner(V2.rows(),2) = V2;
+  V.col(2) = V.col(1).array() * V.col(0).array() + V.col(1).array();
+  Eigen::SparseMatrix<double> L;
+  igl::cotmatrix(V,F,L);
+  Eigen::SparseMatrix<double> G;
+  igl::grad(V,F,G);
+  Eigen::VectorXd dblA;
+  igl::doublearea(V,F,dblA);
+  Eigen::SparseMatrix<double> GTAG = 
+    G.transpose() * (dblA.replicate(3,1).asDiagonal()) * G;
+  test_common::assert_near(
+    Eigen::MatrixXd(L),Eigen::MatrixXd(-0.5*GTAG),igl::EPS<double>());
+}
+

tests/include/igl/grad_intrinsic.cpp

@@ -0,0 +1,29 @@
+#include <test_common.h>
+#include <igl/grad_intrinsic.h>
+#include <igl/triangulated_grid.h>
+#include <igl/edge_lengths.h>
+#include <igl/cotmatrix.h>
+#include <igl/doublearea.h>
+#include <igl/EPS.h>
+
+TEST(grad_intrinsic,laplace_grid)
+{
+  Eigen::MatrixXd V2;
+  Eigen::MatrixXi F;
+  igl::triangulated_grid(3,3,V2,F);
+  Eigen::MatrixXd V = Eigen::MatrixXd::Zero(V2.rows(),3);
+  V.topLeftCorner(V2.rows(),2) = V2;
+  V.col(2) = V.col(1).array() * V.col(0).array() + V.col(1).array();
+  Eigen::SparseMatrix<double> L;
+  igl::cotmatrix(V,F,L);
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V,F,l);
+  Eigen::SparseMatrix<double> G;
+  igl::grad_intrinsic(l,F,G);
+  Eigen::VectorXd dblA;
+  igl::doublearea(l,0,dblA);
+  Eigen::SparseMatrix<double> GTAG = 
+    G.transpose() * (dblA.replicate(2,1).asDiagonal()) * G;
+  test_common::assert_near(
+    Eigen::MatrixXd(L),Eigen::MatrixXd(-0.5*GTAG),igl::EPS<double>());
+}

tests/include/igl/grid.cpp

@@ -0,0 +1,59 @@
+#include <test_common.h>
+#include <igl/grid.h>
+#include <igl/matlab_format.h>
+
+TEST(grid,3d)
+{
+  Eigen::Vector3i res(3,3,3);
+  Eigen::MatrixXd GV;
+  igl::grid(res,GV);
+  const Eigen::MatrixXd GVgt = 
+    (Eigen::MatrixXd(27,3)<<
+      0,  0,  0,
+    0.5,  0,  0,
+      1,  0,  0,
+      0,0.5,  0,
+    0.5,0.5,  0,
+      1,0.5,  0,
+      0,  1,  0,
+    0.5,  1,  0,
+      1,  1,  0,
+      0,  0,0.5,
+    0.5,  0,0.5,
+      1,  0,0.5,
+      0,0.5,0.5,
+    0.5,0.5,0.5,
+      1,0.5,0.5,
+      0,  1,0.5,
+    0.5,  1,0.5,
+      1,  1,0.5,
+      0,  0,  1,
+    0.5,  0,  1,
+      1,  0,  1,
+      0,0.5,  1,
+    0.5,0.5,  1,
+      1,0.5,  1,
+      0,  1,  1,
+    0.5,  1,  1,
+      1,  1,  1).finished();
+  test_common::assert_eq(GV,GVgt);
+}
+
+TEST(grid,2d)
+{
+  Eigen::Vector2i res(3,3);
+  Eigen::MatrixXd GV;
+  igl::grid(res,GV);
+  const Eigen::MatrixXd GVgt = 
+    (Eigen::MatrixXd(9,2)<<
+      0,  0,
+    0.5,  0,
+      1,  0,
+      0,0.5,
+    0.5,0.5,
+      1,0.5,
+      0,  1,
+    0.5,  1,
+      1,  1).finished();
+  test_common::assert_eq(GV,GVgt);
+}

tests/include/igl/intrinsic_delaunay_cotmatrix.cpp

@@ -0,0 +1,73 @@
+#include <test_common.h>
+#include <igl/intrinsic_delaunay_cotmatrix.h>
+#include <igl/EPS.h>
+#include <igl/triangulated_grid.h>
+#include <igl/is_border_vertex.h>
+
+class intrinsic_delaunay_cotmatrix : public ::testing::TestWithParam<std::string> {};
+
+TEST(intrinsic_delaunay_cotmatrix,skewed_grid)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  const int s = 7;
+  igl::triangulated_grid(s,s,V,F);
+  // Skew against diagonal direction
+  V.col(0) -= 1.1*V.col(1);
+  Eigen::SparseMatrix<double> L;
+  Eigen::MatrixXi F_intrinsic;
+  Eigen::MatrixXd l_intrinsic;
+  igl::intrinsic_delaunay_cotmatrix(V,F,L,l_intrinsic,F_intrinsic);
+  Eigen::VectorXi LI,LJ;
+  Eigen::VectorXd LV;
+  igl::find(L,LI,LJ,LV);
+  const auto is_boundary_edge = [](const int i, const int j, const int s)->bool
+  {
+    const auto is_boundary_vertex = [](const int i, const int s)->bool
+    {
+      return (i<s) || (i>=s*s-s) || (i%s == 0) || (i%s == s-1);
+    };
+    return is_boundary_vertex(i,s) && is_boundary_vertex(j,s);
+  };
+  // Off diagonals should be all non-positive
+  for(int k = 0;k<LI.size();k++)
+  {
+    if(LI(k) != LJ(k) && !is_boundary_edge(LI(k),LJ(k),s))
+    {
+      ASSERT_GT(LV(k),-igl::EPS<double>());
+    }
+  }
+}
+
+TEST_P(intrinsic_delaunay_cotmatrix,manifold_meshes)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  test_common::load_mesh(GetParam(), V, F);
+  Eigen::SparseMatrix<double> L;
+  Eigen::MatrixXi F_intrinsic;
+  Eigen::MatrixXd l_intrinsic;
+  igl::intrinsic_delaunay_cotmatrix(V,F,L,l_intrinsic,F_intrinsic);
+  Eigen::VectorXi LI,LJ;
+  Eigen::VectorXd LV;
+  igl::find(L,LI,LJ,LV);
+
+  const std::vector<bool> is_boundary_vertex = igl::is_border_vertex(F);
+  // Off diagonals should be all non-positive
+  for(int k = 0;k<LI.size();k++)
+  {
+    if(LI(k) != LJ(k) && 
+      !(is_boundary_vertex[LI(k)] && is_boundary_vertex[LJ(k)]))
+    {
+      ASSERT_GT(LV(k),-igl::EPS<double>());
+    }
+  }
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+ manifold_meshes,
+ intrinsic_delaunay_cotmatrix,
+ ::testing::ValuesIn(test_common::manifold_meshes()),
+ test_common::string_test_name
+);

tests/include/igl/intrinsic_delaunay_triangulation.cpp

@@ -0,0 +1,123 @@
+#include <test_common.h>
+#include <igl/intrinsic_delaunay_triangulation.h>
+#include <igl/edge_lengths.h>
+#include <igl/triangulated_grid.h>
+#include <igl/is_delaunay.h>
+#include <igl/is_intrinsic_delaunay.h>
+#include <igl/is_edge_manifold.h>
+#include <igl/unique_simplices.h>
+#include <igl/get_seconds.h>
+#include <igl/matlab_format.h>
+
+TEST(intrinsic_delaunay_triangulation, two_triangles)
+{
+  const Eigen::MatrixXd V = 
+    (Eigen::MatrixXd(4,2)<<
+     0,12,
+     1,0,
+     1,20,
+     2,9).finished();
+  const Eigen::MatrixXi FN = 
+    (Eigen::MatrixXi(2,3)<<
+     0,1,2,
+     2,1,3).finished();
+  Eigen::MatrixXd lN;
+  igl::edge_lengths(V,FN,lN);
+  Eigen::MatrixXd l;
+  Eigen::MatrixXi F;
+  igl::intrinsic_delaunay_triangulation( lN, FN, l, F);
+  Eigen::MatrixXd lext;
+  igl::edge_lengths(V,F,lext);
+  test_common::assert_near(l,lext,1e-10);
+}
+
+TEST(intrinsic_delaunay_triangulation, skewed_grid)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F_in;
+  igl::triangulated_grid(4,4,V,F_in);
+  // Skew against diagonal direction
+  V.col(0) -= 1.5*V.col(1);
+  Eigen::MatrixXd l_in;
+  igl::edge_lengths(V,F_in,l_in);
+  Eigen::MatrixXd l;
+  Eigen::MatrixXi F;
+  igl::intrinsic_delaunay_triangulation( l_in, F_in, l, F);
+  Eigen::MatrixXd lext;
+  igl::edge_lengths(V,F,lext);
+  test_common::assert_near(l,lext,1e-10);
+  Eigen::Matrix<bool,Eigen::Dynamic,3> D;
+  igl::is_delaunay(V,F,D);
+  const Eigen::Matrix<bool,Eigen::Dynamic,3> Dtrue = 
+    Eigen::Matrix<bool,Eigen::Dynamic,3>::Constant(F.rows(),F.cols(),true);
+  test_common::assert_eq(D,Dtrue);
+}
+
+TEST(intrinsic_delaunay_triangulation, peaks)
+{
+  Eigen::MatrixXd V2;
+  Eigen::MatrixXi F_in;
+  igl::triangulated_grid(20,20,V2,F_in);
+  Eigen::MatrixXd V(V2.rows(),3);
+  for(int v=0;v<V.rows();v++)
+  {
+    const auto x = (V2(v,0)-0.5)*6.0;
+    const auto y = (V2(v,1)-0.5)*6.0;
+    // peaks.m
+    const auto z = 3.*(1.-x)*(1.-x)*std::exp(-(x*x) - (y+1.)*(y+1.)) +
+      - 10.*(x/5. - x*x*x - y*y*y*y*y)*std::exp(-x*x-y*y) +
+      - 1./3.*std::exp(-(x+1.)*(x+1.) - y*y);
+    V(v,0) = x;
+    V(v,1) = y;
+    V(v,2) = z;
+  }
+  Eigen::MatrixXd l_in;
+  igl::edge_lengths(V,F_in,l_in);
+  Eigen::MatrixXd l;
+  Eigen::MatrixXi F;
+  Eigen::MatrixXi E,uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<int> > uE2E;
+  igl::intrinsic_delaunay_triangulation(l_in,F_in,l,F,E,uE,EMAP,uE2E);
+  Eigen::Matrix<bool,Eigen::Dynamic,3> D;
+  const Eigen::Matrix<bool,Eigen::Dynamic,3> D_gt = 
+    Eigen::Matrix<bool,Eigen::Dynamic,3>::Constant(F.rows(),F.cols(),true);
+  igl::is_intrinsic_delaunay(l,F,uE2E,D);
+  test_common::assert_eq(D,D_gt);
+}
+
+TEST(intrinsic_delaunay_triangulation,tet)
+{
+  const Eigen::MatrixXd V = (Eigen::MatrixXd(4,3)<<
+    10, 4,7,
+     5, 9,0,
+     8, 8,8,
+     1,10,9).finished();
+  const Eigen::MatrixXi F_in = (Eigen::MatrixXi(4,3)<<
+     0,1,2,
+     0,2,3,
+     0,3,1,
+     1,3,2).finished();
+  const Eigen::Matrix<bool,Eigen::Dynamic,3> D_before = 
+    (Eigen::Matrix<bool,Eigen::Dynamic,3>(4,3)<<
+     1,1,1,
+     1,0,1,
+     1,1,0,
+     1,1,1).finished();
+  Eigen::Matrix<bool,Eigen::Dynamic,3> D;
+  Eigen::MatrixXd l_in;
+  igl::edge_lengths(V,F_in,l_in);
+  igl::is_intrinsic_delaunay(l_in,F_in,D);
+  test_common::assert_eq(D,D_before);
+  Eigen::MatrixXd l;
+  Eigen::MatrixXi F;
+  Eigen::MatrixXi E,uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<int> > uE2E;
+  igl::intrinsic_delaunay_triangulation(l_in,F_in,l,F,E,uE,EMAP,uE2E);
+
+  const Eigen::Matrix<bool,Eigen::Dynamic,3> D_after = 
+    Eigen::Matrix<bool,Eigen::Dynamic,3>::Constant(F.rows(),F.cols(),true);
+  igl::is_intrinsic_delaunay(l,F,uE2E,D);
+  test_common::assert_eq(D,D_after);
+}

tests/include/igl/is_delaunay.cpp

@@ -1,6 +1,5 @@
 #include <test_common.h>
 #include <igl/is_delaunay.h>
-#include <igl/matlab_format.h>
 
 TEST(is_delaunay, two_triangles)
 {

tests/include/igl/is_intrinsic_delaunay.cpp

@@ -0,0 +1,37 @@
+#include <test_common.h>
+#include <igl/is_intrinsic_delaunay.h>
+#include <igl/edge_lengths.h>
+
+TEST(is_intrinsic_delaunay, two_triangles)
+{
+  const Eigen::MatrixXd V = 
+    (Eigen::MatrixXd(4,2)<<
+     0,10,
+     1,0,
+     1,20,
+     2,10).finished();
+  const Eigen::MatrixXi FD = 
+    (Eigen::MatrixXi(2,3)<<
+     0,1,3,
+     0,3,2).finished();
+  Eigen::Matrix<bool,Eigen::Dynamic,Eigen::Dynamic> DD,DN;
+  Eigen::MatrixXd lD;
+  igl::edge_lengths(V,FD,lD);
+  igl::is_intrinsic_delaunay(lD,FD,DD);
+  for(int f=0;f<DD.rows();f++)
+  {
+    for(int c=0;c<DD.cols();c++)
+    {
+      ASSERT_TRUE(DD(f,c));
+    }
+  }
+  const Eigen::MatrixXi FN = 
+    (Eigen::MatrixXi(2,3)<<
+     0,1,2,
+     2,1,3).finished();
+  Eigen::MatrixXd lN;
+  igl::edge_lengths(V,FN,lN);
+  igl::is_intrinsic_delaunay(lN,FN,DN);
+  ASSERT_FALSE(DN(0,0));
+  ASSERT_FALSE(DN(1,2));
+}

tests/include/igl/triangulated_grid.cpp

@@ -0,0 +1,23 @@
+#include <test_common.h>
+#include <igl/triangulated_grid.h>
+#include <igl/doublearea.h>
+
+TEST(triangulated_grid,area)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  const int nx = 4;
+  const int ny = 7;
+  igl::triangulated_grid(nx,ny,V,F);
+  ASSERT_EQ(V.rows(),nx*ny);
+  ASSERT_EQ(F.rows(),2*(nx-1)*(ny-1));
+  Eigen::VectorXd dblA;
+  igl::doublearea(V,F,dblA);
+  ASSERT_NEAR(dblA.array().sum(),2.0,1e-10);
+  const Eigen::VectorXd dblAgt = 
+    Eigen::VectorXd::Constant(
+      2*(nx-1)*(ny-1),
+      1,
+      2.0/double(2*(nx-1)*(ny-1)));
+  test_common::assert_near(dblA,dblAgt,1e-10);
+}

tutorial/716_HeatGeodesics/CMakeLists.txt

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

tutorial/716_HeatGeodesics/main.cpp

@@ -0,0 +1,227 @@
+#include "tutorial_shared_path.h"
+#include <igl/read_triangle_mesh.h>
+#include <igl/triangulated_grid.h>
+#include <igl/heat_geodesics.h>
+#include <igl/unproject_onto_mesh.h>
+#include <igl/avg_edge_length.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/opengl/create_shader_program.h>
+#include <igl/opengl/destroy_shader_program.h>
+#include <iostream>
+
+int main(int argc, char *argv[])
+{
+  // Create the peak height field
+  Eigen::MatrixXi F;
+  Eigen::MatrixXd V;
+  igl::read_triangle_mesh( argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F);
+  
+  // Precomputation
+  igl::HeatGeodesicsData<double> data;
+  double t = std::pow(igl::avg_edge_length(V,F),2);
+  const auto precompute = [&]()
+  {
+    if(!igl::heat_geodesics_precompute(V,F,t,data))
+    {
+      std::cerr<<"Error: heat_geodesics_precompute failed."<<std::endl;
+      exit(EXIT_FAILURE);
+    };
+  };
+  precompute();
+
+  // Initialize white
+  Eigen::MatrixXd C = Eigen::MatrixXd::Constant(V.rows(),3,1);
+  igl::opengl::glfw::Viewer viewer;
+  bool down_on_mesh = false;
+  const auto update = [&]()->bool
+  {
+    int fid;
+    Eigen::Vector3f bc;
+    // Cast a ray in the view direction starting from the mouse position
+    double x = viewer.current_mouse_x;
+    double y = viewer.core.viewport(3) - viewer.current_mouse_y;
+    if(igl::unproject_onto_mesh(Eigen::Vector2f(x,y), viewer.core.view,
+      viewer.core.proj, viewer.core.viewport, V, F, fid, bc))
+    {
+      // 3d position of hit
+      const Eigen::RowVector3d m3 = 
+        V.row(F(fid,0))*bc(0) + V.row(F(fid,1))*bc(1) + V.row(F(fid,2))*bc(2);
+      int cid = 0;
+      Eigen::Vector3d(
+          (V.row(F(fid,0))-m3).squaredNorm(),
+          (V.row(F(fid,1))-m3).squaredNorm(),
+          (V.row(F(fid,2))-m3).squaredNorm()).minCoeff(&cid);
+      const int vid = F(fid,cid);
+      C.row(vid)<<1,0,0;
+      Eigen::VectorXd D = Eigen::VectorXd::Zero(data.Grad.cols(),1);
+      D(vid) = 1;
+      igl::heat_geodesics_solve(data,(Eigen::VectorXi(1,1)<<vid).finished(),D);
+      viewer.data().set_colors((D/D.maxCoeff()).replicate(1,3));
+      return true;
+    }
+    return false;
+  };
+  viewer.callback_mouse_down =
+    [&](igl::opengl::glfw::Viewer& viewer, int, int)->bool
+  {
+    if(update())
+    {
+      down_on_mesh = true;
+      return true;
+    }
+    return false;
+  };
+  viewer.callback_mouse_move = 
+    [&](igl::opengl::glfw::Viewer& viewer, int, int)->bool
+    {
+      if(down_on_mesh)
+      {
+        update();
+        return true;
+      }
+      return false;
+    };
+  viewer.callback_mouse_up =
+    [&down_on_mesh](igl::opengl::glfw::Viewer& viewer, int, int)->bool
+  {
+    down_on_mesh = false;
+    return false;
+  };
+  std::cout<<R"(Usage:
+  [click]  Click on shape to pick new geodesic distance source
+  ,/.      Decrease/increase t by factor of 10.0
+  D,d      Toggle using intrinsic Delaunay discrete differential operators
+
+)";
+
+  viewer.callback_key_pressed = 
+    [&](igl::opengl::glfw::Viewer& /*viewer*/, unsigned int key, int mod)->bool
+  {
+    switch(key)
+    {
+    default:
+      return false;
+    case 'D':
+    case 'd':
+      data.use_intrinsic_delaunay = !data.use_intrinsic_delaunay;
+      std::cout<<(data.use_intrinsic_delaunay?"":"not ")<<
+        "using intrinsic delaunay..."<<std::endl;
+      precompute();
+      update();
+      break;
+    case '.':
+    case ',':
+      t *= (key=='.'?10.0:0.1);
+      precompute();
+      update();
+      std::cout<<"t: "<<t<<std::endl;
+      break;
+    }
+    return true;
+  };
+
+  // Show mesh
+  viewer.data().set_mesh(V, F);
+  viewer.data().set_colors(C);
+  viewer.data().show_lines = false;
+  viewer.launch_init(true,false);
+
+  viewer.data().meshgl.init();
+  igl::opengl::destroy_shader_program(viewer.data().meshgl.shader_mesh);
+
+  {
+    std::string mesh_vertex_shader_string =
+R"(#version 150
+uniform mat4 view;
+uniform mat4 proj;
+uniform mat4 normal_matrix;
+in vec3 position;
+in vec3 normal;
+out vec3 position_eye;
+out vec3 normal_eye;
+in vec4 Ka;
+in vec4 Kd;
+in vec4 Ks;
+in vec2 texcoord;
+out vec2 texcoordi;
+out vec4 Kai;
+out vec4 Kdi;
+out vec4 Ksi;
+
+void main()
+{
+  position_eye = vec3 (view * vec4 (position, 1.0));
+  normal_eye = vec3 (normal_matrix * vec4 (normal, 0.0));
+  normal_eye = normalize(normal_eye);
+  gl_Position = proj * vec4 (position_eye, 1.0); //proj * view * vec4(position, 1.0);
+  Kai = Ka;
+  Kdi = Kd;
+  Ksi = Ks;
+  texcoordi = texcoord;
+})";
+
+    std::string mesh_fragment_shader_string =
+R"(#version 150
+uniform mat4 view;
+uniform mat4 proj;
+uniform vec4 fixed_color;
+in vec3 position_eye;
+in vec3 normal_eye;
+uniform vec3 light_position_eye;
+vec3 Ls = vec3 (1, 1, 1);
+vec3 Ld = vec3 (1, 1, 1);
+vec3 La = vec3 (1, 1, 1);
+in vec4 Ksi;
+in vec4 Kdi;
+in vec4 Kai;
+in vec2 texcoordi;
+uniform sampler2D tex;
+uniform float specular_exponent;
+uniform float lighting_factor;
+uniform float texture_factor;
+out vec4 outColor;
+void main()
+{
+vec3 Ia = La * vec3(Kai);    // ambient intensity
+float ni = 30.0; // number of intervals
+float t = 1.0-round(ni*Kdi.r)/ni; // quantize and reverse 
+vec3 Kdiq = clamp(vec3(2.*t,2.*t-1.,6.*t-5.),0,1); // heat map
+
+vec3 vector_to_light_eye = light_position_eye - position_eye;
+vec3 direction_to_light_eye = normalize (vector_to_light_eye);
+float dot_prod = dot (direction_to_light_eye, normalize(normal_eye));
+float clamped_dot_prod = max (dot_prod, 0.0);
+vec3 Id = Ld * Kdiq * clamped_dot_prod;    // Diffuse intensity
+
+vec3 reflection_eye = reflect (-direction_to_light_eye, normalize(normal_eye));
+vec3 surface_to_viewer_eye = normalize (-position_eye);
+float dot_prod_specular = dot (reflection_eye, surface_to_viewer_eye);
+dot_prod_specular = float(abs(dot_prod)==dot_prod) * max (dot_prod_specular, 0.0);
+float specular_factor = pow (dot_prod_specular, specular_exponent);
+vec3 Kfi = 0.5*vec3(Ksi);
+vec3 Lf = Ls;
+float fresnel_exponent = 2*specular_exponent;
+float fresnel_factor = 0;
+{
+  float NE = max( 0., dot( normalize(normal_eye), surface_to_viewer_eye));
+  fresnel_factor = pow (max(sqrt(1. - NE*NE),0.0), fresnel_exponent);
+}
+vec3 Is = Ls * vec3(Ksi) * specular_factor;    // specular intensity
+vec3 If = Lf * vec3(Kfi) * fresnel_factor;     // fresnel intensity
+vec4 color = vec4(lighting_factor * (If + Is + Id) + Ia +
+  (1.0-lighting_factor) * Kdiq,(Kai.a+Ksi.a+Kdi.a)/3);
+outColor = mix(vec4(1,1,1,1), texture(tex, texcoordi), texture_factor) * color;
+if (fixed_color != vec4(0.0)) outColor = fixed_color;
+})";
+
+    igl::opengl::create_shader_program(
+      mesh_vertex_shader_string,
+      mesh_fragment_shader_string,
+      {},
+      viewer.data().meshgl.shader_mesh);
+  }
+
+  viewer.launch_rendering(true);
+  viewer.launch_shut();
+
+}

tutorial/CMakeLists.txt

@@ -150,4 +150,5 @@ if(TUTORIALS_CHAPTER7)
     add_subdirectory("714_MarchingTets")
   endif()
   add_subdirectory("715_MeshImplicitFunction")
+  add_subdirectory("716_HeatGeodesics")
 endif()