Update predicates and avoid heap allocations.

cmake/LibiglDownloadExternal.cmake

@@ -156,7 +156,7 @@ endfunction()
 function(igl_download_predicates)
 	igl_download_project(predicates
 		GIT_REPOSITORY https://github.com/libigl/libigl-predicates
-		GIT_TAG        55c9251878f54cb3c3b8f25f93a7ce4dd60d9317
+		GIT_TAG        ce494af8c90c2cc699146b61c33de102fa96635b
 	)
 endfunction()
 

include/igl/predicates/predicates.cpp

@@ -7,27 +7,30 @@
 // 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 = double;
+using REAL = IGL_PREDICATES_REAL;
 
 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);
+    static_assert(
+            std::is_same<typename Vector2D::Scalar, double>::value ||
+            std::is_same<typename Vector2D::Scalar, float>::value,
+            "Shewchuk's exact predicates only support float and double");
 
-    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>();
+    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]};
 
-    auto r = ::orient2d(
-            const_cast<REAL*>(a.data()),
-            const_cast<REAL*>(b.data()),
-            const_cast<REAL*>(c.data()));
+    const auto r = ::orient2d(a.data(), b.data(), c.data());
 
     if (r > 0) return Orientation::POSITIVE;
     else if (r < 0) return Orientation::NEGATIVE;
@@ -40,18 +43,19 @@ IGL_INLINE Orientation orient3d(
         const Eigen::MatrixBase<Vector3D>& pb,
         const Eigen::MatrixBase<Vector3D>& pc,
         const Eigen::MatrixBase<Vector3D>& pd) {
+    EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Vector3D, 3);
+    static_assert(
+            std::is_same<typename Vector3D::Scalar, double>::value ||
+            std::is_same<typename Vector3D::Scalar, float>::value,
+            "Shewchuk's exact predicates only support float and double");
 
-    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>();
+    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]};
 
-    auto r = ::orient3d(
-            const_cast<REAL*>(a.data()),
-            const_cast<REAL*>(b.data()),
-            const_cast<REAL*>(c.data()),
-            const_cast<REAL*>(d.data()));
+    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;
@@ -64,18 +68,19 @@ IGL_INLINE Orientation incircle(
         const Eigen::MatrixBase<Vector2D>& pb,
         const Eigen::MatrixBase<Vector2D>& pc,
         const Eigen::MatrixBase<Vector2D>& pd) {
+    EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Vector2D, 2);
+    static_assert(
+            std::is_same<typename Vector2D::Scalar, double>::value ||
+            std::is_same<typename Vector2D::Scalar, float>::value,
+            "Shewchuk's exact predicates only support float and double");
 
-    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>();
+    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]};
 
-    auto r = ::incircle(
-            const_cast<REAL*>(a.data()),
-            const_cast<REAL*>(b.data()),
-            const_cast<REAL*>(c.data()),
-            const_cast<REAL*>(d.data()));
+    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;
@@ -89,20 +94,20 @@ IGL_INLINE Orientation insphere(
         const Eigen::MatrixBase<Vector3D>& pc,
         const Eigen::MatrixBase<Vector3D>& pd,
         const Eigen::MatrixBase<Vector3D>& pe) {
+    EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Vector3D, 3);
+    static_assert(
+            std::is_same<typename Vector3D::Scalar, double>::value ||
+            std::is_same<typename Vector3D::Scalar, float>::value,
+            "Shewchuk's exact predicates only support float and double");
 
-    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>();
+    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]};
 
-    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()));
+    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;
@@ -114,155 +119,55 @@ IGL_INLINE Orientation insphere(
 
 #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::orient2d<Eigen::Matrix<double, -1, -1, 0, -1, -1>>
-(
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&
-);
-template
-igl::predicates::Orientation
-igl::predicates::orient2d<Eigen::Matrix<float, -1, -1, 0, -1, -1>>
-(
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -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::orient3d<Eigen::Matrix<double, -1, -1, 0, -1, -1>>
-(
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&
-);
-template
-igl::predicates::Orientation
-igl::predicates::orient3d<Eigen::Matrix<float, -1, -1, 0, -1, -1>>
-(
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -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::incircle<Eigen::Matrix<double, -1, -1, 0, -1, -1>>
-(
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&
-);
-template
-igl::predicates::Orientation
-igl::predicates::incircle<Eigen::Matrix<float, -1, -1, 0, -1, -1>>
-(
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -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>>&
-);
-template
-igl::predicates::Orientation
-igl::predicates::insphere<Eigen::Matrix<double, -1, -1, 0, -1, -1>>
-(
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&
-);
-template
-igl::predicates::Orientation
-igl::predicates::insphere<Eigen::Matrix<float, -1, -1, 0, -1, -1>>
-(
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&,
-const Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>>&
-);
+#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_Vector2D Eigen::Matrix<float, 1, 2>
+IGL_ORIENT2D(IGL_Vector2D);
+IGL_INCIRCLE(IGL_Vector2D);
+#undef IGL_Vector2D
+
+#define IGL_Vector2D Eigen::Matrix<double, 1, 2>
+IGL_ORIENT2D(IGL_Vector2D);
+IGL_INCIRCLE(IGL_Vector2D);
+#undef IGL_Vector2D
+
+#define IGL_Vector2D Eigen::Matrix<float, 2, 1>
+IGL_ORIENT2D(IGL_Vector2D);
+IGL_INCIRCLE(IGL_Vector2D);
+#undef IGL_Vector2D
+
+#define IGL_Vector2D Eigen::Matrix<double, 2, 1>
+IGL_ORIENT2D(IGL_Vector2D);
+IGL_INCIRCLE(IGL_Vector2D);
+#undef IGL_Vector2D
+
+#define IGL_Vector3D Eigen::Matrix<float, 1, 3>
+IGL_ORIENT3D(IGL_Vector3D);
+IGL_INSPHERE(IGL_Vector3D);
+#undef IGL_Vector3D
+
+#define IGL_Vector3D Eigen::Matrix<float, 3, 1>
+IGL_ORIENT3D(IGL_Vector3D);
+IGL_INSPHERE(IGL_Vector3D);
+#undef IGL_Vector3D
+
+#define IGL_Vector3D Eigen::Matrix<double, 1, 3>
+IGL_ORIENT3D(IGL_Vector3D);
+IGL_INSPHERE(IGL_Vector3D);
+#undef IGL_Vector3D
+
+#define IGL_Vector3D Eigen::Matrix<double, 3, 1>
+IGL_ORIENT3D(IGL_Vector3D);
+IGL_INSPHERE(IGL_Vector3D);
+#undef IGL_Vector3D
+
+#undef IGL_ORIENT2D
+#undef IGL_ORIENT3D
+#undef IGL_INCIRCLE
+#undef IGL_INSPHERE
 
 #endif

include/igl/predicates/predicates.h

@@ -26,7 +26,7 @@ namespace igl {
     //   pa, pb, pc  2D points.
     //
     // Output:
-    //   Return POSITIVE if pa, pb, pc is counterclockwise oriented.
+    //   Return POSITIVE if pa, pb, pc are counterclockwise oriented.
     //          NEGATIVE if they are clockwise oriented.
     //          COLLINEAR if they are collinear.
     template<typename Vector2D>

tests/include/igl/predicates/predicates.cpp

@@ -7,7 +7,7 @@ TEST_CASE("predicates", "[igl][predicates]") {
     using Scalar = double;
 
     SECTION("2D") {
-        using Point = Eigen::Matrix<Scalar, -1, -1>;
+        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;
@@ -30,7 +30,7 @@ TEST_CASE("predicates", "[igl][predicates]") {
     }
 
     SECTION("3D") {
-        using Point = Eigen::Matrix<Scalar, -1, -1>;
+        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;