refactor and added extra while loop

include/igl/intrinsic_delaunay_triangulation.cpp

@@ -27,219 +27,250 @@ IGL_INLINE void igl::intrinsic_delaunay_triangulation(
   Eigen::PlainObjectBase<Derivedl> & l,
   Eigen::PlainObjectBase<DerivedF> & F)
 {
-  // We're going to work in place
-  l = l_in;
-  F = F_in;
-
   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 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();
 
-  std::vector<Index> face_queue;
-  face_queue.reserve(32);
-  std::vector<Index> pushed;
-  // 32 is faster than 8
-  pushed.reserve(32);
+  bool any_flips = true;
+  while(any_flips)
+  {
+    any_flips = false;
 
-  // Does edge (a,b) exist in the edges of all faces incident on
-  // existing unique edge uei.
-  //
-  // Inputs:
-  //   a  1st end-point of query edge
-  //   b  2nd end-point of query edge
-  //   uei  index into uE/uE2E of unique edge
-  //   uE2E  map from unique edges to half-edges (see unique_edge_map)
-  //   E  #F*3 by 2 list of half-edges
-  //
-  const auto edge_exists_near = 
-    [&](const Index & a,const Index & b,const Index & uei)->bool
-    {
-      face_queue.clear();
-      pushed.clear();
-      assert(a!=b);
-      // Not handling case where (a,b) is edge of face incident on uei
-      // since this can't happen for edge-flipping.
-      assert(a!=uE(uei,0));
-      assert(a!=uE(uei,1));
-      assert(b!=uE(uei,0));
-      assert(b!=uE(uei,1));
-      // 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
-      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())
+    std::vector<Index> face_queue;
+    face_queue.reserve(32);
+    std::vector<Index> pushed;
+    // 32 is faster than 8
+    pushed.reserve(32);
+
+    // Does edge (a,b) exist in the edges of all faces incident on
+    // existing unique edge uei.
+    //
+    // Inputs:
+    //   a  1st end-point of query edge
+    //   b  2nd end-point of query edge
+    //   uei  index into uE/uE2E of unique edge
+    //   uE2E  map from unique edges to half-edges (see unique_edge_map)
+    //   E  #F*3 by 2 list of half-edges
+    //
+    const auto edge_exists_near = 
+      [&](const Index & a,const Index & b,const Index & uei)->bool
       {
-        const Index f = face_queue.back();
-        face_queue.pop_back();
-        // consider each edge of this face
-        for(int c = 0;c<3;c++)
+        face_queue.clear();
+        pushed.clear();
+        assert(a!=b);
+        // Not handling case where (a,b) is edge of face incident on uei
+        // since this can't happen for edge-flipping.
+        assert(a!=uE(uei,0));
+        assert(a!=uE(uei,1));
+        assert(b!=uE(uei,0));
+        assert(b!=uE(uei,1));
+        // 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
+        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())
         {
-          // 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)
+          const Index f = face_queue.back();
+          face_queue.pop_back();
+          // consider each edge of this face
+          for(int c = 0;c<3;c++)
           {
-            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)
+            // 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)
             {
-              if(fp == fhe)
-              {
-                already_pushed = true;
-                break;
-              }
+              return true;
             }
-            if(!already_pushed)
+            // not incident on either?
+            if(!ona && !onb)
             {
-              pushed.push_back(fhe);
-              face_queue.push_back(fhe);
+              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;
-    };
+        return false;
+      };
 
-  // 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);
-  }
+    // 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())
-  {
-    const Index uei = Q.back();
-    Q.pop_back();
-    if (uE2E[uei].size() == 2) 
+    while(!Q.empty())
     {
-      if(!is_intrinsic_delaunay(l,F,uE2E,uei)) 
+      const Index uei = Q.back();
+      Q.pop_back();
+      if (uE2E[uei].size() == 2) 
       {
-        // 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);
-        // From gptoolbox/mesh/flip_edge.m
-        // "If edge-after-flip already  exists then this will create a non-manifold
-        // edge"
-        // Yes, this can happen: e.g., an edge of a tetrahedron."
-        // "If two edges will be the same edge after flip then this will create a
-        // non-manifold edge."
-        // I dont' think this can happen if we flip one at a time. gptoolbox
-        // flips in parallel.
-
-        // Over 50% of the time is spent doing this check...
-        bool flippable = !edge_exists_near(v3,v4,uei);
-        if(flippable)
+        if(!is_intrinsic_delaunay(l,uE2E,num_faces,uei)) 
         {
-          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;
-          flip_edge(F, E, uE, EMAP, uE2E, uei);
-          // append neighbors to back
-          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);
-          Q.push_back(ue_24);
-          Q.push_back(ue_41);
-          Q.push_back(ue_13);
-          Q.push_back(ue_32);
+          // 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);
+          // From gptoolbox/mesh/flip_edge.m
+          // "If edge-after-flip already  exists then this will create a non-manifold
+          // edge"
+          // Yes, this can happen: e.g., an edge of a tetrahedron."
+          // "If two edges will be the same edge after flip then this will create a
+          // non-manifold edge."
+          // I dont' think this can happen if we flip one at a time. gptoolbox
+          // flips in parallel.
+
+          //// Over 50% of the time is spent doing this check...
+          //bool flippable = !edge_exists_near(v3,v4,uei);
+          //if(flippable)
+          if(true)
+          {
+            any_flips = true;
+            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;
+            flip_edge(F, E, uE, EMAP, uE2E, uei);
+            // append neighbors to back
+            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);
+            Q.push_back(ue_24);
+            Q.push_back(ue_41);
+            Q.push_back(ue_13);
+            Q.push_back(ue_32);
+          }
         }
       }
     }
@@ -249,5 +280,7 @@ IGL_INLINE void igl::intrinsic_delaunay_triangulation(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // 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

@@ -34,6 +34,33 @@ namespace igl
     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/igl/is_intrinsic_delaunay.cpp

@@ -9,6 +9,9 @@
 #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,
@@ -25,35 +28,47 @@ IGL_INLINE void igl::is_intrinsic_delaunay(
   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,F,uE2E,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 DerivedF,
   typename uE2EType,
-  typename ueiType>
+  typename Index>
 IGL_INLINE bool igl::is_intrinsic_delaunay(
   const Eigen::MatrixBase<Derivedl> & l,
-  const Eigen::MatrixBase<DerivedF> & F,
   const std::vector<std::vector<uE2EType> > & uE2E,
-  const ueiType uei)
+  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;
-  typedef typename DerivedF::Scalar Index;
 
   const auto cot_alpha = [](
     const Scalar & a,
@@ -65,19 +80,20 @@ IGL_INLINE bool igl::is_intrinsic_delaunay(
     return (1.0-t*t)/(2*t);
   };
 
-  //      .        //
+  //      2        //
   //     /|\       //
   //   a/ | \d     //
   //   /  e  \     //
   //  /   |   \    //
-  // .α---|-f-β.   //
+  //0.α---|-f-β.3  //
   //  \   |   /    //
   //   \  |  /     //
   //   b\ | /c     //
   //     \|/       //
   //      .        //
+  //      1
 
-  const Index num_faces = F.rows();
+  // 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];
@@ -91,16 +107,39 @@ IGL_INLINE bool igl::is_intrinsic_delaunay(
   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 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> >&);
+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 bool igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, 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&, unsigned long);
 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

@@ -29,26 +29,37 @@ namespace igl
     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
-  //   F  #F by 3 list of triangles indices
   //   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 DerivedF,
     typename uE2EType,
-    typename ueiType>
+    typename Index>
   IGL_INLINE bool is_intrinsic_delaunay(
     const Eigen::MatrixBase<Derivedl> & l,
-    const Eigen::MatrixBase<DerivedF> & F,
     const std::vector<std::vector<uE2EType> > & uE2E,
-    const ueiType uei);
+    const Index num_faces,
+    const Index uei);
 
 }
 #ifndef IGL_STATIC_LIBRARY

tests/include/igl/intrinsic_delaunay_triangulation.cpp

@@ -6,6 +6,7 @@
 #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)
@@ -36,8 +37,8 @@ TEST(intrinsic_delaunay_triangulation, skewed_grid)
   Eigen::MatrixXd V;
   Eigen::MatrixXi F_in;
   igl::triangulated_grid(4,4,V,F_in);
-  // Skew
-  V.col(0) += 1.1*V.col(1);
+  // 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;
@@ -57,7 +58,7 @@ TEST(intrinsic_delaunay_triangulation, peaks)
 {
   Eigen::MatrixXd V2;
   Eigen::MatrixXi F_in;
-  igl::triangulated_grid(1000,1000,V2,F_in);
+  igl::triangulated_grid(6,6,V2,F_in);
   Eigen::MatrixXd V(V2.rows(),3);
   for(int v=0;v<V.rows();v++)
   {
@@ -75,52 +76,64 @@ TEST(intrinsic_delaunay_triangulation, peaks)
   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::MatrixXi Fu;
-    Eigen::VectorXi IA,IC;
-    igl::unique_simplices(F,Fu,IA,IC);
-    ASSERT_EQ(F.rows(),Fu.rows());
-  }
-  // Input better have started manifold, otherwise this test doesn't make sense
-  assert(igl::is_edge_manifold(F_in));
-  ASSERT_TRUE(igl::is_edge_manifold(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::MatrixXd lext;
   //igl::edge_lengths(V,F,lext);
   //std::cout<<igl::matlab_format(V,"V")<<std::endl;
   //std::cout<<igl::matlab_format(F_in.array()+1,"F_in")<<std::endl;
   //std::cout<<igl::matlab_format(F.array()+1,"F")<<std::endl;
   //std::cout<<igl::matlab_format(l,"l")<<std::endl;
-}
-
-TEST(intrinsic_delaunay_triangulation,unflippable_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> Dgt = 
-    (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,Dgt);
-  Eigen::MatrixXd l;
-  Eigen::MatrixXi F;
-  igl::intrinsic_delaunay_triangulation( l_in, F_in, l, F);
-  // Nothing should have changed: no edges are flippable.
-  test_common::assert_eq(F,F_in);
-  test_common::assert_eq(l,l_in);
+  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);
+
+  //{
+  //  Eigen::MatrixXi Fu;
+  //  Eigen::VectorXi IA,IC;
+  //  igl::unique_simplices(F,Fu,IA,IC);
+  //  ASSERT_EQ(F.rows(),Fu.rows());
+  //}
+  // Input better have started manifold, otherwise this test doesn't make sense
+  //assert(igl::is_edge_manifold(F_in));
+  //ASSERT_TRUE(igl::is_edge_manifold(F));
+
 }
+
+//// Not sure if this is a good test... Even though the edge will exist twice
+//the intrinsic triangles on either edge are different...
+//TEST(intrinsic_delaunay_triangulation,unflippable_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> Dgt = 
+//    (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,Dgt);
+//  Eigen::MatrixXd l;
+//  Eigen::MatrixXi F;
+//  igl::intrinsic_delaunay_triangulation( l_in, F_in, l, F);
+//  // Nothing should have changed: no edges are flippable.
+//  test_common::assert_eq(F,F_in);
+//  test_common::assert_eq(l,l_in);
+//}