removed regularity check

include/igl/intrinsic_delaunay_triangulation.cpp

@@ -62,95 +62,6 @@ IGL_INLINE void igl::intrinsic_delaunay_triangulation(
   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);
-
-  // 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())
-      {
-        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;
-    };
-
   // Vector is faster than queue...
   std::vector<Index> Q;
   Q.reserve(uE2E.size());
@@ -233,58 +144,43 @@ IGL_INLINE void igl::intrinsic_delaunay_triangulation(
         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)
-        {
-          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);
-        }
+        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);
       }
     }
   }

include/igl/intrinsic_delaunay_triangulation.h

@@ -21,7 +21,15 @@ namespace igl
   //   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
+  //   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 <