intersection.hpp

Go to the documentation of this file.

#pragma once
 
#include "implementation/atxy.hpp"

 
#include <algorithm>
#include <cassert>
#include <map>
#include <set>
#include <type_traits>
 
 
namespace pgl {
 
namespace detail {

template <class ResultNumber, class ResultLabel, PointConcept APoint, PointConcept BPoint>
constexpr Point<ResultNumber, ResultLabel> carrierCrossing(
    const APoint& a1, const APoint& a2, const BPoint& b1, const BPoint& b2) {
    using InputNumber = std::common_type_t<typename APoint::NumberType,
                                           typename BPoint::NumberType>;
    using Coordinate =
        std::conditional_t<std::floating_point<InputNumber>, ResultNumber,
                           promoted_number_t<promoted_number_t<InputNumber>>>;
 
    const auto wide = [](const auto& value) -> decltype(auto) {
        return detail::asNumber<Coordinate>(value);
    };
    const Coordinate rx = wide(a2.x()) - wide(a1.x());
    const Coordinate ry = wide(a2.y()) - wide(a1.y());
    const Coordinate sx = wide(b2.x()) - wide(b1.x());
    const Coordinate sy = wide(b2.y()) - wide(b1.y());
    const Coordinate ox = wide(b1.x()) - wide(a1.x());
    const Coordinate oy = wide(b1.y()) - wide(a1.y());
 
    const Coordinate determinant = rx * sy - ry * sx;
    const Coordinate along = ox * sy - oy * sx;
 
    // An integral result divides in the wide type, where the division is exact
    // for a crossing that lands on the grid; every other result type divides in
    // itself, which is where its own exactness lives.
    const auto ratio = [&determinant](const Coordinate& numerator) {
        if constexpr (std::integral<ResultNumber>) {
            return static_cast<ResultNumber>(numerator / determinant);
        } else {
            return detail::asNumber<ResultNumber>(numerator) /
                   detail::asNumber<ResultNumber>(determinant);
        }
    };
 
    return Point<ResultNumber, ResultLabel>(
        detail::asNumber<ResultNumber>(a1.x()) + ratio(along * rx),
        detail::asNumber<ResultNumber>(a1.y()) + ratio(along * ry));
}
 
}  // namespace detail
 
// -----------------------------------------------------------------------------
// Point
 
template <class Number, class Label>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, Label>>
Point<Number, Label>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, Label>(*this);
    }
    return {};
}
 
// -----------------------------------------------------------------------------
// Segment
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
Segment<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Segment<PointType, LabelType>::intersection(const OtherSegment& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
 
    if (!boundingBoxesOverlap(other)) {
        return {};
    }
 
    const auto d1 = orientationSign(min(), max(), other.min());
    const auto d2 = orientationSign(min(), max(), other.max());
 
    if (d1 == 0 && d2 == 0) {
        // Both segments are collinear: compare endpoints in ResultPoint so the
        // ternary and ordering stay well-typed when PointType and the other
        // segment's point type
        // differ (e.g. Convex chord clipped from int into Rational).
        const ResultPoint a_min(min());
        const ResultPoint a_max(max());
        const ResultPoint b_min(other.min());
        const ResultPoint b_max(other.max());
        const ResultPoint pminmax = a_max < b_max ? a_max : b_max;
        const ResultPoint pmaxmin = a_min < b_min ? b_min : a_min;
 
        if (pminmax < pmaxmin) {
            return {};
        }
        if (pminmax == pmaxmin) {
            return pminmax;
        }
 
        return ResultSegment(pminmax, pmaxmin);
    }
 
    if (d1 == 0 && containsCollinear(other.min())) {
        return Point<ResultNumber>(other.min());
    }
    if (d2 == 0 && containsCollinear(other.max())) {
        return Point<ResultNumber>(other.max());
    }
 
    const auto d3 = orientationSign(other.min(), other.max(), min());
    if (d3 == 0 && other.containsCollinear(min())) {
        return Point<ResultNumber>(min());
    }
 
    const auto d4 = orientationSign(other.min(), other.max(), max());
    if (d4 == 0 && other.containsCollinear(max())) {
        return Point<ResultNumber>(max());
    }
 
    if (d1 == 0 || d2 == 0 || d3 == 0 || d4 == 0) {
        return {};
    }
 
    if (d1 != d2 && d3 != d4) {
        // The four orientations above have established that the segments cross
        // properly, so their carriers are not parallel.
        return detail::carrierCrossing<ResultNumber, typename PointType::LabelType>(
            min(), max(), other.min(), other.max());
    }
 
    return {};
}
 
// -----------------------------------------------------------------------------
// OrientedSegment
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
OrientedSegment<PointType, LabelType>::intersection(const OtherPoint& other) const {
    return static_cast<Segment<PointType>>(*this).template intersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto OrientedSegment<PointType, LabelType>::intersection(const OtherSegment& other) const {
    return static_cast<Segment<PointType>>(*this).template intersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto OrientedSegment<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return static_cast<Segment<PointType>>(*this).template intersection<ResultNumber>(
        static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
// -----------------------------------------------------------------------------
// Line
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
Line<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Line<Point<ResultNumber, typename PointType::LabelType>>>>
Line<PointType, LabelType>::intersection(const OtherLine& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultLine = Line<ResultPoint>;
 
    if (isDegenerate()) {
        if (other.contains(min())) {
            return ResultPoint(min());
        }
        return {};
    }
 
    if (other.isDegenerate()) {
        if (contains(other.min())) {
            return ResultPoint(other.min());
        }
        return {};
    }
 
    if (contains(other)) {
        return ResultLine(ResultPoint(min()), ResultPoint(max()));
    }
 
    if (parallel(other)) {
        return {};
    }
 
    // Neither line is degenerate and they are not parallel, so they cross once.
    return detail::carrierCrossing<ResultNumber, typename PointType::LabelType>(
        min(), max(), other.min(), other.max());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Line<PointType, LabelType>::intersection(const OtherSegment& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using ResultType = std::optional<std::variant<ResultPoint, ResultSegment>>;
 
    const auto line_intersection = intersection<ResultNumber>(Line<typename OtherSegment::PointType>(other.min(), other.max()));
    if (!line_intersection) {
        return ResultType{};
    }
    if (std::holds_alternative<ResultPoint>(*line_intersection)) {
        const auto& point = std::get<ResultPoint>(*line_intersection);
        if (other.contains(point)) {
            return ResultType(point);
        }
        return ResultType{};
    }
    return ResultType(ResultSegment(ResultPoint(other.min()), ResultPoint(other.max())));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Line<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return intersection<ResultNumber>(static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
// -----------------------------------------------------------------------------
// OrientedLine
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
OrientedLine<PointType, LabelType>::intersection(const OtherPoint& other) const {
    return this->asLine().template intersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Line<Point<ResultNumber, typename PointType::LabelType>>>>
OrientedLine<PointType, LabelType>::intersection(const OtherLine& other) const {
    return this->asLine().template intersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Line<Point<ResultNumber, typename PointType::LabelType>>>>
OrientedLine<PointType, LabelType>::intersection(const OtherOrientedLine& other) const {
    return this->asLine().template intersection<ResultNumber>(
        other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto OrientedLine<PointType, LabelType>::intersection(const OtherSegment& other) const {
    return this->asLine().template intersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto OrientedLine<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return this->asLine().template intersection<ResultNumber>(other);
}
 
// -----------------------------------------------------------------------------
// Ray
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
Ray<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Ray<Point<ResultNumber, typename PointType::LabelType>>>>
Ray<PointType, LabelType>::intersection(const OtherLine& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    if (isDegenerate()) {
        if (other.contains(source())) {
            return ResultPoint(source());
        }
        return {};
    }
 
    const auto line_intersection = this->asLine().template intersection<ResultNumber>(other);
    if (!line_intersection) {
        return {};
    }
 
    if (std::holds_alternative<Line<ResultPoint>>(*line_intersection)) {
        return Ray<ResultPoint>(ResultPoint(source()), ResultPoint(target()));
    }
 
    const auto& point = std::get<ResultPoint>(*line_intersection);
    if (contains(point)) {
        return point;
    }
 
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Ray<Point<ResultNumber, typename PointType::LabelType>>>>
Ray<PointType, LabelType>::intersection(const OtherOrientedLine& other) const {
    return intersection<ResultNumber>(other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Ray<PointType, LabelType>::intersection(const OtherSegment& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    if (isDegenerate()) {
        if (other.contains(source())) {
            return ResultPoint(source());
        }
        return {};
    }
 
    const auto line_intersection =
        this->asLine().template intersection<ResultNumber>(
            Line<typename OtherSegment::PointType>(other.min(), other.max()));
    if (!line_intersection) {
        return {};
    }
 
    if (std::holds_alternative<ResultPoint>(*line_intersection)) {
        const auto& point = std::get<ResultPoint>(*line_intersection);
        if (contains(point) && other.contains(point)) {
            return point;
        }
        return {};
    }
 
    const ResultPoint ray_source(source());
    const bool min_on_ray = contains(other.min());
    const bool max_on_ray = contains(other.max());
 
    if (min_on_ray && max_on_ray) {
        return ResultSegment(ResultPoint(other.min()), ResultPoint(other.max()));
    }
 
    if (min_on_ray) {
        const auto point = ResultPoint(other.min());
        if (other.contains(ray_source) && point != ray_source) {
            return ResultSegment(point, ray_source);
        }
        return point;
    }
 
    if (max_on_ray) {
        const auto point = ResultPoint(other.max());
        if (other.contains(ray_source) && point != ray_source) {
            return ResultSegment(point, ray_source);
        }
        return point;
    }
 
    if (other.contains(ray_source)) {
        return ray_source;
    }
 
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Ray<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return intersection<ResultNumber>(static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Segment<Point<ResultNumber, typename PointType::LabelType>>,
        Ray<Point<ResultNumber, typename PointType::LabelType>>>>
Ray<PointType, LabelType>::intersection(const OtherRay& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using ResultRay = Ray<ResultPoint>;
    if (isDegenerate()) {
        if (other.contains(source())) {
            return ResultPoint(source());
        }
        return {};
    }
 
    if (other.isDegenerate()) {
        if (contains(other.source())) {
            return ResultPoint(other.source());
        }
        return {};
    }
 
    const auto line_intersection =
        this->asLine().template intersection<ResultNumber>(
            other.asLine());
    if (!line_intersection) {
        return {};
    }
 
    if (std::holds_alternative<ResultPoint>(*line_intersection)) {
        const auto& point = std::get<ResultPoint>(*line_intersection);
        if (contains(point) && other.contains(point)) {
            return point;
        }
        return {};
    }
 
    if (source() == other.source()) {
        if (contains(other.target())) {
            return ResultRay(ResultPoint(source()), ResultPoint(target()));
        }
        return ResultPoint(source());
    }
 
    const bool this_contains_other_source = contains(other.source());
    const bool other_contains_this_source = other.contains(source());
 
    if (this_contains_other_source && other_contains_this_source) {
        return ResultSegment(ResultPoint(source()), ResultPoint(other.source()));
    }
 
    if (this_contains_other_source) {
        return ResultRay(ResultPoint(other.source()), ResultPoint(other.target()));
    }
 
    if (other_contains_this_source) {
        return ResultRay(ResultPoint(source()), ResultPoint(target()));
    }
 
    return {};
}
 
// -----------------------------------------------------------------------------
// Halfplane
 
namespace detail {
 
template <class ResultPoint>
constexpr ResultPoint extendRayAlongLine(const ResultPoint& intersection, const ResultPoint& first, const ResultPoint& second) {
    return intersection + (second - first);
}
 
}  // namespace detail
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
Halfplane<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Line<Point<ResultNumber, typename PointType::LabelType>>,
        Ray<Point<ResultNumber, typename PointType::LabelType>>>>
Halfplane<PointType, LabelType>::intersection(const OtherLine& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultLine = Line<ResultPoint>;
    using ResultRay = Ray<ResultPoint>;
 
    if (isDegenerate()) {
        if (other.contains(source())) {
            return ResultPoint(source());
        }
        return {};
    }
 
    if (other.isDegenerate()) {
        if (contains(other.min())) {
            return ResultPoint(other.min());
        }
        return {};
    }
 
    const auto boundary_intersection =
        this->asLine().template intersection<ResultNumber>(other);
 
    if (!boundary_intersection) {
        if (contains(other.min())) {
            return ResultLine(ResultPoint(other.min()), ResultPoint(other.max()));
        }
        return {};
    }
 
    if (std::holds_alternative<ResultPoint>(*boundary_intersection)) {
        const auto& point = std::get<ResultPoint>(*boundary_intersection);
        const ResultPoint first(other.min());
        const ResultPoint second(other.max());
        const auto direction_side =
            orientationDeterminant(source(), target(), other.max()) -
            orientationDeterminant(source(), target(), other.min());
        const auto zero = decltype(direction_side){};
 
        if (zero < direction_side) {
            return ResultRay(point, detail::extendRayAlongLine(point, first, second));
        }
 
        return ResultRay(point, detail::extendRayAlongLine(point, second, first));
    }
 
    return ResultLine(ResultPoint(other.min()), ResultPoint(other.max()));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Line<Point<ResultNumber, typename PointType::LabelType>>,
        Ray<Point<ResultNumber, typename PointType::LabelType>>>>
Halfplane<PointType, LabelType>::intersection(const OtherOrientedLine& other) const {
    return intersection<ResultNumber>(other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Halfplane<PointType, LabelType>::intersection(const OtherSegment& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using ResultLine = Line<ResultPoint>;
    using ResultRay = Ray<ResultPoint>;
 
    if (isDegenerate()) {
        if (other.contains(source())) {
            return ResultPoint(source());
        }
        return {};
    }
 
    if (other.isDegenerate()) {
        if (contains(other.min())) {
            return ResultPoint(other.min());
        }
        return {};
    }
 
    const auto supporting_intersection =
        intersection<ResultNumber>(Line<typename OtherSegment::PointType>(other.min(), other.max()));
    if (!supporting_intersection) {
        return {};
    }
 
    if (std::holds_alternative<ResultPoint>(*supporting_intersection)) {
        const auto& point = std::get<ResultPoint>(*supporting_intersection);
        if (other.contains(point)) {
            return point;
        }
        return {};
    }
 
    if (std::holds_alternative<ResultLine>(*supporting_intersection)) {
        return ResultSegment(ResultPoint(other.min()), ResultPoint(other.max()));
    }
 
    return std::get<ResultRay>(*supporting_intersection).template intersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Halfplane<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return intersection<ResultNumber>(static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr std::optional<
    std::variant<
        Point<ResultNumber, typename PointType::LabelType>,
        Segment<Point<ResultNumber, typename PointType::LabelType>>,
        Ray<Point<ResultNumber, typename PointType::LabelType>>>>
Halfplane<PointType, LabelType>::intersection(const OtherRay& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultRay = Ray<ResultPoint>;
 
    if (isDegenerate()) {
        if (other.contains(source())) {
            return ResultPoint(source());
        }
        return {};
    }
 
    if (other.isDegenerate()) {
        if (contains(other.source())) {
            return ResultPoint(other.source());
        }
        return {};
    }
 
    const auto supporting_intersection =
        intersection<ResultNumber>(other.asLine());
    if (!supporting_intersection) {
        return {};
    }
 
    if (std::holds_alternative<ResultPoint>(*supporting_intersection)) {
        const auto& point = std::get<ResultPoint>(*supporting_intersection);
        if (other.contains(point)) {
            return point;
        }
        return {};
    }
 
    if (std::holds_alternative<Line<ResultPoint>>(*supporting_intersection)) {
        return ResultRay(ResultPoint(other.source()), ResultPoint(other.target()));
    }
 
    return std::get<ResultRay>(*supporting_intersection).template intersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
Halfplane<PointType, LabelType>::intersection(const OtherHalfplane& other) const {
    // Exact and division-free: the region simply stores the two half-planes,
    // discarding one when it is redundant and flagging emptiness when they
    // contradict each other.
    HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>> result(*this);
    result.insert(other);
    return result;
}
 
// -----------------------------------------------------------------------------
// Rectangle
 
namespace detail {
 
template <class PointType, std::size_t Capacity>
struct UniqueIntersectionPoints {
    constexpr void add(const PointType& point) {
        for (std::size_t index = 0; index < size; ++index) {
            if (points[index] == point) {
                return;
            }
        }
 
        points[size] = point;
        ++size;
    }
 
    constexpr const PointType& operator[](std::size_t index) const {
        return points[index];
    }
 
    constexpr PointType minPoint() const {
        PointType result = points[0];
        for (std::size_t index = 1; index < size; ++index) {
            if (points[index] < result) {
                result = points[index];
            }
        }
        return result;
    }
 
    constexpr PointType maxPoint() const {
        PointType result = points[0];
        for (std::size_t index = 1; index < size; ++index) {
            if (result < points[index]) {
                result = points[index];
            }
        }
        return result;
    }
 
    std::array<PointType, Capacity> points{};
    std::size_t size = 0;
};
 
template <class SegmentType, class PointSet>
constexpr void addSegmentEndpoints(const SegmentType& segment, PointSet& points) {
    points.add(segment.min());
    points.add(segment.max());
}
 
template <class PointSet>
constexpr auto pointsToSegmentIntersection(const PointSet& points) {
    using PointType = std::remove_cvref_t<decltype(points[0])>;
    using SegmentType = Segment<PointType>;
    using ResultType = std::optional<std::variant<PointType, SegmentType>>;
 
    if (points.size == 0) {
        return ResultType{};
    }
 
    if (points.size == 1) {
        return ResultType(points[0]);
    }
 
    const auto first = points.minPoint();
    const auto second = points.maxPoint();
    if (first == second) {
        return ResultType(first);
    }
 
    return ResultType(SegmentType(first, second));
}
 
template <class ResultNumber, class LabelType, class RectangleType, class LineType>
constexpr auto rectangleLineIntersection(const RectangleType& rectangle, const LineType& line) {
    using ResultPoint = Point<ResultNumber, LabelType>;
    UniqueIntersectionPoints<ResultPoint, 8> points;
 
    if (line.isDegenerate()) {
        if (rectangle.contains(line.min())) {
            points.add(ResultPoint(line.min()));
        }
        return pointsToSegmentIntersection(points);
    }
 
    const auto rectangle_edges = rectangle.edges();
    for (const auto& edge : rectangle_edges) {
        const Line<typename RectangleType::PointType> edge_line(edge.min(), edge.max());
        const auto intersection = edge_line.template intersection<ResultNumber>(line);
        if (!intersection) {
            continue;
        }
 
        if (std::holds_alternative<ResultPoint>(*intersection)) {
            const auto& point = std::get<ResultPoint>(*intersection);
            if (edge.contains(point)) {
                points.add(point);
            }
        } else {
            const auto& overlap = std::get<Line<ResultPoint>>(*intersection);
            static_cast<void>(overlap);
            addSegmentEndpoints(Segment<ResultPoint>(ResultPoint(edge.min()), ResultPoint(edge.max())), points);
        }
    }
 
    return pointsToSegmentIntersection(points);
}
 
template <class ResultNumber, class LabelType, class RectangleType, class SegmentType>
constexpr auto rectangleSegmentIntersection(const RectangleType& rectangle, const SegmentType& segment) {
    using ResultPoint = Point<ResultNumber, LabelType>;
    UniqueIntersectionPoints<ResultPoint, 10> points;
 
    if (rectangle.contains(segment.min())) {
        points.add(ResultPoint(segment.min()));
    }
    if (rectangle.contains(segment.max())) {
        points.add(ResultPoint(segment.max()));
    }
 
    const auto rectangle_edges = rectangle.edges();
    for (const auto& edge : rectangle_edges) {
        const auto intersection = edge.template intersection<ResultNumber>(segment);
        if (!intersection) {
            continue;
        }
 
        if (std::holds_alternative<ResultPoint>(*intersection)) {
            points.add(std::get<ResultPoint>(*intersection));
        } else {
            addSegmentEndpoints(std::get<Segment<ResultPoint>>(*intersection), points);
        }
    }
 
    return pointsToSegmentIntersection(points);
}
 
template <class ResultNumber, class LabelType, class RectangleType, class RayType>
constexpr auto rectangleRayIntersection(const RectangleType& rectangle, const RayType& ray) {
    using ResultPoint = Point<ResultNumber, LabelType>;
    UniqueIntersectionPoints<ResultPoint, 10> points;
 
    if (rectangle.contains(ray.source())) {
        points.add(ResultPoint(ray.source()));
    }
 
    if (ray.isDegenerate()) {
        return pointsToSegmentIntersection(points);
    }
 
    const Line<Point<ResultNumber, LabelType>> ray_line(ResultPoint(ray.source()), ResultPoint(ray.target()));
    const auto rectangle_edges = rectangle.edges();
    for (const auto& edge : rectangle_edges) {
        const Line<typename RectangleType::PointType> edge_line(edge.min(), edge.max());
        const auto intersection = edge_line.template intersection<ResultNumber>(ray_line);
        if (!intersection) {
            continue;
        }
 
        if (std::holds_alternative<ResultPoint>(*intersection)) {
            const auto& point = std::get<ResultPoint>(*intersection);
            if (edge.contains(point) && ray.contains(point)) {
                points.add(point);
            }
        } else {
            if (ray.contains(edge.min())) {
                points.add(ResultPoint(edge.min()));
            }
            if (ray.contains(edge.max())) {
                points.add(ResultPoint(edge.max()));
            }
        }
    }
 
    return pointsToSegmentIntersection(points);
}
 
template <class ResultNumber, class LabelType, class TriangleType, class LineType>
constexpr auto triangleLineIntersection(const TriangleType& triangle, const LineType& line) {
    using ResultPoint = Point<ResultNumber, LabelType>;
    UniqueIntersectionPoints<ResultPoint, 8> points;
 
    if (line.isDegenerate()) {
        if (triangle.contains(line.min())) {
            points.add(ResultPoint(line.min()));
        }
        return pointsToSegmentIntersection(points);
    }
 
    const auto triangle_edges = triangle.edges();
    for (const auto& edge : triangle_edges) {
        const Line<typename TriangleType::PointType> edge_line(edge.min(), edge.max());
        const auto intersection = edge_line.template intersection<ResultNumber>(line);
        if (!intersection) {
            continue;
        }
 
        if (std::holds_alternative<ResultPoint>(*intersection)) {
            const auto& point = std::get<ResultPoint>(*intersection);
            if (edge.contains(point)) {
                points.add(point);
            }
        } else {
            addSegmentEndpoints(Segment<ResultPoint>(ResultPoint(edge.min()), ResultPoint(edge.max())), points);
        }
    }
 
    return pointsToSegmentIntersection(points);
}
 
template <class ResultNumber, class LabelType, class TriangleType, class SegmentType>
constexpr auto triangleSegmentIntersection(const TriangleType& triangle, const SegmentType& segment) {
    using ResultPoint = Point<ResultNumber, LabelType>;
    UniqueIntersectionPoints<ResultPoint, 10> points;
 
    if (triangle.contains(segment.min())) {
        points.add(ResultPoint(segment.min()));
    }
    if (triangle.contains(segment.max())) {
        points.add(ResultPoint(segment.max()));
    }
 
    const auto triangle_edges = triangle.edges();
    for (const auto& edge : triangle_edges) {
        const auto intersection = edge.template intersection<ResultNumber>(segment);
        if (!intersection) {
            continue;
        }
 
        if (std::holds_alternative<ResultPoint>(*intersection)) {
            points.add(std::get<ResultPoint>(*intersection));
        } else {
            addSegmentEndpoints(std::get<Segment<ResultPoint>>(*intersection), points);
        }
    }
 
    return pointsToSegmentIntersection(points);
}
 
template <class ResultNumber, class LabelType, class TriangleType, class RayType>
constexpr auto triangleRayIntersection(const TriangleType& triangle, const RayType& ray) {
    using ResultPoint = Point<ResultNumber, LabelType>;
    UniqueIntersectionPoints<ResultPoint, 10> points;
 
    if (triangle.contains(ray.source())) {
        points.add(ResultPoint(ray.source()));
    }
 
    if (ray.isDegenerate()) {
        return pointsToSegmentIntersection(points);
    }
 
    const Line<Point<ResultNumber, LabelType>> ray_line(ResultPoint(ray.source()), ResultPoint(ray.target()));
    const auto triangle_edges = triangle.edges();
    for (const auto& edge : triangle_edges) {
        const Line<typename TriangleType::PointType> edge_line(edge.min(), edge.max());
        const auto intersection = edge_line.template intersection<ResultNumber>(ray_line);
        if (!intersection) {
            continue;
        }
 
        if (std::holds_alternative<ResultPoint>(*intersection)) {
            const auto& point = std::get<ResultPoint>(*intersection);
            if (edge.contains(point) && ray.contains(point)) {
                points.add(point);
            }
        } else {
            if (ray.contains(edge.min())) {
                points.add(ResultPoint(edge.min()));
            }
            if (ray.contains(edge.max())) {
                points.add(ResultPoint(edge.max()));
            }
        }
    }
 
    return pointsToSegmentIntersection(points);
}
 
}  // namespace detail
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
Rectangle<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr std::optional<Rectangle<Point<ResultNumber, typename PointType::LabelType>>>
Rectangle<PointType, LabelType>::intersection(const OtherRectangle& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultRectangle = Rectangle<ResultPoint>;
 
    if (!intersects(other)) {
        return {};
    }
 
    const auto min_x = min().x() < other.min().x() ? other.min().x() : min().x();
    const auto min_y = min().y() < other.min().y() ? other.min().y() : min().y();
    const auto max_x = max().x() < other.max().x() ? max().x() : other.max().x();
    const auto max_y = max().y() < other.max().y() ? max().y() : other.max().y();
 
    return ResultRectangle(
        ResultPoint(detail::asNumber<ResultNumber>(min_x), detail::asNumber<ResultNumber>(min_y)),
        ResultPoint(detail::asNumber<ResultNumber>(max_x), detail::asNumber<ResultNumber>(max_y)));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Rectangle<PointType, LabelType>::intersection(const OtherLine& other) const {
    if (empty()) {
        // Nothing meets the empty set.
        return {};
    }
    return detail::rectangleLineIntersection<ResultNumber, typename PointType::LabelType>(*this, other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Rectangle<PointType, LabelType>::intersection(const OtherOrientedLine& other) const {
    return intersection<ResultNumber>(other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Rectangle<PointType, LabelType>::intersection(const OtherSegment& other) const {
    if (empty()) {
        // Nothing meets the empty set.
        return {};
    }
    return detail::rectangleSegmentIntersection<ResultNumber, typename PointType::LabelType>(*this, other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Rectangle<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return intersection<ResultNumber>(static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Rectangle<PointType, LabelType>::intersection(const OtherRay& other) const {
    if (empty()) {
        // Nothing meets the empty set.
        return {};
    }
    return detail::rectangleRayIntersection<ResultNumber, typename PointType::LabelType>(*this, other);
}
 
// A half-plane clips a rectangle to the corners it keeps plus the at most two
// points where the boundary line crosses the rectangle. Walking the four
// corners counterclockwise emits that polygon in order, so no Convex is built
// for the rectangle and the clipped output needs no hull pass, which is what
// delegating to the generic convex clip would cost.
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr auto Rectangle<PointType, LabelType>::intersection(const OtherHalfplane& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using ResultConvex = Convex<ResultPoint>;
    using ResultType = std::optional<std::variant<ResultPoint, ResultSegment, ResultConvex>>;
 
    if (empty()) {
        // Nothing meets the empty set.
        return ResultType{};
    }
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using ResultConvex = Convex<ResultPoint>;
    using ResultType = std::optional<std::variant<ResultPoint, ResultSegment, ResultConvex>>;
 
    const std::array<PointType, 4> corners{min(), bottomRight(), max(), topLeft()};
    using Determinant = decltype(orientationDeterminant(other.source(), other.target(), corners[0]));
    std::array<Determinant, 4> sides{};
    for (std::size_t i = 0; i < 4; ++i) {
        sides[i] = orientationDeterminant(other.source(), other.target(), corners[i]);
    }
 
    // Kept corners and crossings, counterclockwise. A repeated point can only
    // come from a degenerate rectangle, where dropping it is exactly what makes
    // the result collapse to a segment or a point.
    std::array<ResultPoint, 6> clipped{};
    std::size_t count = 0;
    const auto keep = [&clipped, &count](const ResultPoint& point) {
        if (count == 0 || clipped[count - 1] != point) {
            clipped[count++] = point;
        }
    };
 
    const auto boundary = other.asLine();
    for (std::size_t i = 0; i < 4; ++i) {
        const std::size_t next = (i + 1) % 4;
        if (sides[i] >= 0) {
            keep(ResultPoint(corners[i]));
        }
        if (!((sides[i] > 0 && sides[next] < 0) || (sides[i] < 0 && sides[next] > 0))) {
            continue;   // the edge stays on one side; a corner on the line is already kept
        }
        // Edges are axis-parallel and the boundary line strictly crosses this
        // one, so it is not parallel to it and one division gives the crossing,
        // leaving the coordinate along the edge exact.
        if (i % 2 == 0) {   // bottom and top edges are horizontal
            keep(ResultPoint(*boundary.template xAtY<ResultNumber>(corners[i].y()),
                             detail::asNumber<ResultNumber>(corners[i].y())));
        } else {            // right and left edges are vertical
            keep(ResultPoint(detail::asNumber<ResultNumber>(corners[i].x()),
                             *boundary.template yAtX<ResultNumber>(corners[i].x())));
        }
    }
    while (count > 1 && clipped[count - 1] == clipped[0]) {
        --count;   // the walk closed back onto its first point
    }
 
    if (count == 0) {
        return ResultType{};
    }
    if (count == 1) {
        return ResultType(clipped[0]);
    }
    if (count == 2) {
        return ResultType(ResultSegment(clipped[0], clipped[1]));
    }
 
    // A Convex holds its vertices counterclockwise from the lexicographically
    // smallest one, which is where this walk starts only when the corner min()
    // survived the clip.
    std::size_t first = 0;
    for (std::size_t i = 1; i < count; ++i) {
        if (clipped[i] < clipped[first]) {
            first = i;
        }
    }
    std::vector<ResultPoint> vertices;
    vertices.reserve(count);
    for (std::size_t i = 0; i < count; ++i) {
        vertices.push_back(clipped[(first + i) % count]);
    }
    return ResultType(ResultConvex(std::move(vertices), /*trusted=*/true));
}
 
// -----------------------------------------------------------------------------
// Triangle
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
Triangle<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Triangle<PointType, LabelType>::intersection(const OtherLine& other) const {
    return detail::triangleLineIntersection<ResultNumber, typename PointType::LabelType>(*this, other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Triangle<PointType, LabelType>::intersection(const OtherOrientedLine& other) const {
    return intersection<ResultNumber>(other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Triangle<PointType, LabelType>::intersection(const OtherSegment& other) const {
    return detail::triangleSegmentIntersection<ResultNumber, typename PointType::LabelType>(*this, other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Triangle<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return intersection<ResultNumber>(static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Triangle<PointType, LabelType>::intersection(const OtherRay& other) const {
    return detail::triangleRayIntersection<ResultNumber, typename PointType::LabelType>(*this, other);
}
 
// A triangle is convex, so a half-plane (or rectangle) clip is exactly the
// convex-polygon clip; delegating keeps the area overlap (a Convex) instead of
// dropping it, and reuses the exact O(log n + k) routine.
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr auto Triangle<PointType, LabelType>::intersection(const OtherHalfplane& other) const {
    return asConvex().template intersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto Triangle<PointType, LabelType>::intersection(const OtherRectangle& other) const {
    return asConvex().template intersection<ResultNumber>(other);
}
 
// Two triangles are convex, so their intersection is the convex-polygon clip of
// one against the other: an area overlap (Convex), a shared boundary segment, a
// single touch point, or nothing. Delegating to Convex::intersection keeps this
// exact and avoids re-deriving the clip here.
template <class PointType, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr auto Triangle<PointType, LabelType>::intersection(const OtherTriangle& other) const {
    return asConvex().template intersection<ResultNumber>(other.asConvex());
}
 
 
// ---------------------------------------------------------------------------
// Convex
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
Convex<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>> Convex<PointType, LabelType>::intersection(const OtherSegment& other) const {
    auto isec = this->template intersection<ResultNumber>(Line<typename OtherSegment::PointType>(other));
    if (!isec) {
        return {};
    }
    size_t index = isec->index();
    if (index == 0) {
        const auto& p = std::get<0>(*isec);
        if (other.containsCollinear(p)) {
            return p;
        }
        return {};
    }
    const auto& seg = std::get<1>(*isec);
    return seg.template intersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>> Convex<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return intersection<ResultNumber>(Segment<typename OtherOrientedSegment::PointType>(other[0], other[1]));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>> Convex<PointType, LabelType>::intersection(const OtherLine& other) const {
    if (points_.empty()) {
        return {};
    }
 
    using CommonNumberType = std::common_type_t<NumberType, typename OtherLine::NumberType>;
    using CommonPoint = Point<CommonNumberType, typename PointType::LabelType>;
 
    const CommonPoint local_p0(detail::asNumber<CommonNumberType>(other[0].x()) - detail::asNumber<CommonNumberType>(translation_.x()),
                               detail::asNumber<CommonNumberType>(other[0].y()) - detail::asNumber<CommonNumberType>(translation_.y()));
    const CommonPoint local_p1(detail::asNumber<CommonNumberType>(other[1].x()) - detail::asNumber<CommonNumberType>(translation_.x()),
                               detail::asNumber<CommonNumberType>(other[1].y()) - detail::asNumber<CommonNumberType>(translation_.y()));
 
    if (local_p0 == local_p1) {
        return {};
    }
 
    const Line<CommonPoint> localLine(local_p0, local_p1);
 
    auto orientationValue = [&](size_t index) {
        return orientationDeterminant(localLine[0], localLine[1], points_[index]);
    };
 
    if (points_.size() == 1) {
        if (orientationValue(0) == CommonNumberType(0)) {
            Point<ResultNumber, typename PointType::LabelType> result =
                static_cast<Point<ResultNumber, typename PointType::LabelType>>(points_[0]);
            return result + static_cast<Point<ResultNumber, typename PointType::LabelType>>(translation_);
        }
        return {};
    }
 
    if (points_.size() == 2) {
        auto local_segment = Segment<PointType>(points_[0], points_[1]);
        auto isec = local_segment.template intersection<ResultNumber>(localLine);
        if (!isec) {
            return {};
        }
        const auto shift = static_cast<Point<ResultNumber, typename PointType::LabelType>>(translation_);
        if (std::holds_alternative<Point<ResultNumber, typename PointType::LabelType>>(*isec)) {
            return std::get<Point<ResultNumber, typename PointType::LabelType>>(*isec) + shift;
        }
        auto seg = std::get<Segment<Point<ResultNumber, typename PointType::LabelType>>>(*isec);
        return Segment<Point<ResultNumber, typename PointType::LabelType>>(seg[0] + shift, seg[1] + shift);
    }
 
    auto max_it = detail::cyclicMax(points_.begin(), points_.end(),
                                    [&](const PointType& a) {
                                        return orientationDeterminant(localLine[0], localLine[1], a);
                                    });
    auto min_it = detail::cyclicMax(points_.begin(), points_.end(),
                                    [&](const PointType& a) {
                                        return orientationDeterminant(localLine[1], localLine[0], a);
                                    });
 
    const size_t n = points_.size();
    const size_t i_max = static_cast<size_t>(std::distance(points_.begin(), max_it));
    const size_t i_min = static_cast<size_t>(std::distance(points_.begin(), min_it));
 
    const CommonNumberType max_value = orientationValue(i_max);
    const CommonNumberType min_value = orientationValue(i_min);
 
    if (max_value < CommonNumberType(0) || min_value > CommonNumberType(0)) {
        return {};
    }
 
    auto translatePoint = [&](const Point<ResultNumber, typename PointType::LabelType>& point) {
        return point + static_cast<Point<ResultNumber, typename PointType::LabelType>>(translation_);
    };
 
    auto translateSegment = [&](const Segment<Point<ResultNumber, typename PointType::LabelType>>& segment) {
        return Segment<Point<ResultNumber, typename PointType::LabelType>>(
            translatePoint(segment[0]), translatePoint(segment[1]));
    };
 
    if (max_value == CommonNumberType(0) && min_value == CommonNumberType(0)) {
        const CommonPoint direction(localLine[1].x() - localLine[0].x(),
                                    localLine[1].y() - localLine[0].y());
 
        auto extremal_high = detail::cyclicMax(points_.begin(), points_.end(),
            [&](const PointType& a) {
                return detail::asNumber<CommonNumberType>(a.x()) * direction.x() +
                       detail::asNumber<CommonNumberType>(a.y()) * direction.y();
            });
 
        auto extremal_low = detail::cyclicMax(points_.begin(), points_.end(),
            [&](const PointType& a) {
                return -(detail::asNumber<CommonNumberType>(a.x()) * direction.x() +
                         detail::asNumber<CommonNumberType>(a.y()) * direction.y());
            });
 
        Segment<PointType> support(*extremal_low, *extremal_high);
        auto isec = support.template intersection<ResultNumber>(localLine);
        if (!isec) {
            return {};
        }
        if (std::holds_alternative<Point<ResultNumber, typename PointType::LabelType>>(*isec)) {
            return translatePoint(std::get<Point<ResultNumber, typename PointType::LabelType>>(*isec));
        }
        return translateSegment(std::get<Segment<Point<ResultNumber, typename PointType::LabelType>>>(*isec));
    }
 
    auto forwardDist = [&](size_t from, size_t to) {
        return (to + n - from) % n;
    };
 
    auto findBoundaryForward = [&](size_t start, size_t length) -> size_t {
        size_t lo = 1;
        size_t hi = length;
        size_t result = length + 1;
        while (lo <= hi) {
            size_t mid = lo + (hi - lo) / 2;
            if (orientationValue((start + mid) % n) <= CommonNumberType(0)) {
                result = mid;
                hi = mid - 1;
            } else {
                lo = mid + 1;
            }
        }
        return result;
    };
 
    auto findBoundaryBackward = [&](size_t start, size_t length) -> size_t {
        size_t lo = 1;
        size_t hi = length;
        size_t result = length + 1;
        while (lo <= hi) {
            size_t mid = lo + (hi - lo) / 2;
            size_t index = (start + n - mid) % n;
            if (orientationValue(index) <= CommonNumberType(0)) {
                result = mid;
                hi = mid - 1;
            } else {
                lo = mid + 1;
            }
        }
        return result;
    };
 
    const size_t forward_len = forwardDist(i_max, i_min);
    const size_t backward_len = forwardDist(i_min, i_max);
 
    const size_t boundary_forward = findBoundaryForward(i_max, forward_len);
    const size_t boundary_backward = findBoundaryBackward(i_max, backward_len);
 
    if (boundary_forward > forward_len || boundary_backward > backward_len) {
        return {};
    }
 
    const size_t b1 = (i_max + boundary_forward) % n;
    const size_t b1_prev = (b1 + n - 1) % n;
    const size_t b2 = (i_max + n - boundary_backward) % n;
    const size_t b2_next = (b2 + 1) % n;
 
    auto collectPoints = [&](const Segment<PointType>& edge,
                             std::vector<Point<ResultNumber, typename PointType::LabelType>>& points) {
        auto edge_isec = edge.template intersection<ResultNumber>(localLine);
        if (!edge_isec) {
            return;
        }
        if (std::holds_alternative<Point<ResultNumber, typename PointType::LabelType>>(*edge_isec)) {
            points.push_back(std::get<Point<ResultNumber, typename PointType::LabelType>>(*edge_isec));
        } else {
            auto seg = std::get<Segment<Point<ResultNumber, typename PointType::LabelType>>>(*edge_isec);
            points.push_back(seg[0]);
            points.push_back(seg[1]);
        }
    };
 
    std::vector<Point<ResultNumber, typename PointType::LabelType>> points;
    points.reserve(4);
    collectPoints(Segment<PointType>(points_[b1_prev], points_[b1]), points);
    collectPoints(Segment<PointType>(points_[b2], points_[b2_next]), points);
 
    if (points.empty()) {
        return {};
    }
 
    const CommonPoint direction(localLine[1].x() - localLine[0].x(),
                                localLine[1].y() - localLine[0].y());
 
    auto project = [&](const Point<ResultNumber, typename PointType::LabelType>& point) {
        return (detail::asNumber<CommonNumberType>(point.x()) - localLine[0].x()) * direction.x() +
               (detail::asNumber<CommonNumberType>(point.y()) - localLine[0].y()) * direction.y();
    };
 
    std::sort(points.begin(), points.end(), [&](auto const& a, auto const& b) {
        return project(a) < project(b);
    });
    points.erase(std::unique(points.begin(), points.end()), points.end());
 
    if (points.empty()) {
        return {};
    }
    if (points.size() == 1) {
        return translatePoint(points[0]);
    }
 
    return Segment<Point<ResultNumber, typename PointType::LabelType>>(translatePoint(points.front()),
                                                                      translatePoint(points.back()));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>> Convex<PointType, LabelType>::intersection(const OtherOrientedLine& other) const {
    return intersection<ResultNumber>(Line<typename OtherOrientedLine::PointType>(other[0], other[1]));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>> Convex<PointType, LabelType>::intersection(const OtherRay& other) const {
    auto line_isec =  intersection<ResultNumber>(Line<typename OtherRay::PointType>(other[0], other[1]));
    if (!line_isec) {
        return {};
    }
    size_t index = line_isec->index();
    if (index == 0) {
        if (other.containsCollinear(std::get<0>(*line_isec))) {
            return std::get<0>(*line_isec);
        }
        return {};
    }
    // index == 1
    auto seg = std::get<1>(*line_isec);
    return seg.template intersection<ResultNumber>(other);
}
 
// Clip the convex polygon to a closed half-plane. The vertices inside the
// half-plane form one contiguous arc (orientation is unimodal along a convex
// hull), so finding one inside vertex and walking outward to the two boundary
// crossings yields the clipped region in O(log n + k) for an output of size k.
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>, Convex<Point<ResultNumber, typename PointType::LabelType>>>> Convex<PointType, LabelType>::intersection(const OtherHalfplane& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using ResultConvex = Convex<ResultPoint>;
 
    const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(points_.size());
    if (n == 0) {
        return {};
    }
 
    // Any vertex inside the closed half-plane anchors the inside arc; the deepest
    // one is found in O(log n) and is positive when one is strictly inside.
    const auto it = detail::cyclicMaxOrPositive(points_.begin(), points_.end(),
        [&other, this](const PointType& a) {
            return orientationDeterminant(other.source(), other.target(), a + translation_);
        });
    const std::ptrdiff_t start = it - points_.begin();
 
    if (!other.contains(get(start))) {
        return {};   // every vertex is strictly outside the half-plane
    }
 
    // Walk CCW (forward) over the contiguous arc of vertices in the half-plane;
    // get() indexes the vertices cyclically, so f and b may run past the ends.
    std::ptrdiff_t f = start, fSteps = 0;
    while (fSteps + 1 < n && other.contains(get(f + 1))) {
        ++f;
        ++fSteps;
    }
    if (fSteps + 1 == n) {
        // The whole convex lies in the closed half-plane; return it unchanged.
        ResultConvex whole(*this, /*trusted=*/true);
        if (whole.size() == 1) return whole[0];
        if (whole.size() == 2) return ResultSegment(whole[0], whole[1]);
        return whole;
    }
 
    // Walk CW (backward) over the same arc.
    std::ptrdiff_t b = start, bSteps = 0;
    while (bSteps + fSteps + 1 < n && other.contains(get(b - 1))) {
        --b;
        ++bSteps;
    }
 
    // The arc spans b .. f (CCW); the edges leaving it at each end cross the
    // boundary line. A crossing edge has one endpoint inside and one strictly
    // outside, so it meets the boundary line at exactly one point.
    auto crossing = [&](std::ptrdiff_t insideIdx, std::ptrdiff_t outsideIdx) {
        const ResultSegment edge(static_cast<ResultPoint>(get(insideIdx)),
                                 static_cast<ResultPoint>(get(outsideIdx)));
        return std::get<ResultPoint>(*edge.template intersection<ResultNumber>(other.asLine()));
    };
 
    std::vector<ResultPoint> result;
    result.reserve(static_cast<std::size_t>(fSteps + bSteps + 3));
    result.push_back(crossing(b, b - 1));   // entering crossing
    for (std::ptrdiff_t i = b; ; ++i) {
        result.push_back(static_cast<ResultPoint>(get(i)));
        if (i == f) break;
    }
    result.push_back(crossing(f, f + 1));   // leaving crossing
 
    ResultConvex convex(std::move(result));
    if (convex.size() == 1) return convex[0];
    if (convex.size() == 2) return ResultSegment(convex[0], convex[1]);
    return convex;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>, Convex<Point<ResultNumber, typename PointType::LabelType>>>> Convex<PointType, LabelType>::intersection(const OtherRectangle& other) const {
    return intersection<ResultNumber>(other.asConvex());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>, Convex<Point<ResultNumber, typename PointType::LabelType>>>> Convex<PointType, LabelType>::intersection(const OtherTriangle& other) const {
    return intersection<ResultNumber>(other.asConvex());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>, Convex<Point<ResultNumber, typename PointType::LabelType>>>> Convex<PointType, LabelType>::intersection(const OtherConvex& other) const {
    if (size() == 0 || other.size() == 0 || !intersects(other)) {
        return {};
    }
 
    using ResultPointType = Point<ResultNumber, typename PointType::LabelType>;
 
    if (size() == 1) {
        ResultPointType point = static_cast<ResultPointType>(points_[0] + translation_);
        if (other.contains(point)) {
            return point;
        }
        return {};
    }
 
    if (other.size() == 1) {
        ResultPointType point = static_cast<ResultPointType>(other[0]);
        if (contains(point)) {
            return point;
        }
        return {};
    }
 
    std::vector<ResultPointType> isecPoints;
 
    for (auto &edge : other.edges()) {
        auto isec = intersection<ResultNumber>(edge);
        if (isec) {
            if (std::holds_alternative<Point<ResultNumber, typename PointType::LabelType>>(*isec)) {
                isecPoints.push_back(std::get<Point<ResultNumber, typename PointType::LabelType>>(*isec));
            } else {
                auto seg = std::get<Segment<Point<ResultNumber, typename PointType::LabelType>>>(*isec);
                isecPoints.push_back(seg[0]);
                isecPoints.push_back(seg[1]);
            }
        }
    }
    for (auto &edge : edges()) {
        auto isec = other.template intersection<ResultNumber>(edge);
        if (isec) {
            if (std::holds_alternative<ResultPointType>(*isec)) {
                isecPoints.push_back(std::get<ResultPointType>(*isec));
            } else {
                auto seg = std::get<Segment<ResultPointType>>(*isec);
                isecPoints.push_back(seg[0]);
                isecPoints.push_back(seg[1]);
            }
        }
    }
 
    if (isecPoints.empty()) {
        return {};
    }
 
    Convex<ResultPointType> convex(std::move(isecPoints));
 
    if (convex.size() == 1) {
        return convex[0];
    }
    if (convex.size() == 2) {
        return Segment<ResultPointType>(convex[0], convex[1]);
    }
 
    return convex;
}
 
// ---------------------------------------------------------------------------
// Clipping a one-dimensional shape to an area bounded by rings
//
// A polygon and a region differ here only in how many rings bound them, so the
// walk below takes the boundary edges and nothing else: the even-odd rule reads
// a ring system exactly as it reads a single ring, and a point inside a hole is
// inside two rings and hence outside the area. The three shape-specific
// wrappers add only the parameter window their shape occupies on its supporting
// line, which leaves Polygon and PolygonWithHoles a one-line overload each.
 
namespace detail {

template <class ResultPoint>
using LinePieces = std::vector<std::variant<ResultPoint, Segment<ResultPoint>>>;

template <class ResultNumber, class ResultPoint>
struct LineSpan {
    ResultNumber lo, hi;
    ResultPoint plo, phi;
};

template <class ResultNumber, class ResultPoint, class EdgeRange>
constexpr std::vector<LineSpan<ResultNumber, ResultPoint>>
lineAreaSpans(const ResultPoint& a, const ResultPoint& b, const EdgeRange& edges) {
    using Span = LineSpan<ResultNumber, ResultPoint>;
 
    const ResultPoint direction = b - a;
    const Line<ResultPoint> line(a, b);
    auto tOf = [&](const ResultPoint& p) -> ResultNumber { return (p - a) * direction; };
 
    // Signed side of a vertex w.r.t. the line (CCW positive, 0 on the line).
    auto sideOf = [&](const ResultPoint& v) -> int {
        const auto o = orientationSign(a, b, v);
        if (o > 0) return 1;
        if (o < 0) return -1;
        return 0;
    };
 
    auto makeSpan = [](ResultNumber t1, ResultNumber t2, const ResultPoint& p1, const ResultPoint& p2) -> Span {
        return (t1 <= t2) ? Span{t1, t2, p1, p2} : Span{t2, t1, p2, p1};
    };
    std::vector<Span> spans;
 
    // Boundary crossings of the line, used to recover inside/outside by ray
    // parity. A vertex on the line is treated as if perturbed to the -1 side,
    // which counts vertex touches and collinear edges consistently: each edge
    // whose perturbed endpoint signs differ contributes one crossing. Each
    // vertex is the source of exactly one directed edge of its ring, so which
    // way the ring is wound changes nothing here — and hole rings are wound the
    // other way from the outer one.
    std::vector<std::pair<ResultNumber, ResultPoint>> crossings;
 
    for (const auto& edge : edges) {
        using EdgePoint = typename std::remove_cvref_t<decltype(edge)>::PointType;
        const ResultPoint u = static_cast<ResultPoint>(edge.source());
        const ResultPoint w = static_cast<ResultPoint>(edge.target());
        const int su = sideOf(u);
        const int sw = sideOf(w);
 
        if (su == 0 && sw == 0) {
            // Edge collinear with the line: a boundary overlap, always closed.
            spans.push_back(makeSpan(tOf(u), tOf(w), u, w));
            continue;
        }
        if (su == 0) {
            // Vertex on the line: a boundary touch point (recorded once, here,
            // as the source vertex of its edge).
            spans.push_back({tOf(u), tOf(u), u, u});
        }
 
        if (su != 0 && sw != 0 && su != sw) {
            // Transversal crossing through the edge interior.
            const auto isec = line.template intersection<ResultNumber>(
                Segment<EdgePoint>(edge.source(), edge.target()));
            if (isec && std::holds_alternative<ResultPoint>(*isec)) {
                const ResultPoint c = std::get<ResultPoint>(*isec);
                spans.push_back({tOf(c), tOf(c), c, c});
                crossings.emplace_back(tOf(c), c);
            }
        } else {
            // Perturbed crossing located at whichever endpoint is on the line.
            const int eu = (su != 0) ? su : -1;
            const int ew = (sw != 0) ? sw : -1;
            if (eu != ew) {
                const ResultPoint& onLine = (su == 0) ? u : w;
                crossings.emplace_back(tOf(onLine), onLine);
            }
        }
    }
 
    // Walk the crossings in order; the open cell to the right of a crossing is
    // inside the area exactly when an odd number of crossings lie to its left.
    // Each maximal inside cell becomes a closed interval (its endpoints are
    // boundary crossings, hence in the closed area too).
    std::sort(crossings.begin(), crossings.end(),
              [](const auto& x, const auto& y) { return x.first < y.first; });
    std::size_t idx = 0;
    int parity = 0;
    while (idx < crossings.size()) {
        const ResultNumber tcur = crossings[idx].first;
        const ResultPoint pcur = crossings[idx].second;
        std::size_t next = idx;
        while (next < crossings.size() && crossings[next].first == tcur) {
            ++next;
        }
        parity += static_cast<int>(next - idx);
        if ((parity & 1) && next < crossings.size()) {
            spans.push_back({tcur, crossings[next].first, pcur, crossings[next].second});
        }
        idx = next;
    }
 
    // Union the spans: sort by lo, then merge touching or overlapping ones.
    std::sort(spans.begin(), spans.end(),
              [](const Span& x, const Span& y) { return x.lo != y.lo ? x.lo < y.lo : x.hi < y.hi; });
    std::vector<Span> merged;
    for (const Span& s : spans) {
        if (merged.empty() || s.lo > merged.back().hi) {
            merged.push_back(s);
        } else if (s.hi > merged.back().hi) {
            merged.back().hi = s.hi;
            merged.back().phi = s.phi;
        }
    }
    return merged;
}

template <class ResultPoint, class Area, class EdgeRange, class OtherSegment>
constexpr LinePieces<ResultPoint>
areaSegmentIntersection(const Area& area, const EdgeRange& edges, const OtherSegment& other) {
    using ResultNumber = typename ResultPoint::NumberType;
    using ResultSegment = Segment<ResultPoint>;
 
    LinePieces<ResultPoint> result;
    const ResultPoint a(other.min());
    const ResultPoint b(other.max());
 
    // A degenerate segment is just its single point.
    if (other.isDegenerate()) {
        if (area.contains(a)) {
            result.emplace_back(a);
        }
        return result;
    }
 
    // The segment is the parameter window [0, T] of its supporting line.
    const ResultNumber T = (b - a) * (b - a);
    for (const auto& span : lineAreaSpans<ResultNumber>(a, b, edges)) {
        if (span.hi < ResultNumber(0) || span.lo > T) {
            continue;
        }
        const ResultPoint lo = (span.lo < ResultNumber(0)) ? a : span.plo;
        const ResultPoint hi = (span.hi > T) ? b : span.phi;
        if (lo == hi) {
            result.emplace_back(lo);
        } else {
            result.emplace_back(ResultSegment(lo, hi));
        }
    }
    return result;
}

template <class ResultPoint, class Area, class EdgeRange, class OtherLine>
constexpr LinePieces<ResultPoint>
areaLineIntersection(const Area& area, const EdgeRange& edges, const OtherLine& other) {
    using ResultNumber = typename ResultPoint::NumberType;
    using ResultSegment = Segment<ResultPoint>;
 
    LinePieces<ResultPoint> result;
    const ResultPoint a(other.min());
    const ResultPoint b(other.max());
 
    // A degenerate line is a single point.
    if (other.isDegenerate()) {
        if (area.contains(a)) {
            result.emplace_back(a);
        }
        return result;
    }
 
    // The line is its whole parametrization, so nothing is clipped away.
    for (const auto& span : lineAreaSpans<ResultNumber>(a, b, edges)) {
        if (span.lo == span.hi) {
            result.emplace_back(span.plo);
        } else {
            result.emplace_back(ResultSegment(span.plo, span.phi));
        }
    }
    return result;
}

template <class ResultPoint, class Area, class EdgeRange, class OtherRay>
constexpr LinePieces<ResultPoint>
areaRayIntersection(const Area& area, const EdgeRange& edges, const OtherRay& other) {
    using ResultNumber = typename ResultPoint::NumberType;
    using ResultSegment = Segment<ResultPoint>;
 
    LinePieces<ResultPoint> result;
    const ResultPoint a(other.source());
    const ResultPoint b(other.target());
 
    // A degenerate ray is a single point.
    if (other.isDegenerate()) {
        if (area.contains(a)) {
            result.emplace_back(a);
        }
        return result;
    }
 
    // The source is the parameter origin, so the ray is the window t >= 0.
    for (const auto& span : lineAreaSpans<ResultNumber>(a, b, edges)) {
        if (span.hi < ResultNumber(0)) {
            continue;
        }
        const ResultPoint lo = (span.lo < ResultNumber(0)) ? a : span.plo;
        if (lo == span.phi) {
            result.emplace_back(lo);
        } else {
            result.emplace_back(ResultSegment(lo, span.phi));
        }
    }
    return result;
}
 
}  // namespace detail
 
// ---------------------------------------------------------------------------
// Polygon
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
Polygon<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polygon<PointType, LabelType>::intersection(const OtherSegment& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    return detail::areaSegmentIntersection<ResultPoint>(*this, orientedEdgesView(), other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polygon<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return intersection<ResultNumber>(Segment<typename OtherOrientedSegment::PointType>(other[0], other[1]));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polygon<PointType, LabelType>::intersection(const OtherLine& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    return detail::areaLineIntersection<ResultPoint>(*this, orientedEdgesView(), other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polygon<PointType, LabelType>::intersection(const OtherOrientedLine& other) const {
    return intersection<ResultNumber>(other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polygon<PointType, LabelType>::intersection(const OtherRay& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    return detail::areaRayIntersection<ResultPoint>(*this, orientedEdgesView(), other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PolygonConcept OtherPolygon>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>, Polyline<Point<ResultNumber, typename PointType::LabelType>>, Polygon<Point<ResultNumber, typename PointType::LabelType>>>>
Polygon<PointType, LabelType>::intersection(const OtherPolygon& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using ResultPolyline = Polyline<ResultPoint>;
    using ResultPolygon = Polygon<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultPolyline, ResultPolygon>;
 
    // The boundary of A ∩ B is (∂A ∩ B) ∪ (∂B ∩ A): clip every edge of each
    // polygon against the other and collect the resulting boundary pieces. The
    // sets deduplicate shared boundary so each segment/point is stored once.
    std::set<ResultSegment> segments;
    std::set<ResultPoint> touchPoints;
 
    auto clipEdgesAgainst = [&](const auto& edgePolygon, const auto& clipPolygon) {
        for (const auto& edge : edgePolygon.edges()) {
            for (const auto& piece : clipPolygon.template intersection<ResultNumber>(edge)) {
                if (std::holds_alternative<ResultPoint>(piece)) {
                    touchPoints.insert(std::get<ResultPoint>(piece));
                } else {
                    segments.insert(std::get<ResultSegment>(piece));
                }
            }
        }
    };
    clipEdgesAgainst(*this, other);
    clipEdgesAgainst(other, *this);
 
    // Make the collected boundary a consistent planar graph: split any segment at
    // another segment's endpoint lying in its interior, then dedup, so partially
    // overlapping collinear (shared) boundary collapses onto common edges. Only
    // genuine pinch points keep degree > 2 afterwards.
    auto planarize = [](const std::set<ResultSegment>& input) {
        std::vector<ResultPoint> vertices;
        for (const auto& s : input) {
            vertices.push_back(s.min());
            vertices.push_back(s.max());
        }
        std::sort(vertices.begin(), vertices.end());
        vertices.erase(std::unique(vertices.begin(), vertices.end()), vertices.end());
 
        std::set<ResultSegment> out;
        for (const auto& s : input) {
            const ResultPoint a = s.min();
            const ResultPoint b = s.max();
            std::vector<ResultPoint> cut{a, b};
            for (const auto& v : vertices) {
                if (v == a || v == b || !collinear(a, b, v)) {
                    continue;
                }
                if (v.x() < std::min(a.x(), b.x()) || v.x() > std::max(a.x(), b.x()) ||
                    v.y() < std::min(a.y(), b.y()) || v.y() > std::max(a.y(), b.y())) {
                    continue;  // collinear but outside the segment
                }
                cut.push_back(v);
            }
            std::sort(cut.begin(), cut.end());
            cut.erase(std::unique(cut.begin(), cut.end()), cut.end());
            for (std::size_t i = 0; i + 1 < cut.size(); ++i) {
                out.insert(ResultSegment(cut[i], cut[i + 1]));
            }
        }
        return out;
    };
    segments = planarize(segments);
 
    // Build the undirected graph: nodes are endpoints, edges are the segments.
    std::map<ResultPoint, std::vector<ResultPoint>> adjacency;
    for (const auto& segment : segments) {
        adjacency[segment.min()].push_back(segment.max());
        adjacency[segment.max()].push_back(segment.min());
    }
 
    std::vector<Piece> result;
 
    // A touch point that no segment uses is an isolated intersection point.
    for (const auto& point : touchPoints) {
        if (adjacency.find(point) == adjacency.end()) {
            result.emplace_back(point);
        }
    }
 
    using Number = ResultNumber;
    auto degree = [&adjacency](const ResultPoint& p) {
        const auto it = adjacency.find(p);
        return it == adjacency.end() ? std::size_t(0) : it->second.size();
    };
    auto removeEdge = [&adjacency](const ResultPoint& u, const ResultPoint& v) {
        auto& au = adjacency[u];
        au.erase(std::find(au.begin(), au.end(), v));
        auto& av = adjacency[v];
        av.erase(std::find(av.begin(), av.end(), u));
    };
 
    // Peel the open boundary (degree-1 chains, and whiskers hanging off a cycle)
    // into polylines, removing their edges. What remains is a union of closed
    // cycles that may share articulation (pinch) vertices.
    while (true) {
        const ResultPoint* leaf = nullptr;
        for (const auto& [p, neighbors] : adjacency) {
            if (neighbors.size() == 1) {
                leaf = &p;
                break;
            }
        }
        if (!leaf) {
            break;
        }
        std::vector<ResultPoint> path;
        ResultPoint current = *leaf;
        path.push_back(current);
        while (degree(current) == 1) {
            const ResultPoint next = adjacency.at(current)[0];
            removeEdge(current, next);
            path.push_back(next);
            current = next;
        }
        result.emplace_back(ResultPolyline(std::move(path)));
    }
 
    // Trace the remaining cycles as faces of the planar graph: following, at each
    // vertex, the neighbour immediately clockwise from the edge we arrived on
    // visits every face once and hugs its interior. Counterclockwise faces are
    // the filled intersection pieces; the clockwise outer face is dropped. This
    // resolves the articulation (degree > 2) vertices into their separate cycles.
    auto rotationalNext = [&adjacency](const ResultPoint& at, const ResultPoint& from) {
        const auto& neighbors = adjacency.at(at);
        auto dx = [&](const ResultPoint& p) { return p.x() - at.x(); };
        auto dy = [&](const ResultPoint& p) { return p.y() - at.y(); };
        // 0 for directions in [0, 180) degrees, 1 for [180, 360); orders the circle.
        auto half = [&](const ResultPoint& p) {
            const Number y = dy(p);
            if (y > Number(0)) return 0;
            if (y < Number(0)) return 1;
            return dx(p) >= Number(0) ? 0 : 1;
        };
        auto ccwBefore = [&](const ResultPoint& u, const ResultPoint& w) {  // u strictly CCW-before w
            const int hu = half(u), hw = half(w);
            if (hu != hw) {
                return hu < hw;
            }
            return dx(u) * dy(w) - dy(u) * dx(w) > Number(0);
        };
        // The neighbour immediately clockwise from `from`: the CCW-largest one
        // that is CCW-before `from`, wrapping to the CCW-largest overall.
        const ResultPoint* best = nullptr;
        for (const auto& n : neighbors) {
            if (n == from || !ccwBefore(n, from)) {
                continue;
            }
            if (!best || ccwBefore(*best, n)) {
                best = &n;
            }
        }
        if (!best) {
            for (const auto& n : neighbors) {
                if (n != from && (!best || ccwBefore(*best, n))) {
                    best = &n;
                }
            }
        }
        return *best;
    };
 
    std::set<std::pair<ResultPoint, ResultPoint>> usedDart;
    for (const auto& [u, neighbors] : adjacency) {
        for (const auto& v : neighbors) {
            if (usedDart.count({u, v})) {
                continue;
            }
            std::vector<ResultPoint> cycle;
            ResultPoint a = u;
            ResultPoint b = v;
            do {
                usedDart.insert({a, b});
                cycle.push_back(a);
                const ResultPoint c = rotationalNext(b, a);
                a = b;
                b = c;
            } while (a != u || b != v);
 
            // Shoelace twice-area; keep counterclockwise (filled) faces only.
            Number twiceArea(0);
            for (std::size_t i = 0; i < cycle.size(); ++i) {
                const ResultPoint& p = cycle[i];
                const ResultPoint& q = cycle[(i + 1) % cycle.size()];
                twiceArea += p.x() * q.y() - q.x() * p.y();
            }
            if (twiceArea > Number(0)) {
                result.emplace_back(ResultPolygon(cycle));
            }
        }
    }
 
    return result;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>, Polyline<Point<ResultNumber, typename PointType::LabelType>>, Polygon<Point<ResultNumber, typename PointType::LabelType>>>>
Polygon<PointType, LabelType>::intersection(const OtherConvex& other) const {
    return this->template intersection<ResultNumber>(other.asPolygon());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>, Polyline<Point<ResultNumber, typename PointType::LabelType>>, Polygon<Point<ResultNumber, typename PointType::LabelType>>>>
Polygon<PointType, LabelType>::intersection(const OtherTriangle& other) const {
    return this->template intersection<ResultNumber>(other.asConvex());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>, Polyline<Point<ResultNumber, typename PointType::LabelType>>, Polygon<Point<ResultNumber, typename PointType::LabelType>>>>
Polygon<PointType, LabelType>::intersection(const OtherRectangle& other) const {
    return this->template intersection<ResultNumber>(other.asConvex());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>, Segment<Point<ResultNumber, typename PointType::LabelType>>, Polygon<Point<ResultNumber, typename PointType::LabelType>>>>
Polygon<PointType, LabelType>::intersection(const OtherHalfplane& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using ResultPolygon = Polygon<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultSegment, ResultPolygon>;
 
    std::vector<Piece> result;
    if (size() == 0) {
        return result;
    }
 
    // A degenerate half-plane is a single point.
    if (other.isDegenerate()) {
        const ResultPoint p(other.source());
        if (contains(p)) {
            result.emplace_back(p);
        }
        return result;
    }
 
    // The boundary of P ∩ H is (∂P ∩ H) ∪ (∂H ∩ P): clip each polygon edge to
    // the closed half-plane, and clip the half-plane's boundary line to the
    // polygon. The sets deduplicate shared boundary.
    std::set<ResultSegment> segments;
    std::set<ResultPoint> touchPoints;
 
    auto collect = [&](const std::optional<std::variant<ResultPoint, ResultSegment>>& piece) {
        if (!piece) {
            return;
        }
        if (std::holds_alternative<ResultPoint>(*piece)) {
            touchPoints.insert(std::get<ResultPoint>(*piece));
        } else {
            segments.insert(std::get<ResultSegment>(*piece));
        }
    };
 
    for (const auto& edge : edges()) {
        collect(other.template intersection<ResultNumber>(edge));
    }
    for (const auto& piece : intersection<ResultNumber>(Line<typename OtherHalfplane::PointType>(other.source(), other.target()))) {
        if (std::holds_alternative<ResultPoint>(piece)) {
            touchPoints.insert(std::get<ResultPoint>(piece));
        } else {
            segments.insert(std::get<ResultSegment>(piece));
        }
    }
 
    // Make the collected boundary a consistent planar graph (see the polygon
    // overload): split each segment at any other segment's endpoint in its
    // interior, then dedup, so overlapping collinear shared boundary collapses.
    auto planarize = [](const std::set<ResultSegment>& input) {
        std::vector<ResultPoint> vertices;
        for (const auto& s : input) {
            vertices.push_back(s.min());
            vertices.push_back(s.max());
        }
        std::sort(vertices.begin(), vertices.end());
        vertices.erase(std::unique(vertices.begin(), vertices.end()), vertices.end());
 
        std::set<ResultSegment> out;
        for (const auto& s : input) {
            const ResultPoint a = s.min();
            const ResultPoint b = s.max();
            std::vector<ResultPoint> cut{a, b};
            for (const auto& v : vertices) {
                if (v == a || v == b || !collinear(a, b, v)) {
                    continue;
                }
                if (v.x() < std::min(a.x(), b.x()) || v.x() > std::max(a.x(), b.x()) ||
                    v.y() < std::min(a.y(), b.y()) || v.y() > std::max(a.y(), b.y())) {
                    continue;
                }
                cut.push_back(v);
            }
            std::sort(cut.begin(), cut.end());
            cut.erase(std::unique(cut.begin(), cut.end()), cut.end());
            for (std::size_t i = 0; i + 1 < cut.size(); ++i) {
                out.insert(ResultSegment(cut[i], cut[i + 1]));
            }
        }
        return out;
    };
    segments = planarize(segments);
 
    // Build the undirected graph: nodes are endpoints, edges are the segments.
    std::map<ResultPoint, std::vector<ResultPoint>> adjacency;
    for (const auto& segment : segments) {
        adjacency[segment.min()].push_back(segment.max());
        adjacency[segment.max()].push_back(segment.min());
    }
 
    // A touch point that no segment uses is an isolated intersection point.
    for (const auto& point : touchPoints) {
        if (adjacency.find(point) == adjacency.end()) {
            result.emplace_back(point);
        }
    }
 
    using Number = ResultNumber;
    auto degree = [&adjacency](const ResultPoint& p) {
        const auto it = adjacency.find(p);
        return it == adjacency.end() ? std::size_t(0) : it->second.size();
    };
    auto removeEdge = [&adjacency](const ResultPoint& u, const ResultPoint& v) {
        auto& au = adjacency[u];
        au.erase(std::find(au.begin(), au.end(), v));
        auto& av = adjacency[v];
        av.erase(std::find(av.begin(), av.end(), u));
    };
 
    // Peel the open boundary into segments. Every 1D piece lies on the
    // half-plane's straight boundary, so each peeled chain is collinear and is
    // returned as the single segment spanning its two ends.
    while (true) {
        const ResultPoint* leaf = nullptr;
        for (const auto& [p, neighbors] : adjacency) {
            if (neighbors.size() == 1) {
                leaf = &p;
                break;
            }
        }
        if (!leaf) {
            break;
        }
        const ResultPoint start = *leaf;
        ResultPoint current = start;
        while (degree(current) == 1) {
            const ResultPoint next = adjacency.at(current)[0];
            removeEdge(current, next);
            current = next;
        }
        result.emplace_back(ResultSegment(start, current));
    }
 
    // Trace the remaining cycles as faces; keep the counterclockwise (filled)
    // ones. The clockwise-next rule resolves any articulation (pinch) vertices.
    auto rotationalNext = [&adjacency](const ResultPoint& at, const ResultPoint& from) {
        const auto& neighbors = adjacency.at(at);
        auto dx = [&](const ResultPoint& p) { return p.x() - at.x(); };
        auto dy = [&](const ResultPoint& p) { return p.y() - at.y(); };
        auto half = [&](const ResultPoint& p) {
            const Number y = dy(p);
            if (y > Number(0)) return 0;
            if (y < Number(0)) return 1;
            return dx(p) >= Number(0) ? 0 : 1;
        };
        auto ccwBefore = [&](const ResultPoint& u, const ResultPoint& w) {
            const int hu = half(u), hw = half(w);
            if (hu != hw) {
                return hu < hw;
            }
            return dx(u) * dy(w) - dy(u) * dx(w) > Number(0);
        };
        const ResultPoint* best = nullptr;
        for (const auto& n : neighbors) {
            if (n == from || !ccwBefore(n, from)) {
                continue;
            }
            if (!best || ccwBefore(*best, n)) {
                best = &n;
            }
        }
        if (!best) {
            for (const auto& n : neighbors) {
                if (n != from && (!best || ccwBefore(*best, n))) {
                    best = &n;
                }
            }
        }
        return *best;
    };
 
    std::set<std::pair<ResultPoint, ResultPoint>> usedDart;
    for (const auto& [u, neighbors] : adjacency) {
        for (const auto& v : neighbors) {
            if (usedDart.count({u, v})) {
                continue;
            }
            std::vector<ResultPoint> cycle;
            ResultPoint a = u;
            ResultPoint b = v;
            do {
                usedDart.insert({a, b});
                cycle.push_back(a);
                const ResultPoint c = rotationalNext(b, a);
                a = b;
                b = c;
            } while (a != u || b != v);
 
            Number twiceArea(0);
            for (std::size_t i = 0; i < cycle.size(); ++i) {
                const ResultPoint& p = cycle[i];
                const ResultPoint& q = cycle[(i + 1) % cycle.size()];
                twiceArea += p.x() * q.y() - q.x() * p.y();
            }
            if (twiceArea > Number(0)) {
                result.emplace_back(ResultPolygon(cycle));
            }
        }
    }
 
    return result;
}
 
// ---------------------------------------------------------------------------
// Disk
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
Disk<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
// ---------------------------------------------------------------------------
// MonotoneChain
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, MonotoneChainConcept OtherChain>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherChain& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultSegment>;
 
    std::vector<Piece> pieces;
    if (empty() || other.empty()) {
        return pieces;
    }
    if (size() == 1) {
        if (other.contains((*this)[0])) {
            pieces.emplace_back(ResultPoint((*this)[0]));
        }
        return pieces;
    }
    if (other.size() == 1) {
        if (contains(other[0])) {
            pieces.emplace_back(ResultPoint(other[0]));
        }
        return pieces;
    }
 
    // Merge sweep, mirroring intersects(): every edge pair with overlapping
    // x-ranges delegates to the segment-segment intersection. Pairs skipped by
    // a tied advance could only contribute the shared right endpoint, which
    // the pair tested at the tie already reports.
    std::size_t i = 0;
    std::size_t j = 0;
    const std::size_t iEnd = size() - 1;
    const std::size_t jEnd = other.size() - 1;
    while (i < iEnd && j < jEnd) {
        const Segment<PointType> mine((*this)[i], (*this)[i + 1]);
        const Segment<typename OtherChain::PointType> theirs(other[j], other[j + 1]);
        if (!(mine.max().x() < theirs.min().x() || theirs.max().x() < mine.min().x())) {
            if (auto piece = mine.template intersection<ResultNumber>(theirs)) {
                pieces.push_back(std::move(*piece));
            }
        }
        const auto order = mine.max() <=> theirs.max();
        if (order <= 0) {
            ++i;
        }
        if (order >= 0) {
            ++j;
        }
    }
 
    return coalescePieces<ResultNumber>(std::move(pieces));
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::coalescePieces(
    std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                             Segment<Point<ResultNumber, typename PointType::LabelType>>>> pieces) {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultSegment>;
 
    // The fold can report the same touch point through several edge pairs and
    // splits a long collinear overlap at every vertex it crosses; sort the
    // pieces by their lexicographic minimum (segments before points, so a
    // following point is absorbed by the segment sharing its position), then
    // coalesce in one pass.
    std::sort(pieces.begin(), pieces.end(), [](const Piece& lhs, const Piece& rhs) {
        const auto pieceMin = [](const Piece& piece) -> ResultPoint {
            if (const auto* point = std::get_if<ResultPoint>(&piece)) {
                return *point;
            }
            return std::get<ResultSegment>(piece).min();
        };
        const ResultPoint leftMin = pieceMin(lhs);
        const ResultPoint rightMin = pieceMin(rhs);
        if (leftMin != rightMin) {
            return leftMin < rightMin;
        }
        const bool leftIsPoint = std::holds_alternative<ResultPoint>(lhs);
        const bool rightIsPoint = std::holds_alternative<ResultPoint>(rhs);
        if (leftIsPoint != rightIsPoint) {
            return !leftIsPoint;
        }
        if (leftIsPoint) {
            return false;
        }
        return std::get<ResultSegment>(lhs).max() < std::get<ResultSegment>(rhs).max();
    });
 
    std::vector<Piece> result;
    for (const Piece& piece : pieces) {
        if (const auto* point = std::get_if<ResultPoint>(&piece)) {
            if (!result.empty()) {
                // Sorting guarantees any piece that can cover this point has
                // already been placed, and (both pieces lying on this chain's
                // monotone arc) the covering piece is then the previous one.
                const Piece& last = result.back();
                if (const auto* lastPoint = std::get_if<ResultPoint>(&last)) {
                    if (*lastPoint == *point) {
                        continue;
                    }
                } else if (std::get<ResultSegment>(last).contains(*point)) {
                    continue;
                }
            }
            result.push_back(piece);
            continue;
        }
        const auto& segment = std::get<ResultSegment>(piece);
        if (!result.empty()) {
            if (auto* lastSegment = std::get_if<ResultSegment>(&result.back());
                lastSegment != nullptr && !(segment.min() < lastSegment->min()) &&
                !(lastSegment->max() < segment.min()) &&
                collinear(lastSegment->min(), lastSegment->max(), segment.min()) &&
                collinear(lastSegment->min(), lastSegment->max(), segment.max())) {
                // Collinear continuation (chains overlapping across a shared
                // vertex without a bend): extend instead of splitting.
                if (lastSegment->max() < segment.max()) {
                    *lastSegment = ResultSegment(lastSegment->min(), segment.max());
                }
                continue;
            }
        }
        result.push_back(piece);
    }
    return result;
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, class OtherShape>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::edgeFoldIntersection(const OtherShape& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultSegment>;
 
    std::vector<Piece> pieces;
    if (empty()) {
        return pieces;
    }
    if (size() == 1) {
        if (other.contains((*this)[0])) {
            pieces.emplace_back(ResultPoint((*this)[0]));
        }
        return pieces;
    }
    for (std::size_t i = 0; i + 1 < size(); ++i) {
        const auto piece =
            this->template boundaryAt<false>(i).template intersection<ResultNumber>(other);
        if (piece) {
            // Re-wrap in this chain's result types: the delegated intersection
            // labels its points with the other shape's label type.
            pieces.push_back(std::visit(
                [](const auto& value) -> Piece {
                    if constexpr (detail::is_point_v<std::remove_cvref_t<decltype(value)>>) {
                        return Piece(ResultPoint(value));
                    } else {
                        return Piece(ResultSegment(value));
                    }
                },
                *piece));
        }
    }
    return coalescePieces<ResultNumber>(std::move(pieces));
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherSegment& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherOrientedSegment& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherLine& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherOrientedLine& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, RayConcept OtherRay>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherRay& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherHalfplane& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherRectangle& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherTriangle& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
MonotoneChain<PointType, LabelType, Storage>::intersection(const OtherConvex& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
// ---------------------------------------------------------------------------
// Polyline
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
Polyline<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::coalescePieces(
    std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                             Segment<Point<ResultNumber, typename PointType::LabelType>>>> pieces) {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultSegment>;
 
    // Merge the segment pieces first. The kept segments never collinearly
    // touch each other (every insertion absorbs all kept segments it touches),
    // so any kept segment connected to a new one through the growing union
    // must touch the new segment itself: one scan over the kept list per
    // insertion suffices, in any processing order.
    std::vector<ResultSegment> segments;
    for (const Piece& piece : pieces) {
        const auto* segment = std::get_if<ResultSegment>(&piece);
        if (segment == nullptr) {
            continue;
        }
        ResultSegment merged = *segment;
        for (std::size_t i = segments.size(); i-- > 0;) {
            const ResultSegment& kept = segments[i];
            if (collinear(merged.min(), merged.max(), kept.min()) &&
                collinear(merged.min(), merged.max(), kept.max()) &&
                merged.intersects(kept)) {
                // Collinear touching segments cover exactly the span between
                // their lexicographically extreme endpoints.
                merged = ResultSegment(std::min(merged.min(), kept.min()),
                                       std::max(merged.max(), kept.max()));
                segments.erase(segments.begin() + static_cast<std::ptrdiff_t>(i));
            }
        }
        segments.push_back(std::move(merged));
    }
 
    // A touch point can be reported by several edge pairs and may lie on a
    // reported segment; keep one copy of each point no segment covers.
    std::vector<ResultPoint> points;
    for (const Piece& piece : pieces) {
        if (const auto* point = std::get_if<ResultPoint>(&piece)) {
            points.push_back(*point);
        }
    }
    std::sort(points.begin(), points.end());
    points.erase(std::unique(points.begin(), points.end()), points.end());
 
    std::vector<Piece> result;
    result.reserve(segments.size() + points.size());
    for (const ResultPoint& point : points) {
        if (std::none_of(segments.begin(), segments.end(),
                         [&point](const ResultSegment& segment) { return segment.contains(point); })) {
            result.emplace_back(point);
        }
    }
    for (ResultSegment& segment : segments) {
        result.emplace_back(std::move(segment));
    }
 
    // Same output order as MonotoneChain::coalescePieces: lexicographic
    // minimum, segments before points sharing their minimum, then by maximum.
    std::sort(result.begin(), result.end(), [](const Piece& lhs, const Piece& rhs) {
        const auto pieceMin = [](const Piece& piece) -> ResultPoint {
            if (const auto* point = std::get_if<ResultPoint>(&piece)) {
                return *point;
            }
            return std::get<ResultSegment>(piece).min();
        };
        const ResultPoint leftMin = pieceMin(lhs);
        const ResultPoint rightMin = pieceMin(rhs);
        if (leftMin != rightMin) {
            return leftMin < rightMin;
        }
        const bool leftIsPoint = std::holds_alternative<ResultPoint>(lhs);
        const bool rightIsPoint = std::holds_alternative<ResultPoint>(rhs);
        if (leftIsPoint != rightIsPoint) {
            return !leftIsPoint;
        }
        if (leftIsPoint) {
            return false;
        }
        return std::get<ResultSegment>(lhs).max() < std::get<ResultSegment>(rhs).max();
    });
    return result;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, class OtherShape>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::edgeFoldIntersection(const OtherShape& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultSegment>;
 
    std::vector<Piece> pieces;
    if (empty()) {
        return pieces;
    }
    if (size() == 1) {
        if (other.contains((*this)[0])) {
            pieces.emplace_back(ResultPoint((*this)[0]));
        }
        return pieces;
    }
    for (std::size_t i = 0; i + 1 < size(); ++i) {
        const auto piece =
            this->template boundaryAt<false>(i).template intersection<ResultNumber>(other);
        if (piece) {
            // Re-wrap in this polyline's result types: the delegated
            // intersection labels its points with the other shape's label type.
            pieces.push_back(std::visit(
                [](const auto& value) -> Piece {
                    if constexpr (detail::is_point_v<std::remove_cvref_t<decltype(value)>>) {
                        return Piece(ResultPoint(value));
                    } else {
                        return Piece(ResultSegment(value));
                    }
                },
                *piece));
        }
    }
    return coalescePieces<ResultNumber>(std::move(pieces));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherSegment& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherLine& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherOrientedLine& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherRay& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherHalfplane& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherRectangle& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherTriangle& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherConvex& other) const {
    return this->template edgeFoldIntersection<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, MonotoneChainConcept OtherChain>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherChain& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultSegment>;
 
    std::vector<Piece> pieces;
    if (empty() || other.empty()) {
        return pieces;
    }
    if (size() == 1) {
        if (other.contains((*this)[0])) {
            pieces.emplace_back(ResultPoint((*this)[0]));
        }
        return pieces;
    }
    if (other.size() == 1) {
        if (contains(other[0])) {
            pieces.emplace_back(ResultPoint(other[0]));
        }
        return pieces;
    }
    if (!bbox().intersects(other.bbox())) {
        return pieces;
    }
 
    // All-pairs edge test: a self-intersecting polyline has no monotone
    // structure to drive a merge sweep.
    for (std::size_t i = 0; i + 1 < size(); ++i) {
        const Segment<PointType> mine((*this)[i], (*this)[i + 1]);
        for (std::size_t j = 0; j + 1 < other.size(); ++j) {
            const Segment<typename OtherChain::PointType> theirs(other[j], other[j + 1]);
            if (auto piece = mine.template intersection<ResultNumber>(theirs)) {
                pieces.push_back(std::move(*piece));
            }
        }
    }
    return coalescePieces<ResultNumber>(std::move(pieces));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PolylineConcept OtherPolyline>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::intersection(const OtherPolyline& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultSegment>;
 
    std::vector<Piece> pieces;
    if (empty() || other.empty()) {
        return pieces;
    }
    if (size() == 1) {
        if (other.contains((*this)[0])) {
            pieces.emplace_back(ResultPoint((*this)[0]));
        }
        return pieces;
    }
    if (other.size() == 1) {
        if (contains(other[0])) {
            pieces.emplace_back(ResultPoint(other[0]));
        }
        return pieces;
    }
    if (!bbox().intersects(other.bbox())) {
        return pieces;
    }
 
    // All-pairs edge test: neither polyline has a monotone structure to prune
    // the scan.
    for (std::size_t i = 0; i + 1 < size(); ++i) {
        const Segment<PointType> mine((*this)[i], (*this)[i + 1]);
        for (std::size_t j = 0; j + 1 < other.size(); ++j) {
            const Segment<typename OtherPolyline::PointType> theirs(other[j], other[j + 1]);
            if (auto piece = mine.template intersection<ResultNumber>(theirs)) {
                pieces.push_back(std::move(*piece));
            }
        }
    }
    return coalescePieces<ResultNumber>(std::move(pieces));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, class OtherArea>
    requires(PolygonConcept<OtherArea> || PolygonWithHolesConcept<OtherArea>)
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType::LabelType>>>>
Polyline<PointType, LabelType>::polygonIntersection(const OtherArea& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultSegment>;
 
    std::vector<Piece> pieces;
    if (empty()) {
        return pieces;
    }
    if (size() == 1) {
        if (other.contains((*this)[0])) {
            pieces.emplace_back(ResultPoint((*this)[0]));
        }
        return pieces;
    }
    // Unlike the convex edgeFoldIntersection, a non-convex area can split a
    // single edge into several disjoint pieces, so each edge yields a vector.
    for (std::size_t i = 0; i + 1 < size(); ++i) {
        const auto edgePieces =
            this->template boundaryAt<false>(i).template intersection<ResultNumber>(other);
        for (const auto& piece : edgePieces) {
            // Re-wrap in this polyline's result types: the delegated
            // intersection labels its points with the area's label type.
            pieces.push_back(std::visit(
                [](const auto& value) -> Piece {
                    if constexpr (detail::is_point_v<std::remove_cvref_t<decltype(value)>>) {
                        return Piece(ResultPoint(value));
                    } else {
                        return Piece(ResultSegment(value));
                    }
                },
                piece));
        }
    }
    return coalescePieces<ResultNumber>(std::move(pieces));
}
 
// Polygon outranks Polyline and MonotoneChain, so it owns these pairs; the
// lower-ranked chains reach them by forwarding up. The polyline-vs-polygon clip
// itself lives on Polyline (reusing its coalescing), so both overloads delegate
// there -- a monotone chain first views itself as a polyline.
template <class PointType, class LabelType>
template <class ResultNumber, PolylineConcept OtherPolyline>
constexpr auto Polygon<PointType, LabelType>::intersection(const OtherPolyline& other) const {
    return other.template polygonIntersection<ResultNumber>(*this);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, MonotoneChainConcept OtherChain>
constexpr auto Polygon<PointType, LabelType>::intersection(const OtherChain& other) const {
    return other.asPolyline().template polygonIntersection<ResultNumber>(*this);
}
 
 
// ---------------------------------------------------------------------------
// PolygonWithHoles
//
// The one-dimensional operands, which the region-valued booleans of
// implementation/booleans.hpp leave alone: `closure(A° ∩ B°)` is empty for every
// one of them, so they need the plain, unregularized intersection instead. Each
// clip runs over the boundary edges of every ring at once -- see the helpers
// above -- which is all that separates these from Polygon's.
 
template <class PointType_, class TLabel>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherSegment& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    return detail::areaSegmentIntersection<ResultPoint>(*this, orientedEdges(), other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherOrientedSegment& other) const {
    return intersection<ResultNumber>(
        Segment<typename OtherOrientedSegment::PointType>(other[0], other[1]));
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, LineConcept OtherLine>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherLine& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    return detail::areaLineIntersection<ResultPoint>(*this, orientedEdges(), other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherOrientedLine& other) const {
    return intersection<ResultNumber>(other.asLine());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RayConcept OtherRay>
constexpr std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                                   Segment<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherRay& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    return detail::areaRayIntersection<ResultPoint>(*this, orientedEdges(), other);
}
 
// A region outranks Polyline and MonotoneChain, so it owns these pairs; as with
// Polygon, the clip itself lives on Polyline (reusing its coalescing) and a
// monotone chain first views itself as a polyline.
template <class PointType_, class TLabel>
template <class ResultNumber, PolylineConcept OtherPolyline>
constexpr auto PolygonWithHoles<PointType_, TLabel>::intersection(const OtherPolyline& other) const {
    return other.template polygonIntersection<ResultNumber>(*this);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, MonotoneChainConcept OtherChain>
constexpr auto PolygonWithHoles<PointType_, TLabel>::intersection(const OtherChain& other) const {
    return other.asPolyline().template polygonIntersection<ResultNumber>(*this);
}
 
 
// ---------------------------------------------------------------------------
// HalfplaneIntersection
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherPoint& other) const {
    if (contains(other)) {
        return Point<ResultNumber, typename PointType::LabelType>(other);
    }
    return {};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
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>>>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherLine& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    if (empty()) {
        return {};
    }
    if (other.isDegenerate()) {
        if (contains(other[0])) {
            return ResultPoint(other[0]);
        }
        return {};
    }
    if (halfplanes_.empty()) {
        return Line<ResultPoint>(ResultPoint(other[0]), ResultPoint(other[1]));
    }
    const Halfplane<typename OtherLine::PointType> along(other[0], other[1]);
    const auto clip = clipLine(along);
    if (clip.empty) {
        return {};
    }
    const Line<typename OtherLine::PointType> supporting(other[0], other[1]);
    const auto crossing = [&](std::ptrdiff_t idx) {
        const auto isec = halfplanes_[static_cast<std::size_t>(idx)].asLine()
                              .template intersection<ResultNumber>(supporting);
        assert(isec && isec->index() == 0);
        return std::get<0>(*isec);
    };
    if (clip.entry < 0 && clip.exit < 0) {
        return Line<ResultPoint>(ResultPoint(other[0]), ResultPoint(other[1]));
    }
    const ResultPoint source(other[0]);
    const ResultPoint target(other[1]);
    const ResultNumber dx = target.x() - source.x();
    const ResultNumber dy = target.y() - source.y();
    if (clip.entry < 0) {
        const ResultPoint finish = crossing(clip.exit);
        return Ray<ResultPoint>(finish, ResultPoint(finish.x() - dx, finish.y() - dy));
    }
    if (clip.exit < 0) {
        const ResultPoint start = crossing(clip.entry);
        return Ray<ResultPoint>(start, ResultPoint(start.x() + dx, start.y() + dy));
    }
    const auto det = detail::boundaryLinesDeterminant(
        halfplanes_[static_cast<std::size_t>(clip.entry)],
        halfplanes_[static_cast<std::size_t>(clip.exit)], along);
    if (det == decltype(det){}) {
        return crossing(clip.entry);
    }
    return Segment<ResultPoint>(crossing(clip.entry), crossing(clip.exit));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
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>>>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherOrientedLine& other) const {
    return intersection<ResultNumber>(other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                     Segment<Point<ResultNumber, typename PointType::LabelType>>>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherSegment& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    if (empty()) {
        return {};
    }
    if (other.isDegenerate()) {
        if (contains(other.min())) {
            return ResultPoint(other.min());
        }
        return {};
    }
    if (halfplanes_.empty()) {
        return Segment<ResultPoint>(ResultPoint(other.min()), ResultPoint(other.max()));
    }
    const Halfplane<typename OtherSegment::PointType> along(other.min(), other.max());
    const auto clip = clipLine(along);
    if (clip.empty) {
        return {};
    }
    const Line<typename OtherSegment::PointType> supporting(other.min(), other.max());
    const auto crossing = [&](std::ptrdiff_t idx) {
        const auto isec = halfplanes_[static_cast<std::size_t>(idx)].asLine()
                              .template intersection<ResultNumber>(supporting);
        assert(isec && isec->index() == 0);
        return std::get<0>(*isec);
    };
    // Clamp the clip interval to the segment's [0, 1] parameter window; the
    // comparisons against the window edges are point-side tests.
    const bool startsAtMin = clip.entry < 0 ||
        constraintSide(static_cast<std::size_t>(clip.entry), other.min()) >= 0;
    if (!startsAtMin && constraintSide(static_cast<std::size_t>(clip.entry), other.max()) < 0) {
        return {};  // the region begins after the segment ends
    }
    const bool endsAtMax = clip.exit < 0 ||
        constraintSide(static_cast<std::size_t>(clip.exit), other.max()) >= 0;
    if (!endsAtMax && constraintSide(static_cast<std::size_t>(clip.exit), other.min()) < 0) {
        return {};  // the region ends before the segment starts
    }
    if (startsAtMin && endsAtMax) {
        return Segment<ResultPoint>(ResultPoint(other.min()), ResultPoint(other.max()));
    }
    if (startsAtMin) {
        if (constraintSide(static_cast<std::size_t>(clip.exit), other.min()) == 0) {
            return ResultPoint(other.min());
        }
        return Segment<ResultPoint>(ResultPoint(other.min()), crossing(clip.exit));
    }
    if (endsAtMax) {
        if (constraintSide(static_cast<std::size_t>(clip.entry), other.max()) == 0) {
            return ResultPoint(other.max());
        }
        return Segment<ResultPoint>(crossing(clip.entry), ResultPoint(other.max()));
    }
    const auto det = detail::boundaryLinesDeterminant(
        halfplanes_[static_cast<std::size_t>(clip.entry)],
        halfplanes_[static_cast<std::size_t>(clip.exit)], along);
    if (det == decltype(det){}) {
        return crossing(clip.entry);
    }
    return Segment<ResultPoint>(crossing(clip.entry), crossing(clip.exit));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                     Segment<Point<ResultNumber, typename PointType::LabelType>>>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherOrientedSegment& other) const {
    return intersection<ResultNumber>(Segment<typename OtherOrientedSegment::PointType>(other[0], other[1]));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr std::optional<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                     Segment<Point<ResultNumber, typename PointType::LabelType>>,
                                     Ray<Point<ResultNumber, typename PointType::LabelType>>>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherRay& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    if (empty()) {
        return {};
    }
    if (halfplanes_.empty()) {
        return Ray<ResultPoint>(ResultPoint(other.source()), ResultPoint(other.target()));
    }
    const Halfplane<typename OtherRay::PointType> along(other.source(), other.target());
    const auto clip = clipLine(along);
    if (clip.empty) {
        return {};
    }
    const Line<typename OtherRay::PointType> supporting(other.source(), other.target());
    const auto crossing = [&](std::ptrdiff_t idx) {
        const auto isec = halfplanes_[static_cast<std::size_t>(idx)].asLine()
                              .template intersection<ResultNumber>(supporting);
        assert(isec && isec->index() == 0);
        return std::get<0>(*isec);
    };
    // Clamp the clip interval to the ray's [0, +inf) parameter window.
    const bool startsAtSource = clip.entry < 0 ||
        constraintSide(static_cast<std::size_t>(clip.entry), other.source()) >= 0;
    if (clip.exit < 0) {
        if (startsAtSource) {
            return Ray<ResultPoint>(ResultPoint(other.source()), ResultPoint(other.target()));
        }
        const ResultPoint start = crossing(clip.entry);
        const ResultPoint source(other.source());
        const ResultPoint target(other.target());
        return Ray<ResultPoint>(start, ResultPoint(start.x() + (target.x() - source.x()),
                                                   start.y() + (target.y() - source.y())));
    }
    const auto exitSide = constraintSide(static_cast<std::size_t>(clip.exit), other.source());
    if (exitSide < 0) {
        return {};  // the region ends before the ray starts
    }
    if (startsAtSource) {
        if (exitSide == 0) {
            return ResultPoint(other.source());
        }
        return Segment<ResultPoint>(ResultPoint(other.source()), crossing(clip.exit));
    }
    const auto det = detail::boundaryLinesDeterminant(
        halfplanes_[static_cast<std::size_t>(clip.entry)],
        halfplanes_[static_cast<std::size_t>(clip.exit)], along);
    if (det == decltype(det){}) {
        return crossing(clip.entry);
    }
    return Segment<ResultPoint>(crossing(clip.entry), crossing(clip.exit));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherHalfplane& other) const {
    // Half-plane intersections are closed under intersecting with one more
    // half-plane, and no coordinate division is involved, so the result is
    // exact whenever ResultNumber represents the inputs exactly.
    HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>> result(*this);
    result.insert(other);
    return result;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherRectangle& other) const {
    // Intersecting with the rectangle's four edge half-planes stays closed and
    // needs no division, so the result is exact whenever ResultNumber
    // represents the inputs exactly.
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultHalfplane = typename HalfplaneIntersection<ResultPoint>::HalfplaneType;
    HalfplaneIntersection<ResultPoint> result(*this);
    if (other.empty()) {
        // The empty rectangle carries no edge, so the loop below would insert no
        // constraint and hand back this whole region. Force emptiness with two
        // contradictory parallel constraints ({x <= 0} and {x >= 1}), as the
        // Convex and half-plane-intersection overloads do for their own empty
        // operand.
        result.insert(ResultHalfplane(ResultPoint(0, 0), ResultPoint(0, 1)));
        result.insert(ResultHalfplane(ResultPoint(1, 1), ResultPoint(1, 0)));
        return result;
    }
    for (const auto& halfplane : HalfplaneIntersection<ResultPoint>(other)) {
        result.insert(halfplane);
    }
    return result;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherTriangle& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    HalfplaneIntersection<ResultPoint> result(*this);
    for (const auto& halfplane : HalfplaneIntersection<ResultPoint>(other)) {
        result.insert(halfplane);
    }
    return result;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherConvex& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultHalfplane = typename HalfplaneIntersection<ResultPoint>::HalfplaneType;
    HalfplaneIntersection<ResultPoint> result(*this);
    if (other.size() == 0) {
        // The empty convex polygon is the empty set: force emptiness with two
        // contradictory parallel constraints ({x <= 0} and {x >= 1}).
        result.insert(ResultHalfplane(ResultPoint(0, 0), ResultPoint(0, 1)));
        result.insert(ResultHalfplane(ResultPoint(1, 1), ResultPoint(1, 0)));
        return result;
    }
    for (const auto& halfplane : HalfplaneIntersection<ResultPoint>(other)) {
        result.insert(halfplane);
    }
    return result;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneIntersectionConcept OtherRegion>
constexpr HalfplaneIntersection<Point<ResultNumber, typename PointType::LabelType>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherRegion& other) const {
    // Half-plane intersections are closed under intersection and no division
    // is involved, so the result is exact whenever ResultNumber represents
    // the inputs exactly.
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultHalfplane = typename HalfplaneIntersection<ResultPoint>::HalfplaneType;
    HalfplaneIntersection<ResultPoint> result(*this);
    if (other.empty()) {
        // Force emptiness with two contradictory parallel constraints
        // ({x <= 0} and {x >= 1}).
        result.insert(ResultHalfplane(ResultPoint(0, 0), ResultPoint(0, 1)));
        result.insert(ResultHalfplane(ResultPoint(1, 1), ResultPoint(1, 0)));
        return result;
    }
    for (const auto& halfplane : other) {
        result.insert(ResultHalfplane(halfplane));
    }
    return result;
}
 
// A polygon is bounded, so only the part of the region near it matters: the
// region clipped to the polygon's own bounding rectangle has the same
// intersection with it — P ⊆ bbox(P), so A ∩ P = (A ∩ bbox) ∩ P — and, being
// bounded and convex, is a convex polygon, which Polygon::intersection already
// handles. That is where all the work happens; everything here is the
// reduction to it.
//
// The clip is done over exact rationals, and against bbox(P) itself rather
// than an inflated box: an inflated one would move the clip vertices off the
// crossings the answer already carries and deepen them for nothing, which for
// an inexact ResultNumber is the difference between an answer that closes up
// and one that does not.
template <class PointType, class LabelType>
template <class ResultNumber, PolygonConcept OtherPolygon>
std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType::LabelType>>,
                         Polygon<Point<ResultNumber, typename PointType::LabelType>>>>
HalfplaneIntersection<PointType, LabelType>::intersection(const OtherPolygon& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    using ResultSegment = Segment<ResultPoint>;
    using ResultPolyline = Polyline<ResultPoint>;
    using ResultPolygon = Polygon<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultPolyline, ResultPolygon>;
    using ExactNumber = detail::region_exact_number_t<NumberType>;
    using ExactPoint = Point<ExactNumber, typename PointType::LabelType>;
    using ExactRegion = HalfplaneIntersection<ExactPoint>;
 
    std::vector<Piece> result;
    if (empty() || other.size() == 0) {
        return result;
    }
    // The polygon is relabelled into the region's label type so the components
    // come back carrying it, as every other overload here does.
    const Polygon<ExactPoint> exact(other);
 
    // A region with empty interior is its carrier — a point, a segment, a ray
    // or a line — and the polygon clips those itself. Their pieces are points
    // and segments, and a segment is the two-vertex polyline of this contract.
    const auto carrierPieces = [&exact, &result](const auto& region) {
        std::visit(
            [&exact, &result](const auto& carrier) {
                using Carrier = std::remove_cvref_t<decltype(carrier)>;
                if constexpr (PointConcept<Carrier>) {
                    if (exact.contains(carrier)) {
                        result.emplace_back(ResultPoint(carrier));
                    }
                } else {
                    for (const auto& piece : exact.template intersection<ResultNumber>(carrier)) {
                        if (const auto* point = std::get_if<ResultPoint>(&piece)) {
                            result.emplace_back(*point);
                        } else {
                            const ResultSegment& chord = std::get<ResultSegment>(piece);
                            result.emplace_back(
                                ResultPolyline(std::vector<ResultPoint>{chord.min(), chord.max()}));
                        }
                    }
                }
            },
            detail::degenerateRegionCarrier(region));
    };
 
    if (isDegenerate()) {
        carrierPieces(*this);  // unclipped: the polygon handles a ray or a line itself
        return result;
    }
 
    ExactRegion clipped(*this);
    for (const auto& halfplane : ExactRegion(exact.bbox())) {
        clipped.insert(halfplane);
    }
    if (clipped.empty()) {
        return result;  // the region misses the polygon's bounding rectangle
    }
    // A full-dimensional region can still meet the bounding rectangle in a
    // segment or a point, and then so does everything below it.
    if (clipped.isDegenerate()) {
        carrierPieces(clipped);
        return result;
    }
    return exact.template intersection<ResultNumber>(clipped.template asConvex<ExactNumber>());
}
 
}  // namespace pgl