intrinsic delaunay and corresponding test

include/igl/is_delaunay.cpp

@@ -11,20 +11,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 +19,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 +28,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 +62,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];

include/igl/is_intrinsic_delaunay.cpp

@@ -0,0 +1,123 @@
+// 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 <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);
+  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);
+    // 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>
+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)
+{
+  if(uE2E[uei].size() == 1) return true;
+  if(uE2E[uei].size() > 2) return false;
+  typedef typename Derivedl::Scalar Scalar;
+  typedef typename DerivedF::Scalar Index;
+
+  //      .
+  //     /|
+  //   c/ |
+  //   /  |
+  //  /   |
+  // .α   | a
+  //  \   |
+  //   \  |
+  //   b\ |  
+  //     \| 
+  //
+  // tan(α/2)
+  const auto tan_alpha_over_2 = [](
+    const Scalar & a,
+    const Scalar & b,
+    const Scalar & c)->Scalar
+  {
+    // Fisher 2007
+    return sqrt(((a-b+c)*(a+b-c))/((a+b+c)*(-a+b+c)));
+  };
+
+  const auto cot_alpha = [&tan_alpha_over_2](
+    const Scalar & a,
+    const Scalar & b,
+    const Scalar & c)->Scalar
+  {
+    // Fisher 2007
+    const Scalar t = tan_alpha_over_2(a,b,c);
+    return (1.0-t*t)/(2*t);
+  };
+
+  //      .        //
+  //     /|\       //
+  //   a/ | \d     //
+  //   /  e  \     //
+  //  /   |   \    //
+  // .α---|-f-β.   //
+  //  \   |   /    //
+  //   \  |  /     //
+  //   b\ | /c     //
+  //     \|/       //
+  //      .        //
+
+  const Index num_faces = F.rows();
+  assert(uE2E[uei].size() == 2 && "edge should have 2 incident faces");
+  const Index he_left = uE2E[uei][0];
+  const Index he_right = uE2E[uei][1];
+  const Index f_left = he_left%num_faces;
+  const Index c_left = he_left/num_faces;
+  const Index f_right = he_right%num_faces;
+  const Index c_right = he_right/num_faces;
+
+  assert( std::abs(l(f_left,c_left)-l(f_right,c_right) < igl::EPS<Scalar>()) );
+  const Scalar e = l(f_left,c_left);
+  const Scalar a = l(f_left,(c_left+1)%3);
+  const Scalar b = l(f_left,(c_left+2)%3);
+  const Scalar c = l(f_right,(c_right+1)%3);
+  const Scalar d = l(f_right,(c_right+2)%3);
+
+  const Scalar w = cot_alpha(e,a,b) + cot_alpha(e,c,d);
+  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, -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> >&);
+#endif

include/igl/is_intrinsic_delaunay.h

@@ -0,0 +1,58 @@
+// 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);
+  // 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)
+  //   uei  index into uE2E of edge to check
+  // Returns true iff edge is Delaunay
+  template <
+    typename Derivedl,
+    typename DerivedF,
+    typename uE2EType,
+    typename ueiType>
+  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);
+
+}
+#ifndef IGL_STATIC_LIBRARY
+#include "is_intrinsic_delaunay.cpp"
+#endif
+#endif
+

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));
+}