measures.hpp

Go to the documentation of this file.

#pragma once
 
#include <cmath>
#include <optional>
#include <type_traits>
 
#include "implementation/minkowski.hpp"
#include "implementation/orientation.hpp"

 
 
namespace pgl {
 
// -----------------------------------------------------------------------------
// Point
 
template <class LeftNumber, class LeftLabel, class RightNumber, class RightLabel>
constexpr auto operator*(const Point<LeftNumber, LeftLabel>& left, const Point<RightNumber, RightLabel>& right) {
    return left.x() * right.x() + left.y() * right.y();
}
 
// -----------------------------------------------------------------------------
// Segment
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber Segment<PointType, LabelType>::area() const {
    return ResultNumber{};
}
 
template <class PointType, class LabelType>
constexpr typename Segment<PointType, LabelType>::NumberType Segment<PointType, LabelType>::twiceArea() const {
    return NumberType{};
}
 
template <class PointType, class LabelType>
constexpr auto Segment<PointType, LabelType>::squaredLength() const {
    return min().template squaredDistance<NumberType>(max());
}
 
template <class PointType, class LabelType>
template <class ApproximateNumber>
ApproximateNumber Segment<PointType, LabelType>::length() const {
    return min().template distance<ApproximateNumber>(max());
}
 
template <class PointType, class LabelType>
constexpr auto Segment<PointType, LabelType>::lengthL1() const {
    return min().distanceL1(max());
}
 
template <class PointType, class LabelType>
constexpr auto Segment<PointType, LabelType>::lengthLInf() const {
    return min().distanceLInf(max());
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber Segment<PointType, LabelType>::slope() const {
    const auto dy = detail::asNumber<ResultNumber>(max().y()) - detail::asNumber<ResultNumber>(min().y());
    const auto dx = detail::asNumber<ResultNumber>(max().x()) - detail::asNumber<ResultNumber>(min().x());
    return dy / dx;
}
 
template <class PointType, class LabelType>
constexpr Segment<PointType, LabelType> Segment<PointType, LabelType>::diameter() const {
    return *this;
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Segment<PointType, LabelType>::midpoint() const {
    return Point<ResultNumber>(
        (detail::asNumber<ResultNumber>(min().x()) + detail::asNumber<ResultNumber>(max().x())) / static_cast<ResultNumber>(2),
        (detail::asNumber<ResultNumber>(min().y()) + detail::asNumber<ResultNumber>(max().y())) / static_cast<ResultNumber>(2)
    );
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Segment<PointType, LabelType>::pointInside() const {
    return midpoint<ResultNumber>();
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool Segment<PointType, LabelType>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    const auto witness = pointInside<NumberType>();
    if (interiorContains(witness)) {
        return shape.interiorContains(witness);
    }
    // Integer truncation rounded the midpoint onto an endpoint; scaling by 2
    // makes it exact without changing the containment relation.
    return (shape * 2).interiorContains((*this * 2).template pointInside<NumberType>());
}
 
// -----------------------------------------------------------------------------
// OrientedSegment
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber OrientedSegment<PointType, LabelType>::area() const {
    return ResultNumber{};
}
 
template <class PointType, class LabelType>
constexpr typename OrientedSegment<PointType, LabelType>::NumberType OrientedSegment<PointType, LabelType>::twiceArea() const {
    return NumberType{};
}
 
template <class PointType, class LabelType>
constexpr auto OrientedSegment<PointType, LabelType>::squaredLength() const {
    return source().template squaredDistance<NumberType>(target());
}
 
template <class PointType, class LabelType>
template <class ApproximateNumber>
ApproximateNumber OrientedSegment<PointType, LabelType>::length() const {
    return source().template distance<ApproximateNumber>(target());
}
 
template <class PointType, class LabelType>
constexpr auto OrientedSegment<PointType, LabelType>::lengthL1() const {
    return source().distanceL1(target());
}
 
template <class PointType, class LabelType>
constexpr auto OrientedSegment<PointType, LabelType>::lengthLInf() const {
    return source().distanceLInf(target());
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber OrientedSegment<PointType, LabelType>::slope() const {
    const auto dy = detail::asNumber<ResultNumber>(target().y()) - detail::asNumber<ResultNumber>(source().y());
    const auto dx = detail::asNumber<ResultNumber>(target().x()) - detail::asNumber<ResultNumber>(source().x());
    return dy / dx;
}
 
template <class PointType, class LabelType>
constexpr Segment<PointType> OrientedSegment<PointType, LabelType>::diameter() const {
    return static_cast<Segment<PointType>>(*this);
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> OrientedSegment<PointType, LabelType>::midpoint() const {
    return Point<ResultNumber>(
        (detail::asNumber<ResultNumber>(source().x()) + detail::asNumber<ResultNumber>(target().x())) / static_cast<ResultNumber>(2),
        (detail::asNumber<ResultNumber>(source().y()) + detail::asNumber<ResultNumber>(target().y())) / static_cast<ResultNumber>(2));
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> OrientedSegment<PointType, LabelType>::pointInside() const {
    return midpoint<ResultNumber>();
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool OrientedSegment<PointType, LabelType>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    // Orientation does not affect the point set, so defer to the segment.
    return asSegment().pointInsideInteriorContainedIn(shape);
}
 
// -----------------------------------------------------------------------------
// Line
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber Line<PointType, LabelType>::area() const {
    return ResultNumber{};
}
 
template <class PointType, class LabelType>
constexpr typename Line<PointType, LabelType>::NumberType Line<PointType, LabelType>::twiceArea() const {
    return NumberType{};
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber Line<PointType, LabelType>::slope() const {
    const auto dy = detail::asNumber<ResultNumber>(max().y()) - detail::asNumber<ResultNumber>(min().y());
    const auto dx = detail::asNumber<ResultNumber>(max().x()) - detail::asNumber<ResultNumber>(min().x());
    return dy / dx;
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Line<PointType, LabelType>::pointInside() const {
    return Point<ResultNumber>(min());
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool Line<PointType, LabelType>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    // pointInside() lies exactly on the line, so it is a valid interior witness.
    return shape.interiorContains(pointInside<NumberType>());
}
 
// -----------------------------------------------------------------------------
// OrientedLine
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber OrientedLine<PointType, LabelType>::area() const {
    return ResultNumber{};
}
 
template <class PointType, class LabelType>
constexpr typename OrientedLine<PointType, LabelType>::NumberType OrientedLine<PointType, LabelType>::twiceArea() const {
    return NumberType{};
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber OrientedLine<PointType, LabelType>::slope() const {
    const auto dy = detail::asNumber<ResultNumber>(target().y()) - detail::asNumber<ResultNumber>(source().y());
    const auto dx = detail::asNumber<ResultNumber>(target().x()) - detail::asNumber<ResultNumber>(source().x());
    return dy / dx;
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> OrientedLine<PointType, LabelType>::pointInside() const {
    return Point<ResultNumber>(source());
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool OrientedLine<PointType, LabelType>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    // Orientation does not affect the point set, so defer to the line.
    return asLine().pointInsideInteriorContainedIn(shape);
}
 
// -----------------------------------------------------------------------------
// Ray
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber Ray<PointType, LabelType>::area() const {
    return ResultNumber{};
}
 
template <class PointType, class LabelType>
constexpr typename Ray<PointType, LabelType>::NumberType Ray<PointType, LabelType>::twiceArea() const {
    return NumberType{};
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber Ray<PointType, LabelType>::slope() const {
    const auto dy = detail::asNumber<ResultNumber>(target().y()) - detail::asNumber<ResultNumber>(source().y());
    const auto dx = detail::asNumber<ResultNumber>(target().x()) - detail::asNumber<ResultNumber>(source().x());
    return dy / dx;
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Ray<PointType, LabelType>::pointInside() const {
    return Point<ResultNumber>(target());
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool Ray<PointType, LabelType>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    // pointInside() is an exact ray point past the source, so it is a valid
    // interior witness.
    return shape.interiorContains(pointInside<NumberType>());
}
 
// -----------------------------------------------------------------------------
// Rectangle
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber Rectangle<PointType, LabelType>::area() const {
    // width() and height() already report 0 for the empty rectangle, so its
    // inverted corners never reach the product.
    return static_cast<ResultNumber>(width()) * static_cast<ResultNumber>(height());
}
 
template <class PointType, class LabelType>
constexpr auto Rectangle<PointType, LabelType>::twiceArea() const {
    const auto rectangle_area = area();
    return rectangle_area + rectangle_area;
}
 
template <class PointType, class LabelType>
constexpr Segment<PointType> Rectangle<PointType, LabelType>::diameter() const {
    assert(!empty());
    return Segment<PointType>(min(), max());
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Rectangle<PointType, LabelType>::midpoint() const {
    assert(!empty());
    return Point<ResultNumber>(
        (detail::asNumber<ResultNumber>(min().x()) + detail::asNumber<ResultNumber>(max().x())) / static_cast<ResultNumber>(2),
        (detail::asNumber<ResultNumber>(min().y()) + detail::asNumber<ResultNumber>(max().y())) / static_cast<ResultNumber>(2));
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Rectangle<PointType, LabelType>::centroid() const {
    return midpoint<ResultNumber>();
}
 
template <class PointType, class LabelType>
constexpr Disk<PointType, NoLabel> Rectangle<PointType, LabelType>::circumcircle() const {
    assert(!empty());
    return Disk<PointType, NoLabel>(min(), bottomRight(), max()); //Choosen arbitrarly
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Rectangle<PointType, LabelType>::center() const {
    return midpoint<ResultNumber>();
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Rectangle<PointType, LabelType>::pointInside() const {
    return midpoint<ResultNumber>();
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool Rectangle<PointType, LabelType>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    const auto witness = pointInside<NumberType>();
    if (interiorContains(witness)) {
        return shape.interiorContains(witness);
    }
    // Integer truncation rounded the midpoint onto the boundary; scaling by 2
    // makes it exact without changing the containment relation.
    return (shape * 2).interiorContains((*this * 2).template pointInside<NumberType>());
}
 
// -----------------------------------------------------------------------------
// Triangle
 
template <class PointType, class LabelType>
constexpr typename Triangle<PointType, LabelType>::NumberType Triangle<PointType, LabelType>::twiceArea() const {
    return static_cast<NumberType>(detail::abs(orientationDeterminant(a(), b(), c())));
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber Triangle<PointType, LabelType>::area() const {
    const auto area2 = twiceArea();
    return detail::asNumber<ResultNumber>(area2) / static_cast<ResultNumber>(2);
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Triangle<PointType, LabelType>::centroid() const {
    const ResultNumber three = static_cast<ResultNumber>(3);
    return Point<ResultNumber>(
        (detail::asNumber<ResultNumber>(a().x()) +
         detail::asNumber<ResultNumber>(b().x()) +
         detail::asNumber<ResultNumber>(c().x())) / three,
        (detail::asNumber<ResultNumber>(a().y()) +
         detail::asNumber<ResultNumber>(b().y()) +
         detail::asNumber<ResultNumber>(c().y())) / three);
}
 
template <class PointType, class LabelType>
constexpr Disk<PointType, NoLabel> Triangle<PointType, LabelType>::circumcircle() const {
    return Disk<PointType, NoLabel>(a(), b(), c());
}
 
template <class PointType, class LabelType>
constexpr Segment<PointType> Triangle<PointType, LabelType>::diameter() const {
    const auto ab = a().template squaredDistance<NumberType>(b());
    const auto bc = b().template squaredDistance<NumberType>(c());
    const auto ca = c().template squaredDistance<NumberType>(a());
 
    if (ab < bc) {
        if (ca < bc) {
            return Segment<PointType>(b(), c());
        }
        return Segment<PointType>(c(), a());
    }
 
    if (ca < ab) {
        return Segment<PointType>(a(), b());
    }
    return Segment<PointType>(c(), a());
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Triangle<PointType, LabelType>::pointInside() const {
    Point<ResultNumber> p = points_[0] + points_[1];
    p = (p/2 + points_[2])/2;
    return p;
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool Triangle<PointType, LabelType>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    const auto witness = pointInside<NumberType>();
    if (interiorContains(witness)) {
        return shape.interiorContains(witness);
    }
    // pointInside() divides by 4; integer truncation rounded it onto the
    // boundary, so scaling by 4 makes it exact without changing containment.
    return (shape * 4).interiorContains((*this * 4).template pointInside<NumberType>());
}
 
template <class PointType, class LabelType>
constexpr bool Triangle<PointType, LabelType>::isRectangle() const {
    if (isDegenerate()) {
        return false;
    }
 
    const auto ab = b() - a();
    const auto ac = c() - a();
    const auto bc = c() - b();
 
    // dotSign promotes the coordinates before multiplying; the plain dot
    // product operator* does not, and overflows the coordinate type well
    // inside its range (a dot product of 2^32 wraps an int to zero, reporting
    // a right angle where there is none).
    //
    // Three difference vectors carry all three angles: the vectors this test
    // used to build for b and c are just negations of these, and negating one
    // operand of a dot product only flips its sign. ba·bc = -(ab·bc) and
    // ca·cb = ac·bc, and a test against zero cannot tell a sign flip apart.
    return dotSign(ab, ac) == 0 || dotSign(ab, bc) == 0 || dotSign(ac, bc) == 0;
}
 
template <class PointType, class LabelType>
constexpr bool Triangle<PointType, LabelType>::isObtuse() const {
    if (isDegenerate()) {
        return false;
    }
 
    const auto ab = b() - a();
    const auto ac = c() - a();
    const auto bc = c() - b();
 
    // Promoted like isRectangle above, and sharing its three difference
    // vectors. Here the sign does matter: ba·bc = -(ab·bc), so the obtuse
    // angle at b shows up as a *positive* ab·bc, while ca·cb = ac·bc keeps
    // its sign.
    return dotSign(ab, ac) < 0 || dotSign(ab, bc) > 0 || dotSign(ac, bc) < 0;
}
 
template <class PointType, class LabelType>
constexpr bool Triangle<PointType, LabelType>::isIsosceles() const {
    // Squared lengths need the promoted coordinate type for the same reason
    // the dot products above do: |ab|^2 overflows an int coordinate at around
    // 46341, and two sides that merely agree modulo the wrap would be reported
    // as equal. The caller never names this type, so the choice is ours to get
    // right rather than a precision the user asked for.
    using Squared = detail::promoted_number_t<NumberType>;
 
    const auto ab = a().template squaredDistance<Squared>(b());
    const auto bc = b().template squaredDistance<Squared>(c());
    const auto ca = c().template squaredDistance<Squared>(a());
    return ab == bc || bc == ca || ca == ab;
}
 
// -----------------------------------------------------------------------------
// Halfplane
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber Halfplane<PointType, LabelType>::slope() const {
    const auto dy = detail::asNumber<ResultNumber>(target().y()) - detail::asNumber<ResultNumber>(source().y());
    const auto dx = detail::asNumber<ResultNumber>(target().x()) - detail::asNumber<ResultNumber>(source().x());
    return dy / dx;
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Halfplane<PointType, LabelType>::pointInside() const {
    const ResultNumber x = static_cast<ResultNumber>(source().x());
    const ResultNumber y = static_cast<ResultNumber>(source().y());
    const ResultNumber one = static_cast<ResultNumber>(1);
 
    if (source().x() < target().x()) {
        return {x, y + one};
    }
    if (source().x() > target().x()) {
        return {x, y - one};
    }
    if (source().y() < target().y()) {
        return {x - one, y};
    }
    return {x + one, y - one};
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool Halfplane<PointType, LabelType>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    // pointInside() is an exact lattice point strictly inside the halfplane, so
    // it is a valid interior witness.
    return shape.interiorContains(pointInside<NumberType>());
}
 
 
// ---------------------------------------------------------------------------
// Convex
 
template <class PointType, class LabelType>
constexpr auto Convex<PointType, LabelType>::twiceArea() const {
    if (points_.size() < 3) {
        return NumberType(0);
    }
    NumberType sum = 0;
    for (std::size_t i = 0; i < points_.size(); ++i) {
        const auto& p1 = points_[i];
        const auto& p2 = points_[(i + 1) % points_.size()];
        sum += p1.x() * p2.y() - p2.x() * p1.y();
    }
    return pgl::detail::abs(sum);
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr auto Convex<PointType, LabelType>::area() const {
    ResultNumber result = static_cast<ResultNumber>(twiceArea());
    return result / ResultNumber(2);
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Convex<PointType, LabelType>::centroid() const {
    if (points_.empty()) {
        return Point<ResultNumber>();
    }
    if (points_.size() == 1) {
        return Point<ResultNumber>(points_[0]);
    }
    auto area_twice = twiceArea();
    if (points_.size() == 2 || area_twice == NumberType(0)) {
        const  Point<ResultNumber> p1 = static_cast<Point<ResultNumber>>(points_[0]);
        const  Point<ResultNumber> p2 = static_cast<Point<ResultNumber>>(points_[maxIndex()]);
        return (p1 + p2) / ResultNumber(2) + static_cast<Point<ResultNumber>>(translation_);
    }
    ResultNumber cx = 0;
    ResultNumber cy = 0;
    for (std::size_t i = 0; i < points_.size(); ++i) {
        const auto& p1 = points_[i];
        const auto& p2 = points_[(i + 1) % points_.size()];
        const auto cross = p1.x() * p2.y() - p2.x() * p1.y();
        cx += (p1.x() + p2.x()) * cross;
        cy += (p1.y() + p2.y()) * cross;
    }
    return Point<ResultNumber>(cx / (ResultNumber(3) * area_twice), cy / (ResultNumber(3) * area_twice))
           + static_cast<Point<ResultNumber>>(translation_);
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Convex<PointType, LabelType>::verticesCentroid() const {
    if (points_.empty()) {
        return Point<ResultNumber>();
    }
    ResultNumber cx = 0;
    ResultNumber cy = 0;
    for (const auto& vertex : points_) {
        cx += detail::asNumber<ResultNumber>(vertex.x());
        cy += detail::asNumber<ResultNumber>(vertex.y());
    }
    return Point<ResultNumber>(cx / static_cast<ResultNumber>(points_.size()), cy / static_cast<ResultNumber>(points_.size()))
           + static_cast<Point<ResultNumber>>(translation_);
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Convex<PointType, LabelType>::pointInside() const {
    if (points_.empty()) {
        return Point<ResultNumber>();
    }
    if (points_.size() == 1) {
        return Point<ResultNumber>(points_[0]) + static_cast<Point<ResultNumber>>(translation_);
    }
    if (points_.size() == 2) {
        const Point<ResultNumber> p1(points_[0]);
        const Point<ResultNumber> p2(points_[1]);
        return (p1 + p2) / ResultNumber(2) + static_cast<Point<ResultNumber>>(translation_);
    }
    Triangle<PointType> triangle(points_[0], points_[1], points_[2]);
    // Compute the inner triangle's interior point in ResultNumber too; using
    // the (now integer) default would truncate and can fall outside the polygon.
    return Point<ResultNumber>(triangle.template pointInside<ResultNumber>())
           + static_cast<Point<ResultNumber>>(translation_);
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool Convex<PointType, LabelType>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    const auto witness = pointInside<NumberType>();
    if (interiorContains(witness)) {
        return shape.interiorContains(witness);
    }
    // pointInside() divides by up to 4; integer truncation rounded it onto the
    // boundary, so scaling by 4 makes it exact without changing containment.
    return (shape * 4).interiorContains((*this * 4).template pointInside<NumberType>());
}
 
template <class PointType, class LabelType>
constexpr std::vector<std::pair<std::size_t, std::size_t>>
Convex<PointType, LabelType>::antipodalPairs() const {
    const std::size_t n = size();
    std::vector<std::pair<std::size_t, std::size_t>> pairs;
    if (n < 2) {
        return pairs;
    }
    if (n == 2) {
        pairs.emplace_back(0, 1);
        return pairs;
    }
 
    const auto next = [n](std::size_t v) { return v + 1 == n ? std::size_t{0} : v + 1; };
 
    // Cross product of edge(i) = p[i+1]-p[i] with edge(j), in exact promoted
    // arithmetic. As the caliper direction rotates CCW, the i-side caliper
    // reaches the next edge normal before the j-side one exactly when this is
    // negative; a zero means the two edges are parallel (a shared supporting
    // direction touching both edges at once).
    using Coord = detail::promoted_number_t<NumberType>;
    const auto edgeCross = [this, n](std::size_t i, std::size_t j) {
        const auto pi = (*this)[i];
        const auto pi1 = (*this)[i + 1 == n ? std::size_t{0} : i + 1];
        const auto pj = (*this)[j];
        const auto pj1 = (*this)[j + 1 == n ? std::size_t{0} : j + 1];
        const Coord eix = detail::asNumber<Coord>(pi1.x()) - detail::asNumber<Coord>(pi.x());
        const Coord eiy = detail::asNumber<Coord>(pi1.y()) - detail::asNumber<Coord>(pi.y());
        const Coord ejx = detail::asNumber<Coord>(pj1.x()) - detail::asNumber<Coord>(pj.x());
        const Coord ejy = detail::asNumber<Coord>(pj1.y()) - detail::asNumber<Coord>(pj.y());
        return eix * ejy - eiy * ejx;
    };
 
    // Start the calipers at the lexicographically smallest vertex (index 0, the
    // min-x support) and the max-x support; these admit vertical parallel
    // supporting lines, so they form an antipodal pair.
    const std::size_t i0 = maxIndex();
    const std::size_t j0 = 0;
    std::size_t i = i0;
    std::size_t j = j0;
 
    const auto isStart = [i0](std::size_t a, std::size_t b) {
        return (a == i0 && b == j0) || (a == j0 && b == i0);
    };
    const auto emit = [&pairs, &isStart](std::size_t a, std::size_t b) {
        // The starting pair is emitted once up front; skip any later reappearance
        // (it would otherwise close the sweep with a duplicate).
        if (isStart(a, b)) {
            return;
        }
        pairs.emplace_back(std::min(a, b), std::max(a, b));
    };
 
    pairs.emplace_back(std::min(i0, j0), std::max(i0, j0));
 
    // Sweep the caliper direction through half a turn. Advancing whichever side
    // reaches its next edge normal first walks every antipodal pair once.
    while (true) {
        const auto c = edgeCross(i, j);
        if (c < 0) {
            i = next(i);
        } else if (c > 0) {
            j = next(j);
        } else {
            // Parallel edges: vertices i, i+1 are co-extreme with j, j+1, so the
            // three remaining cross pairs of this edge/edge contact are antipodal.
            const std::size_t i2 = next(i);
            const std::size_t j2 = next(j);
            emit(i, j2);
            emit(i2, j);
            i = i2;
            j = j2;
        }
        if (i == j0 && j == i0) {
            break;  // returned to the swapped start: sweep complete.
        }
        emit(i, j);
    }
 
    return pairs;
}
 
template <class PointType, class LabelType>
constexpr Segment<PointType> Convex<PointType, LabelType>::diameter() const {
    const std::size_t n = size();
    if (n == 0) {
        return Segment<PointType>();
    }
    if (n == 1) {
        return Segment<PointType>((*this)[0], (*this)[0]);
    }
 
    // The farthest vertex pair is antipodal, so the diameter is the longest of
    // the O(n) antipodal segments. Compare squared lengths for exactness.
    const auto pairs = antipodalPairs();
    Segment<PointType> best((*this)[pairs.front().first], (*this)[pairs.front().second]);
    auto bestSquared = (*this)[pairs.front().first].template squaredDistance<NumberType>(
        (*this)[pairs.front().second]);
    for (const auto& [i, j] : pairs) {
        const auto pi = (*this)[i];
        const auto pj = (*this)[j];
        const auto squared = pi.template squaredDistance<NumberType>(pj);
        if (bestSquared < squared) {
            bestSquared = squared;
            best = Segment<PointType>(pi, pj);
        }
    }
    return best;
}
 
template <class PointType, class LabelType>
constexpr HalfplaneIntersection<PointType>
Convex<PointType, LabelType>::smallestEnclosingRectangle() const {
    const std::size_t n = size();
    if (n < 3) {
        // No area to enclose: the polygon (empty, a point, or a segment) is its
        // own smallest enclosing rectangle.
        return asHalfplaneIntersection();
    }
 
    // Every quantity below is a difference between vertices, so the promoted
    // coordinate type only has to cover the extent of the polygon.
    using Coord = detail::promoted_number_t<NumberType>;
    // Comparing two candidate areas is one degree past that, so it runs in a
    // type that grows to hold its values: their products are integers for
    // integral coordinates and fractions for rational ones. Floating-point
    // coordinates keep computing in the promoted floating-point type.
    using Wide = std::conditional_t<
        std::floating_point<NumberType>, Coord,
        std::conditional_t<RationalConcept<NumberType>, ERational, BigInt>>;
    struct Vec {
        Coord x, y;
    };
 
    const auto next = [n](std::size_t v) { return v + 1 == n ? std::size_t{0} : v + 1; };
    const auto difference = [this](std::size_t a, std::size_t b) {
        const auto pa = (*this)[a];
        const auto pb = (*this)[b];
        return Vec{detail::asNumber<Coord>(pa.x()) - detail::asNumber<Coord>(pb.x()),
                   detail::asNumber<Coord>(pa.y()) - detail::asNumber<Coord>(pb.y())};
    };
    const auto edgeVector = [&](std::size_t i) { return difference(next(i), i); };
    const auto dot = [](const Vec& u, const Vec& w) { return u.x * w.x + u.y * w.y; };
    const auto cross = [](const Vec& u, const Vec& w) { return u.x * w.y - u.y * w.x; };
 
    // The rectangle flush with edge i rests on three further supports: the
    // extreme vertex along the edge direction, the farthest vertex from the edge
    // line, and the extreme vertex against the edge direction. Around a CCW
    // boundary they appear in exactly that order after i, and as the edge
    // direction turns they only move forward, which is what makes the sweep
    // linear. Each objective is cyclically unimodal along the boundary, so
    // walking forward while the next vertex is strictly better lands on it.
    std::size_t right = 0;
    std::size_t top = 0;
    std::size_t left = 0;
    const auto advance = [&](std::size_t& support, auto better) {
        while (better(difference(next(support), support))) {
            support = next(support);
        }
    };
    const auto advanceRight = [&](const Vec& u) {
        advance(right, [&](const Vec& step) { return dot(u, step) > 0; });
    };
    const auto advanceTop = [&](const Vec& u) {
        advance(top, [&](const Vec& step) { return cross(u, step) > 0; });
    };
    const auto advanceLeft = [&](const Vec& u) {
        advance(left, [&](const Vec& step) { return dot(u, step) < 0; });
    };
 
    // Seed the three supports for edge 0, each starting where the previous one
    // stopped, since that is where its own ascent begins.
    {
        const Vec u = edgeVector(0);
        advanceRight(u);
        top = right;
        advanceTop(u);
        left = top;
        advanceLeft(u);
    }
 
    // A minimum-area enclosing rectangle has a side flush with a polygon edge,
    // so the best of these n candidates is the answer. Measured in the frame of
    // edge i, the rectangle spans [low, high] along u and [0, height] along u
    // turned 90 degrees, all three scaled by |u|; its area is therefore
    // (high - low) * height / (u * u), and two of those fractions are compared
    // by cross-multiplication.
    std::size_t bestEdge = 0, bestRight = 0, bestTop = 0, bestLeft = 0;
    Wide bestWidth(0), bestHeight(0), bestSquaredLength(1);
    for (std::size_t i = 0; i < n; ++i) {
        const Vec u = edgeVector(i);
        advanceRight(u);
        advanceTop(u);
        advanceLeft(u);
 
        const Coord low = dot(u, difference(left, i));
        const Coord high = dot(u, difference(right, i));
        const Coord width = high - low;
        const Coord height = cross(u, difference(top, i));
        const Coord squaredLength = dot(u, u);
 
        const Wide wideWidth = detail::asNumber<Wide>(width);
        const Wide wideHeight = detail::asNumber<Wide>(height);
        const Wide wideSquaredLength = detail::asNumber<Wide>(squaredLength);
        if (i == 0 || wideWidth * wideHeight * bestSquaredLength <
                          bestWidth * bestHeight * wideSquaredLength) {
            bestEdge = i;
            bestRight = right;
            bestTop = top;
            bestLeft = left;
            bestWidth = wideWidth;
            bestHeight = wideHeight;
            bestSquaredLength = wideSquaredLength;
        }
    }
 
    // The four supporting lines, each through the support it touches and along
    // the flush edge or that edge turned 90 degrees, oriented so that the
    // rectangle lies to the left. Nothing here divides, which is why the region
    // is exact in NumberType while its corners are not.
    const PointType base = (*this)[bestEdge];
    const PointType tip = (*this)[next(bestEdge)];
    const PointType along = tip - base;
    const PointType across = along.rotated90();
    const PointType topSupport = (*this)[bestTop];
    const PointType rightSupport = (*this)[bestRight];
    const PointType leftSupport = (*this)[bestLeft];
    return HalfplaneIntersection<PointType>({
        Halfplane<PointType>(base, tip),
        Halfplane<PointType>(rightSupport, rightSupport + across),
        Halfplane<PointType>(topSupport, topSupport - along),
        Halfplane<PointType>(leftSupport, leftSupport - across),
    });
}
 
namespace detail {

template <class PointType, class LabelType>
constexpr std::pair<std::size_t, std::size_t>
minimumWidthSupport(const Convex<PointType, LabelType>& convex) {
    using NumberType = typename PointType::NumberType;
    const std::size_t n = convex.size();
 
    // Every quantity below is a difference between vertices, so the promoted
    // coordinate type only has to cover the extent of the polygon.
    using Coord = promoted_number_t<NumberType>;
    // Comparing two candidate widths squares the supporting distance and
    // cross-multiplies by the squared edge lengths, which is three degrees past
    // that, so it runs in a type that grows to hold its values: their products
    // are integers for integral coordinates and fractions for rational ones.
    // Floating-point coordinates keep computing in the promoted floating-point
    // type.
    using Wide = std::conditional_t<
        std::floating_point<NumberType>, Coord,
        std::conditional_t<RationalConcept<NumberType>, ERational, BigInt>>;
    struct Vec {
        Coord x, y;
    };
 
    const auto next = [n](std::size_t v) { return v + 1 == n ? std::size_t{0} : v + 1; };
    const auto difference = [&convex](std::size_t a, std::size_t b) {
        const auto pa = convex[a];
        const auto pb = convex[b];
        return Vec{asNumber<Coord>(pa.x()) - asNumber<Coord>(pb.x()),
                   asNumber<Coord>(pa.y()) - asNumber<Coord>(pb.y())};
    };
    const auto edgeVector = [&](std::size_t i) { return difference(next(i), i); };
    const auto dot = [](const Vec& u, const Vec& w) { return u.x * w.x + u.y * w.y; };
    const auto cross = [](const Vec& u, const Vec& w) { return u.x * w.y - u.y * w.x; };
 
    std::size_t top = 0;
    const auto advance = [&](auto better) {
        while (better(difference(next(top), top))) {
            top = next(top);
        }
    };
    const auto advanceTop = [&](const Vec& u) {
        advance([&](const Vec& step) { return cross(u, step) > Coord(0); });
    };
 
    // Seed the support for edge 0. The supporting distance stays flat along the
    // edge itself before it starts to rise, so the walk has to pass the extreme
    // vertex along the edge direction first -- that is where its own ascent
    // begins -- or it stalls on the flat step and never leaves vertex 0.
    {
        const Vec u = edgeVector(0);
        advance([&](const Vec& step) { return dot(u, step) > Coord(0); });
        advanceTop(u);
    }
 
    // A minimum-width slab has one supporting line flush with a polygon edge,
    // so the best of these n candidates is the answer. Measured in the frame of
    // edge i, the slab spans [0, height] along u turned 90 degrees, scaled by
    // |u|; its width is therefore height / (u * u) raised to the half, and two
    // of those are compared squared, by cross-multiplication.
    std::size_t bestEdge = 0, bestTop = 0;
    Wide bestSquaredHeight(0), bestSquaredLength(1);
    bool found = false;
    for (std::size_t i = 0; i < n; ++i) {
        const Vec u = edgeVector(i);
        advanceTop(u);
 
        const Coord squaredLength = dot(u, u);
        if (squaredLength == Coord(0)) {
            // A repeated vertex leaves no direction to measure across, and its
            // slab would divide by zero. Canonical hull vertices have no such
            // edge; one reaching here came from a trusted construction.
            continue;
        }
        const Coord height = cross(u, difference(top, i));
 
        const Wide wideHeight = asNumber<Wide>(height);
        const Wide wideSquaredHeight = wideHeight * wideHeight;
        const Wide wideSquaredLength = asNumber<Wide>(squaredLength);
        if (!found || wideSquaredHeight * bestSquaredLength <
                          bestSquaredHeight * wideSquaredLength) {
            found = true;
            bestEdge = i;
            bestTop = top;
            bestSquaredHeight = wideSquaredHeight;
            bestSquaredLength = wideSquaredLength;
        }
    }
    return {bestEdge, bestTop};
}

template <class Number, class PointType, class LabelType>
constexpr std::pair<Number, Number>
minimumWidthFraction(const Convex<PointType, LabelType>& convex, std::size_t edge, std::size_t top) {
    const std::size_t n = convex.size();
    const auto base = convex[edge];
    const auto tip = convex[edge + 1 == n ? std::size_t{0} : edge + 1];
    const auto support = convex[top];
 
    const Number ux = asNumber<Number>(tip.x()) - asNumber<Number>(base.x());
    const Number uy = asNumber<Number>(tip.y()) - asNumber<Number>(base.y());
    const Number vx = asNumber<Number>(support.x()) - asNumber<Number>(base.x());
    const Number vy = asNumber<Number>(support.y()) - asNumber<Number>(base.y());
 
    const Number squaredLength = ux * ux + uy * uy;
    if (squaredLength == Number(0)) {
        // Every edge of the polygon is a repeated vertex, which only a trusted
        // construction produces. It covers one point, whose width is zero;
        // reporting the fraction as 0/1 keeps that the answer instead of a
        // division by zero.
        return {Number(0), Number(1)};
    }
    return {ux * vy - uy * vx, squaredLength};
}
 
}  // namespace detail
 
template <class PointType, class LabelType>
constexpr HalfplaneIntersection<PointType>
Convex<PointType, LabelType>::smallestEnclosingSlab() const {
    const std::size_t n = size();
    if (n < 3) {
        // No width to minimize: the polygon (empty, a point, or a segment) is
        // its own narrowest slab, of width zero.
        return asHalfplaneIntersection();
    }
 
    const auto [bestEdge, bestTop] = detail::minimumWidthSupport(*this);
 
    // The two parallel supporting lines, one along the flush edge and one
    // through the opposite support along that same edge vector, oriented so
    // that the slab lies to the left of both. Nothing here divides, which is
    // why the region is exact in NumberType while its width is not.
    const PointType base = (*this)[bestEdge];
    const PointType tip = (*this)[bestEdge + 1 == n ? std::size_t{0} : bestEdge + 1];
    const PointType along = tip - base;
    const PointType topSupport = (*this)[bestTop];
    return HalfplaneIntersection<PointType>({
        Halfplane<PointType>(base, tip),
        Halfplane<PointType>(topSupport, topSupport - along),
    });
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr ResultNumber Convex<PointType, LabelType>::squaredMinimumWidth() const {
    if (size() < 3) {
        return ResultNumber{};
    }
 
    const auto [bestEdge, bestTop] = detail::minimumWidthSupport(*this);
    const auto [height, squaredLength] =
        detail::minimumWidthFraction<ResultNumber>(*this, bestEdge, bestTop);
    return height * height / squaredLength;
}
 
template <class PointType, class LabelType>
template <class ApproximateNumber>
ApproximateNumber Convex<PointType, LabelType>::minimumWidth() const {
    if (size() < 3) {
        return ApproximateNumber{};
    }
 
    const auto [bestEdge, bestTop] = detail::minimumWidthSupport(*this);
    const auto [height, squaredLength] =
        detail::minimumWidthFraction<ApproximateNumber>(*this, bestEdge, bestTop);
    // One square root, of the edge length rather than of the squared width, so
    // the division's operands are the exactly computed ones.
    return height / std::sqrt(squaredLength);
}
 
// -----------------------------------------------------------------------------
// Polygon
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Polygon<PointType, LabelType>::pointInside() const {
    const std::size_t n = size();
    if (n < 3) {
        // Degenerate polygon (UB per the library contract); fall back to a
        // representative point rather than reading past the end.
        return verticesCentroid<ResultNumber>();
    }
 
    // In canonical form the lexicographically smallest vertex is stored first,
    // so p0 is a convex vertex; a and b are its boundary neighbours, and the
    // triangle (a, p0, b) is counterclockwise and lies locally inside.
    const PointType p0 = (*this)[0];
    const Triangle<PointType> t(p0,get(1),get(-1));
 
    if (t.isDegenerate()) {
        // Degenerate input polygon: all its vertices have the same x coordinate
        return p0;
    }
 
    // Leftmost (lexicographically smallest) vertex inside the closed triangle
    // (a, p0, b), skipping the three corners at indices 0, 1 and n-1. Because
    // p0 is the global lex-min vertex, the diagonal from p0 to this vertex
    // stays inside the polygon.
    std::optional<PointType> q;
    for (std::size_t i = 2; i + 1 < n; ++i) {
        const PointType v = (*this)[i];
        if (t.interiorContains(v)) {
            if (!q.has_value() || v < *q) {
                q = v;
            }
        }
    }
 
    if (q.has_value()) {
        // p0 q is a diagonal; its midpoint is strictly interior (divides by 2).
        return (Point<ResultNumber>(p0) + Point<ResultNumber>(*q)) / ResultNumber(2);
    }
 
    // The triangle is an ear: return its interior point (divides by 4).
    return t.template pointInside<ResultNumber>();
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool Polygon<PointType, LabelType>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    const auto witness = pointInside<NumberType>();
    if (interiorContains(witness)) {
        return shape.interiorContains(witness);
    }
    // pointInside() divides by up to 4; integer truncation rounded it onto the
    // boundary, so scaling by 4 makes it exact without changing containment.
    return (shape * 4).interiorContains((*this * 4).template pointInside<NumberType>());
}
 
// -----------------------------------------------------------------------------
// PolygonWithHoles
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> PolygonWithHoles<PointType, LabelType>::verticesCentroid() const {
    const std::size_t n = vertexCount();
    if (n == 0) {
        return Point<ResultNumber>();
    }
    ResultNumber cx{};
    ResultNumber cy{};
    const auto accumulate = [&](const PolygonType& ring) {
        for (const auto& vertex : ring) {
            cx += detail::asNumber<ResultNumber>(vertex.x());
            cy += detail::asNumber<ResultNumber>(vertex.y());
        }
    };
    accumulate(outer_);
    for (const auto& hole : holes_) {
        accumulate(hole);
    }
    return Point<ResultNumber>(cx / detail::asNumber<ResultNumber>(n), cy / detail::asNumber<ResultNumber>(n));
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> PolygonWithHoles<PointType, LabelType>::centroid() const {
    if (empty()) {
        return Point<ResultNumber>();
    }
    const NumberType netTwiceArea = twiceArea();
    if (netTwiceArea == NumberType(0)) {
        // No area to weight by (a collapsed outer ring, or holes cancelling it
        // out); fall back to the vertex centroid as Polygon does.
        return verticesCentroid<ResultNumber>();
    }
 
    // Each ring contributes its own area-weighted centroid, holes negatively.
    // Weighting by twice the area keeps the single division to the very end.
    ResultNumber cx{};
    ResultNumber cy{};
    const auto accumulate = [&](const PolygonType& ring, int sign) {
        const auto ringTwiceArea = static_cast<ResultNumber>(ring.twiceArea());
        const auto ringCentroid = ring.template centroid<ResultNumber>();
        const ResultNumber weight = ringTwiceArea * detail::asNumber<ResultNumber>(sign);
        cx += ringCentroid.x() * weight;
        cy += ringCentroid.y() * weight;
    };
    accumulate(outer_, 1);
    for (const auto& hole : holes_) {
        accumulate(hole, -1);
    }
    const auto denominator = static_cast<ResultNumber>(netTwiceArea);
    return Point<ResultNumber>(cx / denominator, cy / denominator);
}
 
// -----------------------------------------------------------------------------
// PolygonSet
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> PolygonSet<PointType, LabelType>::verticesCentroid() const {
    const std::size_t n = vertexCount();
    if (n == 0) {
        return Point<ResultNumber>();
    }
    ResultNumber cx{};
    ResultNumber cy{};
    for (const auto& component : components_) {
        const auto componentCentroid = component.template verticesCentroid<ResultNumber>();
        const auto weight = static_cast<ResultNumber>(component.vertexCount());
        cx += componentCentroid.x() * weight;
        cy += componentCentroid.y() * weight;
    }
    return Point<ResultNumber>(cx / detail::asNumber<ResultNumber>(n), cy / detail::asNumber<ResultNumber>(n));
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> PolygonSet<PointType, LabelType>::centroid() const {
    if (empty()) {
        return Point<ResultNumber>();
    }
    const NumberType netTwiceArea = twiceArea();
    if (netTwiceArea == NumberType(0)) {
        // No area to weight by; fall back to the vertex centroid as Polygon and
        // PolygonWithHoles do.
        return verticesCentroid<ResultNumber>();
    }
 
    // The components have disjoint interiors, so the set's centroid is their
    // area-weighted mean. Weighting by twice the area keeps the single division
    // to the very end.
    ResultNumber cx{};
    ResultNumber cy{};
    for (const auto& component : components_) {
        const auto weight = static_cast<ResultNumber>(component.twiceArea());
        const auto componentCentroid = component.template centroid<ResultNumber>();
        cx += componentCentroid.x() * weight;
        cy += componentCentroid.y() * weight;
    }
    const auto denominator = static_cast<ResultNumber>(netTwiceArea);
    return Point<ResultNumber>(cx / denominator, cy / denominator);
}
 
// -----------------------------------------------------------------------------
// MonotoneChain
 
template <class PointType, class LabelType, class Storage>
template <class ApproximateNumber>
ApproximateNumber MonotoneChain<PointType, LabelType, Storage>::length() const {
    ApproximateNumber total{};
    for (std::size_t i = 1; i < points_.size(); ++i) {
        total += points_[i - 1].template distance<ApproximateNumber>(points_[i]);
    }
    return total;
}
 
template <class PointType, class LabelType, class Storage>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::lengthL1() const {
    decltype(std::declval<PointType>().distanceL1(std::declval<PointType>())) total{};
    for (std::size_t i = 1; i < points_.size(); ++i) {
        total += points_[i - 1].distanceL1(points_[i]);
    }
    return total;
}
 
template <class PointType, class LabelType, class Storage>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::lengthLInf() const {
    decltype(std::declval<PointType>().distanceLInf(std::declval<PointType>())) total{};
    for (std::size_t i = 1; i < points_.size(); ++i) {
        total += points_[i - 1].distanceLInf(points_[i]);
    }
    return total;
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber>
constexpr Point<ResultNumber> MonotoneChain<PointType, LabelType, Storage>::pointInside() const {
    if (size() < 2) {
        // Degenerate case
        return size() == 0 ? PointType() : (*this)[0];
    }
 
    return Segment<PointType>((*this)[0], (*this)[1]).template pointInside<ResultNumber>();
}
 
// -----------------------------------------------------------------------------
// Polyline
 
template <class PointType, class LabelType>
template <class ApproximateNumber>
ApproximateNumber Polyline<PointType, LabelType>::length() const {
    ApproximateNumber total{};
    for (std::size_t i = 1; i < points_.size(); ++i) {
        total += points_[i - 1].template distance<ApproximateNumber>(points_[i]);
    }
    return total;
}
 
template <class PointType, class LabelType>
constexpr auto Polyline<PointType, LabelType>::lengthL1() const {
    decltype(std::declval<PointType>().distanceL1(std::declval<PointType>())) total{};
    for (std::size_t i = 1; i < points_.size(); ++i) {
        total += points_[i - 1].distanceL1(points_[i]);
    }
    return total;
}
 
template <class PointType, class LabelType>
constexpr auto Polyline<PointType, LabelType>::lengthLInf() const {
    decltype(std::declval<PointType>().distanceLInf(std::declval<PointType>())) total{};
    for (std::size_t i = 1; i < points_.size(); ++i) {
        total += points_[i - 1].distanceLInf(points_[i]);
    }
    return total;
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> Polyline<PointType, LabelType>::pointInside() const {
    if (size() < 2) {
        // Degenerate case
        return size() == 0 ? PointType() : (*this)[0];
    }
 
    return Segment<PointType>((*this)[0], (*this)[1]).template pointInside<ResultNumber>();
}
 
template <class PointType, class LabelType, class Storage>
template <class OtherShape>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    const auto witness = pointInside<NumberType>();
    if (interiorContains(witness)) {
        return shape.interiorContains(witness);
    }
    // pointInside() is the first edge's midpoint (divides by 2); integer
    // truncation rounded it onto a vertex, so scaling by 2 makes it exact
    // without changing the containment relation.
    return (shape * 2).interiorContains((*this * 2).template pointInside<NumberType>());
}
 
 
// -----------------------------------------------------------------------------
// HalfplaneIntersection
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr auto HalfplaneIntersection<PointType, LabelType>::twiceArea() const {
    if (empty_) {
        return ResultNumber{};
    }
    if (!isBounded()) {
        throw std::logic_error("HalfplaneIntersection::twiceArea requires a bounded region");
    }
    return static_cast<ResultNumber>(asConvex<ResultNumber>().twiceArea());
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr auto HalfplaneIntersection<PointType, LabelType>::area() const {
    if (empty_) {
        return ResultNumber{};
    }
    if (!isBounded()) {
        throw std::logic_error("HalfplaneIntersection::area requires a bounded region");
    }
    return asConvex<ResultNumber>().template area<ResultNumber>();
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> HalfplaneIntersection<PointType, LabelType>::centroid() const {
    if (empty_) {
        return Point<ResultNumber>();
    }
    if (!isBounded()) {
        throw std::logic_error("HalfplaneIntersection::centroid requires a bounded region");
    }
    const auto c = asConvex<ResultNumber>().template centroid<ResultNumber>();
    return Point<ResultNumber>(c.x(), c.y());
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Point<ResultNumber> HalfplaneIntersection<PointType, LabelType>::pointInside() const {
    if (empty_ || halfplanes_.empty()) {
        return Point<ResultNumber>();  // the whole plane's representative is the origin
    }
    if (halfplanes_.size() == 1) {
        return halfplanes_[0].template pointInside<ResultNumber>();
    }
    if (isBounded()) {
        return asConvex<ResultNumber>().template pointInside<ResultNumber>();
    }
    // Unbounded: clip a box around a finite anchor point of the region, then
    // take that bounded piece's interior representative. A vertex, or (for a
    // slab) any point of the first boundary line, anchors the box on the
    // region.
    NumberType anchorX{};
    NumberType anchorY{};
    bool anchored = false;
    for (std::size_t i = 0; i < size(); ++i) {
        if (vertexExists(i)) {
            const auto v = vertex<double>(i);
            anchorX = static_cast<NumberType>(v.x() < 0 ? -std::floor(-v.x()) : std::floor(v.x()));
            anchorY = static_cast<NumberType>(v.y() < 0 ? -std::floor(-v.y()) : std::floor(v.y()));
            anchored = true;
            break;
        }
    }
    if (!anchored) {
        anchorX = static_cast<NumberType>(halfplanes_[0].source().x());
        anchorY = static_cast<NumberType>(halfplanes_[0].source().y());
    }
    const NumberType margin(4);
    const PointType lo(anchorX - margin, anchorY - margin);
    const PointType hi(anchorX + margin, anchorY + margin);
    const PointType lohi(lo.x(), hi.y());
    const PointType hilo(hi.x(), lo.y());
    HalfplaneIntersection clipped(*this);
    clipped.insert(HalfplaneType(lo, hilo));
    clipped.insert(HalfplaneType(hilo, hi));
    clipped.insert(HalfplaneType(hi, lohi));
    clipped.insert(HalfplaneType(lohi, lo));
    if (clipped.empty()) {
        // The anchor rounding missed the region; fall back to a boundary point.
        return Point<ResultNumber>(detail::asNumber<ResultNumber>(halfplanes_[0].source().x()),
                                   detail::asNumber<ResultNumber>(halfplanes_[0].source().y()));
    }
    return clipped.template asConvex<ResultNumber>().template pointInside<ResultNumber>();
}
 
template <class PointType, class LabelType>
template <class OtherShape>
constexpr bool HalfplaneIntersection<PointType, LabelType>::pointInsideInteriorContainedIn(
    const OtherShape& shape) const {
    // Exact witness: rational unless the coordinates are already floating point
    // (region_exact_number_t is not yet visible at this include position).
    using Exact = std::conditional_t<std::is_floating_point_v<NumberType>, double, Rational<BigInt>>;
    return shape.interiorContains(pointInside<Exact>());
}
 
}  // namespace pgl