halfplaneintersection.hpp

Go to the documentation of this file.

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

 
#include <algorithm>
#include <cassert>
#include <compare>
#include <concepts>
#include <cstddef>
#include <functional>
#include <optional>
#include <ostream>
#include <ranges>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <variant>
#include <vector>
 
 
namespace pgl::detail {

template <HalfplaneConcept H>
constexpr int pseudoAngleHalf(const H& h) {
    const auto& s = h.source();
    const auto& t = h.target();
    if (t.y() != s.y()) {
        return t.y() > s.y() ? 0 : 1;
    }
    return t.x() > s.x() ? 0 : 1;
}

template <HalfplaneConcept HA, HalfplaneConcept HB>
constexpr std::partial_ordering directionCross(const HA& a, const HB& b) {
    return crossSign(a.source(), a.target(), b.source(), b.target());
}

template <HalfplaneConcept HA, HalfplaneConcept HB>
constexpr bool directionLess(const HA& a, const HB& b) {
    const int halfA = pseudoAngleHalf(a);
    const int halfB = pseudoAngleHalf(b);
    if (halfA != halfB) {
        return halfA < halfB;
    }
    return directionCross(a, b) > 0;
}

template <HalfplaneConcept HA, HalfplaneConcept HB>
constexpr std::partial_ordering directionDot(const HA& a, const HB& b) {
    return dotSign(a.source(), a.target(), b.source(), b.target());
}

template <HalfplaneConcept HF, HalfplaneConcept HT, HalfplaneConcept HX>
constexpr bool arcContainsDirection(const HF& from, const HT& to, const HX& x) {
    const auto crossFrom = directionCross(from, x);
    if (crossFrom < 0) {
        return false;
    }
    if (crossFrom == 0 && directionDot(from, x) <= 0) {
        return false;  // opposite to the arc start, not on it
    }
    const auto crossTo = directionCross(x, to);
    if (crossTo < 0) {
        return false;
    }
    if (crossTo == 0 && directionDot(to, x) <= 0) {
        return false;  // opposite to the arc end
    }
    return true;
}

template <HalfplaneConcept HA, HalfplaneConcept HB>
constexpr bool directionEqual(const HA& a, const HB& b) {
    return pseudoAngleHalf(a) == pseudoAngleHalf(b) && directionCross(a, b) == 0;
}

template <HalfplaneConcept H1, HalfplaneConcept H2, HalfplaneConcept H3>
constexpr auto boundaryLinesDeterminant(const H1& h1, const H2& h2, const H3& h3) {
    using Common = std::common_type_t<typename H1::NumberType, typename H2::NumberType, typename H3::NumberType>;
    using C = promoted_number_t<promoted_number_t<Common>>;
    const auto row = [](const auto& h) {
        const C sx = static_cast<C>(h.source().x());
        const C sy = static_cast<C>(h.source().y());
        const C tx = static_cast<C>(h.target().x());
        const C ty = static_cast<C>(h.target().y());
        struct Row { C ax, ay, b; };
        return Row{ty - sy, sx - tx, sx * ty - sy * tx};
    };
    const auto r1 = row(h1);
    const auto r2 = row(h2);
    const auto r3 = row(h3);
    return r1.ax * (r2.ay * r3.b - r3.ay * r2.b)
         - r1.ay * (r2.ax * r3.b - r3.ax * r2.b)
         + r1.b * (r2.ax * r3.ay - r3.ax * r2.ay);
}

template <HalfplaneConcept H1, HalfplaneConcept H2, HalfplaneConcept H3>
constexpr int vertexSide(const H1& h1, const H2& h2, const H3& h3) {
    const auto det = boundaryLinesDeterminant(h1, h2, h3);
    const auto zero = decltype(det){};
    const int detSign = det > zero ? 1 : det < zero ? -1 : 0;
    return detSign * signOf(directionCross(h1, h2));
}
 
}  // namespace pgl::detail
 
 
namespace pgl {
 
template <class PointType = Point<>, class Label>
struct HalfplaneIntersection;
 
HalfplaneIntersection() -> HalfplaneIntersection<Point<>, NoLabel>;
 
template <HalfplaneConcept H>
HalfplaneIntersection(const H&) -> HalfplaneIntersection<typename H::PointType, NoLabel>;
 
template <std::ranges::input_range Range>
    requires detail::is_halfplane_v<std::ranges::range_value_t<Range>>
HalfplaneIntersection(Range&&) -> HalfplaneIntersection<typename std::remove_cvref_t<std::ranges::range_value_t<Range>>::PointType, NoLabel>;
 
template <std::ranges::input_range Range>
    requires detail::is_halfplane_v<std::ranges::range_value_t<Range>>
HalfplaneIntersection(Range&&, bool) -> HalfplaneIntersection<typename std::remove_cvref_t<std::ranges::range_value_t<Range>>::PointType, NoLabel>;
 
template <RectangleConcept R>
HalfplaneIntersection(const R&) -> HalfplaneIntersection<typename R::PointType, NoLabel>;
 
template <TriangleConcept T>
HalfplaneIntersection(const T&) -> HalfplaneIntersection<typename T::PointType, NoLabel>;
 
template <ConvexConcept C>
HalfplaneIntersection(const C&) -> HalfplaneIntersection<typename C::PointType, NoLabel>;
 
template <PointConcept P>
HalfplaneIntersection(const P&) -> HalfplaneIntersection<Point<typename P::NumberType, typename P::LabelType>, NoLabel>;
 
template <SegmentConcept S>
HalfplaneIntersection(const S&) -> HalfplaneIntersection<typename S::PointType, NoLabel>;
 
template <LineConcept L>
HalfplaneIntersection(const L&) -> HalfplaneIntersection<typename L::PointType, NoLabel>;

template <class PointType_, class TLabel>
struct HalfplaneIntersection {
    using PointType = PointType_;
    using NumberType = PointType::NumberType;
    using LabelType = TLabel;
    using HalfplaneType = Halfplane<PointType>;
 
    static_assert(detail::is_point_v<PointType>, "HalfplaneIntersection requires pgl::Point defining points");

    constexpr HalfplaneIntersection() = default;

    template <HalfplaneConcept OtherHalfplane>
    constexpr explicit HalfplaneIntersection(const OtherHalfplane& halfplane) {
        insert(halfplane);
    }

    template <std::ranges::input_range Range = std::initializer_list<HalfplaneType>>
        requires std::convertible_to<std::ranges::range_value_t<Range>, HalfplaneType> &&
                 (!detail::is_halfplane_intersection_v<Range>)
    constexpr explicit HalfplaneIntersection(Range&& halfplanes, bool trusted = false) {
        if (trusted) {
            for (const auto& h : halfplanes) {
                halfplanes_.push_back(h);
            }
        } else {
            for (const auto& h : halfplanes) {
                insert(h);
            }
        }
    }

    template <RectangleConcept OtherRectangle>
    constexpr explicit HalfplaneIntersection(const OtherRectangle& rectangle) {
        if (rectangle.empty()) {
            empty_ = true;
            return;
        }
        const PointType lo(rectangle.min());
        const PointType hi(rectangle.max());
        if (lo == hi) {
            // A single point: the four edges would all be degenerate.
            buildPointSlab(lo);
            canonicalizeSorted();
            return;
        }
        if (lo.x() == hi.x() || lo.y() == hi.y()) {
            // A zero-width or zero-height rectangle is the segment lo-hi; two of
            // the four edges would be degenerate, so clamp to the segment.
            buildSegmentSlab(lo, hi);
            canonicalizeSorted();
            return;
        }
        const PointType lohi(lo.x(), hi.y());
        const PointType hilo(hi.x(), lo.y());
        // Directions 0, 90, 180, 270 degrees: already in sorted order.
        halfplanes_.push_back(HalfplaneType(lo, hilo));
        halfplanes_.push_back(HalfplaneType(hilo, hi));
        halfplanes_.push_back(HalfplaneType(hi, lohi));
        halfplanes_.push_back(HalfplaneType(lohi, lo));
    }

    template <TriangleConcept OtherTriangle>
    constexpr explicit HalfplaneIntersection(const OtherTriangle& triangle) {
        if (const auto vertex = triangle.getIfPoint()) {
            // A single point: all three edges would be degenerate.
            buildPointSlab(PointType(*vertex));
            canonicalizeSorted();
            return;
        }
        if (const auto carrier = triangle.getIfSegment()) {
            // Collinear vertices: the edges span the carrier segment, and at
            // least one of them would be degenerate.
            buildSegmentSlab(PointType(carrier->min()), PointType(carrier->max()));
            canonicalizeSorted();
            return;
        }
        // Triangle vertices are CCW when non-degenerate, so each edge's
        // half-plane has the interior on its left.
        for (std::size_t i = 0; i < 3; ++i) {
            halfplanes_.push_back(HalfplaneType(PointType(triangle[i]), PointType(triangle[(i + 1) % 3])));
        }
        canonicalizeSorted();
    }

    template <ConvexConcept OtherConvex>
    constexpr explicit HalfplaneIntersection(const OtherConvex& convex) {
        const std::size_t n = convex.size();
        if (n == 0) {
            empty_ = true;
            return;
        }
        if (convex.isDegenerate()) {
            // A point or a segment: clamp both slab directions. The convex
            // polygon normalizes collinear input down to its two extremes.
            const PointType a(convex[0]);
            const PointType b(convex[n - 1]);
            if (a == b) {
                buildPointSlab(a);
            } else {
                buildSegmentSlab(a, b);
            }
            canonicalizeSorted();
            return;
        }
        for (std::size_t i = 0; i < n; ++i) {
            halfplanes_.push_back(HalfplaneType(PointType(convex[i]), PointType(convex[(i + 1) % n])));
        }
        canonicalizeSorted();
    }

    template <PointConcept OtherPoint>
    constexpr explicit HalfplaneIntersection(const OtherPoint& point) {
        buildPointSlab(PointType(point));
        canonicalizeSorted();
    }

    template <SegmentConcept OtherSegment>
    constexpr explicit HalfplaneIntersection(const OtherSegment& segment) {
        const PointType a(segment[0]);
        const PointType b(segment[1]);
        if (a == b) {
            buildPointSlab(a);
        } else {
            buildSegmentSlab(a, b);
        }
        canonicalizeSorted();
    }

    template <LineConcept OtherLine>
    constexpr explicit HalfplaneIntersection(const OtherLine& line) {
        const PointType a(line[0]);
        const PointType b(line[1]);
        assert(a != b);
        halfplanes_.push_back(HalfplaneType(a, b));
        halfplanes_.push_back(HalfplaneType(b, a));
        degenerate_ = true;
        canonicalizeSorted();
    }

    template <PointConcept OtherPointType, class OtherLabelType>
        requires(std::constructible_from<PointType, const OtherPointType&>)
    constexpr HalfplaneIntersection(const HalfplaneIntersection<OtherPointType, OtherLabelType>& other)
        : empty_(other.empty()), degenerate_(other.empty() ? false : other.isDegenerate()) {
        halfplanes_.reserve(other.size());
        for (const auto& h : other) {
            halfplanes_.push_back(HalfplaneType(h));
        }
        label_ = detail::copyLabel<LabelType>(other);
    }

    template <PointConcept OtherPointType, class OtherLabelType>
        requires(std::constructible_from<PointType, const OtherPointType&>)
    constexpr HalfplaneIntersection& operator=(const HalfplaneIntersection<OtherPointType, OtherLabelType>& other) {
        halfplanes_.clear();
        halfplanes_.reserve(other.size());
        for (const auto& h : other) {
            halfplanes_.push_back(HalfplaneType(h));
        }
        empty_ = other.empty();
        degenerate_ = empty_ ? false : other.isDegenerate();
        label_ = detail::copyLabel<LabelType>(other);
        resetCache();
        return *this;
    }

    template <class A = LabelType>
        requires(detail::has_label_v<A>)
    constexpr A& label() const {
        return label_;
    }

    template <HalfplaneConcept OtherHalfplane>
    constexpr bool insert(const OtherHalfplane& other) {
        if (empty_) {
            return false;
        }
        const HalfplaneType h(PointType(other.source()), PointType(other.target()));
        if (h.isUndefined()) {
            return false;
        }
        if (halfplanes_.empty()) {
            halfplanes_.push_back(h);
            resetCache();
            return true;
        }
        // sup a·p <= b over the region means the region is already inside h.
        const SupStatus supremum = supStatus(h);
        if (supremum == SupStatus::below || supremum == SupStatus::on) {
            return false;
        }
        // The infimum side: inf a·p > b means no point of the region satisfies
        // h; inf == b (attained) means the region collapses onto the minimum
        // face, which has empty interior.
        const SupStatus infimum = supStatus(h.opposite());
        if (infimum == SupStatus::below) {
            halfplanes_.clear();
            empty_ = true;
            degenerate_ = false;
            resetCache();
            return true;
        }
        if (infimum == SupStatus::on) {
            degenerate_ = true;
        }
        std::size_t pos = linearLowerBound(h);
        if (pos < halfplanes_.size() && detail::directionEqual(halfplanes_[pos], h)) {
            // Same direction: h is strictly more restrictive, otherwise the
            // supremum test above would have discarded it.
            halfplanes_[pos] = h;
        } else {
            if (pos > halfplanes_.size()) {
                pos = halfplanes_.size();
            }
            halfplanes_.insert(halfplanes_.begin() + static_cast<std::ptrdiff_t>(pos), h);
        }
        // Cascade: a stored half-plane g between angular neighbors pred and
        // succ is redundant exactly when the wedge pred ∩ succ fits inside g,
        // which (for an angular gap below pi) reduces to the wedge apex lying
        // weakly inside g. Deletions are paid for once, so the cascade is
        // amortized O(1) tests per insertion.
        while (halfplanes_.size() >= 3) {
            const std::size_t s1 = nextIndex(pos);
            const std::size_t s2 = nextIndex(s1);
            if (s2 == pos) {
                break;
            }
            if (detail::directionCross(halfplanes_[pos], halfplanes_[s2]) > 0 &&
                detail::vertexSide(halfplanes_[pos], halfplanes_[s2], halfplanes_[s1]) >= 0) {
                halfplanes_.erase(halfplanes_.begin() + static_cast<std::ptrdiff_t>(s1));
                if (s1 < pos) {
                    --pos;
                }
            } else {
                break;
            }
        }
        while (halfplanes_.size() >= 3) {
            const std::size_t p1 = prevIndex(pos);
            const std::size_t p2 = prevIndex(p1);
            if (p2 == pos) {
                break;
            }
            if (detail::directionCross(halfplanes_[p2], halfplanes_[pos]) > 0 &&
                detail::vertexSide(halfplanes_[p2], halfplanes_[pos], halfplanes_[p1]) >= 0) {
                halfplanes_.erase(halfplanes_.begin() + static_cast<std::ptrdiff_t>(p1));
                if (p1 < pos) {
                    --pos;
                }
            } else {
                break;
            }
        }
        resetCache();
        return true;
    }

    constexpr std::size_t size() const {
        return halfplanes_.size();
    }

    constexpr const HalfplaneType& operator[](std::size_t index) const {
        assert(index < size());
        return halfplanes_[index];
    }

    constexpr const HalfplaneType& get(std::ptrdiff_t index) const {
        const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(size());
        return (*this)[static_cast<std::size_t>(((index % n) + n) % n)];
    }

    constexpr std::ptrdiff_t index(const HalfplaneType& halfplane) const {
        const std::size_t pos = linearLowerBound(halfplane);
        if (pos < size() && halfplanes_[pos] == halfplane) {
            return static_cast<std::ptrdiff_t>(pos);
        }
        return -1;
    }

    constexpr std::vector<HalfplaneType> halfplanes() const {
        return halfplanes_;
    }

    constexpr auto begin() const { return halfplanes_.cbegin(); }
    constexpr auto cbegin() const { return halfplanes_.cbegin(); }
    constexpr auto end() const { return halfplanes_.cend(); }
    constexpr auto cend() const { return halfplanes_.cend(); }

    constexpr bool empty() const {
        return empty_;
    }

    constexpr bool isPlane() const {
        return !empty_ && halfplanes_.empty();
    }

    constexpr bool isDegenerate() const {
        return empty_ || degenerate_;
    }

    constexpr bool isUndefined() const {
        return false;
    }

    [[nodiscard]] constexpr bool isHalfplane() const;

    [[nodiscard]] constexpr std::optional<HalfplaneType> getIfHalfplane() const;

    [[nodiscard]] constexpr bool isLine() const;

    [[nodiscard]] constexpr std::optional<Line<PointType>> getIfLine() const;

    [[nodiscard]] constexpr bool isRay() const;

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr std::optional<Ray<Point<ResultNumber, typename PointType::LabelType>>>
    getIfRay() const;

    [[nodiscard]] constexpr bool isPoint() const;

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
    getIfPoint() const;

    [[nodiscard]] constexpr bool isSegment() const;

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr std::optional<Segment<Point<ResultNumber, typename PointType::LabelType>>>
    getIfSegment() const;

    constexpr bool isBounded() const {
        if (empty_) {
            return true;
        }
        const std::size_t n = size();
        if (n < 3) {
            return false;
        }
        for (std::size_t i = 0; i < n; ++i) {
            if (!(detail::directionCross(halfplanes_[i], halfplanes_[nextIndex(i)]) > 0)) {
                return false;
            }
        }
        return true;
    }

    constexpr std::size_t vertexCount() const {
        if (empty_) {
            return 0;
        }
        const std::size_t n = size();
        if (n < 2) {
            return 0;
        }
        std::size_t count = 0;
        for (std::size_t i = 0; i < n; ++i) {
            if (detail::directionCross(halfplanes_[i], halfplanes_[nextIndex(i)]) > 0) {
                ++count;
            }
        }
        return count;
    }

    constexpr bool vertexExists(std::size_t i) const {
        if (empty_ || size() < 2) {
            return false;
        }
        return detail::directionCross(halfplanes_[i], halfplanes_[nextIndex(i)]) > 0;
    }

    template <class ResultNumber = division_result_t<NumberType>>
    constexpr Point<ResultNumber, typename PointType::LabelType> vertex(std::size_t i) const {
        assert(vertexExists(i));
        const auto isec = halfplanes_[i].asLine().template intersection<ResultNumber>(
            halfplanes_[nextIndex(i)].asLine());
        assert(isec && isec->index() == 0);
        return std::get<0>(*isec);
    }

    template <class ResultNumber = division_result_t<NumberType>>
    constexpr std::vector<Point<ResultNumber, typename PointType::LabelType>> vertices() const {
        std::vector<Point<ResultNumber, typename PointType::LabelType>> result;
        if (empty_ || size() < 2) {
            return result;
        }
        for (std::size_t i = 0; i < size(); ++i) {
            if (vertexExists(i)) {
                result.push_back(vertex<ResultNumber>(i));
            }
        }
        return result;
    }

    template <class ResultNumber = division_result_t<NumberType>>
    constexpr std::variant<Segment<Point<ResultNumber, typename PointType::LabelType>>,
                           Ray<Point<ResultNumber, typename PointType::LabelType>>,
                           Line<Point<ResultNumber, typename PointType::LabelType>>>
    edge(std::size_t i) const {
        using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
        assert(!empty_ && i < size());
        const std::size_t prev = prevIndex(i);
        const bool hasStart = size() >= 2 && detail::directionCross(halfplanes_[prev], halfplanes_[i]) > 0;
        const bool hasEnd = vertexExists(i);
        const auto& h = halfplanes_[i];
        if (hasStart && hasEnd) {
            return Segment<ResultPoint>(vertex<ResultNumber>(prev), vertex<ResultNumber>(i));
        }
        const ResultPoint source(h.source());
        const ResultPoint target(h.target());
        const ResultNumber dx = target.x() - source.x();
        const ResultNumber dy = target.y() - source.y();
        if (hasStart) {
            const ResultPoint start = vertex<ResultNumber>(prev);
            return Ray<ResultPoint>(start, ResultPoint(start.x() + dx, start.y() + dy));
        }
        if (hasEnd) {
            const ResultPoint finish = vertex<ResultNumber>(i);
            return Ray<ResultPoint>(finish, ResultPoint(finish.x() - dx, finish.y() - dy));
        }
        return Line<ResultPoint>(source, target);
    }

    template <class ResultNumber = division_result_t<NumberType>>
    constexpr Convex<Point<ResultNumber, typename PointType::LabelType>> asConvex() const {
        if (empty_) {
            return {};
        }
        if (!isBounded()) {
            throw std::logic_error("HalfplaneIntersection::asConvex requires a bounded region");
        }
        return Convex<Point<ResultNumber, typename PointType::LabelType>>(vertices<ResultNumber>());
    }

    template <class ResultNumber = division_result_t<NumberType>>
    constexpr Convex<Point<ResultNumber, typename PointType::LabelType>> convexHull() const {
        return asConvex<ResultNumber>();
    }

    [[nodiscard]] constexpr bool operator==(const HalfplaneIntersection& other) const {
        if (empty_ != other.empty_) {
            return false;
        }
        if (empty_) {
            return true;
        }
        if (halfplanes_.size() != other.halfplanes_.size()) {
            return false;
        }
        for (std::size_t i = 0; i < halfplanes_.size(); ++i) {
            if (!(halfplanes_[i] == other.halfplanes_[i])) {
                return false;
            }
        }
        return true;
    }

    template<AnyShapeConcept OtherShape>
    [[nodiscard]] constexpr bool samePointSet(const OtherShape& other) const;

    [[nodiscard]] constexpr auto operator<=>(const HalfplaneIntersection& other) const {
        if (auto cmp = empty_ <=> other.empty_; cmp != 0) {
            return cmp;
        }
        if (empty_) {
            return std::strong_ordering::equal;
        }
        if (auto cmp = halfplanes_.size() <=> other.halfplanes_.size(); cmp != 0) {
            return cmp;
        }
        for (std::size_t i = 0; i < halfplanes_.size(); ++i) {
            if (auto cmp = halfplanes_[i] <=> other.halfplanes_[i]; cmp != 0) {
                return cmp;
            }
        }
        return std::strong_ordering::equal;
    }

    template <class ResultNumber = division_result_t<NumberType>>
    constexpr Rectangle<Point<ResultNumber, typename PointType::LabelType>> bbox() const;

    template <class ResultNumber = grid_number_t<typename PointType_::NumberType>>
        requires(detail::extended_integral<ResultNumber> || std::same_as<ResultNumber, BigInt>)
    [[nodiscard]] std::vector<Point<ResultNumber, typename PointType::LabelType>>
    latticePoints() const;

    template <std::floating_point ResultNumber = double>
    constexpr Rectangle<Point<ResultNumber>> fbox() const;

    template <class OtherShape>
        requires MinkowskiSummableConcept<HalfplaneIntersection<PointType_, TLabel>, OtherShape>
    [[nodiscard]] constexpr auto minkowskiSum(const OtherShape& other) const;

    template <class OtherShape>
        requires MinkowskiSummableConcept<HalfplaneIntersection<PointType_, TLabel>, OtherShape>
    [[nodiscard]] constexpr auto minkowskiErosion(const OtherShape& other) const;

    template <PointConcept OtherPoint>
    constexpr HalfplaneIntersection& operator+=(const OtherPoint& translation);

    template <PointConcept OtherPoint>
    constexpr HalfplaneIntersection& operator-=(const OtherPoint& translation);

    template <class Scalar>
        requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
    constexpr HalfplaneIntersection& operator*=(const Scalar& scalar);

    template <class Scalar>
        requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
    constexpr HalfplaneIntersection& operator/=(const Scalar& scalar);

    [[nodiscard]] constexpr HalfplaneIntersection rotated90(int k = 1) const;

    constexpr void rotate90(int k = 1);

    template <class OtherNumber>
    [[nodiscard]] constexpr HalfplaneIntersection scaledUpX(const OtherNumber scalar) const;

    template <class OtherNumber>
    constexpr void scaleUpX(const OtherNumber scalar);

    template <class OtherNumber>
    [[nodiscard]] constexpr HalfplaneIntersection scaledUpY(const OtherNumber scalar) const;

    template <class OtherNumber>
    constexpr void scaleUpY(const OtherNumber scalar);

    template <class OtherNumber>
    [[nodiscard]] constexpr HalfplaneIntersection scaledDownX(const OtherNumber scalar) const;

    template <class OtherNumber>
    constexpr void scaleDownX(const OtherNumber scalar);

    template <class OtherNumber>
    [[nodiscard]] constexpr HalfplaneIntersection scaledDownY(const OtherNumber scalar) const;

    template <class OtherNumber>
    constexpr void scaleDownY(const OtherNumber scalar);
 
    // --- predicates (defined in the implementation layer) ---

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool contains(const OtherPoint& point) const;

    template <SegmentConcept OtherSegment>
    [[nodiscard]] constexpr bool contains(const OtherSegment& other) const;

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool contains(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool contains(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool contains(const OtherOrientedLine& other) const;

    template <RayConcept OtherRay>
    [[nodiscard]] constexpr bool contains(const OtherRay& other) const;

    template <HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool contains(const OtherHalfplane& other) const;

    template <RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool contains(const OtherRectangle& other) const;

    template <TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool contains(const OtherTriangle& other) const;

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool contains(const OtherDisk& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool contains(const OtherConvex& other) const;

    template <MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool contains(const OtherChain& other) const;

    template <PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool contains(const OtherPolyline& other) const;

    template <PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool contains(const OtherPolygon& other) const;

    template <HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool contains(const OtherRegion& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool boundaryContains(const OtherPoint& point) const;

    template <SegmentConcept OtherSegment>
    [[nodiscard]] constexpr bool boundaryContains(const OtherSegment& other) const;

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool boundaryContains(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool boundaryContains(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool boundaryContains(const OtherOrientedLine& other) const;

    template <RayConcept OtherRay>
    [[nodiscard]] constexpr bool boundaryContains(const OtherRay& other) const;

    template <HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool boundaryContains(const OtherHalfplane& other) const;

    template <RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool boundaryContains(const OtherRectangle& other) const;

    template <TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool boundaryContains(const OtherTriangle& other) const;

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool boundaryContains(const OtherDisk& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool boundaryContains(const OtherConvex& other) const;

    template <MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool boundaryContains(const OtherChain& other) const;

    template <PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool boundaryContains(const OtherPolyline& other) const;

    template <PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool boundaryContains(const OtherPolygon& other) const;

    template <HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool boundaryContains(const OtherRegion& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool interiorContains(const OtherPoint& point) const;

    template <SegmentConcept OtherSegment>
    [[nodiscard]] constexpr bool interiorContains(const OtherSegment& other) const;

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool interiorContains(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool interiorContains(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool interiorContains(const OtherOrientedLine& other) const;

    template <RayConcept OtherRay>
    [[nodiscard]] constexpr bool interiorContains(const OtherRay& other) const;

    template <HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool interiorContains(const OtherHalfplane& other) const;

    template <RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool interiorContains(const OtherRectangle& other) const;

    template <TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool interiorContains(const OtherTriangle& other) const;

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool interiorContains(const OtherDisk& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool interiorContains(const OtherConvex& other) const;

    template <MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool interiorContains(const OtherChain& other) const;

    template <PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool interiorContains(const OtherPolyline& other) const;

    template <PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool interiorContains(const OtherPolygon& other) const;

    template <HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool interiorContains(const OtherRegion& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool intersects(const OtherPoint& other) const;

    template <SegmentConcept OtherSegment>
    [[nodiscard]] constexpr bool intersects(const OtherSegment& other) const;

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool intersects(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool intersects(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool intersects(const OtherOrientedLine& other) const;

    template <RayConcept OtherRay>
    [[nodiscard]] constexpr bool intersects(const OtherRay& other) const;

    template <HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool intersects(const OtherHalfplane& other) const;

    template <RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool intersects(const OtherRectangle& other) const;

    template <TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool intersects(const OtherTriangle& other) const;

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool intersects(const OtherDisk& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool intersects(const OtherConvex& other) const;

    template <MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool intersects(const OtherChain& other) const;

    template <PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool intersects(const OtherPolyline& other) const;

    template <PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool intersects(const OtherPolygon& other) const;

    template <HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool intersects(const OtherRegion& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherPoint& other) const;

    template <SegmentConcept OtherSegment>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherSegment& other) const;

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherOrientedLine& other) const;

    template <RayConcept OtherRay>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherRay& other) const;

    template <HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherHalfplane& other) const;

    template <RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherRectangle& other) const;

    template <TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherTriangle& other) const;

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherDisk& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherConvex& other) const;

    template <MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherChain& other) const;

    template <PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherPolyline& other) const;

    template <PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherPolygon& other) const;

    template <HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherRegion& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool separates(const OtherPoint& other) const;

    template <SegmentConcept OtherSegment>
    [[nodiscard]] constexpr bool separates(const OtherSegment& other) const;

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool separates(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool separates(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool separates(const OtherOrientedLine& other) const;

    template <RayConcept OtherRay>
    [[nodiscard]] constexpr bool separates(const OtherRay& other) const;

    template <HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool separates(const OtherHalfplane& other) const;

    template <RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool separates(const OtherRectangle& other) const;

    template <TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool separates(const OtherTriangle& other) const;

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool separates(const OtherDisk& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool separates(const OtherConvex& other) const;

    template <MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool separates(const OtherChain& other) const;

    template <PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool separates(const OtherPolyline& other) const;

    template <PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool separates(const OtherPolygon& other) const;

    template <HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool separates(const OtherRegion& other) const;

    template <PolygonWithHolesConcept OtherHoledRegion>
    [[nodiscard]] constexpr bool contains(const OtherHoledRegion& other) const;

    template <PolygonWithHolesConcept OtherHoledRegion>
    [[nodiscard]] constexpr bool boundaryContains(const OtherHoledRegion& other) const;

    template <PolygonWithHolesConcept OtherHoledRegion>
    [[nodiscard]] constexpr bool interiorContains(const OtherHoledRegion& other) const;

    template <PolygonWithHolesConcept OtherHoledRegion>
    [[nodiscard]] bool separates(const OtherHoledRegion& other) const;
 
    // -------------------------------------------------------------------------
    // A set of regions
    //
    // It outranks every other shape, so the symmetric relations reach it through
    // the rank-based forwarders and only the asymmetric ones are answered here.
    // A set is the union of its components, so it is contained exactly when
    // every component is — no matter what this shape is.

    template<PolygonSetConcept OtherSet>
    [[nodiscard]] constexpr bool contains(const OtherSet& other) const {
        for (const auto& component : other) {
            if (!contains(component)) {
                return false;
            }
        }
        return true;
    }

    template<PolygonSetConcept OtherSet>
    [[nodiscard]] constexpr bool boundaryContains(const OtherSet& other) const {
        for (const auto& component : other) {
            if (!boundaryContains(component)) {
                return false;
            }
        }
        return true;
    }

    template<PolygonSetConcept OtherSet>
    [[nodiscard]] constexpr bool interiorContains(const OtherSet& other) const {
        for (const auto& component : other) {
            if (!interiorContains(component)) {
                return false;
            }
        }
        return true;
    }

    template<PolygonSetConcept OtherSet>
    [[nodiscard]] bool separates(const OtherSet& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool crosses(const OtherPoint& other) const;

    template <SegmentConcept OtherSegment>
    [[nodiscard]] constexpr bool crosses(const OtherSegment& other) const;

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool crosses(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool crosses(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool crosses(const OtherOrientedLine& other) const;

    template <RayConcept OtherRay>
    [[nodiscard]] constexpr bool crosses(const OtherRay& other) const;

    template <HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool crosses(const OtherHalfplane& other) const;

    template <RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool crosses(const OtherRectangle& other) const;

    template <TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool crosses(const OtherTriangle& other) const;

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool crosses(const OtherDisk& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool crosses(const OtherConvex& other) const;

    template <MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool crosses(const OtherChain& other) const;

    template <PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool crosses(const OtherPolyline& other) const;

    template <PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool crosses(const OtherPolygon& other) const;

    template <HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool crosses(const OtherRegion& other) const;

    template <typename OtherShape>
        requires (!PointConcept<OtherShape> &&
                  detail::shapeRank<OtherShape> > detail::shapeRank<HalfplaneIntersection>)
    [[nodiscard]] constexpr bool intersects(const OtherShape& other) const {
        return other.intersects(*this);
    }

    template <typename OtherShape>
        requires (!PointConcept<OtherShape> &&
                  detail::shapeRank<OtherShape> > detail::shapeRank<HalfplaneIntersection>)
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherShape& other) const {
        return other.interiorsIntersect(*this);
    }

    template <typename OtherShape>
        requires (!PointConcept<OtherShape> &&
                  detail::shapeRank<OtherShape> > detail::shapeRank<HalfplaneIntersection> &&
                  requires(const OtherShape& o, const HalfplaneIntersection& self) { o.crosses(self); })
    [[nodiscard]] constexpr bool crosses(const OtherShape& other) const {
        return other.crosses(*this);
    }
 
    // The empty set is a subset of every shape and disjoint from all of them.
    template <class EmptyPoint>
    [[nodiscard]] constexpr bool contains(const EmptyShape<EmptyPoint>&) const {
        return true;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool boundaryContains(const EmptyShape<EmptyPoint>&) const {
        return true;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool interiorContains(const EmptyShape<EmptyPoint>&) const {
        return true;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool intersects(const EmptyShape<EmptyPoint>&) const {
        return false;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool interiorsIntersect(const EmptyShape<EmptyPoint>&) const {
        return false;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool separates(const EmptyShape<EmptyPoint>&) const {
        return false;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool crosses(const EmptyShape<EmptyPoint>&) const {
        return false;
    }

    template <class ResultNumber = NumberType, class EmptyPoint>
    [[nodiscard]] constexpr EmptyShape<EmptyPoint> intersection(const EmptyShape<EmptyPoint>&) const {
        return {};
    }
 
    // Runtime Shape argument: visit the wrapped alternative and re-dispatch to
    // the matching per-shape overload (defined in the implementation layer).

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool contains(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool boundaryContains(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool interiorContains(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool intersects(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool interiorsIntersect(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool separates(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] constexpr bool crosses(const Shape<OtherPoint>& other) const;

    template <class ResultNumber = NumberType, PointConcept OtherPoint>
    [[nodiscard]] constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
    intersection(const OtherPoint& other) const;

    template <class ResultNumber = division_result_t<NumberType>, LineConcept OtherLine>
    [[nodiscard]] constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                                       Segment<Point<ResultNumber, typename PointType::LabelType>>,
                                                       Ray<Point<ResultNumber, typename PointType::LabelType>>,
                                                       Line<Point<ResultNumber, typename PointType::LabelType>>>>
    intersection(const OtherLine& other) const;

    template <class ResultNumber = division_result_t<NumberType>, OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                                       Segment<Point<ResultNumber, typename PointType::LabelType>>,
                                                       Ray<Point<ResultNumber, typename PointType::LabelType>>,
                                                       Line<Point<ResultNumber, typename PointType::LabelType>>>>
    intersection(const OtherOrientedLine& other) const;

    template <class ResultNumber = division_result_t<NumberType>, SegmentConcept OtherSegment>
    [[nodiscard]] constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                                       Segment<Point<ResultNumber, typename PointType::LabelType>>>>
    intersection(const OtherSegment& other) const;

    template <class ResultNumber = division_result_t<NumberType>, OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                                       Segment<Point<ResultNumber, typename PointType::LabelType>>>>
    intersection(const OtherOrientedSegment& other) const;

    template <class ResultNumber = division_result_t<NumberType>, RayConcept OtherRay>
    [[nodiscard]] constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                                       Segment<Point<ResultNumber, typename PointType::LabelType>>,
                                                       Ray<Point<ResultNumber, typename PointType::LabelType>>>>
    intersection(const OtherRay& other) const;

    template <class ResultNumber = NumberType, HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
    intersection(const OtherHalfplane& other) const;

    template <class ResultNumber = NumberType, RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
    intersection(const OtherRectangle& other) const;

    template <class ResultNumber = NumberType, TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
    intersection(const OtherTriangle& other) const;

    template <class ResultNumber = NumberType, ConvexConcept OtherConvex>
    [[nodiscard]] constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
    intersection(const OtherConvex& other) const;

    template <class ResultNumber = NumberType, HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
    intersection(const OtherRegion& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PolygonConcept OtherPolygon>
    [[nodiscard]] std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                           Polyline<Point<ResultNumber, typename PointType::LabelType>>,
                                           Polygon<Point<ResultNumber, typename PointType::LabelType>>>>
    intersection(const OtherPolygon& other) const;

    template <class ResultNumber = division_result_t<NumberType>, typename OtherShape>
        requires (!PointConcept<OtherShape> &&
                  (detail::shapeRank<OtherShape> > detail::shapeRank<HalfplaneIntersection>) &&
                  requires(const OtherShape& o, const HalfplaneIntersection& self) {
                      o.template intersection<ResultNumber>(self);
                  })
    [[nodiscard]] auto intersection(const OtherShape& other) const {
        return other.template intersection<ResultNumber>(*this);
    }

    template <class ResultNumber = division_result_t<NumberType>, typename OtherShape>
        requires (!PointConcept<OtherShape>
                  && (detail::shapeRank<OtherShape> > detail::shapeRank<HalfplaneIntersection>)
                  && requires(const OtherShape& o, const HalfplaneIntersection& self) {
                         o.template regularizedIntersection<ResultNumber>(self);
                     })
    [[nodiscard]] auto regularizedIntersection(const OtherShape& other) const {
        return other.template regularizedIntersection<ResultNumber>(*this);
    }
 
    // --- distances (defined in the implementation layer) ---
    //
    // Zero when the shapes intersect; otherwise the minimum over the region's
    // boundary edges of the edge-to-shape distance. Complexity: O(n) edge
    // queries for n stored half-planes. The empty region has no distance to
    // anything; querying it is undefined behavior.

    template <class ResultNumber = division_result_t<NumberType>, PointConcept OtherPoint>
    [[nodiscard]] constexpr auto squaredDistance(const OtherPoint& other) const;

    template <class ResultNumber = division_result_t<NumberType>, SegmentConcept OtherSegment>
    [[nodiscard]] constexpr auto squaredDistance(const OtherSegment& other) const;

    template <class ResultNumber = division_result_t<NumberType>, OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr auto squaredDistance(const OtherOrientedSegment& other) const;

    template <class ResultNumber = division_result_t<NumberType>, LineConcept OtherLine>
    [[nodiscard]] constexpr auto squaredDistance(const OtherLine& other) const;

    template <class ResultNumber = division_result_t<NumberType>, OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr auto squaredDistance(const OtherOrientedLine& other) const;

    template <class ResultNumber = division_result_t<NumberType>, RayConcept OtherRay>
    [[nodiscard]] constexpr auto squaredDistance(const OtherRay& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr auto squaredDistance(const OtherHalfplane& other) const;

    template <class ResultNumber = division_result_t<NumberType>, RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr auto squaredDistance(const OtherRectangle& other) const;

    template <class ResultNumber = division_result_t<NumberType>, TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr auto squaredDistance(const OtherTriangle& other) const;

    template <class ResultNumber = division_result_t<NumberType>, ConvexConcept OtherConvex>
    [[nodiscard]] constexpr auto squaredDistance(const OtherConvex& other) const;

    template <class ResultNumber = division_result_t<NumberType>, MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr auto squaredDistance(const OtherChain& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr auto squaredDistance(const OtherPolyline& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr auto squaredDistance(const OtherPolygon& other) const;

    template <class ResultNumber = double, DiskConcept OtherDisk>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherDisk& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr auto squaredDistance(const OtherRegion& other) const;

    template <class ResultNumber = division_result_t<NumberType>, class OtherShape>
        requires detail::ClosestPointsPairConcept<HalfplaneIntersection<PointType_, TLabel>, OtherShape>
    [[nodiscard]] constexpr auto closestPoints(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PointConcept OtherPoint>
    [[nodiscard]] constexpr auto distanceL1(const OtherPoint& other) const;

    template <class ResultNumber = division_result_t<NumberType>, SegmentConcept OtherSegment>
    [[nodiscard]] constexpr auto distanceL1(const OtherSegment& other) const;

    template <class ResultNumber = division_result_t<NumberType>, OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr auto distanceL1(const OtherOrientedSegment& other) const;

    template <class ResultNumber = division_result_t<NumberType>, LineConcept OtherLine>
    [[nodiscard]] constexpr auto distanceL1(const OtherLine& other) const;

    template <class ResultNumber = division_result_t<NumberType>, OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr auto distanceL1(const OtherOrientedLine& other) const;

    template <class ResultNumber = division_result_t<NumberType>, RayConcept OtherRay>
    [[nodiscard]] constexpr auto distanceL1(const OtherRay& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr auto distanceL1(const OtherHalfplane& other) const;

    template <class ResultNumber = division_result_t<NumberType>, RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr auto distanceL1(const OtherRectangle& other) const;

    template <class ResultNumber = division_result_t<NumberType>, TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr auto distanceL1(const OtherTriangle& other) const;

    template <class ResultNumber = division_result_t<NumberType>, ConvexConcept OtherConvex>
    [[nodiscard]] constexpr auto distanceL1(const OtherConvex& other) const;

    template <class ResultNumber = division_result_t<NumberType>, MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr auto distanceL1(const OtherChain& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr auto distanceL1(const OtherPolyline& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr auto distanceL1(const OtherPolygon& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr auto distanceL1(const OtherRegion& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PointConcept OtherPoint>
    [[nodiscard]] constexpr auto intersection(const Shape<OtherPoint>& other) const {
        return other.template intersection<ResultNumber>(*this);
    }

    template <class ResultNumber = division_result_t<NumberType>, PointConcept OtherPoint>
    [[nodiscard]] auto regularizedIntersection(const Shape<OtherPoint>& other) const {
        return other.template regularizedIntersection<ResultNumber>(*this);
    }

    template <class ResultNumber = double, PointConcept OtherPoint>
    [[nodiscard]] constexpr auto distanceL1(const Shape<OtherPoint>& other) const {
        return other.template distanceL1<ResultNumber>(*this);
    }

    template <class ResultNumber = division_result_t<NumberType>, PointConcept OtherPoint>
    [[nodiscard]] constexpr auto distanceLInf(const OtherPoint& other) const;

    template <class ResultNumber = division_result_t<NumberType>, SegmentConcept OtherSegment>
    [[nodiscard]] constexpr auto distanceLInf(const OtherSegment& other) const;

    template <class ResultNumber = division_result_t<NumberType>, OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr auto distanceLInf(const OtherOrientedSegment& other) const;

    template <class ResultNumber = division_result_t<NumberType>, LineConcept OtherLine>
    [[nodiscard]] constexpr auto distanceLInf(const OtherLine& other) const;

    template <class ResultNumber = division_result_t<NumberType>, OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr auto distanceLInf(const OtherOrientedLine& other) const;

    template <class ResultNumber = division_result_t<NumberType>, RayConcept OtherRay>
    [[nodiscard]] constexpr auto distanceLInf(const OtherRay& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr auto distanceLInf(const OtherHalfplane& other) const;

    template <class ResultNumber = division_result_t<NumberType>, RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr auto distanceLInf(const OtherRectangle& other) const;

    template <class ResultNumber = division_result_t<NumberType>, TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr auto distanceLInf(const OtherTriangle& other) const;

    template <class ResultNumber = division_result_t<NumberType>, ConvexConcept OtherConvex>
    [[nodiscard]] constexpr auto distanceLInf(const OtherConvex& other) const;

    template <class ResultNumber = division_result_t<NumberType>, MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr auto distanceLInf(const OtherChain& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr auto distanceLInf(const OtherPolyline& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr auto distanceLInf(const OtherPolygon& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr auto distanceLInf(const OtherRegion& other) const;

    template <class ResultNumber = double, PointConcept OtherPoint>
    [[nodiscard]] constexpr auto distanceLInf(const Shape<OtherPoint>& other) const {
        return other.template distanceLInf<ResultNumber>(*this);
    }
 
    // Distances to a higher-ranked shape, forwarded so that each unordered pair
    // needs its implementation only once, on that shape.

    template <class ResultNumber = division_result_t<NumberType>, typename OtherShape>
        requires ((detail::shapeRank<OtherShape> > detail::shapeRank<HalfplaneIntersection>) &&
                  requires(const OtherShape& o, const HalfplaneIntersection& self) {
                      o.template squaredDistance<ResultNumber>(self);
                  })
    [[nodiscard]] constexpr auto squaredDistance(const OtherShape& other) const {
        return other.template squaredDistance<ResultNumber>(*this);
    }

    template <class ResultNumber = division_result_t<NumberType>, typename OtherShape>
        requires ((detail::shapeRank<OtherShape> > detail::shapeRank<HalfplaneIntersection>) &&
                  requires(const OtherShape& o, const HalfplaneIntersection& self) {
                      o.template distanceL1<ResultNumber>(self);
                  })
    [[nodiscard]] constexpr auto distanceL1(const OtherShape& other) const {
        return other.template distanceL1<ResultNumber>(*this);
    }

    template <class ResultNumber = division_result_t<NumberType>, typename OtherShape>
        requires ((detail::shapeRank<OtherShape> > detail::shapeRank<HalfplaneIntersection>) &&
                  requires(const OtherShape& o, const HalfplaneIntersection& self) {
                      o.template distanceLInf<ResultNumber>(self);
                  })
    [[nodiscard]] constexpr auto distanceLInf(const OtherShape& other) const {
        return other.template distanceLInf<ResultNumber>(*this);
    }
 
    // --- measures (defined in the implementation layer) ---

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr auto twiceArea() const;

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr auto area() const;

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr Point<ResultNumber> centroid() const;

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr Point<ResultNumber> pointInside() const;

    template <class OtherShape>
    [[nodiscard]] constexpr bool pointInsideInteriorContainedIn(const OtherShape& shape) const;
 
  private:
    std::vector<HalfplaneType> halfplanes_{};
    bool empty_ = false;
    bool degenerate_ = false;
    [[no_unique_address]] mutable LabelType label_{};
 
    // Memoized hash, mirroring Convex: hashUnset_ means "not yet computed" and
    // the one true hash colliding with it is remapped to hashUnset_ - 1.
    static constexpr std::size_t hashUnset_ = pgl::detail::numeric_limits<std::size_t>::max();
    mutable std::size_t hash_ = hashUnset_;
    friend struct std::hash<HalfplaneIntersection>;
 
    template <class OtherPointType, class OtherLabelType>
    friend struct HalfplaneIntersection;
 
    constexpr void resetCache() const {
        hash_ = hashUnset_;
    }
 
    constexpr std::size_t nextIndex(std::size_t i) const {
        return i + 1 < halfplanes_.size() ? i + 1 : 0;
    }
 
    constexpr std::size_t prevIndex(std::size_t i) const {
        return i == 0 ? halfplanes_.size() - 1 : i - 1;
    }
 
    // Appends the four half-planes of the axis-aligned slab pinning the single
    // point p and marks the region degenerate. Assumes the region is empty.
    constexpr void buildPointSlab(const PointType& p) {
        const PointType px(p.x() + NumberType(1), p.y());
        const PointType py(p.x(), p.y() + NumberType(1));
        halfplanes_.push_back(HalfplaneType(p, px));
        halfplanes_.push_back(HalfplaneType(px, p));
        halfplanes_.push_back(HalfplaneType(p, py));
        halfplanes_.push_back(HalfplaneType(py, p));
        degenerate_ = true;
    }
 
    // Appends the four half-planes clamping to the segment a-b (a != b): the
    // supporting-line slab plus a perpendicular clamp through each endpoint
    // (rotating the direction by 90 degrees keeps integer coordinates). Marks
    // the region degenerate. Assumes the region is empty.
    constexpr void buildSegmentSlab(const PointType& a, const PointType& b) {
        const NumberType dx = b.x() - a.x();
        const NumberType dy = b.y() - a.y();
        const PointType aPerp(a.x() - dy, a.y() + dx);
        const PointType bPerp(b.x() - dy, b.y() + dx);
        halfplanes_.push_back(HalfplaneType(a, b));
        halfplanes_.push_back(HalfplaneType(b, a));
        halfplanes_.push_back(HalfplaneType(aPerp, a));
        halfplanes_.push_back(HalfplaneType(b, bPerp));
        degenerate_ = true;
    }
 
    // Restores the sorted-by-pseudo-angle invariant after bulk construction or
    // a transformation. The input is assumed pairwise direction-distinct.
    constexpr void canonicalizeSorted() {
        std::sort(halfplanes_.begin(), halfplanes_.end(),
                  [](const HalfplaneType& a, const HalfplaneType& b) { return detail::directionLess(a, b); });
        resetCache();
    }
 
    // First stored index whose boundary direction is >= the query's in the
    // linear pseudo-angle order; may equal size().
    template <HalfplaneConcept Query>
    constexpr std::size_t linearLowerBound(const Query& query) const {
        const auto it = std::lower_bound(
            halfplanes_.begin(), halfplanes_.end(), query,
            [](const HalfplaneType& element, const Query& value) { return detail::directionLess(element, value); });
        return static_cast<std::size_t>(it - halfplanes_.begin());
    }
 
    // First stored index whose boundary direction is > the query's in the
    // linear pseudo-angle order; may equal size().
    template <HalfplaneConcept Query>
    constexpr std::size_t linearUpperBound(const Query& query) const {
        const auto it = std::upper_bound(
            halfplanes_.begin(), halfplanes_.end(), query,
            [](const Query& value, const HalfplaneType& element) { return detail::directionLess(value, element); });
        return static_cast<std::size_t>(it - halfplanes_.begin());
    }
 
    // Index of the stored half-plane with the query's boundary direction, or
    // -1 when there is none.
    template <HalfplaneConcept Query>
    constexpr std::ptrdiff_t sameDirectionIndex(const Query& query) const {
        const std::size_t pos = linearLowerBound(query);
        if (pos < halfplanes_.size() && detail::directionEqual(halfplanes_[pos], query)) {
            return static_cast<std::ptrdiff_t>(pos);
        }
        return -1;
    }
 
    // The stored half-planes whose boundary direction lies strictly inside the
    // open half-circle counterclockwise from the query's direction — exactly
    // those whose outward normal has positive dot product with it, i.e. the
    // constraints imposing an upper bound on the parameter along the query
    // line. Returns {first cyclic index, count}; the range is contiguous in
    // cyclic order. Precondition: size() >= 1.
    template <HalfplaneConcept Query>
    constexpr std::pair<std::size_t, std::size_t> leftArc(const Query& query) const {
        const std::size_t n = halfplanes_.size();
        std::size_t lo = linearUpperBound(query);
        if (lo == n) {
            lo = 0;
        }
        std::size_t hi = linearLowerBound(query.opposite());
        if (hi == n) {
            hi = 0;
        }
        std::size_t count = (hi + n - lo) % n;
        if (count == 0 && detail::directionCross(query, halfplanes_[lo]) > 0) {
            count = n;  // every stored direction lies inside the half-circle
        }
        return {lo, count};
    }
 
    // Among the leftArc constraints, the one whose boundary the directed query
    // line crosses first (the binding upper bound). Returns -1 when the arc is
    // empty (the line never leaves the region going forward).
    // Precondition: size() >= 1 and the query is not degenerate.
    template <HalfplaneConcept Query>
    constexpr std::ptrdiff_t exitConstraint(const Query& query) const {
        const auto [lo, count] = leftArc(query);
        if (count == 0) {
            return -1;
        }
        const std::size_t n = halfplanes_.size();
        // The crossing parameters along the arc first decrease, then increase
        // (the arc's boundary lines wrap a convex chain), so a valley binary
        // search finds the minimum. The consecutive-pair comparison is the
        // vertex side test: the shared vertex strictly outside the query
        // half-plane means the parameters are still decreasing.
        std::size_t first = 0;
        std::size_t last = count - 1;
        while (first < last) {
            const std::size_t mid = first + (last - first) / 2;
            const std::size_t i = (lo + mid) % n;
            const std::size_t j = (lo + mid + 1) % n;
            if (detail::vertexSide(halfplanes_[i], halfplanes_[j], query) < 0) {
                first = mid + 1;
            } else {
                last = mid;
            }
        }
        return static_cast<std::ptrdiff_t>((lo + first) % n);
    }
 
    // Side of a point relative to stored constraint i: positive means strictly
    // inside, zero on the boundary, negative strictly outside.
    template <PointConcept OtherPoint>
    constexpr std::partial_ordering constraintSide(std::size_t i, const OtherPoint& point) const {
        return orientationSign(halfplanes_[i].source(), halfplanes_[i].target(), point);
    }
 
    enum class SupStatus { unbounded, above, on, below };
 
    // Status of sup_{p in region} a·p compared to b, where the query
    // half-plane is {p : a·p <= b}: `below`/`on`/`above` when the supremum is
    // attained and compares so with b, `unbounded` when it is infinite.
    // Querying with the opposite half-plane gives the infimum side: `below`
    // then means inf > b, `on` means inf == b (attained).
    // Precondition: !empty_ and size() >= 1; a degenerate query is UB.
    template <HalfplaneConcept Query>
    constexpr SupStatus supStatus(const Query& query) const {
        const std::size_t n = halfplanes_.size();
        // The supremum in direction of the query's outward normal is attained
        // where the stored boundary directions bracket the query's direction.
        const std::size_t pos = linearUpperBound(query);
        const std::size_t predIdx = (pos == 0 ? n : pos) - 1;
        const auto& pred = halfplanes_[predIdx];
        if (detail::directionEqual(pred, query)) {
            // Attained along pred's whole boundary line, where a·p is
            // constant: compare via any boundary point.
            const auto side = orientationSign(query.source(), query.target(), pred.source());
            if (side > 0) {
                return SupStatus::below;
            }
            return side == 0 ? SupStatus::on : SupStatus::above;
        }
        const std::size_t succIdx = pos == n ? 0 : pos;
        if (succIdx == predIdx) {
            return SupStatus::unbounded;  // single stored half-plane, gap 2*pi
        }
        if (!(detail::directionCross(pred, halfplanes_[succIdx]) > 0)) {
            // The query direction lies strictly inside a gap of at least pi:
            // the region is unbounded toward the query's outward normal.
            return SupStatus::unbounded;
        }
        const int side = detail::vertexSide(pred, halfplanes_[succIdx], query);
        if (side > 0) {
            return SupStatus::below;
        }
        return side == 0 ? SupStatus::on : SupStatus::above;
    }
 
    // Exact clip of the directed query line against the region.
    // entry/exit are constraint indices (-1: unbounded on that side);
    // onParallelBoundary reports a constraint parallel to the line whose
    // boundary contains it (the line runs along a face).
    // Precondition: !empty_ and size() >= 1; a degenerate query is UB.
    struct ClipResult {
        bool empty = false;
        bool onParallelBoundary = false;
        std::ptrdiff_t entry = -1;
        std::ptrdiff_t exit = -1;
    };
 
    template <HalfplaneConcept Query>
    constexpr ClipResult clipLine(const Query& query) const {
        ClipResult result;
        const auto reversed = query.opposite();
        // Constraints parallel to the line hold (or fail) along all of it.
        for (const std::ptrdiff_t parallelIdx : {sameDirectionIndex(query), sameDirectionIndex(reversed)}) {
            if (parallelIdx >= 0) {
                const auto side = constraintSide(static_cast<std::size_t>(parallelIdx), query.source());
                if (side < 0) {
                    result.empty = true;
                    return result;
                }
                if (side == 0) {
                    result.onParallelBoundary = true;
                }
            }
        }
        result.exit = exitConstraint(query);
        result.entry = exitConstraint(reversed);
        if (result.entry >= 0 && result.exit >= 0) {
            // The interval is nonempty exactly when t_entry <= t_exit, which
            // (with the entry denominator negative and the exit denominator
            // positive) is the sign of the three-line determinant.
            const auto det = detail::boundaryLinesDeterminant(
                halfplanes_[static_cast<std::size_t>(result.entry)],
                halfplanes_[static_cast<std::size_t>(result.exit)], query);
            if (det > decltype(det){}) {
                result.empty = true;
            }
        }
        return result;
    }
 
    // Compares the clip parameters of the entry and exit constraints strictly:
    // true when t_entry < t_exit (the clipped set has positive length).
    template <HalfplaneConcept Query>
    constexpr bool clipHasLength(const ClipResult& clip, const Query& query) const {
        if (clip.empty) {
            return false;
        }
        if (clip.entry < 0 || clip.exit < 0) {
            return true;
        }
        const auto det = detail::boundaryLinesDeterminant(
            halfplanes_[static_cast<std::size_t>(clip.entry)],
            halfplanes_[static_cast<std::size_t>(clip.exit)], query);
        return det < decltype(det){};
    }
 
    // -1: the point is outside the region, 0: on its boundary, +1: in its
    // interior. O(log n): only the constraints binding along the vertical
    // line through the point (plus the vertical ones) need testing.
    template <PointConcept OtherPoint>
    constexpr int pointStatus(const OtherPoint& point) const {
        if (empty_) {
            return -1;
        }
        if (halfplanes_.empty()) {
            return 1;
        }
        using OtherNumber = typename OtherPoint::NumberType;
        using QueryPoint = Point<OtherNumber>;
        const QueryPoint base(point.x(), point.y());
        const Halfplane<QueryPoint> upward(base, QueryPoint(point.x(), point.y() + OtherNumber(1)));
        int result = 1;
        const std::ptrdiff_t candidates[4] = {
            exitConstraint(upward), exitConstraint(upward.opposite()),
            sameDirectionIndex(upward), sameDirectionIndex(upward.opposite())};
        for (const std::ptrdiff_t idx : candidates) {
            if (idx < 0) {
                continue;
            }
            const auto side = constraintSide(static_cast<std::size_t>(idx), point);
            if (side < 0) {
                return -1;
            }
            if (side == 0) {
                result = 0;
            }
        }
        return result;
    }
 
    // Whether the direction belongs to the region's recession cone: the
    // region contains a translate of every ray with this direction. `strict`
    // additionally requires the interior recession (a·u < 0 for every stored
    // constraint). Precondition: !empty_.
    template <HalfplaneConcept Query>
    constexpr bool recessionContains(const Query& query, bool strict = false) const {
        if (halfplanes_.empty()) {
            return true;
        }
        // a_i·u <= 0 for all i means no stored boundary direction lies
        // strictly inside the open half-circle to the left of u.
        const auto [lo, count] = leftArc(query);
        (void)lo;
        if (count > 0) {
            return false;
        }
        if (!strict) {
            return true;
        }
        return sameDirectionIndex(query) < 0 && sameDirectionIndex(query.opposite()) < 0;
    }
 
    // The recession cone as up to two closed angular arcs, each reported as a
    // pair of half-plane copies whose boundary directions delimit the arc
    // (from, to). A gap of length g >= pi from direction d_i to d_{i+1}
    // contributes the arc [d_i, -d_{i+1}] of length g - pi. O(n).
    // Precondition: !empty_ and size() >= 2 (a single half-plane's recession
    // half-circle and the whole plane are handled by the callers).
    constexpr std::vector<std::pair<HalfplaneType, HalfplaneType>> recessionArcs() const {
        std::vector<std::pair<HalfplaneType, HalfplaneType>> arcs;
        const std::size_t n = halfplanes_.size();
        for (std::size_t i = 0; i < n; ++i) {
            const std::size_t j = nextIndex(i);
            const auto cross = detail::directionCross(halfplanes_[i], halfplanes_[j]);
            // Antiparallel consecutive directions (cross == 0) leave a gap of
            // exactly pi and count as wide too; the same direction is
            // impossible in the canonical form for distinct consecutive
            // constraints, but a single pair wraps twice.
            const bool wide = cross <= 0;
            if (i == j) {
                break;  // n == 1, excluded by precondition but stay safe
            }
            if (wide) {
                arcs.emplace_back(halfplanes_[i], halfplanes_[j].opposite());
            }
        }
        return arcs;
    }
};

template <class PointType_, class TLabel>
constexpr std::optional<HalfplaneIntersection<PointType_>>
Polygon<PointType_, TLabel>::getStarShapedKernel() const {
    using RegionType = HalfplaneIntersection<PointType>;
    if (const auto vertex = getIfPoint()) {
        return RegionType(*vertex);
    }
    if (const auto carrier = getIfSegment()) {
        return RegionType(*carrier);
    }
    if (empty() || isUndefined()) {
        // The empty polygon has no point to see it from, and an undefined one
        // has no interior side for its edges to bound. Either way the kernel is
        // empty, which is reported as no kernel at all.
        return std::nullopt;
    }
    RegionType kernel;
    const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(size());
    for (std::ptrdiff_t i = 0; i < n; ++i) {
        const PointType source = get(i);
        const PointType target = get(i + 1);
        if (source == target) {
            continue;
        }
        kernel.insert(Halfplane<PointType>(source, target));
        if (kernel.empty()) {
            return std::nullopt;
        }
    }
    return kernel;
}

template <class PointType, class LabelType, class TranslationNumber, class TranslationLabel>
constexpr auto operator-(const HalfplaneIntersection<PointType, LabelType>& region,
                         const Point<TranslationNumber, TranslationLabel>& translation) {
    return region + (-translation);
}

template <class PointType, class LabelType, class Scalar>
    requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
constexpr auto operator*(const HalfplaneIntersection<PointType, LabelType>& region, const Scalar& scalar) {
    using ResultPointType = Point<decltype(std::declval<PointType>().x() * scalar), typename PointType::LabelType>;
    HalfplaneIntersection<ResultPointType, LabelType> result(region);
    result *= scalar;
    if constexpr (detail::has_label_v<LabelType>) {
        result.label() = LabelType{};
    }
    return result;
}

template <class Scalar, class PointType, class LabelType>
    requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
constexpr auto operator*(const Scalar& scalar, const HalfplaneIntersection<PointType, LabelType>& region) {
    return region * scalar;
}

template <class PointType, class LabelType, class Scalar>
    requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
constexpr auto operator/(const HalfplaneIntersection<PointType, LabelType>& region, const Scalar& scalar) {
    using ResultPointType = Point<decltype(std::declval<PointType>().x() / scalar), typename PointType::LabelType>;
    HalfplaneIntersection<ResultPointType, LabelType> result(region);
    result /= scalar;
    if constexpr (detail::has_label_v<LabelType>) {
        result.label() = LabelType{};
    }
    return result;
}
 
template <class PointType, class LabelType>
std::ostream& operator<<(std::ostream& stream, const HalfplaneIntersection<PointType, LabelType>& region);
 
}  // namespace pgl