Add support for exact predicates and lexicographic triangulation.

Former-commit-id: 212df95a162f5ee804d6a1f1d32008d78d20f7f9

include/igl/copyleft/cgal/ExactPredicate.h

@@ -0,0 +1,106 @@
+// 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
+    {
+      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;
+
+          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";
+            }
+          }
+
+          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";
+            }
+          }
+
+          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";
+            }
+          }
+
+          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/lexicographic_triangulation.cpp

@@ -0,0 +1,129 @@
+// 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 "ExactPredicate.h"
+#include "../../sortrows.h"
+
+#include <vector>
+
+template<
+  typename DerivedP,
+  typename DerivedF
+  >
+IGL_INLINE void igl::copyleft::cgal::lexicographic_triangulation(
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    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);
+        std::cout << orientation << " ";
+      }
+      std::cout << std::endl;
+
+      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/copyleft/cgal/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_COPYLEFT_CGAL_LEXICOGRAPHIC_TRIANGULATION_H
+#define IGL_COPYLEFT_CGAL_LEXICOGRAPHIC_TRIANGULATION_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+
+      // Given a set of points in 2D, return a lexicographic triangulation of these
+      // points.
+      //
+      // Inputs:
+      //   P  #P by 2 list of vertex positions
+      //
+      // Outputs:
+      //   F  #F by 3 of faces in lexicographic triangulation.
+      template<
+        typename DerivedP,
+        typename DerivedF
+        >
+      IGL_INLINE void lexicographic_triangulation(
+          const Eigen::PlainObjectBase<DerivedP>& P,
+          Eigen::PlainObjectBase<DerivedF>& F);
+    }
+  }
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "lexicographic_triangulation.cpp"
+#endif
+#endif