orientation.hpp

Go to the documentation of this file.

#pragma once
 
#include "shape/emptyshape.hpp"

 
#include <compare>
#include <type_traits>
#include <vector>
 
 
namespace pgl {
 
namespace detail {

template <class... Numbers>
using sign_coordinate_t = detail::promoted_number_t<std::common_type_t<Numbers...>>;
 
template <class AX, class BX, class CX>
using orientation_coordinate_t = sign_coordinate_t<AX, BX, CX>;
template <class AX, class BX>
using dot_coordinate_t = sign_coordinate_t<AX, BX>;

template <class... Numbers>
using incircle_coordinate_t = detail::promoted_number_t<sign_coordinate_t<Numbers...>>;

template <class Ordering>
constexpr int signOf(const Ordering& order) {
    return order > 0 ? 1 : (order < 0 ? -1 : 0);
}

struct Approximate {
    double value = 0.0;  
    double error = 0.0;  
};

struct ApproximatePoint {
    Approximate x;
    Approximate y;
};

constexpr double approximateAbs(double value) { return value < 0.0 ? -value : value; }

inline constexpr double approximateRoundoff = 0x1p-52;

inline constexpr double approximateNarrowing = 0x1p-52;

inline constexpr double approximateConversion = 0x1p-45;

inline constexpr double approximateMargin = 1.0 + 0x1p-20;

constexpr Approximate operator+(const Approximate& a, const Approximate& b) {
    const double value = a.value + b.value;
    return {value, a.error + b.error + approximateRoundoff * approximateAbs(value)};
}

constexpr Approximate operator-(const Approximate& a, const Approximate& b) {
    const double value = a.value - b.value;
    return {value, a.error + b.error + approximateRoundoff * approximateAbs(value)};
}

constexpr Approximate operator*(const Approximate& a, const Approximate& b) {
    const double value = a.value * b.value;
    return {value,
            approximateAbs(a.value) * b.error + approximateAbs(b.value) * a.error +
                a.error * b.error + approximateRoundoff * approximateAbs(value)};
}

template <class Int>
constexpr Approximate approximate(const pgl::Rational<Int>& value) {
    const double quotient = static_cast<double>(value);
    return {quotient, approximateConversion * approximateAbs(quotient)};
}

template <class Number>
constexpr Approximate approximate(const Number& value) {
    if constexpr (std::is_floating_point_v<Number>) {
        if constexpr (numeric_limits<Number>::digits <= numeric_limits<double>::digits &&
                      numeric_limits<Number>::max_exponent <= numeric_limits<double>::max_exponent) {
            return {static_cast<double>(value), 0.0};  // float, double: exact
        } else {
            const double narrowed = static_cast<double>(value);  // long double
            return {narrowed, approximateNarrowing * approximateAbs(narrowed)};
        }
    } else if constexpr (numeric_limits<Number>::is_integer &&
                         numeric_limits<Number>::digits <= numeric_limits<double>::digits) {
        return {static_cast<double>(value), 0.0};  // int and narrower: exact
    } else if constexpr (extended_integral<Number>) {
        const double narrowed = static_cast<double>(value);  // int64_t, int128
        return {narrowed, approximateNarrowing * approximateAbs(narrowed)};
    } else {
        const double narrowed = static_cast<double>(value);  // BigInt and the like
        return {narrowed, approximateConversion * approximateAbs(narrowed)};
    }
}

template <class PointType>
constexpr ApproximatePoint approximatePoint(const PointType& point) {
    return {approximate(point.x()), approximate(point.y())};
}

constexpr std::partial_ordering approximateSign(const Approximate& quantity) {
    const double bound = quantity.error * approximateMargin + 0x1p-1000;
    if (quantity.value > bound) {
        return std::partial_ordering::greater;
    }
    if (quantity.value < -bound) {
        return std::partial_ordering::less;
    }
    return std::partial_ordering::unordered;
}

template <class Coordinate>
inline constexpr bool filtersSign = arbitraryPrecision<Coordinate>;

constexpr std::partial_ordering crossFilter(
    const Approximate& ax, const Approximate& ay,
    const Approximate& bx, const Approximate& by) {
    return approximateSign(ax * by - ay * bx);
}

constexpr std::partial_ordering orientationFilter(
    const ApproximatePoint& a, const ApproximatePoint& b, const ApproximatePoint& c) {
    return crossFilter(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);
}

template <class ANumber, class ALabel, class BNumber, class BLabel, class CNumber, class CLabel>
constexpr std::partial_ordering exactOrientationSign(
    const Point<ANumber, ALabel>& a,
    const Point<BNumber, BLabel>& b,
    const Point<CNumber, CLabel>& c) {
    using Coordinate = orientation_coordinate_t<ANumber, BNumber, CNumber>;
 
    const auto abx = asNumber<Coordinate>(b.x()) - asNumber<Coordinate>(a.x());
    const auto aby = asNumber<Coordinate>(b.y()) - asNumber<Coordinate>(a.y());
    const auto acx = asNumber<Coordinate>(c.x()) - asNumber<Coordinate>(a.x());
    const auto acy = asNumber<Coordinate>(c.y()) - asNumber<Coordinate>(a.y());
 
    return threeWay(abx * acy, aby * acx);
}

template <class PointType, bool Filters>
struct FilteredPoint {
    const PointType* point;
    ApproximatePoint approximation;
};
 
template <class PointType>
struct FilteredPoint<PointType, false> {
    const PointType* point;
};

template <class Coordinate, class PointType>
constexpr FilteredPoint<PointType, filtersSign<Coordinate>> filtered(const PointType& point) {
    if constexpr (filtersSign<Coordinate>) {
        return {&point, approximatePoint(point)};
    } else {
        return {&point};
    }
}

template <class Coordinate, class PointType>
constexpr FilteredPoint<PointType, filtersSign<Coordinate>> filtered(const PointType&&) = delete;

template <class PointType>
constexpr ApproximatePoint approximationOf(const FilteredPoint<PointType, true>& point) {
    return point.approximation;
}
 
template <class PointType>
constexpr ApproximatePoint approximationOf(const FilteredPoint<PointType, false>& point) {
    return approximatePoint(*point.point);
}

template <class APoint, class BPoint, class CPoint, bool Filters>
class DeferredOrientationSign {
public:
    constexpr DeferredOrientationSign(const APoint& a, const BPoint& b, const CPoint& c,
                                      std::partial_ordering proved)
        : a_(&a), b_(&b), c_(&c), proved_(proved) {}

    [[nodiscard]] constexpr bool decided() const {
        return proved_ != std::partial_ordering::unordered;
    }

    [[nodiscard]] constexpr std::partial_ordering value() const {
        return decided() ? proved_ : exactOrientationSign(*a_, *b_, *c_);
    }
 
private:
    const APoint* a_;
    const BPoint* b_;
    const CPoint* c_;
    std::partial_ordering proved_;
};

template <class APoint, class BPoint, class CPoint>
class DeferredOrientationSign<APoint, BPoint, CPoint, false> {
public:
    constexpr DeferredOrientationSign(const APoint& a, const BPoint& b, const CPoint& c)
        : a_(&a), b_(&b), c_(&c) {}
 
    [[nodiscard]] constexpr bool decided() const { return false; }
 
    [[nodiscard]] constexpr std::partial_ordering value() const {
        return exactOrientationSign(*a_, *b_, *c_);
    }
 
private:
    const APoint* a_;
    const BPoint* b_;
    const CPoint* c_;
};

template <class APoint, bool AFilters, class BPoint, bool BFilters, class CPoint, bool CFilters>
constexpr auto orientationSignOf(
    const FilteredPoint<APoint, AFilters>& a,
    const FilteredPoint<BPoint, BFilters>& b,
    const FilteredPoint<CPoint, CFilters>& c) {
    // The gate is this predicate's own, read at the type its exact fallback
    // would evaluate in — never the operands', which say only whether an
    // approximation was already at hand. Deciding it operand-side would both
    // filter where @ref orientationSign does not (an int64_t coordinate keeps
    // approximations for its in-circle tests, whose determinant needs them,
    // but not for these, whose fits in an int128) and fail to filter where it
    // does (`int` vertices against a rational query point).
    using Coordinate = orientation_coordinate_t<typename APoint::NumberType,
                                                typename BPoint::NumberType,
                                                typename CPoint::NumberType>;
    if constexpr (filtersSign<Coordinate>) {
        return DeferredOrientationSign<APoint, BPoint, CPoint, true>{
            *a.point, *b.point, *c.point,
            orientationFilter(approximationOf(a), approximationOf(b), approximationOf(c))};
    } else {
        return DeferredOrientationSign<APoint, BPoint, CPoint, false>{
            *a.point, *b.point, *c.point};
    }
}

template <class... Signs>
constexpr bool allDecided(const Signs&... signs) {
    return (signs.decided() && ...);
}

template <class Coordinate, class PointType>
constexpr FilteredPoint<PointType, filtersSign<Coordinate>> filtered(
    const PointType& point,
    const std::vector<ApproximatePoint>& approximations,
    std::size_t index) {
    if constexpr (filtersSign<Coordinate>) {
        return {&point, approximations[index]};
    } else {
        return {&point};
    }
}

constexpr std::partial_ordering dotFilter(
    const Approximate& ax, const Approximate& ay,
    const Approximate& bx, const Approximate& by) {
    return approximateSign(ax * bx + ay * by);
}
}  // namespace detail

template <class ANumber, class ALabel, class BNumber, class BLabel, class CNumber, class CLabel>
constexpr auto orientationDeterminant(
    const Point<ANumber, ALabel>& a,
    const Point<BNumber, BLabel>& b,
    const Point<CNumber, CLabel>& c) {
    using Coordinate = detail::orientation_coordinate_t<ANumber, BNumber, CNumber>;
 
    const auto abx = detail::asNumber<Coordinate>(b.x()) - detail::asNumber<Coordinate>(a.x());
    const auto aby = detail::asNumber<Coordinate>(b.y()) - detail::asNumber<Coordinate>(a.y());
    const auto acx = detail::asNumber<Coordinate>(c.x()) - detail::asNumber<Coordinate>(a.x());
    const auto acy = detail::asNumber<Coordinate>(c.y()) - detail::asNumber<Coordinate>(a.y());
 
    return abx * acy - aby * acx;
}

template <class ANumber, class ALabel, class BNumber, class BLabel, class CNumber, class CLabel>
constexpr std::partial_ordering orientationSign(
    const Point<ANumber, ALabel> &a,
    const Point<BNumber, BLabel> &b,
    const Point<CNumber, CLabel> &c) {
    using Coordinate = detail::orientation_coordinate_t<ANumber, BNumber, CNumber>;
 
    // A double evaluation that carries its own error bound settles the sign for
    // all but the near-degenerate inputs, sparing the exact arithmetic below.
    // It only reports a sign it has proved, so this is a shortcut, not an
    // approximation; see @ref detail::approximateSign.
    if constexpr (detail::filtersSign<Coordinate>) {
        const detail::Approximate ax = detail::approximate(a.x());
        const detail::Approximate ay = detail::approximate(a.y());
        const std::partial_ordering filtered = detail::crossFilter(
            detail::approximate(b.x()) - ax, detail::approximate(b.y()) - ay,
            detail::approximate(c.x()) - ax, detail::approximate(c.y()) - ay);
        if (filtered != std::partial_ordering::unordered) {
            return filtered;
        }
    }
 
    return detail::exactOrientationSign(a, b, c);
}

template <class UNumber, class ULabel, class VNumber, class VLabel>
constexpr std::partial_ordering crossSign(
    const Point<UNumber, ULabel>& u,
    const Point<VNumber, VLabel>& v) {
    using Coordinate = detail::sign_coordinate_t<UNumber, VNumber>;
 
    if constexpr (detail::filtersSign<Coordinate>) {  // see orientationSign
        const std::partial_ordering filtered = detail::crossFilter(
            detail::approximate(u.x()), detail::approximate(u.y()),
            detail::approximate(v.x()), detail::approximate(v.y()));
        if (filtered != std::partial_ordering::unordered) {
            return filtered;
        }
    }
 
    return detail::threeWay(detail::asNumber<Coordinate>(u.x()) * detail::asNumber<Coordinate>(v.y()),
                            detail::asNumber<Coordinate>(u.y()) * detail::asNumber<Coordinate>(v.x()));
}

template <class ANumber, class ALabel, class BNumber, class BLabel,
          class PNumber, class PLabel, class QNumber, class QLabel>
constexpr std::partial_ordering crossSign(
    const Point<ANumber, ALabel>& a,
    const Point<BNumber, BLabel>& b,
    const Point<PNumber, PLabel>& p,
    const Point<QNumber, QLabel>& q) {
    using Coordinate = detail::sign_coordinate_t<ANumber, BNumber, PNumber, QNumber>;
 
    if constexpr (detail::filtersSign<Coordinate>) {  // see orientationSign
        const std::partial_ordering filtered = detail::crossFilter(
            detail::approximate(b.x()) - detail::approximate(a.x()),
            detail::approximate(b.y()) - detail::approximate(a.y()),
            detail::approximate(q.x()) - detail::approximate(p.x()),
            detail::approximate(q.y()) - detail::approximate(p.y()));
        if (filtered != std::partial_ordering::unordered) {
            return filtered;
        }
    }
 
    const auto abx = detail::asNumber<Coordinate>(b.x()) - detail::asNumber<Coordinate>(a.x());
    const auto aby = detail::asNumber<Coordinate>(b.y()) - detail::asNumber<Coordinate>(a.y());
    const auto pqx = detail::asNumber<Coordinate>(q.x()) - detail::asNumber<Coordinate>(p.x());
    const auto pqy = detail::asNumber<Coordinate>(q.y()) - detail::asNumber<Coordinate>(p.y());
 
    return detail::threeWay(abx * pqy, aby * pqx);
}

template <class ANumber, class ALabel, class BNumber, class BLabel, class CNumber, class CLabel>
constexpr bool collinear(
    const Point<ANumber, ALabel>& a,
    const Point<BNumber, BLabel>& b,
    const Point<CNumber, CLabel>& c) {
    return orientationSign(a, b, c) == 0;
}

template <class ANumber, class ALabel, class BNumber, class BLabel>
constexpr bool sameDirection(
    const Point<ANumber, ALabel>& a1,
    const Point<ANumber, ALabel>& a2,
    const Point<BNumber, BLabel>& b1,
    const Point<BNumber, BLabel>& b2) {
    return crossSign(a1, a2, b1, b2) == 0;
}

template <class ANumber, class ALabel, class BNumber, class BLabel>
constexpr std::partial_ordering dotSign(
    const Point<ANumber, ALabel>& a,
    const Point<BNumber, BLabel>& b) {
    using Coordinate = detail::sign_coordinate_t<ANumber, BNumber>;
 
    if constexpr (detail::filtersSign<Coordinate>) {  // see orientationSign
        const std::partial_ordering filtered = detail::dotFilter(
            detail::approximate(a.x()), detail::approximate(a.y()),
            detail::approximate(b.x()), detail::approximate(b.y()));
        if (filtered != std::partial_ordering::unordered) {
            return filtered;
        }
    }
 
    const auto x = detail::asNumber<Coordinate>(a.x()) * detail::asNumber<Coordinate>(b.x());
    const auto y = detail::asNumber<Coordinate>(a.y()) * detail::asNumber<Coordinate>(b.y());
 
    return detail::threeWay(x, -y);
}

template <class ANumber, class ALabel, class BNumber, class BLabel,
          class PNumber, class PLabel, class QNumber, class QLabel>
constexpr std::partial_ordering dotSign(
    const Point<ANumber, ALabel>& a,
    const Point<BNumber, BLabel>& b,
    const Point<PNumber, PLabel>& p,
    const Point<QNumber, QLabel>& q) {
    using Coordinate = detail::sign_coordinate_t<ANumber, BNumber, PNumber, QNumber>;
 
    if constexpr (detail::filtersSign<Coordinate>) {  // see orientationSign
        const std::partial_ordering filtered = detail::dotFilter(
            detail::approximate(b.x()) - detail::approximate(a.x()),
            detail::approximate(b.y()) - detail::approximate(a.y()),
            detail::approximate(q.x()) - detail::approximate(p.x()),
            detail::approximate(q.y()) - detail::approximate(p.y()));
        if (filtered != std::partial_ordering::unordered) {
            return filtered;
        }
    }
 
    const auto abx = detail::asNumber<Coordinate>(b.x()) - detail::asNumber<Coordinate>(a.x());
    const auto aby = detail::asNumber<Coordinate>(b.y()) - detail::asNumber<Coordinate>(a.y());
    const auto pqx = detail::asNumber<Coordinate>(q.x()) - detail::asNumber<Coordinate>(p.x());
    const auto pqy = detail::asNumber<Coordinate>(q.y()) - detail::asNumber<Coordinate>(p.y());
 
    return detail::threeWay(abx * pqx, -(aby * pqy));
}
 
namespace detail {

constexpr std::partial_ordering inCircleFilter(
    const ApproximatePoint& a,
    const ApproximatePoint& b,
    const ApproximatePoint& c,
    const ApproximatePoint& d) {
    const Approximate adx = a.x - d.x;
    const Approximate ady = a.y - d.y;
    const Approximate bdx = b.x - d.x;
    const Approximate bdy = b.y - d.y;
    const Approximate cdx = c.x - d.x;
    const Approximate cdy = c.y - d.y;
    const Approximate abdet = adx * bdy - bdx * ady;
    const Approximate bcdet = bdx * cdy - cdx * bdy;
    const Approximate cadet = cdx * ady - adx * cdy;
    const Approximate alift = adx * adx + ady * ady;
    const Approximate blift = bdx * bdx + bdy * bdy;
    const Approximate clift = cdx * cdx + cdy * cdy;
    return approximateSign(alift * bcdet + blift * cadet + clift * abdet);
}

template <class ANumber, class ALabel, class BNumber, class BLabel, class CNumber, class CLabel, class DNumber, class DLabel>
constexpr std::partial_ordering exactInCircleSign(
    const Point<ANumber, ALabel>& a,
    const Point<BNumber, BLabel>& b,
    const Point<CNumber, CLabel>& c,
    const Point<DNumber, DLabel>& d) {
    using Coordinate = incircle_coordinate_t<ANumber, BNumber, CNumber, DNumber>;
 
    const auto adx = asNumber<Coordinate>(a.x()) - asNumber<Coordinate>(d.x());
    const auto ady = asNumber<Coordinate>(a.y()) - asNumber<Coordinate>(d.y());
    const auto bdx = asNumber<Coordinate>(b.x()) - asNumber<Coordinate>(d.x());
    const auto bdy = asNumber<Coordinate>(b.y()) - asNumber<Coordinate>(d.y());
    const auto cdx = asNumber<Coordinate>(c.x()) - asNumber<Coordinate>(d.x());
    const auto cdy = asNumber<Coordinate>(c.y()) - asNumber<Coordinate>(d.y());
    const auto abdet = adx * bdy - bdx * ady;
    const auto bcdet = bdx * cdy - cdx * bdy;
    const auto cadet = cdx * ady - adx * cdy;
    const auto alift = adx * adx + ady * ady;
    const auto blift = bdx * bdx + bdy * bdy;
    const auto clift = cdx * cdx + cdy * cdy;
    return threeWay(alift * bcdet + blift * cadet, -clift * abdet);
}

template <class APoint, bool AFilters, class BPoint, bool BFilters,
          class CPoint, bool CFilters, class DPoint, bool DFilters>
constexpr std::partial_ordering inCircleSignOf(
    const FilteredPoint<APoint, AFilters>& a,
    const FilteredPoint<BPoint, BFilters>& b,
    const FilteredPoint<CPoint, CFilters>& c,
    const FilteredPoint<DPoint, DFilters>& d) {
    using Coordinate = incircle_coordinate_t<typename APoint::NumberType,
                                             typename BPoint::NumberType,
                                             typename CPoint::NumberType,
                                             typename DPoint::NumberType>;
    if constexpr (filtersSign<Coordinate>) {  // see orientationSignOf
        const std::partial_ordering filtered = inCircleFilter(
            approximationOf(a), approximationOf(b), approximationOf(c), approximationOf(d));
        if (filtered != std::partial_ordering::unordered) {
            return filtered;
        }
    }
    return exactInCircleSign(*a.point, *b.point, *c.point, *d.point);
}
 
}  // namespace detail

template <class ANumber, class ALabel, class BNumber, class BLabel, class CNumber, class CLabel, class DNumber, class DLabel>
constexpr auto inCircleDeterminant(
    const Point<ANumber, ALabel>& a,
    const Point<BNumber, BLabel>& b,
    const Point<CNumber, CLabel>& c,
    const Point<DNumber, DLabel>& d) {
    using Coordinate = detail::incircle_coordinate_t<ANumber, BNumber, CNumber, DNumber>;
 
    const auto adx = detail::asNumber<Coordinate>(a.x()) - detail::asNumber<Coordinate>(d.x());
    const auto ady = detail::asNumber<Coordinate>(a.y()) - detail::asNumber<Coordinate>(d.y());
    const auto bdx = detail::asNumber<Coordinate>(b.x()) - detail::asNumber<Coordinate>(d.x());
    const auto bdy = detail::asNumber<Coordinate>(b.y()) - detail::asNumber<Coordinate>(d.y());
    const auto cdx = detail::asNumber<Coordinate>(c.x()) - detail::asNumber<Coordinate>(d.x());
    const auto cdy = detail::asNumber<Coordinate>(c.y()) - detail::asNumber<Coordinate>(d.y());
    const auto abdet = adx * bdy - bdx * ady;
    const auto bcdet = bdx * cdy - cdx * bdy;
    const auto cadet = cdx * ady - adx * cdy;
    const auto alift = adx * adx + ady * ady;
    const auto blift = bdx * bdx + bdy * bdy;
    const auto clift = cdx * cdx + cdy * cdy;
    return alift * bcdet + blift * cadet + clift * abdet;
}

template <class ANumber, class ALabel, class BNumber, class BLabel, class CNumber, class CLabel, class DNumber, class DLabel>
constexpr std::partial_ordering inCircleSign(
    const Point<ANumber, ALabel>& a,
    const Point<BNumber, BLabel>& b,
    const Point<CNumber, CLabel>& c,
    const Point<DNumber, DLabel>& d) {
    using Coordinate = detail::incircle_coordinate_t<ANumber, BNumber, CNumber, DNumber>;
 
    // A double evaluation that carries its own error bound settles the sign for
    // all but the near-degenerate inputs, sparing the exact arithmetic below.
    // It only reports a sign it has proved, so this is a shortcut, not an
    // approximation; see @ref detail::inCircleFilter. A caller running this
    // repeatedly over one point array should keep the conversions — see
    // @ref detail::filtered — and call @ref detail::inCircleSignOf instead.
    if constexpr (detail::filtersSign<Coordinate>) {
        const std::partial_ordering filtered = detail::inCircleFilter(
            detail::approximatePoint(a), detail::approximatePoint(b),
            detail::approximatePoint(c), detail::approximatePoint(d));
        if (filtered != std::partial_ordering::unordered) {
            return filtered;
        }
    }
 
    return detail::exactInCircleSign(a, b, c, d);
}
 
}  // namespace pgl