Refactor exact predicates class into functions.

Former-commit-id: 0eaf21093d619249d6ff481e8a6bfbdd1bd2345b

include/igl/copyleft/cgal/ExactPredicate.h

@@ -1,140 +0,0 @@
-// 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_EXACT_PREDICATE
-#define IGL_COPYLEFT_CGAL_EXACT_PREDICATE
-
-#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-
-namespace igl
-{
-  namespace copyleft
-  {
-    namespace cgal
-    {
-      // This class extracts the exact predicates routines from CGAL kernels.
-      //
-      // Example:
-      //    Eigen::PlainObjectBase<DerivedV> V = ...;
-      //    typedef typename DerivedV::Scalar Scalar;
-      //    typedef igl::copyleft::cgal::ExactPredicate<Scalar> Predicate;
-      //    auto result = Predicate::orientation(
-      //      {V(0, 0), V(0, 1)},
-      //      {V(1, 0), V(1, 1)},
-      //      {V(2, 0), V(2, 1)});
-      template<typename Scalar>
-      class ExactPredicate {
-        public:
-          typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
-          typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
-          typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
-                  Epeck, Epick>::type Kernel;
-
-          // Inputs:
-          //   pa,pb,pc   2D points.
-          // Output:
-          //   1 if pa,pb,pc are counterclockwise oriented.
-          //   0 if pa,pb,pc are counterclockwise oriented.
-          //  -1 if pa,pb,pc are clockwise oriented.
-          static short orientation(const Scalar pa[2], const Scalar pb[2], const Scalar pc[2])
-          {
-            switch(CGAL::orientation(
-                  typename Kernel::Point_2(pa[0], pa[1]),
-                  typename Kernel::Point_2(pb[0], pb[1]),
-                  typename Kernel::Point_2(pc[0], pc[1]))) {
-              case CGAL::LEFT_TURN:
-                return 1;
-              case CGAL::RIGHT_TURN:
-                return -1;
-              case CGAL::COLLINEAR:
-                return 0;
-              default:
-                throw "Invalid orientation";
-            }
-          }
-
-          // Inputs:
-          //   pa,pb,pc,pd  3D points.
-          // Output:
-          //   1 if pa,pb,pc,pd forms a tet of positive volume.
-          //   0 if pa,pb,pc,pd are coplanar.
-          //  -1 if pa,pb,pc,pd forms a tet of negative volume.
-          static short orientation(const Scalar pa[3], const Scalar pb[3], const Scalar pc[3], const Scalar pd[3])
-          {
-            switch(CGAL::orientation(
-                  typename Kernel::Point_3(pa[0], pa[1], pa[2]),
-                  typename Kernel::Point_3(pb[0], pb[1], pb[2]),
-                  typename Kernel::Point_3(pc[0], pc[1], pc[2]),
-                  typename Kernel::Point_3(pd[0], pd[1], pd[2]))) {
-              case CGAL::POSITIVE:
-                return 1;
-              case CGAL::NEGATIVE:
-                return -1;
-              case CGAL::COPLANAR:
-                return 0;
-              default:
-                throw "Invalid orientation";
-            }
-          }
-
-          // Inputs:
-          //   pa,pb,pc,pd  2D points.
-          // Output:
-          //   1 if pd is inside of the oriented circle formed by pa,pb,pc.
-          //   0 if pd is co-circular with pa,pb,pc.
-          //  -1 if pd is outside of the oriented circle formed by pa,pb,pc.
-          static short incircle(const Scalar pa[2], const Scalar pb[2], const Scalar pc[2], const Scalar pd[2])
-          {
-            switch(CGAL::side_of_oriented_circle(
-                  typename Kernel::Point_2(pa[0], pa[1]),
-                  typename Kernel::Point_2(pb[0], pb[1]),
-                  typename Kernel::Point_2(pc[0], pc[1]),
-                  typename Kernel::Point_2(pd[0], pd[1]))) {
-              case CGAL::ON_POSITIVE_SIDE:
-                return 1;
-              case CGAL::ON_NEGATIVE_SIDE:
-                return -1;
-              case CGAL::ON_ORIENTED_BOUNDARY:
-                return 0;
-              default:
-                throw "Invalid incircle result";
-            }
-          }
-
-          // Inputs:
-          //   pa,pb,pc,pd,pe  3D points.
-          // Output:
-          //   1 if pe is inside of the oriented sphere formed by pa,pb,pc,pd.
-          //   0 if pe is co-spherical with pa,pb,pc,pd.
-          //  -1 if pe is outside of the oriented sphere formed by pa,pb,pc,pd.
-          static short insphere(const Scalar pa[3], const Scalar pb[3], const Scalar pc[3], const Scalar pd[3],
-              const Scalar pe[3])
-          {
-            switch(CGAL::side_of_oriented_sphere(
-                  typename Kernel::Point_3(pa[0], pa[1], pa[2]),
-                  typename Kernel::Point_3(pb[0], pb[1], pb[2]),
-                  typename Kernel::Point_3(pc[0], pc[1], pc[2]),
-                  typename Kernel::Point_3(pd[0], pd[1], pd[2]),
-                  typename Kernel::Point_3(pe[0], pe[1], pe[2]))) {
-              case CGAL::ON_POSITIVE_SIDE:
-                return 1;
-              case CGAL::ON_NEGATIVE_SIDE:
-                return -1;
-              case CGAL::ON_ORIENTED_BOUNDARY:
-                return 0;
-              default:
-                throw "Invalid incircle result";
-            }
-          }
-      };
-    }
-  }
-}
-
-#endif

include/igl/copyleft/cgal/delaunay_triangulation.cpp

@@ -7,13 +7,9 @@
 // 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>
+#include "../../delaunay_triangulation.h"
+#include "orient2D.h"
+#include "incircle.h"
 
 template<
   typename DerivedV,
@@ -22,60 +18,7 @@ IGL_INLINE void igl::copyleft::cgal::delaunay_triangulation(
     const Eigen::PlainObjectBase<DerivedV>& V,
     Eigen::PlainObjectBase<DerivedF>& F)
 {
-  assert(V.cols() == 2);
-  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();
-  if (num_faces == 0) {
-    // Input points are degenerate.  No faces will be generated.
-    return;
-  }
-  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);
-        }
-      }
-    }
-  }
+  igl::delaunay_triangulation(V, orient2D<Scalar>, incircle<Scalar>, F);
 }
 

include/igl/copyleft/cgal/incircle.cpp

@@ -0,0 +1,39 @@
+// 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/.
+
+#include "incircle.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::incircle(
+    const Scalar pa[2],
+    const Scalar pb[2],
+    const Scalar pc[2],
+    const Scalar pd[2])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::side_of_oriented_circle(
+        typename Kernel::Point_2(pa[0], pa[1]),
+        typename Kernel::Point_2(pb[0], pb[1]),
+        typename Kernel::Point_2(pc[0], pc[1]),
+        typename Kernel::Point_2(pd[0], pd[1]))) {
+    case CGAL::ON_POSITIVE_SIDE:
+      return 1;
+    case CGAL::ON_NEGATIVE_SIDE:
+      return -1;
+    case CGAL::ON_ORIENTED_BOUNDARY:
+      return 0;
+    default:
+      throw "Invalid incircle result";
+  }
+}

include/igl/copyleft/cgal/incircle.h

@@ -0,0 +1,39 @@
+// 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_INCIRCLE_H
+#define IGL_COPYLEFT_CGAL_INCIRCLE_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc,pd  2D points.
+      // Output:
+      //   1 if pd is inside of the oriented circle formed by pa,pb,pc.
+      //   0 if pd is co-circular with pa,pb,pc.
+      //  -1 if pd is outside of the oriented circle formed by pa,pb,pc.
+      template <typename Scalar>
+      IGL_INLINE short incircle(
+          const Scalar pa[2],
+          const Scalar pb[2],
+          const Scalar pc[2],
+          const Scalar pd[2]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "incircle.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/insphere.cpp

@@ -0,0 +1,41 @@
+// 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/.
+
+#include "insphere.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::insphere(
+    const Scalar pa[3],
+    const Scalar pb[3],
+    const Scalar pc[3],
+    const Scalar pd[3],
+    const Scalar pe[3])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::side_of_oriented_sphere(
+        typename Kernel::Point_3(pa[0], pa[1], pa[2]),
+        typename Kernel::Point_3(pb[0], pb[1], pb[2]),
+        typename Kernel::Point_3(pc[0], pc[1], pc[2]),
+        typename Kernel::Point_3(pd[0], pd[1], pd[2]),
+        typename Kernel::Point_3(pe[0], pe[1], pe[2]))) {
+    case CGAL::ON_POSITIVE_SIDE:
+      return 1;
+    case CGAL::ON_NEGATIVE_SIDE:
+      return -1;
+    case CGAL::ON_ORIENTED_BOUNDARY:
+      return 0;
+    default:
+      throw "Invalid incircle result";
+  }
+}

include/igl/copyleft/cgal/insphere.h

@@ -0,0 +1,40 @@
+// 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_INSPHERE_H
+#define IGL_COPYLEFT_CGAL_INSPHERE_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc,pd,pe  3D points.
+      // Output:
+      //   1 if pe is inside of the oriented sphere formed by pa,pb,pc,pd.
+      //   0 if pe is co-spherical with pa,pb,pc,pd.
+      //  -1 if pe is outside of the oriented sphere formed by pa,pb,pc,pd.
+      template <typename Scalar>
+      IGL_INLINE short insphere(
+          const Scalar pa[3],
+          const Scalar pb[3],
+          const Scalar pc[3],
+          const Scalar pd[3],
+          const Scalar pe[3]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "insphere.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/lexicographic_triangulation.cpp

@@ -1,17 +1,15 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-//
+// 
 // Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
 //                    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
+// 
+// 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 "lexicographic_triangulation.h"
-#include "ExactPredicate.h"
-#include "../../sortrows.h"
-
-#include <vector>
+#include "../../lexicographic_triangulation.h"
+#include "orient2D.h"
 
 template<
   typename DerivedP,
@@ -22,106 +20,5 @@ IGL_INLINE void igl::copyleft::cgal::lexicographic_triangulation(
     Eigen::PlainObjectBase<DerivedF>& F)
 {
   typedef typename DerivedP::Scalar Scalar;
-  typedef typename igl::copyleft::cgal::ExactPredicate<Scalar> Predicate;
-  const size_t num_pts = P.rows();
-  if (num_pts < 3) {
-    throw "At least 3 points are required for triangulation!";
-  }
-
-  // Sort points in lexicographic order.
-  DerivedP ordered_P;
-  Eigen::VectorXi order;
-  igl::sortrows(P, true, ordered_P, order);
-
-  std::vector<Eigen::Vector3i> faces;
-  std::list<size_t> boundary;
-  const Scalar p0[] = {ordered_P(0, 0), ordered_P(0, 1)};
-  const Scalar p1[] = {ordered_P(1, 0), ordered_P(1, 1)};
-  for (size_t i=2; i<num_pts; i++) {
-    const Scalar curr_p[] = {ordered_P(i, 0), ordered_P(i, 1)};
-    if (faces.size() == 0) {
-      // All points processed so far are collinear.
-      // Check if the current point is collinear with every points before it.
-      auto orientation = Predicate::orientation(p0, p1, curr_p);
-      if (orientation != 0) {
-        // Add a fan of triangles eminating from curr_p.
-        if (orientation > 0) {
-          for (size_t j=0; j<=i-2; j++) {
-            faces.push_back({order[j], order[j+1], order[i]});
-          }
-        } else if (orientation < 0) {
-          for (size_t j=0; j<=i-2; j++) {
-            faces.push_back({order[j+1], order[j], order[i]});
-          }
-        }
-        // Initialize current boundary.
-        boundary.insert(boundary.end(), order.data(), order.data()+i+1);
-        if (orientation < 0) {
-          boundary.reverse();
-        }
-      }
-    } else {
-      const size_t bd_size = boundary.size();
-      assert(bd_size >= 3);
-      std::vector<short> orientations;
-      for (auto itr=boundary.begin(); itr!=boundary.end(); itr++) {
-        auto next_itr = std::next(itr, 1);
-        if (next_itr == boundary.end()) {
-          next_itr = boundary.begin();
-        }
-        const Scalar bd_p0[] = {P(*itr, 0), P(*itr, 1)};
-        const Scalar bd_p1[] = {P(*next_itr, 0), P(*next_itr, 1)};
-        auto orientation = Predicate::orientation(bd_p0, bd_p1, curr_p);
-        if (orientation < 0) {
-          faces.push_back({*next_itr, *itr, order[i]});
-        }
-        orientations.push_back(orientation);
-      }
-
-      auto left_itr = boundary.begin();
-      auto right_itr = boundary.begin();
-      auto curr_itr = boundary.begin();
-      for (size_t j=0; j<bd_size; j++, curr_itr++) {
-        size_t prev = (j+bd_size-1) % bd_size;
-        if (orientations[j] >= 0 && orientations[prev] < 0) {
-          right_itr = curr_itr;
-        } else if (orientations[j] < 0 && orientations[prev] >= 0) {
-          left_itr = curr_itr;
-        }
-      }
-      assert(left_itr != right_itr);
-
-      for (auto itr=left_itr; itr!=right_itr; itr++) {
-        if (itr == boundary.end()) itr = boundary.begin();
-        if (itr == right_itr) break;
-        if (itr == left_itr || itr == right_itr) continue;
-        itr = boundary.erase(itr);
-        if (itr == boundary.begin()) {
-            itr = boundary.end();
-        } else {
-            itr--;
-        }
-      }
-
-      if (right_itr == boundary.begin()) {
-        assert(std::next(left_itr, 1) == boundary.end());
-        boundary.insert(boundary.end(), order[i]);
-      } else {
-        assert(std::next(left_itr, 1) == right_itr);
-        boundary.insert(right_itr, order[i]);
-      }
-    }
-  }
-
-  const size_t num_faces = faces.size();
-  if (num_faces == 0) {
-      // All input points are collinear.
-      // Do nothing here.
-  } else {
-      F.resize(num_faces, 3);
-      for (size_t i=0; i<num_faces; i++) {
-          F.row(i) = faces[i];
-      }
-  }
+  igl::lexicographic_triangulation(P, orient2D<Scalar>, F);
 }
-

include/igl/copyleft/cgal/orient2D.cpp

@@ -0,0 +1,37 @@
+// 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/.
+
+#include "orient2D.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::orient2D(
+    const Scalar pa[2],
+    const Scalar pb[2],
+    const Scalar pc[2])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::orientation(
+        typename Kernel::Point_2(pa[0], pa[1]),
+        typename Kernel::Point_2(pb[0], pb[1]),
+        typename Kernel::Point_2(pc[0], pc[1]))) {
+    case CGAL::LEFT_TURN:
+      return 1;
+    case CGAL::RIGHT_TURN:
+      return -1;
+    case CGAL::COLLINEAR:
+      return 0;
+    default:
+      throw "Invalid orientation";
+  }
+}

include/igl/copyleft/cgal/orient2D.h

@@ -0,0 +1,38 @@
+// 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_ORIENT_2D_H
+#define IGL_COPYLEFT_CGAL_ORIENT_2D_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc   2D points.
+      // Output:
+      //   1 if pa,pb,pc are counterclockwise oriented.
+      //   0 if pa,pb,pc are counterclockwise oriented.
+      //  -1 if pa,pb,pc are clockwise oriented.
+      template <typename Scalar>
+      IGL_INLINE short orient2D(
+          const Scalar pa[2],
+          const Scalar pb[2],
+          const Scalar pc[2]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "orient2D.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/orient3D.cpp

@@ -0,0 +1,39 @@
+// 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/.
+
+#include "orient3D.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::orient3D(
+    const Scalar pa[3],
+    const Scalar pb[3],
+    const Scalar pc[3],
+    const Scalar pd[3])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::orientation(
+        typename Kernel::Point_3(pa[0], pa[1], pa[2]),
+        typename Kernel::Point_3(pb[0], pb[1], pb[2]),
+        typename Kernel::Point_3(pc[0], pc[1], pc[2]),
+        typename Kernel::Point_3(pd[0], pd[1], pd[2]))) {
+    case CGAL::POSITIVE:
+      return 1;
+    case CGAL::NEGATIVE:
+      return -1;
+    case CGAL::COPLANAR:
+      return 0;
+    default:
+      throw "Invalid orientation";
+  }
+}

include/igl/copyleft/cgal/orient3D.h

@@ -0,0 +1,39 @@
+// 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_ORIENT_3D_H
+#define IGL_COPYLEFT_CGAL_ORIENT_3D_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc,pd  3D points.
+      // Output:
+      //   1 if pa,pb,pc,pd forms a tet of positive volume.
+      //   0 if pa,pb,pc,pd are coplanar.
+      //  -1 if pa,pb,pc,pd forms a tet of negative volume.
+      template <typename Scalar>
+      IGL_INLINE short orient3D(
+          const Scalar pa[3],
+          const Scalar pb[3],
+          const Scalar pc[3],
+          const Scalar pd[3]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "orient3D.cpp"
+#endif
+#endif

include/igl/delaunay_triangulation.cpp

@@ -0,0 +1,83 @@
+// 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 "flip_edge.h"
+#include "lexicographic_triangulation.h"
+#include "unique_edge_map.h"
+
+#include <vector>
+#include <sstream>
+
+template<
+  typename DerivedV,
+  typename Orient2D,
+  typename InCircle,
+  typename DerivedF>
+IGL_INLINE void igl::delaunay_triangulation(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    Orient2D orient2D,
+    InCircle incircle,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  assert(V.cols() == 2);
+  typedef typename DerivedF::Scalar Index;
+  typedef typename DerivedV::Scalar Scalar;
+  igl::lexicographic_triangulation(V, orient2D, F);
+  const size_t num_faces = F.rows();
+  if (num_faces == 0) {
+    // Input points are degenerate.  No faces will be generated.
+    return;
+  }
+  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,&incircle](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 = 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/delaunay_triangulation.h

@@ -0,0 +1,44 @@
+// 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_DELAUNAY_TRIANGULATION_H
+#define IGL_DELAUNAY_TRIANGULATION_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // 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 Orient2D,
+    typename InCircle,
+    typename DerivedF
+    >
+  IGL_INLINE void delaunay_triangulation(
+      const Eigen::PlainObjectBase<DerivedV>& V,
+      Orient2D orient2D,
+      InCircle incircle,
+      Eigen::PlainObjectBase<DerivedF>& F);
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "delaunay_triangulation.cpp"
+#endif
+#endif

include/igl/lexicographic_triangulation.cpp

@@ -0,0 +1,128 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//                    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/.
+
+#include "lexicographic_triangulation.h"
+#include "sortrows.h"
+
+#include <vector>
+#include <list>
+
+template<
+  typename DerivedP,
+  typename Orient2D,
+  typename DerivedF
+  >
+IGL_INLINE void igl::lexicographic_triangulation(
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    Orient2D orient2D,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  typedef typename DerivedP::Scalar Scalar;
+  const size_t num_pts = P.rows();
+  if (num_pts < 3) {
+    throw "At least 3 points are required for triangulation!";
+  }
+
+  // Sort points in lexicographic order.
+  DerivedP ordered_P;
+  Eigen::VectorXi order;
+  igl::sortrows(P, true, ordered_P, order);
+
+  std::vector<Eigen::Vector3i> faces;
+  std::list<size_t> boundary;
+  const Scalar p0[] = {ordered_P(0, 0), ordered_P(0, 1)};
+  const Scalar p1[] = {ordered_P(1, 0), ordered_P(1, 1)};
+  for (size_t i=2; i<num_pts; i++) {
+    const Scalar curr_p[] = {ordered_P(i, 0), ordered_P(i, 1)};
+    if (faces.size() == 0) {
+      // All points processed so far are collinear.
+      // Check if the current point is collinear with every points before it.
+      auto orientation = orient2D(p0, p1, curr_p);
+      if (orientation != 0) {
+        // Add a fan of triangles eminating from curr_p.
+        if (orientation > 0) {
+          for (size_t j=0; j<=i-2; j++) {
+            faces.push_back({order[j], order[j+1], order[i]});
+          }
+        } else if (orientation < 0) {
+          for (size_t j=0; j<=i-2; j++) {
+            faces.push_back({order[j+1], order[j], order[i]});
+          }
+        }
+        // Initialize current boundary.
+        boundary.insert(boundary.end(), order.data(), order.data()+i+1);
+        if (orientation < 0) {
+          boundary.reverse();
+        }
+      }
+    } else {
+      const size_t bd_size = boundary.size();
+      assert(bd_size >= 3);
+      std::vector<short> orientations;
+      for (auto itr=boundary.begin(); itr!=boundary.end(); itr++) {
+        auto next_itr = std::next(itr, 1);
+        if (next_itr == boundary.end()) {
+          next_itr = boundary.begin();
+        }
+        const Scalar bd_p0[] = {P(*itr, 0), P(*itr, 1)};
+        const Scalar bd_p1[] = {P(*next_itr, 0), P(*next_itr, 1)};
+        auto orientation = orient2D(bd_p0, bd_p1, curr_p);
+        if (orientation < 0) {
+          faces.push_back({*next_itr, *itr, order[i]});
+        }
+        orientations.push_back(orientation);
+      }
+
+      auto left_itr = boundary.begin();
+      auto right_itr = boundary.begin();
+      auto curr_itr = boundary.begin();
+      for (size_t j=0; j<bd_size; j++, curr_itr++) {
+        size_t prev = (j+bd_size-1) % bd_size;
+        if (orientations[j] >= 0 && orientations[prev] < 0) {
+          right_itr = curr_itr;
+        } else if (orientations[j] < 0 && orientations[prev] >= 0) {
+          left_itr = curr_itr;
+        }
+      }
+      assert(left_itr != right_itr);
+
+      for (auto itr=left_itr; itr!=right_itr; itr++) {
+        if (itr == boundary.end()) itr = boundary.begin();
+        if (itr == right_itr) break;
+        if (itr == left_itr || itr == right_itr) continue;
+        itr = boundary.erase(itr);
+        if (itr == boundary.begin()) {
+            itr = boundary.end();
+        } else {
+            itr--;
+        }
+      }
+
+      if (right_itr == boundary.begin()) {
+        assert(std::next(left_itr, 1) == boundary.end());
+        boundary.insert(boundary.end(), order[i]);
+      } else {
+        assert(std::next(left_itr, 1) == right_itr);
+        boundary.insert(right_itr, order[i]);
+      }
+    }
+  }
+
+  const size_t num_faces = faces.size();
+  if (num_faces == 0) {
+      // All input points are collinear.
+      // Do nothing here.
+  } else {
+      F.resize(num_faces, 3);
+      for (size_t i=0; i<num_faces; i++) {
+          F.row(i) = faces[i];
+      }
+  }
+}
+

include/igl/lexicographic_triangulation.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//                    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_LEXICOGRAPHIC_TRIANGULATION_H
+#define IGL_LEXICOGRAPHIC_TRIANGULATION_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Given a set of points in 2D, return a lexicographic triangulation of these
+  // points.
+  //
+  // Inputs:
+  //   P  #P by 2 list of vertex positions
+  //   orient2D  A functor such that orient2D(pa, pb, pc) returns
+  //               1 if pa,pb,pc forms a conterclockwise triangle.
+  //              -1 if pa,pb,pc forms a clockwise triangle.
+  //               0 if pa,pb,pc are collinear.
+  //              where the argument pa,pb,pc are of type Scalar[2].
+  //
+  // Outputs:
+  //   F  #F by 3 of faces in lexicographic triangulation.
+  template<
+    typename DerivedP,
+    typename Orient2D,
+    typename DerivedF
+    >
+  IGL_INLINE void lexicographic_triangulation(
+      const Eigen::PlainObjectBase<DerivedP>& P,
+      Orient2D orient2D,
+      Eigen::PlainObjectBase<DerivedF>& F);
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "lexicographic_triangulation.cpp"
+#endif
+#endif