Add support for Shewchuk's predicates.

CMakeLists.txt

@@ -29,6 +29,7 @@ option(LIBIGL_WITH_OPENGL_GLFW_IMGUI "Use ImGui"                    ON)
 option(LIBIGL_WITH_PNG               "Use PNG"                      ON)
 option(LIBIGL_WITH_TETGEN            "Use Tetgen"                   ON)
 option(LIBIGL_WITH_TRIANGLE          "Use Triangle"                 ON)
+option(LIBIGL_WITH_PREDICATES        "Use exact predicates"         ON)
 option(LIBIGL_WITH_XML               "Use XML"                      ON)
 option(LIBIGL_WITH_PYTHON            "Use Python"                   ${LIBIGL_BUILD_PYTHON})
 ### End

cmake/LibiglDownloadExternal.cmake

@@ -152,6 +152,14 @@ function(igl_download_catch2)
 	)
 endfunction()
 
+## Predicates
+function(igl_download_predicates)
+	igl_download_project(predicates
+		GIT_REPOSITORY https://github.com/libigl/libigl-predicates
+		GIT_TAG        55c9251878f54cb3c3b8f25f93a7ce4dd60d9317
+	)
+endfunction()
+
 ################################################################################
 
 ## Test data

cmake/libigl.cmake

@@ -32,6 +32,7 @@ option(LIBIGL_WITH_OPENGL_GLFW_IMGUI "Use ImGui"                    OFF)
 option(LIBIGL_WITH_PNG               "Use PNG"                      OFF)
 option(LIBIGL_WITH_TETGEN            "Use Tetgen"                   OFF)
 option(LIBIGL_WITH_TRIANGLE          "Use Triangle"                 OFF)
+option(LIBIGL_WITH_PREDICATES        "Use exact predicates"         OFF)
 option(LIBIGL_WITH_XML               "Use XML"                      OFF)
 option(LIBIGL_WITH_PYTHON            "Use Python"                   OFF)
 option(LIBIGL_WITHOUT_COPYLEFT       "Disable Copyleft libraries"   OFF)
@@ -411,6 +412,19 @@ if(LIBIGL_WITH_TRIANGLE)
   target_include_directories(igl_triangle ${IGL_SCOPE} ${TRIANGLE_DIR})
 endif()
 
+################################################################################
+### Compile the predicates part ###
+if(LIBIGL_WITH_PREDICATES)
+  set(PREDICATES_DIR "${LIBIGL_EXTERNAL}/predicates")
+  if(NOT TARGET predicates)
+    igl_download_predicates()
+    add_subdirectory("${PREDICATES_DIR}" "predicates")
+  endif()
+  compile_igl_module("predicates")
+  target_link_libraries(igl_predicates ${IGL_SCOPE} predicates)
+  target_include_directories(igl_predicates ${IGL_SCOPE} ${PREDICATES_DIR})
+endif()
+
 ################################################################################
 ### Compile the xml part ###
 if(LIBIGL_WITH_XML)

include/igl/predicates/predicates.cpp

@@ -0,0 +1,189 @@
+#include <igl/predicates/predicates.h>
+#include <predicates.h>
+
+namespace igl {
+namespace predicates {
+
+using REAL = double;
+
+template<typename Vector2D>
+IGL_INLINE Orientation orient2d(
+        const Eigen::MatrixBase<Vector2D>& pa,
+        const Eigen::MatrixBase<Vector2D>& pb,
+        const Eigen::MatrixBase<Vector2D>& pc) {
+
+    using Point = Eigen::Matrix<REAL, Eigen::Dynamic, Eigen::Dynamic>;
+    const Point& a = pa.template cast<REAL>();
+    const Point& b = pb.template cast<REAL>();
+    const Point& c = pc.template cast<REAL>();
+
+    auto r = ::orient2d(
+            const_cast<REAL*>(a.data()),
+            const_cast<REAL*>(b.data()),
+            const_cast<REAL*>(c.data()));
+
+    if (r > 0) return Orientation::POSITIVE;
+    else if (r < 0) return Orientation::NEGATIVE;
+    else return Orientation::COLINEAR;
+}
+
+template<typename Vector3D>
+IGL_INLINE Orientation orient3d(
+        const Eigen::MatrixBase<Vector3D>& pa,
+        const Eigen::MatrixBase<Vector3D>& pb,
+        const Eigen::MatrixBase<Vector3D>& pc,
+        const Eigen::MatrixBase<Vector3D>& pd) {
+
+    using Point = Eigen::Matrix<REAL, Eigen::Dynamic, Eigen::Dynamic>;
+    const Point& a = pa.template cast<REAL>();
+    const Point& b = pb.template cast<REAL>();
+    const Point& c = pc.template cast<REAL>();
+    const Point& d = pd.template cast<REAL>();
+
+    auto r = ::orient3d(
+            const_cast<REAL*>(a.data()),
+            const_cast<REAL*>(b.data()),
+            const_cast<REAL*>(c.data()),
+            const_cast<REAL*>(d.data()));
+
+    if (r > 0) return Orientation::POSITIVE;
+    else if (r < 0) return Orientation::NEGATIVE;
+    else return Orientation::COPLANAR;
+}
+
+template<typename Vector2D>
+IGL_INLINE Orientation incircle(
+        const Eigen::MatrixBase<Vector2D>& pa,
+        const Eigen::MatrixBase<Vector2D>& pb,
+        const Eigen::MatrixBase<Vector2D>& pc,
+        const Eigen::MatrixBase<Vector2D>& pd) {
+
+    using Point = Eigen::Matrix<REAL, Eigen::Dynamic, Eigen::Dynamic>;
+    const Point& a = pa.template cast<REAL>();
+    const Point& b = pb.template cast<REAL>();
+    const Point& c = pc.template cast<REAL>();
+    const Point& d = pd.template cast<REAL>();
+
+    auto r = ::incircle(
+            const_cast<REAL*>(a.data()),
+            const_cast<REAL*>(b.data()),
+            const_cast<REAL*>(c.data()),
+            const_cast<REAL*>(d.data()));
+
+    if (r > 0) return Orientation::INSIDE;
+    else if (r < 0) return Orientation::OUTSIDE;
+    else return Orientation::COCIRCULAR;
+}
+
+template<typename Vector3D>
+IGL_INLINE Orientation insphere(
+        const Eigen::MatrixBase<Vector3D>& pa,
+        const Eigen::MatrixBase<Vector3D>& pb,
+        const Eigen::MatrixBase<Vector3D>& pc,
+        const Eigen::MatrixBase<Vector3D>& pd,
+        const Eigen::MatrixBase<Vector3D>& pe) {
+
+    using Point = Eigen::Matrix<REAL, Eigen::Dynamic, Eigen::Dynamic>;
+    const Point& a = pa.template cast<REAL>();
+    const Point& b = pb.template cast<REAL>();
+    const Point& c = pc.template cast<REAL>();
+    const Point& d = pd.template cast<REAL>();
+    const Point& e = pe.template cast<REAL>();
+
+    auto r = ::insphere(
+            const_cast<REAL*>(a.data()),
+            const_cast<REAL*>(b.data()),
+            const_cast<REAL*>(c.data()),
+            const_cast<REAL*>(d.data()),
+            const_cast<REAL*>(e.data()));
+
+    if (r > 0) return Orientation::INSIDE;
+    else if (r < 0) return Orientation::OUTSIDE;
+    else return Orientation::COSPHERICAL;
+}
+
+}
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template
+igl::predicates::Orientation
+igl::predicates::orient2d<Eigen::Matrix<double, -1, -1, 1, -1, -1>>
+(
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&
+);
+template
+igl::predicates::Orientation
+igl::predicates::orient2d<Eigen::Matrix<float, -1, -1, 1, -1, -1>>
+(
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&
+);
+
+
+template
+igl::predicates::Orientation
+igl::predicates::orient3d<Eigen::Matrix<double, -1, -1, 1, -1, -1>>
+(
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&
+);
+template
+igl::predicates::Orientation
+igl::predicates::orient3d<Eigen::Matrix<float, -1, -1, 1, -1, -1>>
+(
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&
+);
+
+
+template
+igl::predicates::Orientation
+igl::predicates::incircle<Eigen::Matrix<double, -1, -1, 1, -1, -1>>
+(
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&
+);
+template
+igl::predicates::Orientation
+igl::predicates::incircle<Eigen::Matrix<float, -1, -1, 1, -1, -1>>
+(
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&
+);
+
+
+template
+igl::predicates::Orientation
+igl::predicates::insphere<Eigen::Matrix<double, -1, -1, 1, -1, -1>>
+(
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&
+);
+template
+igl::predicates::Orientation
+igl::predicates::insphere<Eigen::Matrix<float, -1, -1, 1, -1, -1>>
+(
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&,
+const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 1, -1, -1>>&
+);
+
+#endif

include/igl/predicates/predicates.h

@@ -0,0 +1,50 @@
+#pragma once
+#ifndef IGL_PREDICATES_PREDICATES_H
+#define IGL_PREDICATES_PREDICATES_H
+
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+
+namespace igl {
+  namespace predicates {
+    enum class Orientation {
+      POSITIVE=1, INSIDE=1,
+      NEGATIVE=-1, OUTSIDE=-1,
+      COLINEAR=0, COPLANAR=0, COCIRCULAR=0, COSPHERICAL=0, DEGENERATE=0
+    };
+
+    template<typename Vector2D>
+    IGL_INLINE Orientation orient2d(
+        const Eigen::MatrixBase<Vector2D>& pa,
+        const Eigen::MatrixBase<Vector2D>& pb,
+        const Eigen::MatrixBase<Vector2D>& pc);
+
+    template<typename Vector3D>
+    IGL_INLINE Orientation orient3d(
+        const Eigen::MatrixBase<Vector3D>& pa,
+        const Eigen::MatrixBase<Vector3D>& pb,
+        const Eigen::MatrixBase<Vector3D>& pc,
+        const Eigen::MatrixBase<Vector3D>& pd);
+
+    template<typename Vector2D>
+    IGL_INLINE Orientation incircle(
+        const Eigen::MatrixBase<Vector2D>& pa,
+        const Eigen::MatrixBase<Vector2D>& pb,
+        const Eigen::MatrixBase<Vector2D>& pc,
+        const Eigen::MatrixBase<Vector2D>& pd);
+
+    template<typename Vector3D>
+    IGL_INLINE Orientation insphere(
+        const Eigen::MatrixBase<Vector3D>& pa,
+        const Eigen::MatrixBase<Vector3D>& pb,
+        const Eigen::MatrixBase<Vector3D>& pc,
+        const Eigen::MatrixBase<Vector3D>& pd,
+        const Eigen::MatrixBase<Vector3D>& pe);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "predicates.cpp"
+#endif
+
+#endif