split functions into files; add vectorized orient3d (#2479)

include/igl/predicates/IGL_PREDICATES_ASSERT_SCALAR.h

@@ -0,0 +1,15 @@
+#ifndef IGL_PREDICATES_ASSERT_SCALAR_H
+#include <type_traits>
+#ifdef LIBIGL_PREDICATES_USE_FLOAT
+#define IGL_PREDICATES_ASSERT_SCALAR(Vector)                        \
+  static_assert(                                         \
+    std::is_same<typename Vector::Scalar, float>::value, \
+    "Shewchuk's exact predicates only support float")
+#else
+#define IGL_PREDICATES_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
+#endif

include/igl/predicates/Orientation.h

@@ -0,0 +1,30 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// 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_ORIENTATION_H
+#define IGL_PREDICATES_ORIENTATION_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+  namespace predicates {
+    /// Types of orientations and other predicate results.
+    ///
+    /// \fileinfo
+    enum class Orientation {
+      POSITIVE=1, INSIDE=1,
+      NEGATIVE=-1, OUTSIDE=-1,
+      COLLINEAR=0, COPLANAR=0, COCIRCULAR=0, COSPHERICAL=0, DEGENERATE=0
+    };
+  }
+}
+
+
+#endif

include/igl/predicates/delaunay_triangulation.cpp

@@ -8,6 +8,8 @@
 
 #include "delaunay_triangulation.h"
 #include "predicates.h"
+#include "orient2d.h"
+#include "incircle.h"
 #include "../delaunay_triangulation.h"
 
 template<

include/igl/predicates/ear_clipping.cpp

@@ -8,6 +8,7 @@
 
 #include "ear_clipping.h"
 #include "predicates.h"
+#include "orient2d.h"
 #include "segment_segment_intersect.h"
 #include "../turning_number.h"
 

include/igl/predicates/exactinit.cpp

@@ -0,0 +1,24 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// 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 "exactinit.h"
+#include <predicates.h>
+
+IGL_INLINE void igl::predicates::exactinit() {
+  // Thread-safe initialization using Meyers' singleton
+  class MySingleton {
+  public:
+    static MySingleton& instance() {
+      static MySingleton instance;
+      return instance;
+    }
+  private:
+    MySingleton() { ::exactinit(); }
+  };
+  MySingleton::instance();
+}

include/igl/predicates/exactinit.h

@@ -0,0 +1,32 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// 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_EXACTINIT_H
+#define IGL_PREDICATES_EXACTINIT_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+  namespace predicates {
+    /// Initialize internal variable used by predciates. Must be called before
+    /// using exact predicates. It is safe to call this function from multiple
+    /// threads.
+    ///
+    /// \fileinfo
+    IGL_INLINE void exactinit();
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "exactinit.cpp"
+#endif
+
+#endif
+

include/igl/predicates/incircle.cpp

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+//
+// 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 <predicates.h>
+
+namespace igl {
+namespace predicates {
+
+using REAL = IGL_PREDICATES_REAL;
+#include "IGL_PREDICATES_ASSERT_SCALAR.h"
+
+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_PREDICATES_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;
+}
+
+}
+}
+
+#ifdef IGL_STATIC_LIBRARY
+#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_MATRIX(T, R, C) Eigen::Matrix<T, R, C>
+IGL_INCIRCLE(IGL_MATRIX(float, 1, 2));
+IGL_INCIRCLE(IGL_MATRIX(float, 2, 1));
+#ifndef LIBIGL_PREDICATES_USE_FLOAT
+IGL_INCIRCLE(IGL_MATRIX(double, 1, 2));
+IGL_INCIRCLE(IGL_MATRIX(double, 2, 1));
+#endif
+#undef IGL_MATRIX
+#undef IGL_INCIRCLE
+#endif
+

include/igl/predicates/incircle.h

@@ -0,0 +1,43 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+//
+// 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_INCIRCLE_H
+#define IGL_PREDICATES_INCIRCLE_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include "Orientation.h"
+
+namespace igl {
+  namespace predicates {
+    /// Decide whether a point is inside/outside/on a circle.
+    ///
+    /// @param[in] pa  2D point on circle
+    /// @param[in] pb  2D point on circle
+    /// @param[in] pc  2D point on circle
+    /// @param[in] pd  2D point query
+    /// @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.
+    ///
+    /// \fileinfo
+    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);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "incircle.cpp"
+#endif
+
+#endif
+

include/igl/predicates/insphere.cpp

@@ -0,0 +1,59 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+//
+// 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 <predicates.h>
+
+namespace igl {
+namespace predicates {
+
+using REAL = IGL_PREDICATES_REAL;
+#include "IGL_PREDICATES_ASSERT_SCALAR.h"
+
+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_PREDICATES_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;
+}
+
+
+}
+}
+
+#ifdef IGL_STATIC_LIBRARY
+#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_INSPHERE(IGL_MATRIX(float, 1, 3));
+IGL_INSPHERE(IGL_MATRIX(float, 3, 1));
+#ifndef LIBIGL_PREDICATES_USE_FLOAT
+IGL_INSPHERE(IGL_MATRIX(double, 1, 3));
+IGL_INSPHERE(IGL_MATRIX(double, 3, 1));
+#endif
+#undef IGL_MATRIX
+#undef IGL_INSPHERE
+#endif
+
+

include/igl/predicates/insphere.h

@@ -0,0 +1,46 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+//
+// 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_INSPHERE_H
+#define IGL_PREDICATES_INSPHERE_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include "Orientation.h"
+
+namespace igl {
+  namespace predicates {
+    /// Decide whether a point is inside/outside/on a sphere.
+    ///
+    /// @param[in] pa  2D point on sphere
+    /// @param[in] pb  2D point on sphere
+    /// @param[in] pc  2D point on sphere
+    /// @param[in] pd  2D point on sphere
+    /// @param[in] pe  2D point query
+    /// @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.
+    ///
+    /// \fileinfo
+    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 "insphere.cpp"
+#endif
+
+#endif
+
+

include/igl/predicates/lexicographic_triangulation.cpp

@@ -9,6 +9,7 @@
 #include "lexicographic_triangulation.h"
 #include "../lexicographic_triangulation.h"
 #include "predicates.h"
+#include "orient2d.h"
 
 template<
   typename DerivedV,

include/igl/predicates/orient2d.cpp

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+//
+// 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 <predicates.h>
+
+namespace igl {
+namespace predicates {
+
+using REAL = IGL_PREDICATES_REAL;
+#include "IGL_PREDICATES_ASSERT_SCALAR.h"
+
+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_PREDICATES_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;
+}
+
+}
+}
+
+#ifdef IGL_STATIC_LIBRARY
+#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_MATRIX(T, R, C) Eigen::Matrix<T, R, C>
+IGL_ORIENT2D(IGL_MATRIX(float, 1, 2));
+IGL_ORIENT2D(IGL_MATRIX(float, 2, 1));
+#ifndef LIBIGL_PREDICATES_USE_FLOAT
+IGL_ORIENT2D(IGL_MATRIX(double, 1, 2));
+IGL_ORIENT2D(IGL_MATRIX(double, 2, 1));
+#endif
+#undef IGL_MATRIX
+#undef IGL_ORIENT2D
+#endif

include/igl/predicates/orient2d.h

@@ -0,0 +1,40 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+//
+// 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_ORIENT2D_H
+#define IGL_PREDICATES_ORIENT2D_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include "Orientation.h"
+
+namespace igl {
+  namespace predicates {
+    /// Compute the orientation of the triangle formed by pa, pb, pc.
+    ///
+    /// @param[in] pa  2D point on line
+    /// @param[in] pb  2D point on line
+    /// @param[in] pc  2D query point.
+    /// @return POSITIVE if pa, pb, pc are counterclockwise oriented.
+    ///          NEGATIVE if they are clockwise oriented.
+    ///          COLLINEAR if they are collinear.
+    ///
+    /// \fileinfo
+    template<typename Vector2D>
+    IGL_INLINE Orientation orient2d(
+        const Eigen::MatrixBase<Vector2D>& pa,
+        const Eigen::MatrixBase<Vector2D>& pb,
+        const Eigen::MatrixBase<Vector2D>& pc);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "orient2d.cpp"
+#endif
+
+#endif

include/igl/predicates/orient3d.cpp

@@ -0,0 +1,95 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+//
+// 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 "../parallel_for.h"
+#include <predicates.h>
+#include "exactinit.h"
+
+namespace igl {
+namespace predicates {
+
+using REAL = IGL_PREDICATES_REAL;
+#include "IGL_PREDICATES_ASSERT_SCALAR.h"
+
+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_PREDICATES_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 DerivedA,
+   typename DerivedB,
+   typename DerivedC,
+   typename DerivedD,
+   typename DerivedR>
+IGL_INLINE void orient3d(
+    const Eigen::MatrixBase<DerivedA>& A,
+    const Eigen::MatrixBase<DerivedB>& B,
+    const Eigen::MatrixBase<DerivedC>& C,
+    const Eigen::MatrixBase<DerivedD>& D,
+    Eigen::PlainObjectBase<DerivedR>& R)
+{
+  igl::predicates::exactinit();
+  typedef typename DerivedR::Scalar RScalar;
+  typedef typename DerivedA::Scalar Scalar;
+  typedef Eigen::Matrix<Scalar, 1, 3> RowVector3S;
+
+  // max(A.rows(),B.rows(),C.rows(),D.rows()) is the number of points
+  const int np = std::max(
+    std::max(A.rows(), B.rows()),
+    std::max(C.rows(), D.rows()));
+  R.resize(np, 1);
+  igl::parallel_for(np, [&](const int p)
+    {
+      // Not sure if these copies are needed
+      const RowVector3S a = A.row(p % A.rows());
+      const RowVector3S b = B.row(p % B.rows());
+      const RowVector3S c = C.row(p % C.rows());
+      const RowVector3S d = D.row(p % D.rows());
+      // Compute the orientation
+      R(p) = static_cast<RScalar>(igl::predicates::orient3d(a, b, c, d));
+    },1000);
+}
+
+}
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::predicates::orient3d<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(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::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
+
+#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_MATRIX(T, R, C) Eigen::Matrix<T, R, C>
+IGL_ORIENT3D(IGL_MATRIX(float, 1, 3));
+IGL_ORIENT3D(IGL_MATRIX(float, 3, 1));
+#ifndef LIBIGL_PREDICATES_USE_FLOAT
+IGL_ORIENT3D(IGL_MATRIX(double, 1, 3));
+IGL_ORIENT3D(IGL_MATRIX(double, 3, 1));
+#endif
+#undef IGL_MATRIX
+#undef IGL_ORIENT3D
+
+#endif

include/igl/predicates/orient3d.h

@@ -0,0 +1,64 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <[email protected]>
+//
+// 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_ORIENT3D_H
+#define IGL_PREDICATES_ORIENT3D_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include "Orientation.h"
+
+namespace igl {
+  namespace predicates {
+    /// Compute the orientation of the tetrahedron formed by pa, pb, pc, pd.
+    ///
+    /// @param[in] pa  3D point on plane
+    /// @param[in] pb  3D point on plane
+    /// @param[in] pc  3D point on plane
+    /// @param[in] pd  3D query point 
+    ///  @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.
+    ///
+    /// \fileinfo
+    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);
+    /// Compute the orientation of the tetrahedron formed by each 4-tuple of
+    /// points
+    ///
+    /// @param[in] A  #P|1 by 3 matrix of 3D points
+    /// @param[in] B  #P|1 by 3 matrix of 3D points
+    /// @param[in] C  #P|1 by 3 matrix of 3D points
+    /// @param[in] D  #P|1 by 3 matrix of 3D points
+    /// @param[out] R  #P vector of orientations
+    ///
+    template 
+      <typename DerivedA,
+       typename DerivedB,
+       typename DerivedC,
+       typename DerivedD,
+       typename DerivedR>
+    IGL_INLINE void orient3d(
+        const Eigen::MatrixBase<DerivedA>& A,
+        const Eigen::MatrixBase<DerivedB>& B,
+        const Eigen::MatrixBase<DerivedC>& C,
+        const Eigen::MatrixBase<DerivedD>& D,
+        Eigen::PlainObjectBase<DerivedR>& R);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "orient3d.cpp"
+#endif
+
+#endif
+

include/igl/predicates/point_inside_convex_polygon.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 
 #include "point_inside_convex_polygon.h"
+#include "orient2d.h"
 
 template <typename DerivedP, typename DerivedQ>
 IGL_INLINE bool igl::predicates::point_inside_convex_polygon(

include/igl/predicates/predicates.cpp

@@ -1,177 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2019 Qingnan Zhou <[email protected]>
-//
-// 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 "./predicates.h"
-// This is a different file also called predicates.h
-#include <predicates.h>
-#include <type_traits>
-
-namespace igl {
-namespace predicates {
-
-using REAL = IGL_PREDICATES_REAL;
-
-#ifdef LIBIGL_PREDICATES_USE_FLOAT
-#define IGL_PREDICATES_ASSERT_SCALAR(Vector)                        \
-  static_assert(                                         \
-    std::is_same<typename Vector::Scalar, float>::value, \
-    "Shewchuk's exact predicates only support float")
-#else
-#define IGL_PREDICATES_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
-
-IGL_INLINE void exactinit() {
-  // Thread-safe initialization using Meyers' singleton
-  class MySingleton {
-  public:
-    static MySingleton& instance() {
-      static MySingleton instance;
-      return instance;
-    }
-  private:
-    MySingleton() { ::exactinit(); }
-  };
-  MySingleton::instance();
-}
-
-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_PREDICATES_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_PREDICATES_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_PREDICATES_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_PREDICATES_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_PREDICATES_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

@@ -1,108 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2019 Qingnan Zhou <[email protected]>
-//
-// 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_inline.h"
-#include <Eigen/Core>
-
-namespace igl {
-  namespace predicates {
-    /// Types of orientations and other predicate results.
-    ///
-    /// \fileinfo
-    enum class Orientation {
-      POSITIVE=1, INSIDE=1,
-      NEGATIVE=-1, OUTSIDE=-1,
-      COLLINEAR=0, COPLANAR=0, COCIRCULAR=0, COSPHERICAL=0, DEGENERATE=0
-    };
-
-    /// Initialize internal variable used by predciates. Must be called before
-    /// using exact predicates. It is safe to call this function from multiple
-    /// threads.
-    ///
-    /// \fileinfo
-    IGL_INLINE void exactinit();
-
-    /// Compute the orientation of the triangle formed by pa, pb, pc.
-    ///
-    /// @param[in] pa  2D point on line
-    /// @param[in] pb  2D point on line
-    /// @param[in] pc  2D query point.
-    /// @return POSITIVE if pa, pb, pc are counterclockwise oriented.
-    ///          NEGATIVE if they are clockwise oriented.
-    ///          COLLINEAR if they are collinear.
-    ///
-    /// \fileinfo
-    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.
-    ///
-    /// @param[in] pa  2D point on plane
-    /// @param[in] pb  2D point on plane
-    /// @param[in] pc  2D point on plane
-    /// @param[in] pd  2D query point 
-    ///  @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.
-    ///
-    /// \fileinfo
-    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.
-    ///
-    /// @param[in] pa  2D point on circle
-    /// @param[in] pb  2D point on circle
-    /// @param[in] pc  2D point on circle
-    /// @param[in] pd  2D point query
-    /// @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.
-    ///
-    /// \fileinfo
-    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.
-    ///
-    /// @param[in] pa  2D point on sphere
-    /// @param[in] pb  2D point on sphere
-    /// @param[in] pc  2D point on sphere
-    /// @param[in] pd  2D point on sphere
-    /// @param[in] pe  2D point query
-    /// @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.
-    ///
-    /// \fileinfo
-    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

include/igl/predicates/segment_segment_intersect.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 
 #include "segment_segment_intersect.h"
+#include "orient2d.h"
 
 // https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/
 template <typename DerivedP>

include/igl/predicates/triangle_triangle_intersect.cpp

@@ -1,5 +1,8 @@
 #include "triangle_triangle_intersect.h"
 #include "predicates.h"
+#include "exactinit.h"
+#include "orient2d.h"
+#include "orient3d.h"
 #include <Eigen/Core>
 #include <Eigen/Geometry>
 

tests/include/igl/predicates/orient3d.cpp

@@ -0,0 +1,65 @@
+#include <test_common.h>
+#include <igl/predicates/orient2d.h>
+#include <igl/predicates/orient3d.h>
+#include <igl/predicates/incircle.h>
+#include <igl/predicates/insphere.h>
+#include <igl/predicates/exactinit.h>
+#include <limits>
+
+// Didn't have the stamina to break the tests into separate files but they
+// should be
+TEST_CASE("orient3d", "[igl][predicates]") 
+{
+  using namespace igl::predicates;
+  {
+    Eigen::MatrixXd A(3,3);
+    Eigen::MatrixXd B(3,3);
+    Eigen::MatrixXd C(3,3);
+    Eigen::MatrixXd D(3,3);
+    A<<
+      0.0, 0.0, 0.0,
+      0.0, 0.0, 0.0,
+      0.0, 0.0, 0.0;
+    B<<
+      1.0, 0.0, 0.0,
+      1.0, 0.0, 0.0,
+      1.0, 0.0, 0.0;
+    C<<
+      0.0, 1.0, 0.0,
+      0.0, 1.0, 0.0,
+      0.0, 1.0, 0.0;
+    D<<
+      0.0, 0.0, 1.0,
+      0.0, 0.0, -1.0,
+      0.0, 0.0, 0.0;
+    Eigen::VectorXi R;
+    igl::predicates::orient3d(A,B,C,D,R);
+    REQUIRE(R.size() == 3);
+    REQUIRE(R(0) == int(igl::predicates::Orientation::NEGATIVE));
+    REQUIRE(R(1) == int(igl::predicates::Orientation::POSITIVE));
+    REQUIRE(R(2) == int(igl::predicates::Orientation::COPLANAR));
+  }
+  {
+    Eigen::MatrixXd A(1,3);
+    Eigen::MatrixXd B(1,3);
+    Eigen::MatrixXd C(1,3);
+    Eigen::MatrixXd D(3,3);
+    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,
+      0.0, 0.0, -1.0,
+      0.0, 0.0, 0.0;
+    Eigen::VectorXi R;
+    igl::predicates::orient3d(A,B,C,D,R);
+    REQUIRE(R.size() == 3);
+    REQUIRE(R(0) == int(igl::predicates::Orientation::NEGATIVE));
+    REQUIRE(R(1) == int(igl::predicates::Orientation::POSITIVE));
+    REQUIRE(R(2) == int(igl::predicates::Orientation::COPLANAR));
+  }
+}
+

tests/include/igl/predicates/predicates.cpp

@@ -1,7 +1,13 @@
 #include <test_common.h>
-#include <igl/predicates/predicates.h>
+#include <igl/predicates/orient2d.h>
+#include <igl/predicates/orient3d.h>
+#include <igl/predicates/incircle.h>
+#include <igl/predicates/insphere.h>
+#include <igl/predicates/exactinit.h>
 #include <limits>
 
+// Didn't have the stamina to break the tests into separate files but they
+// should be
 TEST_CASE("predicates", "[igl][predicates]") {
     using namespace igl::predicates;
     using Scalar = double;

tests/include/igl/triangle/predicates.cpp

@@ -1,5 +1,5 @@
 #include <test_common.h>
-#include <igl/predicates/predicates.h>
+#include <igl/predicates/exactinit.h>
 #include <igl/triangle/triangulate.h>
 #include <limits>