Delaunay triangulation with naive edge flipping approach.

Former-commit-id: b4c99cb7f71dc58310de57bafa4ac32d89073fe6

include/igl/copyleft/cgal/delaunay_triangulation.cpp

@@ -0,0 +1,76 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Qingnan Zhou <qnzhou@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 "delaunay_triangulation.h"
+#include "../../unique_edge_map.h"
+#include "../../flip_edge.h"
+#include "ExactPredicate.h"
+#include "lexicographic_triangulation.h"
+
+#include <vector>
+#include <sstream>
+
+template<
+  typename DerivedV,
+  typename DerivedF>
+IGL_INLINE void igl::copyleft::cgal::delaunay_triangulation(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  typedef typename DerivedF::Scalar Index;
+  typedef typename DerivedV::Scalar Scalar;
+  typedef typename igl::copyleft::cgal::ExactPredicate<Scalar> Predicate;
+  igl::copyleft::cgal::lexicographic_triangulation(V, F);
+  const size_t num_faces = F.rows();
+  assert(F.cols() == 3);
+
+  Eigen::MatrixXi E;
+  Eigen::MatrixXi uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<Index> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+
+  auto is_delaunay = [&V,&F,&uE2E,num_faces](size_t uei) {
+    auto& half_edges = uE2E[uei];
+    if (half_edges.size() != 2) {
+      throw "Cannot flip non-manifold or boundary edge";
+    }
+
+    const size_t f1 = half_edges[0] % num_faces;
+    const size_t f2 = half_edges[1] % num_faces;
+    const size_t c1 = half_edges[0] / num_faces;
+    const size_t c2 = half_edges[1] / num_faces;
+    assert(c1 < 3);
+    assert(c2 < 3);
+    assert(f1 != f2);
+    const size_t v1 = F(f1, (c1+1)%3);
+    const size_t v2 = F(f1, (c1+2)%3);
+    const size_t v4 = F(f1, c1);
+    const size_t v3 = F(f2, c2);
+    const Scalar p1[] = {V(v1, 0), V(v1, 1)};
+    const Scalar p2[] = {V(v2, 0), V(v2, 1)};
+    const Scalar p3[] = {V(v3, 0), V(v3, 1)};
+    const Scalar p4[] = {V(v4, 0), V(v4, 1)};
+    auto orientation = Predicate::incircle(p1, p2, p4, p3);
+    return orientation <= 0;
+  };
+
+  bool all_delaunay = false;
+  while(!all_delaunay) {
+    all_delaunay = true;
+    for (size_t i=0; i<uE2E.size(); i++) {
+      if (uE2E[i].size() == 2) {
+        if (!is_delaunay(i)) {
+          all_delaunay = false;
+          flip_edge(F, E, uE, EMAP, uE2E, i);
+        }
+      }
+    }
+  }
+}
+

include/igl/copyleft/cgal/delaunay_triangulation.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_DELAUNAY_TRIANGULATION_H
+#define IGL_COPYLEFT_CGAL_DELAUNAY_TRIANGULATION_H
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+
+      // Given a set of points in 2D, return a Delaunay triangulation of these
+      // points.
+      //
+      // Inputs:
+      //   V  #V by 2 list of vertex positions
+      //
+      // Outputs:
+      //   F  #F by 3 of faces in Delaunay triangulation.
+      template<
+        typename DerivedV,
+        typename DerivedF
+        >
+      IGL_INLINE void delaunay_triangulation(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          Eigen::PlainObjectBase<DerivedF>& F);
+    }
+  }
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "delaunay_triangulation.cpp"
+#endif
+#endif

include/igl/flip_edge.cpp

@@ -0,0 +1,148 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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/.
+
+template <
+  typename DerivedF,
+  typename DerivedE,
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename uE2EType>
+IGL_INLINE void igl::flip_edge(
+  Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DeriveduE> & uE,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+  std::vector<std::vector<uE2EType> > & uE2E,
+  const size_t uei)
+{
+  typedef typename DerivedF::Scalar Index;
+  const size_t num_faces = F.rows();
+  if (F.cols() != 3) {
+    throw "Edge flip only works on triangle mesh";
+  }
+  //          v1                 v1
+  //          /|\                / \
+  //         / | \              /f1 \
+  //     v3 /f2|f1\ v4  =>  v3 /_____\ v4
+  //        \  |  /            \ f2  /
+  //         \ | /              \   /
+  //          \|/                \ /
+  //          v2                 v2
+  auto& half_edges = uE2E[uei];
+  if (half_edges.size() != 2) {
+    throw "Cannot flip non-manifold or boundary edge";
+  }
+
+  const size_t f1 = half_edges[0] % num_faces;
+  const size_t f2 = half_edges[1] % num_faces;
+  const size_t c1 = half_edges[0] / num_faces;
+  const size_t c2 = half_edges[1] / num_faces;
+  assert(c1 < 3);
+  assert(c2 < 3);
+
+  assert(f1 != f2);
+  const size_t v1 = F(f1, (c1+1)%3);
+  const size_t v2 = F(f1, (c1+2)%3);
+  const size_t v4 = F(f1, c1);
+  const size_t v3 = F(f2, c2);
+  assert(F(f2, (c2+2)%3) == v1);
+  assert(F(f2, (c2+1)%3) == v2);
+
+  const size_t e_12 = half_edges[0];
+  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_21 = half_edges[1];
+  const size_t e_13 = f2 + ((c2 + 1) % 3) * num_faces;
+  const size_t e_32 = f2 + ((c2 + 2) % 3) * num_faces;
+  assert(E(e_12, 0) == v1);
+  assert(E(e_12, 1) == v2);
+  assert(E(e_24, 0) == v2);
+  assert(E(e_24, 1) == v4);
+  assert(E(e_41, 0) == v4);
+  assert(E(e_41, 1) == v1);
+  assert(E(e_21, 0) == v2);
+  assert(E(e_21, 1) == v1);
+  assert(E(e_13, 0) == v1);
+  assert(E(e_13, 1) == v3);
+  assert(E(e_32, 0) == v3);
+  assert(E(e_32, 1) == v2);
+
+  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];
+
+  F(f1, 0) = v1;
+  F(f1, 1) = v3;
+  F(f1, 2) = v4;
+  F(f2, 0) = v2;
+  F(f2, 1) = v4;
+  F(f2, 2) = v3;
+
+  uE(uei, 0) = v3;
+  uE(uei, 1) = v4;
+
+  const size_t new_e_34 = f1;
+  const size_t new_e_41 = f1 + num_faces;
+  const size_t new_e_13 = f1 + num_faces*2;
+  const size_t new_e_43 = f2;
+  const size_t new_e_32 = f2 + num_faces;
+  const size_t new_e_24 = f2 + num_faces*2;
+
+  E(new_e_34, 0) = v3;
+  E(new_e_34, 1) = v4;
+  E(new_e_41, 0) = v4;
+  E(new_e_41, 1) = v1;
+  E(new_e_13, 0) = v1;
+  E(new_e_13, 1) = v3;
+  E(new_e_43, 0) = v4;
+  E(new_e_43, 1) = v3;
+  E(new_e_32, 0) = v3;
+  E(new_e_32, 1) = v2;
+  E(new_e_24, 0) = v2;
+  E(new_e_24, 1) = v4;
+
+  EMAP[new_e_34] = uei;
+  EMAP[new_e_43] = uei;
+  EMAP[new_e_41] = ue_41;
+  EMAP[new_e_13] = ue_13;
+  EMAP[new_e_32] = ue_32;
+  EMAP[new_e_24] = ue_24;
+
+  auto replace = [](std::vector<Index>& array, Index old_v, Index new_v) {
+    std::replace(array.begin(), array.end(), old_v, new_v);
+  };
+  replace(uE2E[uei], e_12, new_e_34);
+  replace(uE2E[uei], e_21, new_e_43);
+  replace(uE2E[ue_13], e_13, new_e_13);
+  replace(uE2E[ue_32], e_32, new_e_32);
+  replace(uE2E[ue_24], e_24, new_e_24);
+  replace(uE2E[ue_41], e_41, new_e_41);
+
+#ifndef NDEBUG
+  auto sanity_check = [&](size_t ue) {
+    const auto& adj_faces = uE2E[ue];
+    if (adj_faces.size() == 2) {
+      const size_t first_f  = adj_faces[0] % num_faces;
+      const size_t first_c  = adj_faces[0] / num_faces;
+      const size_t second_f = adj_faces[1] % num_faces;
+      const size_t second_c = adj_faces[1] / num_faces;
+      const size_t vertex_0 = F(first_f, (first_c+1) % 3);
+      const size_t vertex_1 = F(first_f, (first_c+2) % 3);
+      assert(vertex_0 == F(second_f, (second_c+2) % 3));
+      assert(vertex_1 == F(second_f, (second_c+1) % 3));
+    }
+  };
+
+  sanity_check(uei);
+  sanity_check(ue_13);
+  sanity_check(ue_32);
+  sanity_check(ue_24);
+  sanity_check(ue_41);
+#endif
+}

include/igl/flip_edge.h

@@ -0,0 +1,51 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_FLIP_EDGE_H
+#define IGL_FLIP_EDGE_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Flip an edge in a triangle mesh.  The edge specified by uei must have
+  // exactly **two** adjacent faces.  Violation will result in exception.
+  // Another warning: edge flipping could convert manifold mesh into
+  // non-manifold.
+  //
+  // Inputs:
+  //   F    #F by 3 list of triangles.
+  //   E    #F*3 by 2 list of all of directed edges
+  //   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
+  //   ue   index into uE the edge to be flipped.
+  //
+  // Output:
+  //   Updated F, E, uE, EMAP and uE2E.
+  template <
+    typename DerivedF,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType>
+  IGL_INLINE void flip_edge(
+    Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DeriveduE> & uE,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    std::vector<std::vector<uE2EType> > & uE2E,
+    const size_t uei);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "flip_edge.cpp"
+#endif
+#endif