Merge pull request #1163 from PyMesh/dev

Shewchuk's predicates support

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        ce494af8c90c2cc699146b61c33de102fa96635b
+	)
+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,20 @@ 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})
+  target_compile_definitions(igl_predicates ${IGL_SCOPE} -DLIBIGL_WITH_PREDICATES)
+endif()
+
 ################################################################################
 ### Compile the xml part ###
 if(LIBIGL_WITH_XML)

include/igl/predicates/predicates.cpp

@@ -0,0 +1,162 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 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 <igl/predicates/predicates.h>
+#include <predicates.h>
+#include <type_traits>
+
+namespace igl {
+namespace predicates {
+
+using REAL = IGL_PREDICATES_REAL;
+
+#ifdef LIBIGL_PREDICATES_USE_FLOAT
+#define IGL_ASSERT_SCALAR(Vector)                        \
+  static_assert(                                         \
+    std::is_same<typename Vector::Scalar, float>::value, \
+    "Shewchuk's exact predicates only support float")
+#else
+#define IGL_ASSERT_SCALAR(Vector)                           \
+  static_assert(                                            \
+    std::is_same<typename Vector::Scalar, double>::value || \
+    std::is_same<typename Vector::Scalar, float>::value,    \
+    "Shewchuk's exact predicates only support float and double")
+#endif
+
+template<typename Vector2D>
+IGL_INLINE Orientation orient2d(
+    const Eigen::MatrixBase<Vector2D>& pa,
+    const Eigen::MatrixBase<Vector2D>& pb,
+    const Eigen::MatrixBase<Vector2D>& pc)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Vector2D, 2);
+  IGL_ASSERT_SCALAR(Vector2D);
+
+  using Point = Eigen::Matrix<REAL, 2, 1>;
+  Point a{pa[0], pa[1]};
+  Point b{pb[0], pb[1]};
+  Point c{pc[0], pc[1]};
+
+  const auto r = ::orient2d(a.data(), b.data(), c.data());
+
+  if (r > 0) return Orientation::POSITIVE;
+  else if (r < 0) return Orientation::NEGATIVE;
+  else return Orientation::COLLINEAR;
+}
+
+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)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Vector3D, 3);
+  IGL_ASSERT_SCALAR(Vector3D);
+
+  using Point = Eigen::Matrix<REAL, 3, 1>;
+  Point a{pa[0], pa[1], pa[2]};
+  Point b{pb[0], pb[1], pb[2]};
+  Point c{pc[0], pc[1], pc[2]};
+  Point d{pd[0], pd[1], pd[2]};
+
+  const auto r = ::orient3d(a.data(), b.data(), c.data(), 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)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Vector2D, 2);
+  IGL_ASSERT_SCALAR(Vector2D);
+
+  using Point = Eigen::Matrix<REAL, 2, 1>;
+  Point a{pa[0], pa[1]};
+  Point b{pb[0], pb[1]};
+  Point c{pc[0], pc[1]};
+  Point d{pd[0], pd[1]};
+
+  const auto r = ::incircle(a.data(), b.data(), c.data(), 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)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Vector3D, 3);
+  IGL_ASSERT_SCALAR(Vector3D);
+
+  using Point = Eigen::Matrix<REAL, 3, 1>;
+  Point a{pa[0], pa[1], pa[2]};
+  Point b{pb[0], pb[1], pb[2]};
+  Point c{pc[0], pc[1], pc[2]};
+  Point d{pd[0], pd[1], pd[2]};
+  Point e{pe[0], pe[1], pe[2]};
+
+  const auto r = ::insphere(a.data(), b.data(), c.data(), d.data(), e.data());
+
+  if (r > 0) return Orientation::INSIDE;
+  else if (r < 0) return Orientation::OUTSIDE;
+  else return Orientation::COSPHERICAL;
+}
+
+}
+}
+
+#undef IGL_ASSERT_SCALAR
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+
+#define IGL_ORIENT2D(Vector) template igl::predicates::Orientation igl::predicates::orient2d<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
+#define IGL_INCIRCLE(Vector) template igl::predicates::Orientation igl::predicates::incircle<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
+#define IGL_ORIENT3D(Vector) template igl::predicates::Orientation igl::predicates::orient3d<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
+#define IGL_INSPHERE(Vector) template igl::predicates::Orientation igl::predicates::insphere<Vector>(const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&, const Eigen::MatrixBase<Vector>&)
+
+#define IGL_MATRIX(T, R, C) Eigen::Matrix<T, R, C>
+IGL_ORIENT2D(IGL_MATRIX(float, 1, 2));
+IGL_INCIRCLE(IGL_MATRIX(float, 1, 2));
+IGL_ORIENT2D(IGL_MATRIX(float, 2, 1));
+IGL_INCIRCLE(IGL_MATRIX(float, 2, 1));
+IGL_ORIENT3D(IGL_MATRIX(float, 1, 3));
+IGL_INSPHERE(IGL_MATRIX(float, 1, 3));
+IGL_ORIENT3D(IGL_MATRIX(float, 3, 1));
+IGL_INSPHERE(IGL_MATRIX(float, 3, 1));
+
+#ifndef LIBIGL_PREDICATES_USE_FLOAT
+IGL_ORIENT2D(IGL_MATRIX(double, 1, 2));
+IGL_INCIRCLE(IGL_MATRIX(double, 1, 2));
+IGL_ORIENT2D(IGL_MATRIX(double, 2, 1));
+IGL_INCIRCLE(IGL_MATRIX(double, 2, 1));
+IGL_ORIENT3D(IGL_MATRIX(double, 1, 3));
+IGL_INSPHERE(IGL_MATRIX(double, 1, 3));
+IGL_ORIENT3D(IGL_MATRIX(double, 3, 1));
+IGL_INSPHERE(IGL_MATRIX(double, 3, 1));
+#endif
+#undef IGL_MATRIX
+
+#undef IGL_ORIENT2D
+#undef IGL_ORIENT3D
+#undef IGL_INCIRCLE
+#undef IGL_INSPHERE
+
+#endif

include/igl/predicates/predicates.h

@@ -0,0 +1,95 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 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/.
+#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,
+      COLLINEAR=0, COPLANAR=0, COCIRCULAR=0, COSPHERICAL=0, DEGENERATE=0
+    };
+
+    // Compute the orientation of the triangle formed by pa, pb, pc.
+    //
+    // Input:
+    //   pa, pb, pc  2D points.
+    //
+    // Output:
+    //   Return POSITIVE if pa, pb, pc are counterclockwise oriented.
+    //          NEGATIVE if they are clockwise oriented.
+    //          COLLINEAR if they are collinear.
+    template<typename Vector2D>
+    IGL_INLINE Orientation orient2d(
+        const Eigen::MatrixBase<Vector2D>& pa,
+        const Eigen::MatrixBase<Vector2D>& pb,
+        const Eigen::MatrixBase<Vector2D>& pc);
+
+    // Compute the orientation of the tetrahedron formed by pa, pb, pc, pd.
+    //
+    // Input:
+    //   pa, pb, pc, pd  3D points.
+    //
+    // Output:
+    //   Return POSITIVE if pd is "below" the oriented plane formed by pa, pb and pc.
+    //          NEGATIVE if pd is "above" the plane.
+    //          COPLANAR if pd is on the plane.
+    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);
+
+    // Decide whether a point is inside/outside/on a circle.
+    //
+    // Input:
+    //   pa, pb, pc  2D points that defines an oriented circle.
+    //   pd          2D query point.
+    //
+    // Output:
+    //   Return INSIDE if pd is inside of the circle defined by pa, pb and pc.
+    //          OUSIDE if pd is outside of the circle.
+    //          COCIRCULAR pd is exactly on the circle.
+    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);
+
+    // Decide whether a point is inside/outside/on a sphere.
+    //
+    // Input:
+    //   pa, pb, pc, pd  3D points that defines an oriented sphere.
+    //   pe              3D query point.
+    //
+    // Output:
+    //   Return INSIDE if pe is inside of the sphere defined by pa, pb, pc and pd.
+    //          OUSIDE if pe is outside of the sphere.
+    //          COSPHERICAL pd is exactly on the sphere.
+    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

tests/CMakeLists.txt

@@ -76,6 +76,14 @@ if(LIBIGL_WITH_EMBREE)
   target_link_libraries(libigl_tests PUBLIC igl::embree)
 endif()
 
+if(LIBIGL_WITH_PREDICATES)
+  file(GLOB TEST_SRC_FILES ./include/igl/predicates/*.cpp)
+  file(GLOB TEST_INC_FILES ./include/igl/predicates/*.h ./include/igl/predicates/*.inl)
+  target_sources(libigl_tests PRIVATE ${TEST_SRC_FILES} ${TEST_INC_FILES})
+
+  target_link_libraries(libigl_tests PUBLIC igl::predicates)
+endif()
+
 file(GLOB TEST_SRC_FILES ./include/igl/*.cpp)
 file(GLOB TEST_INC_FILES ./include/igl/*.h ./include/igl/*.inl)
 target_sources(libigl_tests PRIVATE ${TEST_SRC_FILES} ${TEST_INC_FILES})

tests/include/igl/predicates/predicates.cpp

@@ -0,0 +1,54 @@
+#include <test_common.h>
+#include <igl/predicates/predicates.h>
+#include <limits>
+
+TEST_CASE("predicates", "[igl][predicates]") {
+    using namespace igl::predicates;
+    using Scalar = double;
+
+    SECTION("2D") {
+        using Point = Eigen::Matrix<Scalar, 2, 1>;
+        Point a(2,1),b(2,1),c(2,1),d(2,1),e(2,1),f(2,1);
+        a << 0.0, 0.0;
+        b << 1.0, 0.0;
+        c << 1.0, 1.0;
+        d << 2.0, 0.0;
+        e << 0.0, 1.0;
+        f << 0.5, 0.5;
+
+        REQUIRE(orient2d(a, b, c) == Orientation::POSITIVE);
+        REQUIRE(orient2d(a, c, b) == Orientation::NEGATIVE);
+        REQUIRE(orient2d(a, b, b) == Orientation::COLLINEAR);
+        REQUIRE(orient2d(a, a, a) == Orientation::COLLINEAR);
+        REQUIRE(orient2d(a, b, d) == Orientation::COLLINEAR);
+        REQUIRE(orient2d(a, f, c) == Orientation::COLLINEAR);
+
+        REQUIRE(incircle(a,b,c,e) == Orientation::COCIRCULAR);
+        REQUIRE(incircle(a,b,c,a) == Orientation::COCIRCULAR);
+        REQUIRE(incircle(a,b,c,d) == Orientation::OUTSIDE);
+        REQUIRE(incircle(a,b,c,f) == Orientation::INSIDE);
+    }
+
+    SECTION("3D") {
+        using Point = Eigen::Matrix<Scalar, 3, 1>;
+        Point a(3,1),b(3,1),c(3,1),d(3,1),e(3,1),f(3,1);
+        a << 0.0, 0.0, 0.0;
+        b << 1.0, 0.0, 0.0;
+        c << 0.0, 1.0, 0.0;
+        d << 0.0, 0.0, 1.0;
+        e << 1.0, 1.0, 1.0;
+        f << std::numeric_limits<Scalar>::epsilon(), 0.0, 0.0;
+
+        REQUIRE(orient3d(a, b, c, d) == Orientation::NEGATIVE);
+        REQUIRE(orient3d(a, b, d, c) == Orientation::POSITIVE);
+        REQUIRE(orient3d(a, b, d, d) == Orientation::COPLANAR);
+        REQUIRE(orient3d(a, a, a, a) == Orientation::COPLANAR);
+        REQUIRE(orient3d(a, b, f, c) == Orientation::COPLANAR);
+
+        REQUIRE(insphere(a, b, c, d, e) == Orientation::COSPHERICAL);
+        REQUIRE(insphere(a, b, d, e, c) == Orientation::COSPHERICAL);
+        REQUIRE(insphere(b, c, e, d, ((a+b)*0.5).eval()) == Orientation::INSIDE);
+        REQUIRE(insphere(b, c, e, d, (-f).eval()) == Orientation::OUTSIDE);
+        REQUIRE(insphere(f, b, d, c, e) == Orientation::INSIDE);
+    }
+}