predicates.hpp

#pragma once
 
#include <Eigen/Core>
 
namespace wmtk::utils::predicates {
 
enum class Orientation {
    POSITIVE = 1,
    INSIDE = 1,
    NEGATIVE = -1,
    OUTSIDE = -1,
    COLLINEAR = 0,
    COPLANAR = 0,
    COCIRCULAR = 0,
    COSPHERICAL = 0,
    DEGENERATE = 0
};
 
namespace detail {
 
int orient2d_sign(double ax, double ay, double bx, double by, double cx, double cy);
int orient3d_sign(
    double ax,
    double ay,
    double az,
    double bx,
    double by,
    double bz,
    double cx,
    double cy,
    double cz,
    double dx,
    double dy,
    double dz);
 
inline Orientation from_sign(const int s)
{
    return s > 0 ? Orientation::POSITIVE : (s < 0 ? Orientation::NEGATIVE : Orientation::COLLINEAR);
}
 
} // namespace detail
 
inline void exactinit() {}
 
template <typename D0, typename D1, typename D2>
inline Orientation orient2d(
    const Eigen::MatrixBase<D0>& a,
    const Eigen::MatrixBase<D1>& b,
    const Eigen::MatrixBase<D2>& c)
{
    return detail::from_sign(detail::orient2d_sign(a[0], a[1], b[0], b[1], c[0], c[1]));
}
 
template <typename D0, typename D1, typename D2, typename D3>
inline Orientation orient3d(
    const Eigen::MatrixBase<D0>& a,
    const Eigen::MatrixBase<D1>& b,
    const Eigen::MatrixBase<D2>& c,
    const Eigen::MatrixBase<D3>& d)
{
    return detail::from_sign(
        detail::
            orient3d_sign(a[0], a[1], a[2], b[0], b[1], b[2], c[0], c[1], c[2], d[0], d[1], d[2]));
}
 
template <typename D0, typename D1, typename D2>
inline bool is_degenerate(
    const Eigen::MatrixBase<D0>& v0,
    const Eigen::MatrixBase<D1>& v1,
    const Eigen::MatrixBase<D2>& v2)
{
    if constexpr (Eigen::MatrixBase<D0>::RowsAtCompileTime == 3) {
        // Three points in 3D are collinear only if they are collinear in all three
        // coordinate-plane projections; any one projection can flatten a real triangle.
        for (int dim = 0; dim < 3; ++dim) {
            const Eigen::Vector2d p0(v0[dim], v0[(dim + 1) % 3]);
            const Eigen::Vector2d p1(v1[dim], v1[(dim + 1) % 3]);
            const Eigen::Vector2d p2(v2[dim], v2[(dim + 1) % 3]);
            if (orient2d(p0, p1, p2) != Orientation::COLLINEAR) {
                return false;
            }
        }
        return true;
    } else {
        return orient2d(v0, v1, v2) == Orientation::COLLINEAR;
    }
}
 
template <typename D0, typename D1, typename D2, typename D3>
inline bool is_degenerate(
    const Eigen::MatrixBase<D0>& v0,
    const Eigen::MatrixBase<D1>& v1,
    const Eigen::MatrixBase<D2>& v2,
    const Eigen::MatrixBase<D3>& v3)
{
    return orient3d(v0, v1, v2, v3) == Orientation::COPLANAR;
}
 
inline bool faces_along(
    const Eigen::Vector3d& a,
    const Eigen::Vector3d& b,
    const Eigen::Vector3d& c,
    const Eigen::Vector3d& n)
{
    const Eigen::Vector3d below = a - (1 + a.lpNorm<Eigen::Infinity>()) * n.normalized();
    return orient3d(a, b, c, below) == Orientation::POSITIVE;
}
 
} // namespace wmtk::utils::predicates