intrinsic delaunay (re-)triangulation [Fisher 2007] + tests

include/igl/intrinsic_delaunay_triangulation.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 "intrinsic_delaunay_triangulation.h"
+#include "is_intrinsic_delaunay.h"
+#include "tan_half_angle.h"
+#include "unique_edge_map.h"
+#include "flip_edge.h"
+#include <iostream>
+
+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)
+{
+  // 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);
+  typedef typename DerivedF::Scalar Index;
+  typedef typename Derivedl::Scalar Scalar;
+  const Index num_faces = F.rows();
+
+  // Copied from delaunay_triangulation
+  bool all_delaunay = false;
+  while(!all_delaunay) 
+  {
+    all_delaunay = true;
+    for (size_t uei=0; uei<uE2E.size(); uei++) 
+    {
+      if (uE2E[uei].size() == 2) 
+      {
+        if(!is_intrinsic_delaunay(l,F,uE2E,uei)) 
+        {
+          all_delaunay = false;
+          // update l just before flipping edge
+          //      .        //
+          //     /|\       //
+          //   a/ | \d     //
+          //   /  e  \     //
+          //  /   |   \    //
+          // .----|-f--.   //
+          //  \   |   /    //
+          //   \  |  /     //
+          //   b\α|δ/c     //
+          //     \|/       //
+          //      .        //
+          // Compute intrinsic length of oppposite edge
+          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);
+          // 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);
+          // 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
+          //
+          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);
+        }
+      }
+    }
+  }
+}
+
+#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::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,43 @@
+// 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>
+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
+  //
+  // 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);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "intrinsic_delaunay_triangulation.cpp"
+#endif
+
+#endif

tests/include/igl/intrinsic_delaunay_triangulation.cpp

@@ -0,0 +1,50 @@
+#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>
+
+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
+  V.col(0) += 1.1*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);
+}