arrangement.hpp

Go to the documentation of this file.

#pragma once
 
#include "algorithm/triangulation.hpp"

 
#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <map>
#include <memory>
#include <numeric>
#include <optional>
#include <random>
#include <ranges>
#include <set>
#include <span>
#include <stdexcept>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <variant>
#include <vector>
 
namespace pgl {
 
namespace detail {
 
struct SimpleBoundariesTag {};
inline constexpr SimpleBoundariesTag simpleBoundaries;

template <class ExactPoint>
int ringOrientation(const std::vector<ExactPoint>& ring) {
    if (ring.size() < 3) {
        return 0;
    }
    const std::size_t leftmost =
        static_cast<std::size_t>(std::min_element(ring.begin(), ring.end()) - ring.begin());
    const ExactPoint& corner = ring[leftmost];
    const ExactPoint& ahead = ring[(leftmost + 1) % ring.size()];
    const ExactPoint& behind = ring[(leftmost + ring.size() - 1) % ring.size()];
    const auto turn = orientationSign(corner, ahead, behind);
    if (turn > 0) {
        return 1;
    }
    return turn < 0 ? -1 : 0;
}

template <class ExactPoint>
void splitWalkIntoRings(const std::vector<ExactPoint>& walk,
                        std::vector<std::vector<ExactPoint>>& out) {
    std::vector<ExactPoint> pending;
    std::map<ExactPoint, std::size_t> position;
    for (const ExactPoint& vertex : walk) {
        const auto seen = position.find(vertex);
        if (seen != position.end()) {
            const std::size_t from = seen->second;
            out.emplace_back(pending.begin() + static_cast<std::ptrdiff_t>(from), pending.end());
            for (std::size_t i = from; i < pending.size(); ++i) {
                position.erase(pending[i]);
            }
            pending.resize(from);
        }
        position.emplace(vertex, pending.size());
        pending.push_back(vertex);
    }
    if (!pending.empty()) {
        out.push_back(std::move(pending));
    }
}
 
}  // namespace detail

template <class PointType_, class TLabel>
class Arrangement {
    // Tags that keep the three handle families distinct types, so a vertex
    // handle can never be passed where a halfedge or face one is meant.
    struct VertexTag;
    struct HalfedgeTag;
    struct FaceTag;
    class TrapezoidPointLocation;
 
public:
    using VertexId = detail::Handle<VertexTag>;
    using HalfedgeId = detail::Handle<HalfedgeTag>;
    using FaceId = detail::Handle<FaceTag>;
    using CellId = std::variant<VertexId, HalfedgeId, FaceId>;

    using PointType = PointType_;
    using NumberType = typename PointType::NumberType;
    using LabelType = TLabel;
    using SegmentType = Segment<PointType, TLabel>;
    using LineType = Line<PointType, TLabel>;
    using OrientedSegmentType = OrientedSegment<PointType, TLabel>;
    using OrientedLineType = OrientedLine<PointType, TLabel>;
    using RayType = Ray<PointType, TLabel>;
    using HalfedgeType = std::variant<OrientedSegmentType, OrientedLineType, RayType>;
    using EdgeType = std::variant<SegmentType, LineType, RayType>;
    using IntersectionId = std::variant<HalfedgeId, VertexId>;
 
    static_assert(detail::is_point_v<PointType>, "Arrangement requires pgl::Point vertices");

    Arrangement() {
        buildFaces();
    }

    template <std::ranges::input_range ShapeRange>
    explicit Arrangement(const ShapeRange& shapes) {
        std::vector<InputSegment> segments;
        std::vector<InputCurve> curves;
        std::vector<PointType> isolated;
        collect(shapes, segments, curves, isolated);
        build(segments, curves, isolated, false);
    }
 
    // Internal fast path for boolean operands whose boundary rings are known
    // simple and non-overlapping within each input shape.
    template <std::ranges::input_range ShapeRange>
    Arrangement(const ShapeRange& shapes, detail::SimpleBoundariesTag) {
        std::vector<InputSegment> segments;
        std::vector<InputCurve> curves;
        std::vector<PointType> isolated;
        collect(shapes, segments, curves, isolated);
        build(segments, curves, isolated, true);
    }

    template <std::ranges::input_range ShapeRange, std::ranges::input_range PointRange>
    Arrangement(const ShapeRange& shapes, const PointRange& points) {
        std::vector<InputSegment> segments;
        std::vector<InputCurve> curves;
        std::vector<PointType> isolated;
        collect(shapes, segments, curves, isolated);
        for (const auto& point : points) {
            isolated.emplace_back(point);
        }
        build(segments, curves, isolated, false);
    }
 
    // -------------------------------------------------------------------------
    // Cells

    [[nodiscard]] std::size_t vertexCount() const {
        return points_.size();
    }

    [[nodiscard]] std::size_t halfedgeCount() const {
        return origin_.size();
    }

    [[nodiscard]] std::size_t edgeCount() const {
        return origin_.size() / 2;
    }

    [[nodiscard]] std::size_t faceCount() const {
        return outerCycle_.size();
    }

    [[nodiscard]] const std::vector<PointType>& vertices() const {
        return points_;
    }

    [[nodiscard]] std::vector<SegmentType> boundedEdges() const {
        std::vector<SegmentType> result;
        result.reserve(edgeCount());
        for (std::size_t i = 0; i < edgeGeometry_.size(); ++i) {
            const EdgeGeometry& geometry = edgeGeometry_[i];
            if (geometry.kind != EdgeKind::segment) {
                continue;
            }
            SegmentType edge(geometry.a, geometry.b);
            if constexpr (detail::has_label_v<TLabel>) {
                edge.label() = edgeLabel_[i];
            }
            result.push_back(std::move(edge));
        }
        return result;
    }

    [[nodiscard]] std::vector<EdgeType> edges() const {
        std::vector<EdgeType> result;
        result.reserve(edgeCount());
        for (std::size_t i = 0; i < edgeGeometry_.size(); ++i) {
            const EdgeGeometry& geometry = edgeGeometry_[i];
            if (geometry.kind == EdgeKind::segment) {
                SegmentType edge(geometry.a, geometry.b);
                if constexpr (detail::has_label_v<TLabel>) {
                    edge.label() = edgeLabel_[i];
                }
                result.emplace_back(std::move(edge));
            } else if (geometry.kind == EdgeKind::line) {
                LineType edge(geometry.a, geometry.b);
                if constexpr (detail::has_label_v<TLabel>) {
                    edge.label() = edgeLabel_[i];
                }
                result.emplace_back(std::move(edge));
            } else {
                RayType edge(geometry.a, geometry.b);
                if constexpr (detail::has_label_v<TLabel>) {
                    edge.label() = edgeLabel_[i];
                }
                result.emplace_back(std::move(edge));
            }
        }
        return result;
    }

    [[nodiscard]] const PointType& operator[](VertexId v) const {
        if (!v.valid() || v.index() >= points_.size()) {
            throw std::logic_error("the fictitious arrangement vertex has no finite position");
        }
        return points_[v.index()];
    }

    [[nodiscard]] HalfedgeType operator[](HalfedgeId h) const {
        assert(h.valid() && h.index() < origin_.size());
        const EdgeGeometry& geometry = edgeGeometry_[h.index() / 2];
        const TLabel& label = edgeLabel_[h.index() / 2];
        if (geometry.kind == EdgeKind::segment) {
            OrientedSegmentType segment(points_[origin_[h.index()]],
                                        points_[origin_[h.index() ^ 1]]);
            if constexpr (detail::has_label_v<TLabel>) {
                segment.label() = label;
            }
            return segment;
        }
        if (geometry.kind == EdgeKind::line) {
            OrientedLineType line(h.index() % 2 == 0 ? geometry.a : geometry.b,
                                  h.index() % 2 == 0 ? geometry.b : geometry.a);
            if constexpr (detail::has_label_v<TLabel>) {
                line.label() = label;
            }
            return line;
        }
        RayType ray(geometry.a, geometry.b);
        if constexpr (detail::has_label_v<TLabel>) {
            ray.label() = label;
        }
        return ray;
    }
 
    // -------------------------------------------------------------------------
    // Incidence

    [[nodiscard]] HalfedgeId twin(HalfedgeId h) const {
        assert(h.valid() && h.index() < origin_.size());
        return HalfedgeId(h.index() ^ 1);
    }

    [[nodiscard]] HalfedgeId next(HalfedgeId h) const {
        assert(h.valid() && h.index() < next_.size());
        return HalfedgeId(next_[h.index()]);
    }

    [[nodiscard]] VertexId source(HalfedgeId h) const {
        assert(h.valid() && h.index() < origin_.size());
        return VertexId(origin_[h.index()]);
    }

    [[nodiscard]] VertexId target(HalfedgeId h) const {
        return source(twin(h));
    }

    [[nodiscard]] FaceId face(HalfedgeId h) const {
        assert(h.valid() && h.index() < face_.size());
        return FaceId(face_[h.index()]);
    }

    [[nodiscard]] HalfedgeId outgoing(VertexId v) const {
        assert(v.valid() && v.index() < outgoing_.size());
        return outgoing_[v.index()];
    }

    [[nodiscard]] std::size_t degree(VertexId v) const {
        assert(v.valid() && v.index() < outgoing_.size());
        const HalfedgeId start = outgoing_[v.index()];
        if (!start.valid()) {
            return 0;
        }
        std::size_t count = 0;
        HalfedgeId h = start;
        do {
            ++count;
            h = next(twin(h));
        } while (h != start);
        return count;
    }

    [[nodiscard]] std::vector<HalfedgeId> outgoingHalfedges(VertexId v) const {
        assert(v.valid() && v.index() < outgoing_.size());
        std::vector<HalfedgeId> halfedges;
        const HalfedgeId start = outgoing_[v.index()];
        if (!start.valid()) {
            return halfedges;
        }
        HalfedgeId h = start;
        do {
            halfedges.push_back(h);
            h = next(twin(h));
        } while (h != start);
        return halfedges;
    }

    [[nodiscard]] Graph<VertexId> asGraph() const {
        Graph<VertexId> result;
        for (std::uint32_t v = 0; v < topologicalVertexCount(); ++v) {
            result.addVertex(VertexId(v));
        }
        for (std::size_t h = 0; h < origin_.size(); h += 2) {
            result.addEdge(VertexId(origin_[h]), VertexId(origin_[h + 1]));
        }
        return result;
    }
 
    // -------------------------------------------------------------------------
    // The uniform cell interface

    template <class ResultNumber = NumberType>
    [[nodiscard]] Point<ResultNumber> witness(VertexId v) const {
        return Point<ResultNumber>((*this)[v]);
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] Point<ResultNumber> witness(HalfedgeId h) const {
        assert(h.valid() && h.index() < origin_.size());
        const EdgeGeometry& geometry = edgeGeometry_[h.index() / 2];
        const PointType& a = geometry.a;
        const PointType& b = geometry.b;
        const ResultNumber two = static_cast<ResultNumber>(NumberType(2));
        return Point<ResultNumber>(
            (detail::asNumber<ResultNumber>(a.x()) + detail::asNumber<ResultNumber>(b.x())) / two,
            (detail::asNumber<ResultNumber>(a.y()) + detail::asNumber<ResultNumber>(b.y())) / two);
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] Point<ResultNumber> witness(FaceId f) const {
        assert(f.valid() && !isUnbounded(f));
        if (hasSimpleBoundary(f)) {
            return ringWitness<ResultNumber>(outerCycle_[f.index()]);
        }
        return sweptWitness<ResultNumber>(f);
    }
 
    // -------------------------------------------------------------------------
    // Faces

    [[nodiscard]] bool hasSimpleBoundary(FaceId f) const {
        assert(f.valid() && f.index() < outerCycle_.size());
        if (!outerCycle_[f.index()].valid() || !innerCycles(f).empty()) {
            return false;
        }
        const std::uint32_t start = outerCycle_[f.index()].index();
        std::uint32_t h = start;
        do {
            if (face_[h ^ 1] == f.index()) {
                return false;
            }
            h = next_[h];
        } while (h != start);
        return true;
    }
 
private:
    // The witness of a simple ring, by the argument Polygon::pointInside uses:
    // the ring's leftmost vertex is convex, so either the triangle it makes with
    // its two neighbours is an ear — and its interior point will do — or the
    // nearest vertex inside that triangle cuts a diagonal whose midpoint is
    // interior. Reading the cycle in place keeps this allocation-free.
    template <class ResultNumber>
    [[nodiscard]] Point<ResultNumber> ringWitness(HalfedgeId start) const {
        std::uint32_t leftmost = start.index();
        std::uint32_t beforeLeftmost = start.index();
        std::uint32_t previous = start.index();
        for (std::uint32_t h = next_[start.index()]; h != start.index();
             previous = h, h = next_[h]) {
            if (points_[origin_[h]] < points_[origin_[leftmost]]) {
                leftmost = h;
                beforeLeftmost = previous;
            }
        }
        if (leftmost == start.index()) {
            beforeLeftmost = previous;  // the cycle's last halfedge closes onto it
        }
        const PointType& corner = points_[origin_[leftmost]];
        const PointType& ahead = points_[origin_[next_[leftmost]]];
        const PointType& behind = points_[origin_[beforeLeftmost]];
        const Triangle<PointType> ear(corner, ahead, behind);
        assert(!ear.isDegenerate());
 
        // The leftmost vertex inside the ear, skipping the ear's own corners.
        const PointType* diagonal = nullptr;
        for (std::uint32_t h = next_[next_[leftmost]]; h != leftmost; h = next_[h]) {
            const PointType& vertex = points_[origin_[h]];
            if (ear.interiorContains(vertex) &&
                (diagonal == nullptr || vertex < *diagonal)) {
                diagonal = &vertex;
            }
        }
        if (diagonal != nullptr) {
            return (Point<ResultNumber>(corner) + Point<ResultNumber>(*diagonal)) /
                   static_cast<ResultNumber>(NumberType(2));
        }
        return ear.template pointInside<ResultNumber>();
    }
 
    // The witness of a face whose boundary is anything else.
    template <class ResultNumber>
    [[nodiscard]] Point<ResultNumber> sweptWitness(FaceId f) const {
        const HalfedgeId seed = outerCycle_[f.index()];
        const PointType& a = points_[origin_[seed.index()]];
        const PointType& b = points_[origin_[seed.index() ^ 1]];
 
        // Twice the midpoint of the seed halfedge, and the normal pointing into
        // the face — the face is the one to the left of `seed`, and rotating a
        // direction a quarter turn counterclockwise points left. Keeping the
        // midpoint doubled leaves every quantity below a polynomial in the input
        // coordinates, so the only division is the final one.
        //
        // Those polynomials reach degree four, since ordering two hits along the
        // ray cross-multiplies two quadratic parameters, so they are formed in
        // @ref WideNumber. That holds every one of them for a coordinate type
        // whose promotion grows — and for a fixed-width one, coordinates up to
        // the fourth root of the widened range, which for `int` is about `2^28`.
        const auto wide = [](const NumberType& value) {
            return static_cast<WideNumber>(value);
        };
        const WideNumber midX = wide(a.x()) + wide(b.x());
        const WideNumber midY = wide(a.y()) + wide(b.y());
        const WideNumber normalX = wide(a.y()) - wide(b.y());
        const WideNumber normalY = wide(b.x()) - wide(a.x());
        const WideNumber two(2);
 
        // The first point where the inward ray meets the boundary, as the exact
        // fraction hitNumerator / hitDenominator with a positive denominator.
        WideNumber hitNumerator(0);
        WideNumber hitDenominator(0);
        const auto offer = [&](WideNumber numerator, WideNumber denominator) {
            if (denominator < WideNumber(0)) {
                numerator = -numerator;
                denominator = -denominator;
            }
            if (numerator <= WideNumber(0)) {
                return;  // behind the ray's start, or at it
            }
            if (hitDenominator == WideNumber(0) ||
                numerator * hitDenominator < hitNumerator * denominator) {
                hitNumerator = numerator;
                hitDenominator = denominator;
            }
        };
 
        forEachBoundaryHalfedge(f, [&](HalfedgeId h) {
            if (h.index() / 2 == seed.index() / 2) {
                return;  // the ray leaves its own edge, and meets it nowhere else
            }
            const PointType& p = points_[origin_[h.index()]];
            const PointType& q = points_[origin_[h.index() ^ 1]];
            const WideNumber edgeX = wide(q.x()) - wide(p.x());
            const WideNumber edgeY = wide(q.y()) - wide(p.y());
            // Twice the vector from the ray's start to the edge's first endpoint.
            const WideNumber toEdgeX = two * wide(p.x()) - midX;
            const WideNumber toEdgeY = two * wide(p.y()) - midY;
            const WideNumber denominator = normalX * edgeY - normalY * edgeX;
            const WideNumber alongEdge = toEdgeX * normalY - toEdgeY * normalX;
            if (denominator != WideNumber(0)) {
                // A proper crossing: the ray meets the edge's line at parameter
                // (toEdge x edge) / 2·denominator, inside the edge when the
                // parameter along the edge stays within [0, 1].
                WideNumber along = alongEdge;
                WideNumber scale = two * denominator;
                if (scale < WideNumber(0)) {
                    along = -along;
                    scale = -scale;
                }
                if (along < WideNumber(0) || along > scale) {
                    return;
                }
                offer(toEdgeX * edgeY - toEdgeY * edgeX, two * denominator);
            } else if (alongEdge == WideNumber(0)) {
                // The edge lies along the ray: it blocks it at whichever of its
                // endpoints comes first.
                const WideNumber squaredNormal = normalX * normalX + normalY * normalY;
                for (const PointType& endpoint : {p, q}) {
                    offer((two * wide(endpoint.x()) - midX) * normalX +
                              (two * wide(endpoint.y()) - midY) * normalY,
                          two * squaredNormal);
                }
            }
        });
 
        // A bounded face confines the ray, so it is always stopped.
        assert(hitDenominator != WideNumber(0));
        const ResultNumber scale = static_cast<ResultNumber>(two * hitDenominator);
        return Point<ResultNumber>(
            static_cast<ResultNumber>(midX * hitDenominator + hitNumerator * normalX) / scale,
            static_cast<ResultNumber>(midY * hitDenominator + hitNumerator * normalY) / scale);
    }
 
public:
    // -------------------------------------------------------------------------
    // Faces, continued

    [[nodiscard]] bool isUnbounded() const {
        return infinity_.valid();
    }

    [[nodiscard]] bool isUnbounded(HalfedgeId h) const {
        assert(h.valid() && h.index() < origin_.size());
        return infinity_.valid() &&
               (source(h) == infinity_ || target(h) == infinity_);
    }

    [[nodiscard]] bool isUnbounded(FaceId f) const {
        assert(f.valid() && f.index() < outerCycle_.size());
        return unboundedFace_[f.index()];
    }

    [[nodiscard]] bool isFictitious(VertexId v) const {
        assert(v.valid() && v.index() < topologicalVertexCount());
        return infinity_.valid() && v == infinity_;
    }

    [[nodiscard]] HalfedgeId outerCycle(FaceId f) const {
        assert(f.valid() && f.index() < outerCycle_.size());
        return outerCycle_[f.index()];
    }

    [[nodiscard]] std::span<const HalfedgeId> innerCycles(FaceId f) const {
        assert(f.valid() && f.index() + 1 < innerOffset_.size());
        const std::size_t from = innerOffset_[f.index()];
        const std::size_t to = innerOffset_[f.index() + 1];
        return std::span<const HalfedgeId>(innerCycle_.data() + from, to - from);
    }

    [[nodiscard]] std::vector<HalfedgeId> boundaryOf(FaceId f) const {
        assert(f.valid() && f.index() < outerCycle_.size());
        std::vector<HalfedgeId> boundary;
        const auto walkCycle = [&](HalfedgeId start) {
            HalfedgeId h = start;
            do {
                boundary.push_back(h);
                h = next(h);
            } while (h != start);
        };
        if (outerCycle_[f.index()].valid()) {
            walkCycle(outerCycle_[f.index()]);
        }
        for (HalfedgeId inner : innerCycles(f)) {
            walkCycle(inner);
        }
        return boundary;
    }

    [[nodiscard]] std::vector<HalfedgeId> outerBoundaryOf(FaceId f) const {
        assert(f.valid() && f.index() < outerCycle_.size());
        std::vector<HalfedgeId> boundary;
        const HalfedgeId start = outerCycle_[f.index()];
        if (!start.valid()) {
            return boundary;
        }
        HalfedgeId h = start;
        do {
            boundary.push_back(h);
            h = next(h);
        } while (h != start);
        return boundary;
    }

    [[nodiscard]] std::vector<std::vector<HalfedgeId>> innerBoundariesOf(FaceId f) const {
        assert(f.valid() && f.index() < outerCycle_.size());
        std::vector<std::vector<HalfedgeId>> boundaries;
        boundaries.reserve(innerCycles(f).size());
        for (HalfedgeId start : innerCycles(f)) {
            std::vector<HalfedgeId>& boundary = boundaries.emplace_back();
            HalfedgeId h = start;
            do {
                boundary.push_back(h);
                h = next(h);
            } while (h != start);
        }
        return boundaries;
    }

    template <class ResultNumber = NumberType>
    [[nodiscard]] PolygonWithHoles<Point<ResultNumber>> polygonWithHoles(FaceId f) const {
        if (!f.valid() || f.index() >= outerCycle_.size() || isUnbounded(f)) {
            throw std::logic_error(
                "Arrangement::polygonWithHoles is only defined for bounded faces");
        }
        using ExactPolygon = Polygon<PointType>;
 
        std::vector<ExactPolygon> rings;
        collectRings(cycleRing(outerCycle_[f.index()]), rings);
        assert(!rings.empty());
        // An outer cycle that pinches shut comes apart into several rings, and
        // the one holding the rest — the largest, since they are nested — is the
        // outer boundary.
        const auto largest =
            std::max_element(rings.begin(), rings.end(),
                             [](const ExactPolygon& left, const ExactPolygon& right) {
                                 return left.twiceArea() < right.twiceArea();
                             });
        if (largest != rings.begin()) {
            std::iter_swap(rings.begin(), largest);
        }
        std::vector<ExactPolygon> holes(rings.begin() + 1, rings.end());
        for (HalfedgeId inner : innerCycles(f)) {
            collectRings(cycleRing(inner), holes);
        }
        const PolygonWithHoles<PointType> exact(std::move(rings.front()), std::move(holes));
        return PolygonWithHoles<Point<ResultNumber>>(exact);
    }

    template <class ResultNumber = NumberType>
    [[nodiscard]] HalfplaneIntersection<Point<ResultNumber>>
    halfplaneIntersection(FaceId f) const {
        if (!f.valid() || f.index() >= outerCycle_.size()) {
            throw std::logic_error(
                "Arrangement::halfplaneIntersection requires a valid face");
        }
 
        using ResultPoint = Point<ResultNumber>;
        using ResultHalfplane = Halfplane<ResultPoint>;
        std::vector<ResultHalfplane> halfplanes;
        const auto appendCycle = [&](HalfedgeId start) {
            HalfedgeId h = start;
            do {
                if (face(twin(h)) != f) {
                    const EdgeGeometry& geometry = edgeGeometry_[h.index() / 2];
                    const bool forward = h.index() % 2 == 0;
                    halfplanes.emplace_back(
                        ResultPoint(forward ? geometry.a : geometry.b),
                        ResultPoint(forward ? geometry.b : geometry.a));
                }
                h = next(h);
            } while (h != start);
        };
 
        if (outerCycle_[f.index()].valid()) {
            appendCycle(outerCycle_[f.index()]);
        } else {
            for (HalfedgeId start : innerCycles(f)) {
                bool reachesInfinity = false;
                HalfedgeId h = start;
                do {
                    reachesInfinity = reachesInfinity || isUnbounded(h);
                    h = next(h);
                } while (h != start);
                if (reachesInfinity) {
                    appendCycle(start);
                }
            }
        }
        return HalfplaneIntersection<ResultPoint>(std::move(halfplanes));
    }
 
    // -------------------------------------------------------------------------
    // Labels and history

    [[nodiscard]] const TLabel& label(HalfedgeId h) const {
        assert(h.valid() && h.index() < origin_.size());
        return edgeLabel_[h.index() / 2];
    }

    [[nodiscard]] TLabel& label(HalfedgeId h) {
        assert(h.valid() && h.index() < origin_.size());
        return edgeLabel_[h.index() / 2];
    }

    [[nodiscard]] const TLabel& label(FaceId f) const {
        assert(f.valid() && f.index() < faceLabel_.size());
        return faceLabel_[f.index()];
    }

    [[nodiscard]] TLabel& label(FaceId f) {
        assert(f.valid() && f.index() < faceLabel_.size());
        return faceLabel_[f.index()];
    }

    [[nodiscard]] std::span<const std::uint32_t> originsOf(HalfedgeId h) const {
        assert(h.valid() && h.index() < origin_.size());
        const std::size_t edge = h.index() / 2;
        const std::size_t from = originOffset_[edge];
        const std::size_t to = originOffset_[edge + 1];
        return std::span<const std::uint32_t>(originIndex_.data() + from, to - from);
    }

    [[nodiscard]] std::vector<std::uint32_t> originsOf(VertexId v) const {
        assert(v.valid() && v.index() < outgoing_.size());
        std::vector<std::uint32_t> origins;
        const HalfedgeId start = outgoing_[v.index()];
        if (!start.valid()) {
            return origins;
        }
        HalfedgeId h = start;
        do {
            const std::span<const std::uint32_t> edgeOrigins = originsOf(h);
            origins.insert(origins.end(), edgeOrigins.begin(), edgeOrigins.end());
            h = next(twin(h));
        } while (h != start);
        std::sort(origins.begin(), origins.end());
        origins.erase(std::unique(origins.begin(), origins.end()), origins.end());
        return origins;
    }
 
    // -------------------------------------------------------------------------
    // Location

    void buildPointLocation();

    template <class UniformRandomBitGenerator>
    void buildPointLocation(UniformRandomBitGenerator&& generator);

    void clearPointLocation() noexcept {
        pointLocation_.reset();
    }

    [[nodiscard]] bool hasPointLocation() const noexcept {
        return static_cast<bool>(pointLocation_);
    }

    [[nodiscard]] FaceId locateFace(const PointType& p) const;

    [[nodiscard]] CellId locateCell(const PointType& p) const;
 
    // -------------------------------------------------------------------------
    // Intersection traversal

    template <class Q, class Fn>
        requires (OrientedSegmentConcept<Q> || OrientedLineConcept<Q> ||
                  RayConcept<Q> || MonotoneChainConcept<Q> || PolylineConcept<Q>)
    bool visitIntersecting(const Q& r, Fn fn) const;

    template <class Q>
        requires (OrientedSegmentConcept<Q> || OrientedLineConcept<Q> ||
                  RayConcept<Q> || MonotoneChainConcept<Q> || PolylineConcept<Q>)
    [[nodiscard]] std::vector<IntersectionId> reportIntersecting(const Q& r) const {
        std::vector<IntersectionId> result;
        visitIntersecting(r, [&](const IntersectionId& id) { result.push_back(id); });
        return result;
    }

    template <class Q>
        requires (OrientedSegmentConcept<Q> || OrientedLineConcept<Q> ||
                  RayConcept<Q> || MonotoneChainConcept<Q> || PolylineConcept<Q>)
    [[nodiscard]] std::optional<IntersectionId> firstIntersecting(const Q& r) const {
        std::optional<IntersectionId> result;
        visitIntersecting(r, [&](const IntersectionId& id) {
            result = id;
            return true;
        });
        return result;
    }

    template <class Q>
        requires (OrientedSegmentConcept<Q> || OrientedLineConcept<Q> ||
                  RayConcept<Q> || MonotoneChainConcept<Q> || PolylineConcept<Q>)
    [[nodiscard]] bool emptyIntersecting(const Q& r) const {
        return !visitIntersecting(r, [](const IntersectionId&) { return true; });
    }
 
private:
    // Two promotions above the coordinates, which is what it takes to hold a
    // product of three coordinate differences exactly — one degree more than the
    // sign predicates promote for. It is the width the geometry that goes beyond
    // a sign works in: the crossing abscissae @ref halfedgeLeftOf orders, the
    // ordinate @ref pointAt cuts a carrier at, and the ray parameters
    // @ref sweptWitness compares, that last one quartic and so the one place a
    // fixed-width coordinate type is bounded further. A type that grows to hold
    // its values, or an approximate one, is its own promotion and pays nothing.
    using WideNumber =
        detail::promoted_number_t<detail::promoted_number_t<NumberType>>;
 
    // Point location without the index: the face is the one left of the nearest
    // edge to the west, and a point with nothing to its west lies in the face
    // the boundary at infinity opens onto.
    [[nodiscard]] FaceId locateFaceLinear(const PointType& p) const {
        const HalfedgeId h = halfedgeLeftOf(p);
        if (h.valid()) {
            return face(h);
        }
        if (infinity_.valid()) {
            return face(infinityBoundaryAtWest(p));
        }
        return FaceId(0);
    }
 
    // The same without the index, but naming the cell the point lies on rather
    // than the face around it: a vertex wins over an edge, and an edge over the
    // face, since a point on the boundary belongs to the lower-dimensional cell.
    [[nodiscard]] CellId locateCellLinear(const PointType& p) const {
        for (std::uint32_t v = 0; v < points_.size(); ++v) {
            if (points_[v] == p) {
                return VertexId(v);
            }
        }
        for (std::uint32_t edge = 0; edge < edgeGeometry_.size(); ++edge) {
            const EdgeGeometry& geometry = edgeGeometry_[edge];
            if (orientationSign(geometry.a, geometry.b, p) != 0) {
                continue;
            }
            bool contains = geometry.kind == EdgeKind::line;
            if (geometry.kind == EdgeKind::segment) {
                contains = p.x() >= std::min(geometry.a.x(), geometry.b.x()) &&
                           p.x() <= std::max(geometry.a.x(), geometry.b.x()) &&
                           p.y() >= std::min(geometry.a.y(), geometry.b.y()) &&
                           p.y() <= std::max(geometry.a.y(), geometry.b.y());
            } else if (geometry.kind == EdgeKind::ray) {
                // Ahead of the source along the ray, which on the carrier is the
                // whole of it. The predicate promotes the product of the two
                // differences, which a raw multiplication of coordinates would
                // not.
                contains = dotSign(geometry.a, p, geometry.a, geometry.b) >= 0;
            }
            if (contains) {
                return HalfedgeId(2 * edge);
            }
        }
        return locateFaceLinear(p);
    }
 
    // Vertices as the halfedge structure counts them, which is the geometric
    // vertices plus the single point at infinity when the arrangement has one.
    [[nodiscard]] std::size_t topologicalVertexCount() const {
        return points_.size() + (infinity_.valid() ? 1 : 0);
    }
 
    // How far an edge runs: a segment is bounded at both ends, a ray only at
    // the first, and a line at neither.
    enum class EdgeKind : std::uint8_t { segment, ray, line };
 
    // The defining coordinates of an edge, kept beside the topology so a
    // predicate never has to walk the halfedge structure to rebuild them. For a
    // ray, `a` is the source and `b` only fixes the direction; for a line the
    // two merely span it.
    struct EdgeGeometry {
        EdgeKind kind;
        PointType a;
        PointType b;
    };
 
    // One straight piece of a query, reported in order along the piece. The
    // seen vectors carry the cells a chain has already reported across its
    // pieces, and are null for a query that is a single piece.
    template <class Q, class Fn>
    bool visitStraightIntersecting(const Q& piece, Fn& fn,
                                   std::vector<bool>* seenVertices,
                                   std::vector<bool>* seenEdges) const;
 
    // A query that degenerated to a single point, which meets at most one cell.
    template <PointConcept Q, class Fn>
    bool visitPointIntersecting(const Q& point, Fn& fn,
                                std::vector<bool>* seenVertices,
                                std::vector<bool>* seenEdges) const;
 
    // An input segment, with the position of the shape it came from.
    struct InputSegment {
        Segment<PointType> segment;
        std::uint32_t origin;
        [[no_unique_address]] TLabel label;
    };
 
    // One piece of an input segment after splitting, before twin pieces of
    // overlapping input are merged into a single edge.
    struct Piece {
        PointType a;
        PointType b;
        std::uint32_t origin;
        [[no_unique_address]] TLabel label;
    };
 
    // A segment, ray, or line before collinear inputs are overlaid.
    struct InputCurve {
        EdgeKind kind;
        PointType a;
        PointType b;
        std::uint32_t origin;
        [[no_unique_address]] TLabel label;
    };
 
    // -------------------------------------------------------------------------
    // Input normalization
 
    // Splits a shape range into the segments it contributes and the points it
    // reduces to, keeping each shape's position in the range as its origin.
    template <class ShapeRange>
    static void collect(const ShapeRange& shapes, std::vector<InputSegment>& segments,
                        std::vector<InputCurve>& curves, std::vector<PointType>& isolated) {
        std::uint32_t index = 0;
        for (const auto& shape : shapes) {
            append(shape, index, segments, curves, isolated);
            ++index;
        }
    }
 
    // The shapes the arrangement can be built from directly or through Shape.
    template <class InputShape>
    static constexpr bool isSupported =
        detail::is_empty_shape_v<InputShape> || detail::is_point_v<InputShape> ||
        detail::is_segment_v<InputShape> || detail::is_oriented_segment_v<InputShape> ||
        detail::is_line_v<InputShape> || detail::is_oriented_line_v<InputShape> ||
        detail::is_ray_v<InputShape> ||
        detail::is_polyline_v<InputShape> || detail::is_monotone_chain_v<InputShape> ||
        detail::is_triangle_v<InputShape> || detail::is_rectangle_v<InputShape> ||
        detail::is_convex_v<InputShape> || detail::is_polygon_v<InputShape> ||
        detail::is_polygon_with_holes_v<InputShape>;
 
    // Appends the cut segments — and the isolated points — of one input shape.
    template <class InputShape>
    static void append(const InputShape& shape, std::uint32_t index,
                       std::vector<InputSegment>& segments, std::vector<InputCurve>& curves,
                       std::vector<PointType>& isolated) {
        static_assert(isSupported<InputShape> || detail::is_shape_v<InputShape>,
                      "Arrangement accepts points, segment-bounded shapes, lines, and rays");
        const auto addSegment = [&](const auto& edge) {
            const PointType a(edge[0]);
            const PointType b(edge[1]);
            if (a == b) {
                isolated.push_back(a);
                return;
            }
            InputSegment input{Segment<PointType>(a, b), index, TLabel{}};
            if constexpr (detail::has_label_v<TLabel>) {
                input.label = detail::copyLabel<TLabel>(shape);
            }
            segments.push_back(std::move(input));
            InputCurve curve{EdgeKind::segment, a < b ? a : b, a < b ? b : a, index, TLabel{}};
            if constexpr (detail::has_label_v<TLabel>) {
                curve.label = detail::copyLabel<TLabel>(shape);
            }
            curves.push_back(std::move(curve));
        };
 
        if constexpr (detail::is_shape_v<InputShape>) {
            std::visit([&](const auto& alternative) {
                if constexpr (isSupported<std::remove_cvref_t<decltype(alternative)>>) {
                    append(alternative, index, segments, curves, isolated);
                } else {
                    throw std::invalid_argument(
                        "Arrangement accepts only points, segment-bounded shapes, lines, and rays");
                }
            }, shape.variant());
        } else if constexpr (detail::is_empty_shape_v<InputShape>) {
            (void)shape;
        } else if constexpr (detail::is_point_v<InputShape>) {
            isolated.emplace_back(shape);
        } else if constexpr (detail::is_segment_v<InputShape> ||
                             detail::is_oriented_segment_v<InputShape>) {
            addSegment(shape);
        } else if constexpr (detail::is_line_v<InputShape> ||
                             detail::is_oriented_line_v<InputShape> ||
                             detail::is_ray_v<InputShape>) {
            const PointType a(shape[0]);
            const PointType b(shape[1]);
            if (a == b) {
                isolated.push_back(a);
                return;
            }
            constexpr EdgeKind kind = detail::is_ray_v<InputShape> ? EdgeKind::ray : EdgeKind::line;
            InputCurve curve{kind, a, b, index, TLabel{}};
            if constexpr (detail::has_label_v<TLabel>) {
                curve.label = detail::copyLabel<TLabel>(shape);
            }
            curves.push_back(std::move(curve));
        } else if constexpr (detail::is_polygon_with_holes_v<InputShape>) {
            for (const auto& edge : shape.edges()) {
                addSegment(edge);
            }
        } else if constexpr (detail::is_polyline_v<InputShape> ||
                             detail::is_monotone_chain_v<InputShape>) {
            if (shape.size() == 1) {
                isolated.emplace_back(shape[0]);
            } else {
                for (const auto& edge : shape.edgesView()) {
                    addSegment(edge);
                }
            }
        } else if constexpr (requires { shape.edgesView(); }) {
            for (const auto& edge : shape.edgesView()) {
                addSegment(edge);
            }
        } else {
            for (const auto& edge : shape.edges()) {
                addSegment(edge);
            }
        }
    }
 
    // -------------------------------------------------------------------------
    // Construction
 
    void build(std::vector<InputSegment>& segments, std::vector<InputCurve>& curves,
               std::vector<PointType>& isolated, bool simpleBoundaries) {
        const bool hasUnbounded = std::ranges::any_of(
            curves, [](const InputCurve& curve) { return curve.kind != EdgeKind::segment; });
        if (hasUnbounded) {
            buildUnbounded(curves, isolated);
            return;
        }
        std::vector<Piece> pieces = split(segments, isolated, simpleBoundaries);
        internVertices(pieces, isolated);
        simplifyStoredCoordinates();
        syncVertexApproximations();
        wireHalfedges();
        buildFaces();
    }
 
    // A stretch of one carrier that some input curve covers, as a parameter
    // range along that carrier. An absent bound runs to infinity, and `origin`
    // names the input shape the stretch came from.
    struct CarrierInterval {
        std::optional<NumberType> low;
        std::optional<NumberType> high;
        std::uint32_t origin;
        [[no_unique_address]] TLabel label;
    };
 
    // A supporting line, with every input stretch laid over it. Collinear input
    // shares one carrier, so overlap is settled by merging intervals instead of
    // by intersecting curves. `usesX` picks the coordinate that parameterizes
    // it, which is the abscissa unless the line is vertical. `cuts` holds the
    // parameters where the carrier must be split.
    struct Carrier {
        PointType a;
        PointType b;
        bool usesX;
        std::vector<CarrierInterval> intervals;
        std::vector<NumberType> cuts;
    };
 
    // A piece of a carrier between consecutive cuts, so no other curve crosses
    // its interior. These are the pieces the halfedge structure is wired from,
    // and `origins` records every input shape covering the piece.
    struct AtomicCurve {
        std::uint32_t carrier;
        std::optional<NumberType> low;
        std::optional<NumberType> high;
        std::vector<std::uint32_t> origins;
        [[no_unique_address]] TLabel label;
    };
 
    // Where a point on a carrier falls along it. Vertical carriers are
    // parameterized by ordinate, everything else by abscissa, so the
    // parameter is always strictly monotone along the carrier.
    static NumberType parameterOf(const Carrier& carrier, const PointType& point) {
        return carrier.usesX ? point.x() : point.y();
    }
 
    // Whether an interval already spans the whole of another one, absent bounds
    // reaching infinity. This is what makes a repeated or contained input
    // stretch add nothing.
    static bool covers(const CarrierInterval& interval,
                       const std::optional<NumberType>& low,
                       const std::optional<NumberType>& high) {
        if (!low.has_value()) {
            if (interval.low.has_value()) {
                return false;
            }
        } else if (interval.low.has_value() && *low < *interval.low) {
            return false;
        }
        if (!high.has_value()) {
            if (interval.high.has_value()) {
                return false;
            }
        } else if (interval.high.has_value() && *interval.high < *high) {
            return false;
        }
        return true;
    }
 
    // Whether an interval holds a single parameter, ends included.
    static bool covers(const CarrierInterval& interval, const NumberType& value) {
        return (!interval.low.has_value() || !(value < *interval.low)) &&
               (!interval.high.has_value() || !(*interval.high < value));
    }
 
    // The point at a parameter along a carrier.
    static PointType pointAt(const Carrier& carrier, const NumberType& parameter) {
        PointType point = [&] {
            if (carrier.usesX) {
                // The ordinate reached at that abscissa, as a displacement from
                // the carrier's own: the product of two coordinate differences
                // outgrows a narrow coordinate type long before the quotient
                // does, and the quotient is what an arrangement whose vertices
                // are representable has to land on.
                const auto wide = [](const NumberType& value) {
                    return static_cast<WideNumber>(value);
                };
                const WideNumber dx = wide(carrier.b.x()) - wide(carrier.a.x());
                const WideNumber dy = wide(carrier.b.y()) - wide(carrier.a.y());
                const WideNumber rise =
                    (wide(parameter) - wide(carrier.a.x())) * dy / dx;
                return PointType(parameter,
                                 carrier.a.y() + detail::asNumber<NumberType>(rise));
            }
            return PointType(carrier.a.x(), parameter);
        }();
        // A point cut out of a carrier becomes an interned vertex, so it is
        // hashed and compared the moment it is built and read for the rest of the
        // arrangement's life. Reducing the chain of arithmetic above once, here,
        // is what puts all of that on the normalized fast path.
        if constexpr (pgl::is_Rational_v<NumberType>) {
            point.x().simplify();
            point.y().simplify();
        }
        return point;
    }
 
    // Whether a curve is collinear with a carrier, and so belongs on it rather
    // than on one of its own.
    static bool sameCarrier(const Carrier& carrier, const InputCurve& curve) {
        return orientationSign(carrier.a, carrier.b, curve.a) == 0 &&
               orientationSign(carrier.a, carrier.b, curve.b) == 0;
    }
 
    // General normalization for input containing a ray or a line. Collinear
    // inputs are overlaid as intervals on one exact carrier; intersections of
    // distinct carriers then become additional interval endpoints.
    void buildUnbounded(const std::vector<InputCurve>& curves,
                        const std::vector<PointType>& isolated) {
        std::vector<Carrier> carriers;
        for (const InputCurve& curve : curves) {
            auto found = std::find_if(carriers.begin(), carriers.end(),
                                      [&](const Carrier& carrier) {
                                          return sameCarrier(carrier, curve);
                                      });
            if (found == carriers.end()) {
                PointType a = curve.a;
                PointType b = curve.b;
                if (b < a) {
                    std::swap(a, b);
                }
                carriers.push_back(Carrier{a, b, !(a.x() == b.x()), {}, {}});
                found = std::prev(carriers.end());
            }
            Carrier& carrier = *found;
            const NumberType ta = parameterOf(carrier, curve.a);
            const NumberType tb = parameterOf(carrier, curve.b);
            CarrierInterval interval{{}, {}, curve.origin, curve.label};
            if (curve.kind == EdgeKind::segment) {
                interval.low = std::min(ta, tb);
                interval.high = std::max(ta, tb);
            } else if (curve.kind == EdgeKind::ray) {
                if (ta < tb) {
                    interval.low = ta;
                } else {
                    interval.high = ta;
                }
            }
            if (interval.low.has_value()) {
                carrier.cuts.push_back(*interval.low);
            }
            if (interval.high.has_value()) {
                carrier.cuts.push_back(*interval.high);
            }
            carrier.intervals.push_back(std::move(interval));
        }
 
        for (Carrier& carrier : carriers) {
            for (const PointType& point : isolated) {
                if (orientationSign(carrier.a, carrier.b, point) != 0) {
                    continue;
                }
                const NumberType value = parameterOf(carrier, point);
                if (std::ranges::any_of(carrier.intervals,
                                        [&](const CarrierInterval& interval) {
                                            return covers(interval, value);
                                        })) {
                    carrier.cuts.push_back(value);
                }
            }
        }
 
        for (std::size_t i = 0; i < carriers.size(); ++i) {
            for (std::size_t j = i + 1; j < carriers.size(); ++j) {
                const Line<PointType> first(carriers[i].a, carriers[i].b);
                const Line<PointType> second(carriers[j].a, carriers[j].b);
                const auto intersection = first.template intersection<NumberType>(second);
                if (!intersection || !std::holds_alternative<PointType>(*intersection)) {
                    continue;
                }
                const PointType& point = std::get<PointType>(*intersection);
                const NumberType ti = parameterOf(carriers[i], point);
                const NumberType tj = parameterOf(carriers[j], point);
                const bool onFirst = std::ranges::any_of(
                    carriers[i].intervals,
                    [&](const CarrierInterval& interval) { return covers(interval, ti); });
                const bool onSecond = std::ranges::any_of(
                    carriers[j].intervals,
                    [&](const CarrierInterval& interval) { return covers(interval, tj); });
                if (onFirst && onSecond) {
                    carriers[i].cuts.push_back(ti);
                    carriers[j].cuts.push_back(tj);
                }
            }
        }
 
        std::vector<AtomicCurve> atoms;
        for (std::uint32_t c = 0; c < carriers.size(); ++c) {
            Carrier& carrier = carriers[c];
            // A cut coming from a carrier crossing is an unreduced fraction, and
            // it is read by the sort and unique here, by the interval tests in
            // emit, and again by pointAt once it reaches an atom. One gcd apiece
            // covers all of it; the same reasoning as split's cut lists.
            if constexpr (pgl::is_Rational_v<NumberType>) {
                for (NumberType& cut : carrier.cuts) {
                    cut.simplify();
                }
            }
            std::sort(carrier.cuts.begin(), carrier.cuts.end());
            carrier.cuts.erase(std::unique(carrier.cuts.begin(), carrier.cuts.end()),
                               carrier.cuts.end());
 
            const auto emit = [&](std::optional<NumberType> low,
                                  std::optional<NumberType> high) {
                const CarrierInterval* first = nullptr;
                std::vector<std::uint32_t> origins;
                for (const CarrierInterval& interval : carrier.intervals) {
                    if (!covers(interval, low, high)) {
                        continue;
                    }
                    if (first == nullptr || interval.origin < first->origin) {
                        first = &interval;
                    }
                    origins.push_back(interval.origin);
                }
                if (first == nullptr) {
                    return;
                }
                std::sort(origins.begin(), origins.end());
                origins.erase(std::unique(origins.begin(), origins.end()), origins.end());
                atoms.push_back(AtomicCurve{c, std::move(low), std::move(high),
                                            std::move(origins), first->label});
            };
 
            if (carrier.cuts.empty()) {
                emit({}, {});
                continue;
            }
            emit({}, carrier.cuts.front());
            for (std::size_t i = 0; i + 1 < carrier.cuts.size(); ++i) {
                if (!(carrier.cuts[i] == carrier.cuts[i + 1])) {
                    emit(carrier.cuts[i], carrier.cuts[i + 1]);
                }
            }
            emit(carrier.cuts.back(), {});
        }
 
        std::unordered_map<PointType, std::uint32_t> vertexOf;
        const auto idOf = [&](const PointType& point) {
            const auto found = vertexOf.find(point);
            if (found != vertexOf.end()) {
                return found->second;
            }
            const auto id = static_cast<std::uint32_t>(points_.size());
            points_.push_back(point);
            vertexOf.emplace(point, id);
            return id;
        };
 
        for (const AtomicCurve& atom : atoms) {
            if (atom.low.has_value()) {
                idOf(pointAt(carriers[atom.carrier], *atom.low));
            }
            if (atom.high.has_value()) {
                idOf(pointAt(carriers[atom.carrier], *atom.high));
            }
        }
        for (const PointType& point : isolated) {
            idOf(point);
        }
        infinity_ = VertexId(static_cast<std::uint32_t>(points_.size()));
 
        originOffset_.push_back(0);
        for (const AtomicCurve& atom : atoms) {
            const Carrier& carrier = carriers[atom.carrier];
            if (atom.low.has_value() && atom.high.has_value()) {
                const PointType a = pointAt(carrier, *atom.low);
                const PointType b = pointAt(carrier, *atom.high);
                origin_.push_back(idOf(a));
                origin_.push_back(idOf(b));
                edgeGeometry_.push_back({EdgeKind::segment, a, b});
            } else if (atom.low.has_value() || atom.high.has_value()) {
                const bool increasing = atom.low.has_value();
                const PointType sourcePoint = pointAt(
                    carrier, increasing ? *atom.low : *atom.high);
                const NumberType dx = carrier.b.x() - carrier.a.x();
                const NumberType dy = carrier.b.y() - carrier.a.y();
                const PointType directionPoint(
                    increasing ? sourcePoint.x() + dx : sourcePoint.x() - dx,
                    increasing ? sourcePoint.y() + dy : sourcePoint.y() - dy);
                origin_.push_back(idOf(sourcePoint));
                origin_.push_back(infinity_.index());
                edgeGeometry_.push_back({EdgeKind::ray, sourcePoint, directionPoint});
            } else {
                origin_.push_back(infinity_.index());
                origin_.push_back(infinity_.index());
                edgeGeometry_.push_back({EdgeKind::line, carrier.a, carrier.b});
            }
            edgeLabel_.push_back(atom.label);
            originIndex_.insert(originIndex_.end(), atom.origins.begin(), atom.origins.end());
            originOffset_.push_back(static_cast<std::uint32_t>(originIndex_.size()));
        }
 
        simplifyStoredCoordinates();
        syncVertexApproximations();
        next_.assign(origin_.size(), 0);
        face_.assign(origin_.size(), 0);
        outgoing_.assign(topologicalVertexCount(), HalfedgeId());
        wireHalfedgesUnbounded();
        buildFacesUnbounded();
    }

    static std::vector<Piece> split(std::vector<InputSegment>& segments,
                                    const std::vector<PointType>& isolated,
                                    bool simpleBoundaries) {
        using IntegralPoint = Point<std::int64_t>;
        using IntegralSegment = Segment<IntegralPoint>;
        constexpr bool mayNeedIntegralNarrowing =
            is_Rational_v<NumberType> || std::same_as<NumberType, BigInt>;
 
        // Equal geometry adjacent, and within a group the contributing shapes in
        // the order internVertices expects to see them.
        std::sort(segments.begin(), segments.end(),
                  [](const InputSegment& left, const InputSegment& right) {
                      if (!(left.segment == right.segment)) {
                          return left.segment < right.segment;
                      }
                      return left.origin < right.origin;
                  });
        // One entry per group plus a trailing sentinel, so group g occupies
        // segments[group[g] .. group[g + 1]).
        std::vector<std::size_t> group;
        for (std::size_t i = 0; i < segments.size(); ++i) {
            if (i == 0 || !(segments[i].segment == segments[i - 1].segment)) {
                group.push_back(i);
            }
        }
        group.push_back(segments.size());
        const std::size_t count = group.size() - 1;
 
        // A segment's endpoints are in lexicographic order, so min().x() is the
        // left end of its x-projection and max().x() the right one, with nothing
        // to compute. The y-extent is not ordered, hence the explicit minmax.
        std::vector<NumberType> right, low, high;
        right.reserve(count);
        low.reserve(count);
        high.reserve(count);
        for (std::size_t i = 0; i < count; ++i) {
            const Segment<PointType>& current = segments[group[i]].segment;
            right.push_back(current.max().x());
            const auto [lo, hi] = std::minmax(current.min().y(), current.max().y());
            low.push_back(lo);
            high.push_back(hi);
        }
 
        std::vector<std::uint32_t> order(count);
        for (std::size_t i = 0; i < count; ++i) {
            order[i] = static_cast<std::uint32_t>(i);
        }
        std::sort(order.begin(), order.end(), [&](std::uint32_t a, std::uint32_t b) {
            return segments[group[a]].segment.min().x() < segments[group[b]].segment.min().x();
        });
 
        std::vector<std::vector<PointType>> cuts(count);
        std::vector<std::optional<IntegralSegment>> integral;
        if constexpr (mayNeedIntegralNarrowing) {
            integral.resize(count);
        }
        for (std::size_t i = 0; i < count; ++i) {
            const Segment<PointType>& segment = segments[group[i]].segment;
            cuts[i].push_back(segment.min());
            cuts[i].push_back(segment.max());
 
            // A rational arrangement is commonly fed lattice segments — in
            // particular, every convex-piece Minkowski sum keeps integral
            // vertices until this overlay creates its crossings. Remember that
            // narrower spelling so the many rejected segment pairs use native
            // int64 coordinates and int128 determinants instead of promoting
            // Rational<BigInt> predicates all the way to BigInt. The exact
            // intersection is still constructed in NumberType below.
            const auto asIntegral = [](const PointType& point)
                -> std::optional<IntegralPoint> {
                const auto store = [](const auto& x, const auto& y)
                    -> std::optional<IntegralPoint> {
                    if (!detail::representableAs<std::int64_t>(x) ||
                        !detail::representableAs<std::int64_t>(y)) {
                        return std::nullopt;
                    }
                    return IntegralPoint(detail::narrowTo<std::int64_t>(x),
                                         detail::narrowTo<std::int64_t>(y));
                };
 
                if constexpr (is_Rational_v<NumberType>) {
                    if (!point.x().isInteger() || !point.y().isInteger()) {
                        return std::nullopt;
                    }
                    using Integer = rational_int_t<NumberType>;
                    return store(Integer(static_cast<Integer>(point.x())),
                                 Integer(static_cast<Integer>(point.y())));
                } else if constexpr (detail::extended_integral<NumberType> ||
                                     std::same_as<NumberType, BigInt>) {
                    return store(point.x(), point.y());
                } else {
                    return std::nullopt;
                }
            };
            if constexpr (mayNeedIntegralNarrowing) {
                const auto a = asIntegral(segment.min());
                const auto b = asIntegral(segment.max());
                if (a && b) {
                    integral[i].emplace(*a, *b);
                }
            }
        }
 
        const auto meet = [&](std::size_t a, std::size_t b) {
            const auto add = [&](const auto& piece) {
                if (!piece) {
                    return;
                }
                if (const auto* point = std::get_if<0>(&*piece)) {
                    cuts[a].emplace_back(*point);
                    cuts[b].emplace_back(*point);
                } else {
                    const auto& overlap = std::get<1>(*piece);
                    for (const auto& end : {overlap.min(), overlap.max()}) {
                        cuts[a].emplace_back(end);
                        cuts[b].emplace_back(end);
                    }
                }
            };
 
            if constexpr (mayNeedIntegralNarrowing) {
                if (integral[a] && integral[b]) {
                    add(integral[a]->template intersection<NumberType>(*integral[b]));
                    return;
                }
            }
            add(segments[group[a]].segment.template intersection<NumberType>(
                segments[group[b]].segment));
        };
 
        // For two large simple boundaries the arrangement normally has only a
        // small number of red-blue crossings, while the edges of each jagged
        // ring have long, mutually overlapping x-projections. The box sweep
        // below then examines quadratically many pairs it can reject only with
        // an exact predicate. Bentley--Ottmann follows the actual crossings and
        // makes that case O((n + k) log n). Keep the box sweep for many-shape and
        // dense overlays: its compact vectors and native-integer rejection are
        // substantially cheaper when k itself is large.
        const std::uint32_t originCount = segments.empty()
            ? 0
            : 1 + std::ranges::max(segments, {}, &InputSegment::origin).origin;
        const bool useIntersectionSweep = simpleBoundaries && originCount == 2 && count >= 256 &&
                                          !std::floating_point<NumberType>;
        if (useIntersectionSweep) {
            if constexpr (!std::floating_point<NumberType>) {
                const bool allIntegral = mayNeedIntegralNarrowing &&
                    std::ranges::all_of(
                        integral, [](const auto& segment) { return segment.has_value(); });
                if (allIntegral) {
                    std::vector<IntegralSegment> unique;
                    unique.reserve(count);
                    for (const auto& segment : integral) {
                        unique.push_back(*segment);
                    }
                    detail::BentleyOttmann<NumberType, IntegralSegment> sweep;
                    for (const auto& pair : sweep.findIntersections(unique)) {
                        const auto indexOf = [&](const IntegralSegment& segment) {
                            return static_cast<std::size_t>(
                                std::lower_bound(unique.begin(), unique.end(), segment) -
                                unique.begin());
                        };
                        const std::size_t a = indexOf(pair[0]);
                        const std::size_t b = indexOf(pair[1]);
                        if (a != b) {
                            meet(a, b);
                        }
                    }
                } else {
                    std::vector<Segment<PointType>> unique;
                    unique.reserve(count);
                    for (std::size_t i = 0; i < count; ++i) {
                        unique.push_back(segments[group[i]].segment);
                    }
                    detail::BentleyOttmann<NumberType, Segment<PointType>> sweep;
                    for (const auto& pair : sweep.findIntersections(unique)) {
                        const auto indexOf = [&](const Segment<PointType>& segment) {
                            return static_cast<std::size_t>(
                                std::lower_bound(unique.begin(), unique.end(), segment) -
                                unique.begin());
                        };
                        const std::size_t a = indexOf(pair[0]);
                        const std::size_t b = indexOf(pair[1]);
                        if (a != b) {
                            meet(a, b);
                        }
                    }
                }
            }
        } else {
            // The active list is compacted by the same pass that tests it, so a
            // group is dropped exactly once and expiry costs nothing beyond the
            // comparison the test needed anyway.
            std::vector<std::uint32_t> active;
            for (const std::uint32_t current : order) {
                const NumberType& left = segments[group[current]].segment.min().x();
                std::size_t write = 0;
                for (std::size_t read = 0; read < active.size(); ++read) {
                    const std::uint32_t other = active[read];
                    bool expired;
                    bool missesInY;
                    if constexpr (mayNeedIntegralNarrowing) {
                        if (integral[current] && integral[other]) {
                            const IntegralSegment& currentSegment = *integral[current];
                            const IntegralSegment& otherSegment = *integral[other];
                            expired = otherSegment.max().x() < currentSegment.min().x();
                            const auto [currentLow, currentHigh] =
                                std::minmax(currentSegment.min().y(), currentSegment.max().y());
                            const auto [otherLow, otherHigh] =
                                std::minmax(otherSegment.min().y(), otherSegment.max().y());
                            missesInY = otherHigh < currentLow || currentHigh < otherLow;
                        } else {
                            expired = right[other] < left;
                            missesInY = high[other] < low[current] || high[current] < low[other];
                        }
                    } else {
                        expired = right[other] < left;
                        missesInY = high[other] < low[current] || high[current] < low[other];
                    }
                    if (expired) {
                        continue;  // its projection closed before this one opened
                    }
                    active[write++] = other;
                    if (missesInY) {
                        continue;  // boxes overlap in x but miss in y
                    }
                    meet(other, current);
                }
                active.resize(write);
                active.push_back(current);
            }
        }
 
        std::vector<Piece> pieces;
        for (std::size_t i = 0; i < count; ++i) {
            for (const PointType& point : isolated) {
                if (segments[group[i]].segment.contains(point)) {
                    cuts[i].push_back(point);
                }
            }
            // A crossing arrives from intersection() as a fraction whose
            // normalization this type defers, and every read of one reduces it
            // again and keeps nothing. Each cut is about to be read many times
            // over — by the sort and unique just below, by the piece copies they
            // feed, and by everything internVertices and the wiring passes do
            // with those. Reducing once per distinct cut, before any of that,
            // collapses all of it to a single gcd apiece.
            if constexpr (pgl::is_Rational_v<NumberType>) {
                for (PointType& cut : cuts[i]) {
                    cut.x().simplify();
                    cut.y().simplify();
                }
            }
            // Every cut lies on the segment, so the lexicographic point order is
            // the linear order along it.
            std::sort(cuts[i].begin(), cuts[i].end());
            cuts[i].erase(std::unique(cuts[i].begin(), cuts[i].end()), cuts[i].end());
            for (std::size_t k = 0; k + 1 < cuts[i].size(); ++k) {
                for (std::size_t s = group[i]; s < group[i + 1]; ++s) {
                    pieces.push_back(Piece{cuts[i][k], cuts[i][k + 1], segments[s].origin,
                                           segments[s].label});
                }
            }
        }
        return pieces;
    }
 
    // Turns the split pieces into vertices and edges: equal pieces — the
    // overlapping and duplicated input of the same stretch — become one edge
    // remembering every input shape that produced it.
    void internVertices(std::vector<Piece>& pieces, const std::vector<PointType>& isolated) {
        // The piece endpoints were reduced by @ref split, before the cuts were
        // copied into pieces, so the sort and the hash lookups below already read
        // them on the normalized fast path.
        std::sort(pieces.begin(), pieces.end(), [](const Piece& left, const Piece& right) {
            if (!(left.a == right.a)) {
                return left.a < right.a;
            }
            if (!(left.b == right.b)) {
                return left.b < right.b;
            }
            return left.origin < right.origin;
        });
 
        std::unordered_map<PointType, std::uint32_t> vertexOf;
        const auto idOf = [&](const PointType& point) {
            const auto found = vertexOf.find(point);
            if (found != vertexOf.end()) {
                return found->second;
            }
            const auto id = static_cast<std::uint32_t>(points_.size());
            points_.push_back(point);
            vertexOf.emplace(point, id);
            return id;
        };
 
        originOffset_.push_back(0);
        for (std::size_t i = 0; i < pieces.size();) {
            origin_.push_back(idOf(pieces[i].a));
            origin_.push_back(idOf(pieces[i].b));
            edgeGeometry_.push_back({EdgeKind::segment, pieces[i].a, pieces[i].b});
            edgeLabel_.push_back(pieces[i].label);
            // The pieces of one stretch are adjacent and ordered by their input
            // position, so the same shape is caught by looking one back only.
            std::size_t j = i;
            while (j < pieces.size() && pieces[j].a == pieces[i].a && pieces[j].b == pieces[i].b) {
                if (j == i || pieces[j].origin != pieces[j - 1].origin) {
                    originIndex_.push_back(pieces[j].origin);
                }
                ++j;
            }
            originOffset_.push_back(static_cast<std::uint32_t>(originIndex_.size()));
            i = j;
        }
 
        for (const PointType& point : isolated) {
            idOf(point);
        }
        next_.assign(origin_.size(), 0);
        face_.assign(origin_.size(), 0);
        outgoing_.assign(points_.size(), HalfedgeId());
    }
 
    // Reduces every rational coordinate the arrangement keeps — the interned
    // vertex positions and the defining geometry of each edge — to lowest terms,
    // once, before anything reads them.
    //
    // Two kinds of coordinate arrive here unreduced. A constructed one (a
    // crossing, a point cut out of a carrier) is a fraction whose normalization
    // this type defers, and an input one may have been handed to us unreduced by
    // the caller. Either way a stored coordinate is read over and over — by the
    // wiring and face passes immediately below, by point location, and by every
    // query the caller later makes — and each of those reads recomputes the same
    // gcd and throws it away. This spends one gcd per coordinate and keeps it, so
    // all of that lands on the already-normalized fast path. On a value already
    // in lowest terms it is a flag test.
    //
    // This is the belt to @ref split's braces: it covers the coordinates that
    // never passed through a cut list — isolated input points, the source and
    // direction of a ray, the two points defining a line — and costs nothing on
    // the ones that did.
    void simplifyStoredCoordinates() {
        if constexpr (pgl::is_Rational_v<NumberType>) {
            const auto reduce = [](PointType& point) {
                point.x().simplify();
                point.y().simplify();
            };
            for (PointType& point : points_) {
                reduce(point);
            }
            for (EdgeGeometry& geometry : edgeGeometry_) {
                reduce(geometry.a);
                reduce(geometry.b);
            }
        }
    }
 
    // The coordinate type the vertex approximations are kept for: every sign
    // predicate that reads points_ back is an orientation over three of them.
    using VertexCoordinate =
        detail::orientation_coordinate_t<NumberType, NumberType, NumberType>;
 
    // Converts every interned vertex once, for those predicates to read back.
    //
    // They run inside comparators — the rotational sort around a vertex, the
    // status order of the sweep — so a vertex takes part in a logarithmic number
    // of them and converting per predicate re-derives the same few doubles over
    // and over. Must run after @ref simplifyStoredCoordinates: an unreduced
    // fraction can overflow double in both parts, leaving a quotient the filter
    // can only abstain on. Stays empty, and unindexed, where
    // @ref detail::filtersSign says exact arithmetic is cheap enough already.
    void syncVertexApproximations() {
        if constexpr (detail::filtersSign<VertexCoordinate>) {
            vertexApproximations_.clear();
            vertexApproximations_.reserve(points_.size());
            for (const PointType& point : points_) {
                vertexApproximations_.push_back(detail::approximatePoint(point));
            }
        }
    }
 
    // An interned vertex paired with the approximation kept for it. The two
    // arrays are parallel, and a vertex added without one would be read as some
    // other vertex's approximation — a sign the filter then proves wrong rather
    // than abstains on — so the pairing is asserted rather than trusted.
    [[nodiscard]] auto filteredVertex(std::uint32_t index) const {
        assert(!detail::filtersSign<VertexCoordinate> ||
               vertexApproximations_.size() == points_.size());
        return detail::filtered<VertexCoordinate>(points_[index], vertexApproximations_, index);
    }
 
    // Sorts the halfedges leaving each vertex counterclockwise and links them:
    // arriving at a vertex along one edge, the boundary of the face on the left
    // leaves along the next edge clockwise, which is the previous one in
    // counterclockwise order.
    void wireHalfedges() {
        std::vector<std::vector<std::uint32_t>> fan(points_.size());
        for (std::uint32_t h = 0; h < origin_.size(); ++h) {
            fan[origin_[h]].push_back(h);
        }
        for (std::uint32_t v = 0; v < fan.size(); ++v) {
            std::vector<std::uint32_t>& around = fan[v];
            if (around.empty()) {
                continue;
            }
            const PointType& center = points_[v];
            // Half 0 holds the directions with an angle in [0, pi), half 1 the
            // rest, so the comparison never needs an angle, only a sign.
            const auto half = [&](std::uint32_t h) {
                const PointType& to = points_[origin_[h ^ 1]];
                if (to.y() > center.y()) {
                    return 0;
                }
                if (to.y() < center.y()) {
                    return 1;
                }
                return to.x() > center.x() ? 0 : 1;
            };
            const auto filteredCenter = filteredVertex(v);
            std::sort(around.begin(), around.end(), [&](std::uint32_t left, std::uint32_t right) {
                const int leftHalf = half(left);
                const int rightHalf = half(right);
                if (leftHalf != rightHalf) {
                    return leftHalf < rightHalf;
                }
                return detail::orientationSignOf(filteredCenter,
                                                 filteredVertex(origin_[left ^ 1]),
                                                 filteredVertex(origin_[right ^ 1]))
                           .value() > 0;
            });
            const std::size_t degree = around.size();
            for (std::size_t i = 0; i < degree; ++i) {
                next_[around[i] ^ 1] = around[(i + degree - 1) % degree];
            }
            outgoing_[v] = HalfedgeId(around.front());
        }
    }
 
    // The direction a halfedge leaves its origin in, with the point that fixes
    // where its carrier sits. Ordering these is what puts the halfedges around
    // a vertex in rotational order.
    struct FanDirection {
        NumberType dx;
        NumberType dy;
        PointType anchor;
        std::uint32_t halfedge;
    };
 
    // The direction of a halfedge, taken from the interned endpoints when the
    // edge is a segment so a split edge leaves along its own piece rather than
    // along the whole original. A line's two halfedges escape opposite ways.
    [[nodiscard]] FanDirection fanDirection(std::uint32_t h) const {
        const EdgeGeometry& geometry = edgeGeometry_[h / 2];
        NumberType dx = geometry.b.x() - geometry.a.x();
        NumberType dy = geometry.b.y() - geometry.a.y();
        if (geometry.kind == EdgeKind::segment) {
            const PointType& from = points_[origin_[h]];
            const PointType& to = points_[origin_[h ^ 1]];
            dx = to.x() - from.x();
            dy = to.y() - from.y();
        } else if (geometry.kind == EdgeKind::line && h % 2 == 0) {
            dx = -dx;
            dy = -dy;
        }
        return {dx, dy, geometry.a, h};
    }
 
    // Which half-turn a direction falls in, so the angular comparison below can
    // use a cross product without ever spanning more than half a turn. Due east
    // counts as the upper half and due west as the lower, which puts the cut
    // just below the positive x axis.
    static int directionHalf(const FanDirection& direction) {
        if (direction.dy > NumberType(0)) {
            return 0;
        }
        if (direction.dy < NumberType(0)) {
            return 1;
        }
        return direction.dx > NumberType(0) ? 0 : 1;
    }
 
    // Counterclockwise order of two directions around a shared origin, starting
    // from the cut `directionHalf` sets. Exact throughout: the angle itself is
    // never formed, only compared.
    static bool fanLess(const FanDirection& left, const FanDirection& right) {
        const int leftHalf = directionHalf(left);
        const int rightHalf = directionHalf(right);
        if (leftHalf != rightHalf) {
            return leftHalf < rightHalf;
        }
        using Vector = Point<NumberType>;
        const Vector leftDirection(left.dx, left.dy);
        const auto cross = crossSign(leftDirection, Vector(right.dx, right.dy));
        if (cross != 0) {
            return cross > 0;
        }
        // Parallel ends with the same escape direction meet only at infinity.
        // Their transverse order is the order of their carriers there; it
        // reverses automatically at the opposite end of a line.
        const auto offset = crossSign(
            leftDirection, Vector(right.anchor.x() - left.anchor.x(),
                                  right.anchor.y() - left.anchor.y()));
        if (offset != 0) {
            return offset > 0;
        }
        return left.halfedge < right.halfedge;
    }
 
    // The same order seen from the point at infinity, where the plane is
    // traversed the other way round and every end is entered rather than left.
    static bool infinityFanLess(const FanDirection& left, const FanDirection& right) {
        return fanLess(right, left);
    }
 
    // Links the halfedges into face cycles for an arrangement carrying rays or
    // lines. Sorting each vertex's fan rotationally makes the next halfedge of
    // a cycle the one just clockwise of the twin, which is what walks a face
    // boundary; the vertex at infinity uses the reversed order.
    void wireHalfedgesUnbounded() {
        std::vector<std::vector<std::uint32_t>> fan(topologicalVertexCount());
        for (std::uint32_t h = 0; h < origin_.size(); ++h) {
            fan[origin_[h]].push_back(h);
        }
        for (std::uint32_t v = 0; v < fan.size(); ++v) {
            std::vector<std::uint32_t>& around = fan[v];
            if (around.empty()) {
                continue;
            }
            const bool atInfinity = infinity_.valid() && v == infinity_.index();
            std::sort(around.begin(), around.end(), [&](std::uint32_t left, std::uint32_t right) {
                return atInfinity
                    ? infinityFanLess(fanDirection(left), fanDirection(right))
                    : fanLess(fanDirection(left), fanDirection(right));
            });
            const std::size_t degree = around.size();
            for (std::size_t i = 0; i < degree; ++i) {
                next_[around[i] ^ 1] = around[(i + degree - 1) % degree];
            }
            outgoing_[v] = HalfedgeId(around.front());
            if (infinity_.valid() && v == infinity_.index()) {
                infinityFan_ = around;
            }
        }
    }
 
    // The unbounded face a point with no edge to its west lies in, found by
    // placing a due-west ray from the point into the fan at infinity.
    [[nodiscard]] HalfedgeId infinityBoundaryAtWest(const PointType& point) const {
        assert(!infinityFan_.empty());
        // At an end exactly collinear with the query ray, choose the sector on
        // its +y side. This is the same symbolic perturbation halfedgeLeftOf
        // uses at finite vertices.
        const FanDirection query{NumberType(-1), NumberType(0), point, 0};
        const auto after = std::upper_bound(
            infinityFan_.begin(), infinityFan_.end(), query,
            [&](const FanDirection& value, std::uint32_t halfedge) {
                return infinityFanLess(value, fanDirection(halfedge));
            });
        const std::uint32_t nextDirection =
            after == infinityFan_.end() ? infinityFan_.front() : *after;
        return HalfedgeId(nextDirection ^ 1);
    }
 
    // Turns the halfedge cycles into faces for an arrangement carrying rays or
    // lines. A cycle passing through infinity bounds an unbounded face, and the
    // rest nest as holes the way the bounded builder handles them.
    void buildFacesUnbounded() {
        constexpr std::uint32_t none = ~std::uint32_t{};
        const std::uint32_t halfedges = static_cast<std::uint32_t>(origin_.size());
        std::vector<std::uint32_t> cycleOf(halfedges, none);
        std::vector<std::uint32_t> representative;
        std::vector<bool> reachesInfinity;
        for (std::uint32_t h = 0; h < halfedges; ++h) {
            if (cycleOf[h] != none) {
                continue;
            }
            const auto id = static_cast<std::uint32_t>(representative.size());
            representative.push_back(h);
            bool unbounded = false;
            std::uint32_t walk = h;
            do {
                cycleOf[walk] = id;
                unbounded = unbounded || origin_[walk] == infinity_.index();
                walk = next_[walk];
            } while (walk != h);
            reachesInfinity.push_back(unbounded);
        }
 
        const auto cycles = static_cast<std::uint32_t>(representative.size());
        std::vector<std::uint32_t> parent(cycles);
        for (std::uint32_t i = 0; i < cycles; ++i) {
            parent[i] = i;
        }
        const auto root = [&parent](std::uint32_t x) {
            while (parent[x] != x) {
                parent[x] = parent[parent[x]];
                x = parent[x];
            }
            return x;
        };
 
        std::vector<bool> isOuter(cycles, false);
        for (std::uint32_t id = 0; id < cycles; ++id) {
            if (reachesInfinity[id]) {
                continue;
            }
            const std::uint32_t leftmost = leftmostVertexOf(representative[id]);
            if (turnsLeftEverywhereAt(representative[id], leftmost)) {
                isOuter[id] = true;
                continue;
            }
            const HalfedgeId left = halfedgeLeftOf(points_[leftmost]);
            const HalfedgeId other = left.valid() ? left : infinityBoundaryAtWest(points_[leftmost]);
            parent[root(id)] = root(cycleOf[other.index()]);
        }
 
        std::vector<std::uint32_t> infinityCycles;
        infinityCycles.reserve(cycles);
        for (std::uint32_t id = 0; id < cycles; ++id) {
            if (reachesInfinity[id]) {
                infinityCycles.push_back(id);
            }
        }
        const PointType zero(NumberType(0), NumberType(0));
        const std::uint32_t westCycle = cycleOf[infinityBoundaryAtWest(zero).index()];
        const auto west = std::find(infinityCycles.begin(), infinityCycles.end(), westCycle);
        if (west != infinityCycles.end()) {
            std::iter_swap(infinityCycles.begin(), west);
        }
 
        std::vector<std::uint32_t> faceOfComponent(cycles, none);
        outerCycle_.clear();
        unboundedFace_.clear();
        for (const std::uint32_t id : infinityCycles) {
            const std::uint32_t component = root(id);
            if (faceOfComponent[component] != none) {
                continue;
            }
            faceOfComponent[component] = static_cast<std::uint32_t>(outerCycle_.size());
            outerCycle_.push_back(HalfedgeId());
            unboundedFace_.push_back(true);
        }
        for (std::uint32_t id = 0; id < cycles; ++id) {
            if (!isOuter[id]) {
                continue;
            }
            const std::uint32_t component = root(id);
            assert(faceOfComponent[component] == none);
            faceOfComponent[component] = static_cast<std::uint32_t>(outerCycle_.size());
            outerCycle_.push_back(HalfedgeId(representative[id]));
            unboundedFace_.push_back(false);
        }
 
        std::vector<std::vector<HalfedgeId>> inner(outerCycle_.size());
        for (std::uint32_t id = 0; id < cycles; ++id) {
            const std::uint32_t f = faceOfComponent[root(id)];
            assert(f != none);
            if (!isOuter[id]) {
                inner[f].push_back(HalfedgeId(representative[id]));
            }
            std::uint32_t walk = representative[id];
            do {
                face_[walk] = f;
                walk = next_[walk];
            } while (walk != representative[id]);
        }
 
        innerCycle_.clear();
        innerOffset_.assign(1, 0);
        for (const std::vector<HalfedgeId>& cyclesOfFace : inner) {
            innerCycle_.insert(innerCycle_.end(), cyclesOfFace.begin(), cyclesOfFace.end());
            innerOffset_.push_back(static_cast<std::uint32_t>(innerCycle_.size()));
        }
        faceLabel_.assign(outerCycle_.size(), TLabel{});
    }
 
    // Traces the `next` cycles, tells the outer boundary cycles (counterclockwise)
    // from the inner ones (clockwise, or degenerate), and gathers the cycles that
    // bound the same face.
    void buildFaces() {
        constexpr std::uint32_t none = ~std::uint32_t{};
        const std::uint32_t halfedges = static_cast<std::uint32_t>(origin_.size());
 
        std::vector<std::uint32_t> cycleOf(halfedges, none);
        std::vector<std::uint32_t> representative;
        for (std::uint32_t h = 0; h < halfedges; ++h) {
            if (cycleOf[h] != none) {
                continue;
            }
            const auto id = static_cast<std::uint32_t>(representative.size());
            representative.push_back(h);
            std::uint32_t walk = h;
            do {
                cycleOf[walk] = id;
                walk = next_[walk];
            } while (walk != h);
        }
 
        // Cycle `n` stands for the outside of everything: an inner cycle with no
        // edge to its left bounds the unbounded face.
        const auto cycles = static_cast<std::uint32_t>(representative.size());
        std::vector<std::uint32_t> parent(cycles + 1);
        for (std::uint32_t i = 0; i < parent.size(); ++i) {
            parent[i] = i;
        }
        const auto root = [&parent](std::uint32_t x) {
            while (parent[x] != x) {
                parent[x] = parent[parent[x]];
                x = parent[x];
            }
            return x;
        };
 
        // Which cycles are outer boundaries, and, for the inner ones, the vertex
        // each of them has to ask what holds it from. The questions are collected
        // rather than asked one by one: one sweep answers the whole batch, where
        // a horizontal ray per cycle would scan every edge again for each.
        std::vector<bool> isOuter(cycles, false);
        std::vector<std::uint32_t> asking;
        std::vector<std::uint32_t> askedFrom;
        for (std::uint32_t id = 0; id < cycles; ++id) {
            const std::uint32_t leftmost = leftmostVertexOf(representative[id]);
            if (turnsLeftEverywhereAt(representative[id], leftmost)) {
                isOuter[id] = true;
                continue;
            }
            asking.push_back(id);
            askedFrom.push_back(leftmost);
        }
        // The leftmost vertex of an inner cycle has the face this cycle bounds
        // immediately to its left, so whatever edge is there bounds the same
        // face. Merging them all afterwards builds the same partition merging
        // them one at a time would.
        const std::vector<HalfedgeId> toTheLeft = halfedgesLeftOf(askedFrom);
        for (std::size_t i = 0; i < asking.size(); ++i) {
            const std::uint32_t other =
                toTheLeft[i].valid() ? cycleOf[toTheLeft[i].index()] : cycles;
            parent[root(asking[i])] = root(other);
        }
 
        // The unbounded face comes first, then one face per outer cycle.
        std::vector<std::uint32_t> faceOfComponent(cycles + 1, none);
        faceOfComponent[root(cycles)] = 0;
        outerCycle_.assign(1, HalfedgeId());
        for (std::uint32_t id = 0; id < cycles; ++id) {
            if (!isOuter[id]) {
                continue;
            }
            // A component holds at most one counterclockwise cycle: it is the
            // face's outer boundary, and the rest of the component is its holes.
            assert(faceOfComponent[root(id)] == none);
            faceOfComponent[root(id)] = static_cast<std::uint32_t>(outerCycle_.size());
            outerCycle_.push_back(HalfedgeId(representative[id]));
        }
 
        std::vector<std::vector<HalfedgeId>> inner(outerCycle_.size());
        for (std::uint32_t id = 0; id < cycles; ++id) {
            const std::uint32_t f = faceOfComponent[root(id)];
            assert(f != none);
            if (!isOuter[id]) {
                inner[f].push_back(HalfedgeId(representative[id]));
            }
            std::uint32_t walk = representative[id];
            do {
                face_[walk] = f;
                walk = next_[walk];
            } while (walk != representative[id]);
        }
 
        innerOffset_.assign(1, 0);
        for (const std::vector<HalfedgeId>& cyclesOfFace : inner) {
            innerCycle_.insert(innerCycle_.end(), cyclesOfFace.begin(), cyclesOfFace.end());
            innerOffset_.push_back(static_cast<std::uint32_t>(innerCycle_.size()));
        }
        faceLabel_.assign(outerCycle_.size(), TLabel{});
        unboundedFace_.assign(outerCycle_.size(), false);
        unboundedFace_.front() = true;
    }
 
    // -------------------------------------------------------------------------
    // Geometry helpers

    [[nodiscard]] HalfedgeId halfedgeLeftOf(const PointType& p) const {
        HalfedgeId best;
        WideNumber bestNumerator(0);
        WideNumber bestUpX(0);
        WideNumber bestUpY(1);
        const auto wide = [](const NumberType& value) {
            return static_cast<WideNumber>(value);
        };
        for (std::uint32_t h = 0; h < origin_.size(); h += 2) {
            const EdgeGeometry& geometry = edgeGeometry_[h / 2];
            const PointType& a = geometry.a;
            const PointType& b = geometry.b;
            if (a.y() == b.y()) {
                continue;  // a horizontal edge crosses no horizontal line
            }
            // Whether the edge runs upwards is a comparison, not a subtraction,
            // and it is all the height tests below need.
            const bool upwards = a.y() < b.y();
            if (geometry.kind == EdgeKind::segment) {
                const NumberType& lowY = upwards ? a.y() : b.y();
                const NumberType& highY = upwards ? b.y() : a.y();
                if (p.y() < lowY || !(p.y() < highY)) {
                    continue;
                }
            } else if (geometry.kind == EdgeKind::ray) {
                if (upwards ? p.y() < a.y() : !(p.y() < a.y())) {
                    continue;
                }
            }
            // The edge directed upwards, so that the query point lies strictly
            // right of it exactly when the crossing is strictly left of the
            // query point — which is the orientation predicate, and needs no
            // arithmetic here at all.
            const PointType& low = upwards ? a : b;
            const PointType& high = upwards ? b : a;
            if (!(orientationSign(low, high, p) < 0)) {
                continue;  // to the right of the query point, or through it
            }
            // How far right of the query point the edge crosses its horizontal
            // line, as the fraction numerator / upY — the upward direction's own
            // ordinate is the positive denominator. No division, so exact
            // coordinates stay exact.
            const WideNumber upX = wide(high.x()) - wide(low.x());
            const WideNumber upY = wide(high.y()) - wide(low.y());
            const WideNumber numerator =
                (wide(low.x()) - wide(p.x())) * upY - (wide(low.y()) - wide(p.y())) * upX;
            if (best.valid()) {
                const WideNumber here = numerator * bestUpY;
                const WideNumber there = bestNumerator * upY;
                if (here < there) {
                    continue;
                }
                if (here == there) {
                    // The same crossing point: the edge leaving it clockwise of
                    // the incumbent is the one whose left side holds the query.
                    if (!(upX * bestUpY - upY * bestUpX > WideNumber(0))) {
                        continue;
                    }
                }
            }
            // The halfedge running downwards has the crossing's right-hand side,
            // where the query point lies, on its left.
            best = HalfedgeId(upwards ? h + 1 : h);
            bestNumerator = numerator;
            bestUpX = upX;
            bestUpY = upY;
        }
        return best;
    }

    struct SweepOrder {
        using is_transparent = void;
        const Arrangement* arrangement;
        const std::vector<std::uint32_t>* position;
 
        bool operator()(std::uint32_t left, std::uint32_t right) const {
            const std::vector<std::uint32_t>& origin = arrangement->origin_;
            const std::uint32_t leftLow = origin[left ^ 1];
            const std::uint32_t rightLow = origin[right ^ 1];
            const auto vertex = [this](std::uint32_t index) {
                return arrangement->filteredVertex(index);
            };
            if (leftLow == rightLow) {
                // One vertex, both leaving it upwards: the one leaning further
                // right crosses the line further right.
                return detail::orientationSignOf(vertex(leftLow), vertex(origin[left]),
                                                 vertex(origin[right]))
                           .value() < 0;
            }
            if ((*position)[leftLow] > (*position)[rightLow]) {
                return detail::orientationSignOf(vertex(rightLow), vertex(origin[right]),
                                                 vertex(leftLow))
                           .value() > 0;
            }
            return detail::orientationSignOf(vertex(leftLow), vertex(origin[left]),
                                             vertex(rightLow))
                       .value() < 0;
        }
 
        // The edge is strictly left of the point when the point is strictly to
        // the right of it. A point *on* the edge is not, which is what keeps the
        // edges through a query vertex from answering it.
        bool operator()(std::uint32_t left, const PointType& p) const {
            return detail::orientationSignOf(
                       arrangement->filteredVertex(arrangement->origin_[left ^ 1]),
                       arrangement->filteredVertex(arrangement->origin_[left]),
                       detail::filtered<VertexCoordinate>(p))
                       .value() < 0;
        }
 
        bool operator()(const PointType& p, std::uint32_t right) const {
            return detail::orientationSignOf(
                       arrangement->filteredVertex(arrangement->origin_[right ^ 1]),
                       arrangement->filteredVertex(arrangement->origin_[right]),
                       detail::filtered<VertexCoordinate>(p))
                       .value() > 0;
        }
    };

    [[nodiscard]] std::vector<HalfedgeId> halfedgesLeftOf(
        const std::vector<std::uint32_t>& queries) const {
        const std::uint64_t edges = origin_.size() / 2;
        const std::uint64_t asked = queries.size();
        // log2 of the number of edges, near enough, and without a cast to double.
        const std::uint64_t depth = std::bit_width(edges);
        // What one sweep comparison costs in straddle tests, as a ratio so the
        // calibrated number stays an integer. Measured at about a half on exact
        // rational coordinates: the two predicates cost much the same, and the
        // sweep stays well inside the `(E + Q) log E` its bound allows, its
        // status holding only the edges the line currently crosses.
        constexpr std::uint64_t perComparisonNum = 1;
        constexpr std::uint64_t perComparisonDen = 2;
        if (!infinity_.valid() &&
            perComparisonDen * asked * edges >
                perComparisonNum * (2 * edges + asked) * depth) {
            return sweepHalfedgesLeftOf(queries);
        }
        std::vector<HalfedgeId> answer;
        answer.reserve(queries.size());
        for (const std::uint32_t query : queries) {
            answer.push_back(halfedgeLeftOf(points_[query]));
        }
        return answer;
    }

    [[nodiscard]] std::vector<HalfedgeId> sweepHalfedgesLeftOf(
        const std::vector<std::uint32_t>& queries) const {
        // Vertices by height, then abscissa. A vertex is its position here, so
        // the events sort on an integer and never on a coordinate.
        std::vector<std::uint32_t> byHeight(points_.size());
        for (std::uint32_t v = 0; v < byHeight.size(); ++v) {
            byHeight[v] = v;
        }
        std::sort(byHeight.begin(), byHeight.end(), [this](std::uint32_t a, std::uint32_t b) {
            if (!(points_[a].y() == points_[b].y())) {
                return points_[a].y() < points_[b].y();
            }
            return points_[a].x() < points_[b].x();
        });
        std::vector<std::uint32_t> position(points_.size());
        for (std::uint32_t i = 0; i < byHeight.size(); ++i) {
            position[byHeight[i]] = i;
        }
 
        // What happens at a vertex, in the order it has to happen.
        enum Phase : std::uint8_t { leaves = 0, joins = 1, asks = 2 };
        struct Event {
            std::uint32_t at;       // the vertex, as its position in `byHeight`
            std::uint32_t subject;  // a downward halfedge, or a query's position
            Phase phase;
        };
        std::vector<Event> events;
        events.reserve(origin_.size() + queries.size());
        for (std::uint32_t h = 0; h < origin_.size(); h += 2) {
            if (points_[origin_[h]].y() == points_[origin_[h + 1]].y()) {
                continue;  // horizontal: it crosses no horizontal line
            }
            const std::uint32_t downward =
                points_[origin_[h]].y() > points_[origin_[h + 1]].y() ? h : h + 1;
            events.push_back({position[origin_[downward]], downward, leaves});
            events.push_back({position[origin_[downward ^ 1]], downward, joins});
        }
        for (std::uint32_t q = 0; q < queries.size(); ++q) {
            events.push_back({position[queries[q]], q, asks});
        }
        std::sort(events.begin(), events.end(), [](const Event& left, const Event& right) {
            if (left.at != right.at) {
                return left.at < right.at;
            }
            return left.phase < right.phase;
        });
 
        std::vector<HalfedgeId> answer(queries.size());
        using Status = std::set<std::uint32_t, SweepOrder>;
        Status line(SweepOrder{this, &position});
        // Erasing by key would compare an edge that ends on the line against the
        // edges still crossing it, and an edge that has reached its top no longer
        // has a side of it to be on. Every edge remembers where it sits instead.
        std::vector<typename Status::iterator> seat(origin_.size() / 2);
        for (const Event& event : events) {
            if (event.phase == leaves) {
                line.erase(seat[event.subject / 2]);
            } else if (event.phase == joins) {
                const auto placed = line.insert(event.subject);
                assert(placed.second);
                seat[event.subject / 2] = placed.first;
            } else {
                const auto above = line.lower_bound(points_[queries[event.subject]]);
                if (above != line.begin()) {
                    answer[event.subject] = HalfedgeId(*std::prev(above));
                }
            }
        }
        return answer;
    }
 
    // Calls `fn` on every halfedge bounding the face, outer cycle first.
    template <class Function>
    void forEachBoundaryHalfedge(FaceId f, const Function& fn) const {
        const auto walkCycle = [&](HalfedgeId start) {
            std::uint32_t h = start.index();
            do {
                fn(HalfedgeId(h));
                h = next_[h];
            } while (h != start.index());
        };
        if (outerCycle_[f.index()].valid()) {
            walkCycle(outerCycle_[f.index()]);
        }
        for (HalfedgeId inner : innerCycles(f)) {
            walkCycle(inner);
        }
    }
 
    // The lexicographically smallest vertex a boundary cycle visits.
    [[nodiscard]] std::uint32_t leftmostVertexOf(std::uint32_t start) const {
        std::uint32_t leftmost = origin_[start];
        for (std::uint32_t h = next_[start]; h != start; h = next_[h]) {
            if (points_[origin_[h]] < points_[leftmost]) {
                leftmost = origin_[h];
            }
        }
        return leftmost;
    }

    [[nodiscard]] bool turnsLeftEverywhereAt(std::uint32_t start, std::uint32_t vertex) const {
        std::uint32_t h = start;
        do {
            // (h, ahead) is the pair of halfedges meeting at the vertex between
            // them, so their origins are the previous and the next vertex.
            const std::uint32_t ahead = next_[h];
            if (origin_[ahead] == vertex &&
                !(detail::orientationSignOf(filteredVertex(vertex),
                                            filteredVertex(origin_[ahead ^ 1]),
                                            filteredVertex(origin_[h]))
                      .value() > 0)) {
                return false;
            }
            h = ahead;
        } while (h != start);
        return true;
    }
 
    // The vertices of a boundary cycle, in order.
    [[nodiscard]] std::vector<PointType> cycleRing(HalfedgeId start) const {
        std::vector<PointType> ring;
        std::uint32_t h = start.index();
        do {
            ring.push_back(points_[origin_[h]]);
            h = next_[h];
        } while (h != start.index());
        return ring;
    }
 
    // Turns one boundary cycle into the polygons that describe the same area:
    // the spikes a dangling edge leaves behind are dropped, a cycle pinching
    // shut at a vertex is cut there, and every ring comes out counterclockwise
    // and in canonical form.
    static void collectRings(std::vector<PointType> walk, std::vector<Polygon<PointType>>& out) {
        pruneSpikes(walk);
        if (walk.size() < 3) {
            return;
        }
        std::vector<std::vector<PointType>> rings;
        detail::splitWalkIntoRings(walk, rings);
        for (std::vector<PointType>& ring : rings) {
            pruneSpikes(ring);
            if (ring.size() < 3) {
                continue;
            }
            const int orientation = detail::ringOrientation(ring);
            if (orientation == 0) {
                continue;
            }
            if (orientation < 0) {
                std::reverse(ring.begin(), ring.end());
            }
            // Counterclockwise and rotated onto its smallest vertex is exactly
            // the canonical form, so the polygon needs no normalization.
            std::rotate(ring.begin(), std::min_element(ring.begin(), ring.end()), ring.end());
            out.emplace_back(std::move(ring), /*trusted=*/true);
        }
    }
 
    // Drops the stretches a boundary walk covers twice, once each way: they
    // bound no area, and no ring may repeat a vertex around them.
    static void pruneSpikes(std::vector<PointType>& ring) {
        bool changed = true;
        while (changed && ring.size() >= 3) {
            changed = false;
            for (std::size_t i = 0; i < ring.size() && ring.size() >= 3; ++i) {
                const std::size_t size = ring.size();
                if (!(ring[(i + size - 1) % size] == ring[(i + 1) % size])) {
                    continue;
                }
                const auto at = static_cast<std::ptrdiff_t>(i);
                if (i + 1 < size) {
                    ring.erase(ring.begin() + at, ring.begin() + at + 2);
                } else {
                    ring.erase(ring.begin() + at);
                    ring.erase(ring.begin());
                }
                changed = true;
                break;
            }
        }
        if (ring.size() < 3) {
            ring.clear();
        }
    }

    class TrapezoidPointLocation {
        using WorkNumber = division_result_t<NumberType>;
        using WorkPoint = Point<WorkNumber>;
        using WorkLine = OrientedLine<WorkPoint>;
        using IntegralPoint = Point<std::int64_t>;
        using IntegralLine = OrientedLine<IntegralPoint>;
        using StoredLine = std::variant<WorkLine, IntegralLine>;
        using QueryInteger = rational_int_t<WorkNumber>;
        using WideIntegralPoint = Point<QueryInteger>;
        static constexpr std::uint32_t none = ~std::uint32_t{};
 
        // An abscissa of the sheared frame that a wall stands on. A vertex wall
        // carries the vertex, whose abscissa is `x + eps*y`; the two infinite
        // eps levels are the ends a vertical unbounded edge escapes to, past
        // every vertex on that abscissa yet within an infinitesimal of it.
        struct Abscissa {
            PointType point;
            std::int8_t epsInfinity = 0;
        };
 
        // Sheared order of two wall abscissae, formed without ever building a
        // sheared coordinate: the plane's own abscissa decides, and the
        // ordinate only separates what shares it.
        static bool abscissaLess(const Abscissa& left, const Abscissa& right) {
            if (left.point.x() != right.point.x()) {
                return left.point.x() < right.point.x();
            }
            if (left.epsInfinity != right.epsInfinity) {
                return left.epsInfinity < right.epsInfinity;
            }
            return left.epsInfinity == 0 && left.point.y() < right.point.y();
        }
 
        // An abscissa of the sheared frame that is not a wall: the pair stands
        // for `a + eps*b`. Sampling produces these and nothing stores them.
        struct SampleAbscissa {
            WorkNumber a{};
            WorkNumber b{};
        };
 
        // A point of the sheared plane, `value + eps*epsilon`. A curve meets an
        // abscissa strictly inside an infinitesimally thin slab at such a point
        // and at no ordinary one.
        struct DualPoint {
            WorkPoint value;
            WorkPoint epsilon;
        };
 
        // A vertical bound of a trapezoid. Walls are held in increasing
        // abscissa, so their index orders every finite bound outright.
        struct Bound {
            // -1 and +1 are the two infinities; zero carries a wall.
            std::int8_t infinity = 0;
            std::uint32_t wall = none;
 
            static Bound negativeInfinity() { return Bound{-1, none}; }
            static Bound positiveInfinity() { return Bound{1, none}; }
 
            bool operator==(const Bound& other) const {
                return infinity == other.infinity && wall == other.wall;
            }
        };
 
        // Live trapezoids immediately to the right of one transformed vertex.
        // Entries are append-only; an old entry is ignored once its trapezoid
        // becomes inactive. This is the local adjacency needed to walk a new
        // curve across the map without searching the history DAG again.
        struct Wall {
            std::vector<std::uint32_t> starts;
        };
 
        struct Curve {
            // Original coordinates, oriented left to right after the shear.
            StoredLine queryLine;
            Bound left;
            Bound right;
            HalfedgeId leftToRight;
            WorkNumber originalDx{};
            WorkNumber originalDy{};
        };
 
        struct Query {
            const PointType* point = nullptr;
            std::optional<IntegralPoint> integralPoint;
            std::optional<WideIntegralPoint> wideIntegralPoint;
        };
 
        struct Trapezoid {
            Bound left;
            Bound right;
            std::uint32_t bottom = none;
            std::uint32_t top = none;
            std::uint32_t leaf = none;
            FaceId face;
            bool active = true;
        };
 
        enum class NodeKind : std::uint8_t { leaf, x, curve };
 
        struct Node {
            NodeKind kind = NodeKind::leaf;
            std::uint32_t value = 0;  // trapezoid, curve, or wall
            std::uint32_t low = none;
            std::uint32_t high = none;
        };
 
    public:
        template <class UniformRandomBitGenerator>
        TrapezoidPointLocation(const Arrangement& arrangement,
                               UniformRandomBitGenerator&& generator) {
            makeWalls(arrangement);
            makeVertexIndex(arrangement);
            makeCurves(arrangement);
 
            const std::uint32_t initial = newTrapezoid(
                Bound::negativeInfinity(), Bound::positiveInfinity(), none, none);
            root_ = trapezoids_[initial].leaf;
 
            std::vector<std::uint32_t> order(curves_.size());
            for (std::uint32_t i = 0; i < order.size(); ++i) {
                order[i] = i;
            }
            std::shuffle(order.begin(), order.end(),
                         std::forward<UniformRandomBitGenerator>(generator));
            for (const std::uint32_t curve : order) {
                insertCurve(curve);
            }
            labelTrapezoids(arrangement);
            std::vector<WorkLine>().swap(constructionLines_);
        }
 
        [[nodiscard]] FaceId locateFace(const Arrangement& arrangement,
                                        const PointType& point) const {
            const Query query = makeQuery(point);
            std::uint32_t node = root_;
            while (nodes_[node].kind != NodeKind::leaf) {
                const Node& decision = nodes_[node];
                if (decision.kind == NodeKind::x) {
                    // A query standing on the wall goes left, which is the
                    // existing -x perturbation.
                    node = againstWall(point, decision.value) <= 0 ? decision.low
                                                                  : decision.high;
                    continue;
                }
                const Curve& curve = curves_[decision.value];
                const auto orientation = sideOf(curve, query);
                int side = orientation > 0 ? 1 : (orientation < 0 ? -1 : 0);
                if (side == 0) {
                    // Sign of cross(d, (-eps,+eps^2)): dy is primary and
                    // dx resolves the horizontal case.
                    side = curve.originalDy > WorkNumber(0)
                               ? 1
                               : (curve.originalDy < WorkNumber(0)
                                      ? -1
                                      : (curve.originalDx > WorkNumber(0) ? 1 : -1));
                }
                node = side < 0 ? decision.low : decision.high;
            }
            const FaceId result = trapezoids_[nodes_[node].value].face;
            assert(result.valid() && result.index() < arrangement.faceCount());
            return result;
        }
 
        [[nodiscard]] CellId locateCell(const Arrangement& arrangement,
                                        const PointType& point) const {
            const auto vertex = std::lower_bound(
                vertices_.begin(), vertices_.end(), point,
                [&](VertexId v, const PointType& p) { return arrangement.points_[v.index()] < p; });
            if (vertex != vertices_.end() && arrangement.points_[vertex->index()] == point) {
                return *vertex;
            }
 
            const Query query = makeQuery(point);
            std::uint32_t node = root_;
            while (nodes_[node].kind != NodeKind::leaf) {
                const Node& decision = nodes_[node];
                if (decision.kind == NodeKind::x) {
                    node = againstWall(point, decision.value) < 0 ? decision.low
                                                                 : decision.high;
                    continue;
                }
                const Curve& curve = curves_[decision.value];
                const auto side = sideOf(curve, query);
                if (side == 0) {
                    return HalfedgeId(2 * (curve.leftToRight.index() / 2));
                }
                node = side < 0 ? decision.low : decision.high;
            }
            return trapezoids_[nodes_[node].value].face;
        }
 
        [[nodiscard]] VertexId indexedVertex(std::size_t index) const {
            return vertices_[index];
        }
 
        // Vertices whose abscissa lies within the closed range, in increasing
        // order. An absent bound is unbounded on that side. The index is sorted
        // lexicographically on the original coordinates, so a query confined to
        // a narrow band of abscissae -- a vertical ray above all -- reaches only
        // the vertices that band can hold.
        template <class Number, class Fn>
        void visitCandidateVertices(const Arrangement& arrangement,
                                    const std::optional<Number>& low,
                                    const std::optional<Number>& high,
                                    Fn&& fn) const {
            auto first = vertices_.begin();
            auto last = vertices_.end();
            if (low) {
                first = std::lower_bound(
                    first, last, *low, [&](VertexId v, const Number& x) {
                        return arrangement.points_[v.index()].x() < x;
                    });
            }
            if (high) {
                last = std::upper_bound(
                    first, last, *high, [&](const Number& x, VertexId v) {
                        return x < arrangement.points_[v.index()].x();
                    });
            }
            for (auto it = first; it != last; ++it) {
                fn(*it);
            }
        }
 
        [[nodiscard]] HalfedgeId indexedHalfedge(std::size_t index) const {
            return HalfedgeId(2 * (curves_[index].leftToRight.index() / 2));
        }
 
        // Walks the search DAG carrying the query clipped to the region each
        // node stands for, held as an interval of the query's parameter. Both
        // kinds of decision cut that interval with a straight line -- an x node
        // with the wall it splits on, a curve node with the curve's supporting
        // line -- so a child is entered only where the query really reaches it,
        // and the region a node stands for stays exactly the intersection of
        // the cuts along the path. A curve is offered where its supporting line
        // still crosses the clipped query, which is a superset of the
        // intersecting arrangement edges: the caller deduplicates the offers
        // and applies the exact predicate.
        //
        // Nodes are taken earliest first, by the smallest parameter the node
        // can still hold, so once a node comes out no candidate found later can
        // precede it. That bound goes to `fn` as the frontier, which lets the
        // caller settle its answer for everything behind the frontier and stop
        // the walk there instead of running the query out to its far end. `fn`
        // is called with an invalid halfedge to carry the frontier alone, and
        // stops the walk by returning `true`. An absent frontier is the start
        // of an unbounded query, where nothing is settled yet.
        template <class Parameter, class Q, class Fn>
        bool visitCandidateHalfedges(const Q& query, Fn&& fn) const {
            const Parameter zero(0);
 
            const Parameter x0(query[0].x());
            const Parameter y0(query[0].y());
            const Parameter dx = Parameter(query[1].x()) - x0;
            const Parameter dy = Parameter(query[1].y()) - y0;
 
            // The parameters still in play, each end absent where the query
            // runs off to infinity.
            struct Span {
                std::uint32_t node;
                std::optional<Parameter> low;
                std::optional<Parameter> high;
            };
 
            // Keeps only the parameters where `value + t * rate` holds the
            // requested sign, and answers whether any parameter is left.
            const auto narrow = [&zero](Span& span, const Parameter& value,
                                        const Parameter& rate, bool nonNegative) {
                if (rate == zero) {
                    return nonNegative ? !(value < zero) : !(value > zero);
                }
                const Parameter root = -value / rate;
                if ((rate > zero) == nonNegative) {
                    if (!span.low || *span.low < root) {
                        span.low = root;
                    }
                } else {
                    if (!span.high || root < *span.high) {
                        span.high = root;
                    }
                }
                return !(span.low && span.high && *span.high < *span.low);
            };
 
            // The same, for a quantity of the sheared frame held as the pair
            // `(value, epsValue)` standing for `value + eps*epsValue` and
            // changing at `(rate, epsRate)`. The plane's own level decides
            // unless it vanishes identically, and only then does the eps level
            // speak. Where the plane's level merely reaches zero, keeping that
            // one parameter on both sides costs a candidate and spares the eps
            // level a comparison.
            const auto narrowSheared = [&](Span& span, const Parameter& value,
                                           const Parameter& rate,
                                           const Parameter& epsValue,
                                           const Parameter& epsRate, bool nonNegative) {
                if (rate != zero) {
                    return narrow(span, value, rate, nonNegative);
                }
                if (value != zero) {
                    return (value > zero) == nonNegative;
                }
                return narrow(span, epsValue, epsRate, nonNegative);
            };
 
            // The side of a curve's supporting line, as `value + t * rate`.
            const auto sideAlongQuery = [&](const Curve& curve) {
                return std::visit(
                    [&](const auto& line) {
                        const Parameter sourceX(line.source().x());
                        const Parameter sourceY(line.source().y());
                        const Parameter edgeX = Parameter(line.target().x()) - sourceX;
                        const Parameter edgeY = Parameter(line.target().y()) - sourceY;
                        return std::pair<Parameter, Parameter>(
                            edgeX * (y0 - sourceY) - edgeY * (x0 - sourceX),
                            edgeX * dy - edgeY * dx);
                    },
                    curve.queryLine);
            };
 
            // Orders the heap so the span holding the earliest parameter, an
            // absent bound being earliest of all, comes out first.
            const auto laterFirst = [](const Span& left, const Span& right) {
                if (!left.low) {
                    return false;
                }
                if (!right.low) {
                    return true;
                }
                return *right.low < *left.low;
            };
 
            Span start{root_, std::nullopt, std::nullopt};
            if constexpr (!OrientedLineConcept<Q>) {
                start.low = zero;
                if constexpr (!RayConcept<Q>) {
                    start.high = Parameter(1);
                }
            }
 
            std::vector<Span> pending;
            pending.push_back(std::move(start));
            const auto offer = [&](Span span) {
                pending.push_back(std::move(span));
                std::push_heap(pending.begin(), pending.end(), laterFirst);
            };
 
            while (!pending.empty()) {
                std::pop_heap(pending.begin(), pending.end(), laterFirst);
                const Span span = std::move(pending.back());
                pending.pop_back();
 
                const Parameter* frontier = span.low ? &*span.low : nullptr;
                if (fn(HalfedgeId(), frontier)) {
                    return true;
                }
 
                const Node& node = nodes_[span.node];
                if (node.kind == NodeKind::leaf) {
                    continue;
                }
 
                if (node.kind == NodeKind::x) {
                    // A query sitting exactly on the wall enters both sides.
                    const Abscissa& wall = wallXs_[node.value];
                    const Parameter splitX(wall.point.x());
                    Span low = span;
                    Span high = span;
                    bool below = false;
                    bool above = false;
                    if (wall.epsInfinity != 0) {
                        // The wall stands infinitesimally beside an abscissa of
                        // the plane, so the plane's own level settles the side
                        // wherever it settles anything at all.
                        below = narrow(low, splitX - x0, -dx, true);
                        above = narrow(high, x0 - splitX, dx, true);
                    } else {
                        const Parameter splitY(wall.point.y());
                        below = narrowSheared(low, splitX - x0, -dx, splitY - y0, -dy,
                                              true);
                        above = narrowSheared(high, x0 - splitX, dx, y0 - splitY, dy,
                                              true);
                    }
                    if (below) {
                        low.node = node.low;
                        offer(std::move(low));
                    }
                    if (above) {
                        high.node = node.high;
                        offer(std::move(high));
                    }
                    continue;
                }
 
                const Curve& curve = curves_[node.value];
                const auto [value, rate] = sideAlongQuery(curve);
                Span low = span;
                const bool below = narrow(low, value, rate, false);
                Span high = span;
                const bool above = narrow(high, value, rate, true);
                if (below && above) {
                    // The clipped query reaches both sides, so it meets the
                    // supporting line inside the region this node stands for.
                    const HalfedgeId h(2 * (curve.leftToRight.index() / 2));
                    if (fn(h, frontier)) {
                        return true;
                    }
                }
                if (below) {
                    low.node = node.low;
                    offer(std::move(low));
                }
                if (above) {
                    high.node = node.high;
                    offer(std::move(high));
                }
            }
            return false;
        }
 
    private:
        static bool less(const Bound& left, const Bound& right) {
            if (left.infinity != right.infinity) {
                return left.infinity < right.infinity;
            }
            return left.infinity == 0 && left.wall < right.wall;
        }
 
        static Bound maximum(const Bound& left, const Bound& right) {
            return less(left, right) ? right : left;
        }
 
        static Bound minimum(const Bound& left, const Bound& right) {
            return less(left, right) ? left : right;
        }
 
        [[nodiscard]] Bound wallBound(const PointType& point,
                                      std::int8_t epsInfinity) const {
            const Abscissa key{point, epsInfinity};
            const auto found =
                std::lower_bound(wallXs_.begin(), wallXs_.end(), key, abscissaLess);
            assert(found != wallXs_.end() && !abscissaLess(key, *found));
            return Bound{0, static_cast<std::uint32_t>(found - wallXs_.begin())};
        }
 
        // Where an ordinary query point falls against a wall. Its abscissa is
        // the point read lexicographically, so nothing is computed here.
        [[nodiscard]] int againstWall(const PointType& point, std::uint32_t wall) const {
            const Abscissa& abscissa = wallXs_[wall];
            if (point.x() != abscissa.point.x()) {
                return point.x() < abscissa.point.x() ? -1 : 1;
            }
            if (abscissa.epsInfinity != 0) {
                return abscissa.epsInfinity > 0 ? -1 : 1;
            }
            if (point.y() != abscissa.point.y()) {
                return point.y() < abscissa.point.y() ? -1 : 1;
            }
            return 0;
        }
 
        // The same for a point of the sheared plane: the plane's own abscissa
        // first, then the eps level the shear introduces, then the eps^2 level
        // a sampled point can reach and a wall never does.
        [[nodiscard]] int againstWall(const DualPoint& point, std::uint32_t wall) const {
            const Abscissa& abscissa = wallXs_[wall];
            const WorkNumber x(abscissa.point.x());
            if (point.value.x() != x) {
                return point.value.x() < x ? -1 : 1;
            }
            if (abscissa.epsInfinity != 0) {
                return abscissa.epsInfinity > 0 ? -1 : 1;
            }
            const WorkNumber y(abscissa.point.y());
            const WorkNumber first = point.value.y() + point.epsilon.x();
            if (first != y) {
                return first < y ? -1 : 1;
            }
            if (point.epsilon.y() != WorkNumber(0)) {
                return point.epsilon.y() < WorkNumber(0) ? -1 : 1;
            }
            return 0;
        }
 
        void makeVertexIndex(const Arrangement& arrangement) {
            vertices_.reserve(arrangement.points_.size());
            for (std::uint32_t v = 0; v < arrangement.points_.size(); ++v) {
                vertices_.push_back(VertexId(v));
            }
            std::sort(vertices_.begin(), vertices_.end(), [&](VertexId left, VertexId right) {
                return arrangement.points_[left.index()] < arrangement.points_[right.index()];
            });
        }
 
        void makeWalls(const Arrangement& arrangement) {
            wallXs_.reserve(arrangement.points_.size());
            for (const PointType& point : arrangement.points_) {
                wallXs_.push_back(Abscissa{point, 0});
            }
            // A vertical unbounded edge occupies a single abscissa of the plane
            // and escapes to an end infinitesimally beside it, past every
            // vertex there. Those ends bound trapezoids, so they are walls too.
            for (const EdgeGeometry& geometry : arrangement.edgeGeometry_) {
                if (geometry.kind == EdgeKind::segment ||
                    geometry.a.x() != geometry.b.x()) {
                    continue;
                }
                const bool upward = geometry.a.y() < geometry.b.y();
                if (geometry.kind == EdgeKind::line || upward) {
                    wallXs_.push_back(Abscissa{geometry.a, 1});
                }
                if (geometry.kind == EdgeKind::line || !upward) {
                    wallXs_.push_back(Abscissa{geometry.a, -1});
                }
            }
            std::sort(wallXs_.begin(), wallXs_.end(), abscissaLess);
            wallXs_.erase(std::unique(wallXs_.begin(), wallXs_.end(),
                                      [](const Abscissa& left, const Abscissa& right) {
                                          return !abscissaLess(left, right) &&
                                                 !abscissaLess(right, left);
                                      }),
                          wallXs_.end());
            walls_.resize(wallXs_.size());
        }
 
        void makeCurves(const Arrangement& arrangement) {
            curves_.reserve(arrangement.edgeGeometry_.size());
            constructionLines_.reserve(arrangement.edgeGeometry_.size());
            for (std::uint32_t edge = 0; edge < arrangement.edgeGeometry_.size(); ++edge) {
                const EdgeGeometry& geometry = arrangement.edgeGeometry_[edge];
                WorkNumber originalDx = WorkNumber(geometry.b.x()) - WorkNumber(geometry.a.x());
                WorkNumber originalDy = WorkNumber(geometry.b.y()) - WorkNumber(geometry.a.y());
                // Lexicographic order is the abscissa order of the sheared
                // frame, so no coordinate has to be formed to read it off.
                const bool forward = geometry.a < geometry.b;
                if (!forward) {
                    originalDx = -originalDx;
                    originalDy = -originalDy;
                }
 
                Bound left = Bound::negativeInfinity();
                Bound right = Bound::positiveInfinity();
                if (geometry.kind == EdgeKind::segment) {
                    left = wallBound(forward ? geometry.a : geometry.b, 0);
                    right = wallBound(forward ? geometry.b : geometry.a, 0);
                } else {
                    // A vertical edge spans no abscissa of the plane at all, so
                    // an unbounded one escapes to the infinitesimal side of the
                    // abscissa it stands on rather than to infinity.
                    const bool vertical = geometry.a.x() == geometry.b.x();
                    if (geometry.kind == EdgeKind::ray) {
                        (forward ? left : right) = wallBound(geometry.a, 0);
                        if (vertical) {
                            (forward ? right : left) =
                                wallBound(geometry.a, forward ? 1 : -1);
                        }
                    } else if (vertical) {
                        left = wallBound(geometry.a, -1);
                        right = wallBound(geometry.a, 1);
                    }
                }
 
                WorkPoint originalA(WorkNumber(geometry.a.x()), WorkNumber(geometry.a.y()));
                WorkPoint originalB(WorkNumber(geometry.b.x()), WorkNumber(geometry.b.y()));
                if (!forward) {
                    std::swap(originalA, originalB);
                }
                // The shear has determinant one, so orientation is identical in
                // both frames and the original coordinates serve throughout.
                StoredLine queryLine(std::in_place_type<WorkLine>, originalA, originalB);
                if constexpr (is_Rational_v<WorkNumber>) {
                    if (auto integral = WorkLine(originalA, originalB)
                                            .template integralLine<std::int64_t>()) {
                        queryLine = std::move(*integral);
                    }
                }
                curves_.push_back(Curve{std::move(queryLine), left, right,
                                        HalfedgeId(2 * edge + (forward ? 0 : 1)),
                                        originalDx, originalDy});
                constructionLines_.emplace_back(originalA, originalB);
            }
        }
 
        // Which side of a curve a point of the sheared plane falls on: the
        // ordinary orientation decides, and where that vanishes the point's own
        // infinitesimal offset does.
        [[nodiscard]] std::partial_ordering constructionSideOf(
            std::uint32_t curve, const DualPoint& point) const {
            const std::partial_ordering side = std::visit(
                [&](const auto& line) {
                    return orientationSign(line.source(), line.target(), point.value);
                },
                curves_[curve].queryLine);
            if (side != 0) {
                return side;
            }
            const WorkNumber cross = curves_[curve].originalDx * point.epsilon.y() -
                                     curves_[curve].originalDy * point.epsilon.x();
            return cross <=> WorkNumber(0);
        }
 
        [[nodiscard]] static std::partial_ordering sideOf(const Curve& curve,
                                                           const Query& query) {
            return std::visit(
                [&](const auto& line) {
                    using Line = std::remove_cvref_t<decltype(line)>;
                    if constexpr (std::same_as<Line, IntegralLine>) {
                        if (query.integralPoint) {
                            return orientationSign(line.source(), line.target(),
                                                   *query.integralPoint);
                        }
                        if (query.wideIntegralPoint) {
                            return orientationSign(line.source(), line.target(),
                                                   *query.wideIntegralPoint);
                        }
                    }
                    return orientationSign(line.source(), line.target(), *query.point);
                },
                curve.queryLine);
        }
 
        [[nodiscard]] Query makeQuery(const PointType& point) const {
            Query query{&point, std::nullopt, std::nullopt};
            classifyQuery(query, point);
            return query;
        }
 
        void classifyQuery(Query& query, const PointType& original) const {
            const auto storeIntegral = [&](QueryInteger x, QueryInteger y) {
                if (detail::representableAs<std::int64_t>(x) &&
                    detail::representableAs<std::int64_t>(y)) {
                    query.integralPoint.emplace(detail::narrowTo<std::int64_t>(x),
                                                detail::narrowTo<std::int64_t>(y));
                } else {
                    query.wideIntegralPoint.emplace(std::move(x), std::move(y));
                }
            };
 
            // Curve nodes operate in the original coordinates, and so do x
            // nodes now that the shear is symbolic; an integral query point
            // only ever needs its own coordinates.
            if constexpr (is_Rational_v<NumberType>) {
                if (original.x().isInteger() && original.y().isInteger()) {
                    storeIntegral(
                        QueryInteger(static_cast<rational_int_t<NumberType>>(
                            original.x())),
                        QueryInteger(static_cast<rational_int_t<NumberType>>(
                            original.y())));
                }
            } else if constexpr (detail::extended_integral<NumberType> ||
                                 std::same_as<NumberType, BigInt>) {
                storeIntegral(QueryInteger(original.x()), QueryInteger(original.y()));
            }
        }
 
        // An abscissa strictly between two bounds. Where the bounds share the
        // abscissa of the plane -- vertices stacked above one another -- what
        // separates them lives at the eps level, and the sample goes there.
        [[nodiscard]] SampleAbscissa sample(const Bound& left, const Bound& right) const {
            assert(less(left, right));
            const Abscissa* low = left.infinity == 0 ? &wallXs_[left.wall] : nullptr;
            const Abscissa* high = right.infinity == 0 ? &wallXs_[right.wall] : nullptr;
            if (low == nullptr && high == nullptr) {
                return SampleAbscissa{WorkNumber(0), WorkNumber(0)};
            }
            if (low == nullptr) {
                return SampleAbscissa{WorkNumber(high->point.x()) - WorkNumber(1),
                                      WorkNumber(0)};
            }
            if (high == nullptr) {
                return SampleAbscissa{WorkNumber(low->point.x()) + WorkNumber(1),
                                      WorkNumber(0)};
            }
            const WorkNumber lowX(low->point.x());
            const WorkNumber highX(high->point.x());
            if (lowX != highX) {
                return SampleAbscissa{(lowX + highX) / WorkNumber(2), WorkNumber(0)};
            }
            // Both bounds stand on one abscissa of the plane, so the eps level
            // is what is left to separate, its infinities included.
            if (low->epsInfinity < 0) {
                return SampleAbscissa{lowX, high->epsInfinity > 0
                                                ? WorkNumber(0)
                                                : WorkNumber(high->point.y()) -
                                                      WorkNumber(1)};
            }
            if (high->epsInfinity > 0) {
                return SampleAbscissa{lowX, WorkNumber(low->point.y()) + WorkNumber(1)};
            }
            return SampleAbscissa{lowX, (WorkNumber(low->point.y()) +
                                         WorkNumber(high->point.y())) /
                                            WorkNumber(2)};
        }
 
        // The point of a curve whose sheared abscissa is `x`. Inside a slab of
        // infinitesimal width the curve reaches the abscissa only infinitesimally
        // away from the plane's own points, which is what the eps part carries;
        // a vertical curve is the one that lives there outright.
        [[nodiscard]] DualPoint pointAt(std::uint32_t curve,
                                        const SampleAbscissa& x) const {
            const WorkPoint& source = constructionLines_[curve].source();
            const WorkNumber& dx = curves_[curve].originalDx;
            const WorkNumber& dy = curves_[curve].originalDy;
            const WorkPoint zero(WorkNumber(0), WorkNumber(0));
            if (dx == WorkNumber(0)) {
                // The eps level of the abscissa is the ordinate outright.
                assert(x.a == source.x());
                return DualPoint{WorkPoint(source.x(), x.b), zero};
            }
            const WorkNumber along = (x.a - source.x()) / dx;
            const WorkNumber y = source.y() + along * dy;
            // Sliding along the curve moves the abscissa by the slide itself,
            // so covering the eps level takes the ordinate still missing.
            const WorkNumber slide = x.b - y;
            if (slide == WorkNumber(0)) {
                return DualPoint{WorkPoint(x.a, y), zero};
            }
            return DualPoint{WorkPoint(x.a, y), WorkPoint(slide, slide * dy / dx)};
        }
 
        [[nodiscard]] bool crosses(std::uint32_t curveId,
                                   const Trapezoid& trapezoid) const {
            const Curve& curve = curves_[curveId];
            const Bound left = maximum(curve.left, trapezoid.left);
            const Bound right = minimum(curve.right, trapezoid.right);
            if (!less(left, right)) {
                return false;
            }
            const DualPoint point = pointAt(curveId, sample(left, right));
            if (trapezoid.bottom != none &&
                constructionSideOf(trapezoid.bottom, point) <= 0) {
                return false;
            }
            if (trapezoid.top != none &&
                constructionSideOf(trapezoid.top, point) >= 0) {
                return false;
            }
            return true;
        }
 
        [[nodiscard]] std::uint32_t newLeaf(std::uint32_t trapezoid) {
            const std::uint32_t id = static_cast<std::uint32_t>(nodes_.size());
            nodes_.push_back(Node{NodeKind::leaf, trapezoid, none, none});
            return id;
        }
 
        [[nodiscard]] std::uint32_t newTrapezoid(Bound left, Bound right,
                                                  std::uint32_t bottom,
                                                  std::uint32_t top) {
            if (trapezoids_.size() >= none || nodes_.size() >= none) {
                throw std::length_error("Arrangement point-location index exceeds 32-bit capacity");
            }
            const std::uint32_t id = static_cast<std::uint32_t>(trapezoids_.size());
            trapezoids_.push_back(Trapezoid{std::move(left), std::move(right), bottom, top,
                                            none, FaceId(), true});
            if (trapezoids_.back().left.infinity == 0) {
                walls_[trapezoids_.back().left.wall].starts.push_back(id);
            }
            trapezoids_.back().leaf = newLeaf(id);
            return id;
        }
 
        [[nodiscard]] std::uint32_t newNode(Node node) {
            if (nodes_.size() >= none) {
                throw std::length_error("Arrangement point-location index exceeds 32-bit capacity");
            }
            const std::uint32_t id = static_cast<std::uint32_t>(nodes_.size());
            nodes_.push_back(std::move(node));
            return id;
        }
 
        [[nodiscard]] Node curveNode(std::uint32_t curve, std::uint32_t below,
                                     std::uint32_t above) const {
            return Node{NodeKind::curve, curve, below, above};
        }
 
        [[nodiscard]] Node xNode(const Bound& bound, std::uint32_t left,
                                 std::uint32_t right) const {
            assert(bound.infinity == 0);
            return Node{NodeKind::x, bound.wall, left, right};
        }
 
        // Chooses an abscissa just to the right of a curve's left end, but
        // before the next arrangement wall. No vertical wall of the trapezoidal
        // map can occur inside that open slab.
        [[nodiscard]] SampleAbscissa probeRightOf(const Bound& left,
                                                  const Bound& right) const {
            const std::size_t following =
                left.infinity < 0 ? 0 : static_cast<std::size_t>(left.wall) + 1;
            Bound next = Bound::positiveInfinity();
            if (following < wallXs_.size()) {
                next = Bound{0, static_cast<std::uint32_t>(following)};
            }
            return sample(left, minimum(next, right));
        }
 
        [[nodiscard]] std::uint32_t locateTrapezoid(const DualPoint& point) const {
            std::uint32_t nodeId = root_;
            while (nodes_[nodeId].kind != NodeKind::leaf) {
                const Node& node = nodes_[nodeId];
                if (node.kind == NodeKind::x) {
                    const int side = againstWall(point, node.value);
                    assert(side != 0);
                    nodeId = side < 0 ? node.low : node.high;
                    continue;
                }
                const auto side = constructionSideOf(node.value, point);
                assert(side != 0);
                nodeId = side < 0 ? node.low : node.high;
            }
            const std::uint32_t trapezoid = nodes_[nodeId].value;
            assert(trapezoids_[trapezoid].active);
            return trapezoid;
        }
 
        // One DAG query finds the first trapezoid. Every subsequent one is
        // selected from the live trapezoids beginning at the current right
        // wall, so insertion work follows the curve instead of revisiting its
        // complete search history.
        [[nodiscard]] std::vector<std::uint32_t> crossedByWalk(
            std::uint32_t curveId) const {
            const Curve& curve = curves_[curveId];
            const SampleAbscissa probe = probeRightOf(curve.left, curve.right);
            std::uint32_t current = locateTrapezoid(pointAt(curveId, probe));
            std::vector<std::uint32_t> crossed;
 
            for (;;) {
                const Trapezoid& trapezoid = trapezoids_[current];
                assert(trapezoid.active && crosses(curveId, trapezoid));
                crossed.push_back(current);
 
                if (!less(trapezoid.right, curve.right)) {
                    break;
                }
                assert(trapezoid.right.infinity == 0);
 
                std::uint32_t next = none;
                for (const std::uint32_t candidate :
                     walls_[trapezoid.right.wall].starts) {
                    const Trapezoid& adjacent = trapezoids_[candidate];
                    if (adjacent.active && adjacent.left == trapezoid.right &&
                        crosses(curveId, adjacent)) {
                        assert(next == none);
                        next = candidate;
                    }
                }
                if (next == none) {
                    throw std::logic_error(
                        "arrangement edge lost its adjacent trapezoid");
                }
                current = next;
            }
            return crossed;
        }
 
        void insertCurve(std::uint32_t curveId) {
            const Curve& curve = curves_[curveId];
            const std::vector<std::uint32_t> crossed = crossedByWalk(curveId);
            if (crossed.empty()) {
                throw std::logic_error("arrangement edge did not cross its trapezoidal map");
            }
 
            std::uint32_t previousAbove = none;
            std::uint32_t previousBelow = none;
            for (const std::uint32_t oldId : crossed) {
                const Trapezoid old = trapezoids_[oldId];
                const Bound left = maximum(curve.left, old.left);
                const Bound right = minimum(curve.right, old.right);
 
                std::uint32_t below;
                if (previousBelow != none &&
                    trapezoids_[previousBelow].bottom == old.bottom &&
                    trapezoids_[previousBelow].right == left) {
                    below = previousBelow;
                    trapezoids_[below].right = right;
                } else {
                    below = newTrapezoid(left, right, old.bottom, curveId);
                }
 
                std::uint32_t above;
                if (previousAbove != none && trapezoids_[previousAbove].top == old.top &&
                    trapezoids_[previousAbove].right == left) {
                    above = previousAbove;
                    trapezoids_[above].right = right;
                } else {
                    above = newTrapezoid(left, right, curveId, old.top);
                }
 
                const bool hasLeftCap = less(old.left, left);
                const bool hasRightCap = less(right, old.right);
                std::uint32_t leftCap = none;
                std::uint32_t rightCap = none;
                if (hasLeftCap) {
                    leftCap = newTrapezoid(old.left, left, old.bottom, old.top);
                }
                if (hasRightCap) {
                    rightCap = newTrapezoid(right, old.right, old.bottom, old.top);
                }
 
                Node replacement = curveNode(curveId, trapezoids_[below].leaf,
                                             trapezoids_[above].leaf);
                if (hasRightCap) {
                    const std::uint32_t split = newNode(std::move(replacement));
                    replacement = xNode(right, split, trapezoids_[rightCap].leaf);
                }
                if (hasLeftCap) {
                    const std::uint32_t split = newNode(std::move(replacement));
                    replacement = xNode(left, trapezoids_[leftCap].leaf, split);
                }
                nodes_[old.leaf] = std::move(replacement);
                trapezoids_[oldId].active = false;
                previousBelow = below;
                previousAbove = above;
            }
        }
 
        // The face above a vertical strip that no edge crosses, read off the
        // fan at infinity, where the strip escapes. The strip's own abscissa is
        // part of the query: an end rising straight up is parallel to it, and
        // the fan separates parallel ends by where they stand.
        //
        // A strip standing on a single abscissa can tie with such an end, and
        // then either neighbouring sector answers: an end escaping straight up
        // from an abscissa the map leaves free is a ray or a line that the free
        // stretch itself passes through, so the two sides of it are one face.
        [[nodiscard]] FaceId faceAboveStrip(const Arrangement& arrangement,
                                            const WorkNumber& abscissa) const {
            const auto beforeStrip = [&](std::uint32_t halfedge) {
                const FanDirection end = arrangement.fanDirection(halfedge);
                if (directionHalf(end) != 0) {
                    return false;
                }
                if (end.dx != NumberType(0)) {
                    return end.dx > NumberType(0);
                }
                return abscissa < WorkNumber(end.anchor.x());
            };
            // The fan at infinity runs against the plain fan, so the ends
            // before the strip are exactly its prefix.
            const auto after = std::partition_point(
                arrangement.infinityFan_.begin(), arrangement.infinityFan_.end(),
                [&](std::uint32_t halfedge) { return !beforeStrip(halfedge); });
            const std::uint32_t nextDirection =
                after == arrangement.infinityFan_.end() ? arrangement.infinityFan_.front() : *after;
            return arrangement.face(HalfedgeId(nextDirection ^ 1));
        }
 
        // The face of a trapezoid with no curve above and none below. Nothing
        // meets the vertical line through it, so the whole line above it is
        // free of edges but for the vertices standing on it, and the face is
        // the one that line runs into at the top.
        //
        // Where the shear was a number that face was always the outside of the
        // map. It is still an unbounded face here, but no longer always the
        // same one: an infinitesimal shear leaves a vertical line spanning one
        // abscissa, so the trapezoids on either side of one are free of curves
        // and yet lie in the two faces it separates.
        [[nodiscard]] FaceId emptyFace(const Arrangement& arrangement,
                                       const Trapezoid& trapezoid) const {
            if (!arrangement.infinity_.valid()) {
                return FaceId(0);
            }
            return faceAboveStrip(arrangement, sample(trapezoid.left, trapezoid.right).a);
        }
 
        void labelTrapezoids(const Arrangement& arrangement) {
            for (Trapezoid& trapezoid : trapezoids_) {
                if (!trapezoid.active) {
                    continue;
                }
                if (trapezoid.bottom != none) {
                    trapezoid.face = arrangement.face(curves_[trapezoid.bottom].leftToRight);
                    if (trapezoid.top != none) {
                        assert(trapezoid.face == arrangement.face(
                            arrangement.twin(curves_[trapezoid.top].leftToRight)));
                    }
                } else if (trapezoid.top != none) {
                    trapezoid.face = arrangement.face(
                        arrangement.twin(curves_[trapezoid.top].leftToRight));
                } else {
                    trapezoid.face = emptyFace(arrangement, trapezoid);
                }
            }
        }
 
        std::uint32_t root_ = none;
        std::vector<VertexId> vertices_;
        std::vector<Abscissa> wallXs_;  // increasing
        std::vector<Wall> walls_;       // parallel to wallXs_
        std::vector<Curve> curves_;
        // Used only while constructing the RIC and released by the constructor:
        // the curves again, in the one representation sampling can divide in.
        std::vector<WorkLine> constructionLines_;
        std::vector<Trapezoid> trapezoids_;
        std::vector<Node> nodes_;
    };
 
    std::vector<PointType> points_;
    // Approximations of points_, parallel to it, or empty where the filter would
    // not pay for itself. See @ref syncVertexApproximations.
    std::vector<detail::ApproximatePoint> vertexApproximations_;
    VertexId infinity_;                            // symbolic vertex, absent for bounded input
    std::vector<HalfedgeId> outgoing_;         // one per vertex; invalid when isolated
    std::vector<std::uint32_t> origin_;        // one per halfedge
    std::vector<std::uint32_t> next_;          // one per halfedge
    std::vector<std::uint32_t> face_;          // one per halfedge
    std::vector<TLabel> edgeLabel_;            // one per edge
    std::vector<EdgeGeometry> edgeGeometry_;   // finite defining geometry, one per edge
    std::vector<std::uint32_t> originOffset_;  // one per edge, plus a closing entry
    std::vector<std::uint32_t> originIndex_;   // input positions, grouped by edge
    std::vector<HalfedgeId> outerCycle_;       // one per face; invalid for unbounded faces
    std::vector<std::uint32_t> innerOffset_;   // one per face, plus a closing entry
    std::vector<HalfedgeId> innerCycle_;       // inner cycles, grouped by face
    std::vector<TLabel> faceLabel_;            // one per face
    std::vector<bool> unboundedFace_;           // one per face
    std::vector<std::uint32_t> infinityFan_;    // outgoing ends in counterclockwise order
    std::shared_ptr<const TrapezoidPointLocation> pointLocation_;
};
 
template <class PointType, class TLabel>
template <class Q, class Fn>
bool Arrangement<PointType, TLabel>::visitStraightIntersecting(
    const Q& piece, Fn& fn, std::vector<bool>* seenVertices,
    std::vector<bool>* seenEdges) const {
    using CommonNumber = std::common_type_t<NumberType, typename Q::NumberType>;
    using Parameter = division_result_t<CommonNumber>;
 
    // A cell the piece meets, at the parameter where it first meets it. An edge
    // the piece runs along from infinity has no such parameter, and sorts ahead
    // of everything instead.
    struct Event {
        bool negativeInfinity = false;
        Parameter parameter{};
        IntersectionId id;
    };
 
    const Parameter ax(piece[0].x());
    const Parameter ay(piece[0].y());
    const Parameter dx = Parameter(piece[1].x()) - ax;
    const Parameter dy = Parameter(piece[1].y()) - ay;
    const Parameter squaredLength = dx * dx + dy * dy;
 
    // Where a point falls along the piece, as the parameter that is 0 at its
    // first endpoint and 1 at its second. A degenerate piece puts everything
    // at 0, which is the only parameter it has.
    const auto parameterOf = [&](const auto& point) {
        if (squaredLength == Parameter(0)) {
            return Parameter(0);
        }
        return ((Parameter(point.x()) - ax) * dx +
                (Parameter(point.y()) - ay) * dy) /
               squaredLength;
    };
 
    // The parameter at which the piece first meets an edge. Crossing edges meet
    // at one parameter; a collinear edge is met along a stretch, and the earliest
    // of it is what counts. A query bounded at its start never reports earlier
    // than that start.
    const auto firstOnEdge = [&](const EdgeGeometry& geometry) {
        Event event;
        if (squaredLength == Parameter(0)) {
            return event;
        }
 
        const Parameter ex = Parameter(geometry.b.x()) - Parameter(geometry.a.x());
        const Parameter ey = Parameter(geometry.b.y()) - Parameter(geometry.a.y());
        Parameter denominator = dx * ey - dy * ex;
        if (denominator != Parameter(0)) {
            Parameter numerator =
                (Parameter(geometry.a.x()) - ax) * ey -
                (Parameter(geometry.a.y()) - ay) * ex;
            event.parameter = numerator / denominator;
            return event;
        }
 
        if (geometry.kind == EdgeKind::line) {
            event.negativeInfinity = true;
        } else if (geometry.kind == EdgeKind::ray && ex * dx + ey * dy < Parameter(0)) {
            event.negativeInfinity = true;
        } else {
            event.parameter = parameterOf(geometry.a);
            if (geometry.kind == EdgeKind::segment) {
                event.parameter = std::min(event.parameter, parameterOf(geometry.b));
            }
        }
 
        if constexpr (!OrientedLineConcept<Q>) {
            if (event.negativeInfinity || event.parameter < Parameter(0)) {
                event.negativeInfinity = false;
                event.parameter = Parameter(0);
            }
        }
        return event;
    };
 
    // Order along the piece, breaking a tie towards the vertex, since a cell the
    // piece meets at the same parameter as an edge is the edge's endpoint and
    // stands for the contact. Handle order settles the rest, so the traversal is
    // reproducible.
    const auto eventLess = [](const Event& left, const Event& right) {
        if (left.negativeInfinity != right.negativeInfinity) {
            return left.negativeInfinity;
        }
        if (!left.negativeInfinity && left.parameter != right.parameter) {
            return left.parameter < right.parameter;
        }
        const bool leftVertex = std::holds_alternative<VertexId>(left.id);
        const bool rightVertex = std::holds_alternative<VertexId>(right.id);
        if (leftVertex != rightVertex) {
            return leftVertex;
        }
        return std::visit([](const auto& id) { return id.index(); }, left.id) <
               std::visit([](const auto& id) { return id.index(); }, right.id);
    };
 
    // The abscissae the piece can reach, open on a side the piece runs off to.
    std::optional<Parameter> lowX;
    std::optional<Parameter> highX;
    if (Parameter(piece[0].x()) == Parameter(piece[1].x())) {
        lowX = Parameter(piece[0].x());
        highX = lowX;
    } else if constexpr (RayConcept<Q>) {
        (Parameter(piece[1].x()) > Parameter(piece[0].x()) ? lowX : highX) =
            Parameter(piece[0].x());
    } else if constexpr (!OrientedLineConcept<Q>) {
        lowX = std::min(Parameter(piece[0].x()), Parameter(piece[1].x()));
        highX = std::max(Parameter(piece[0].x()), Parameter(piece[1].x()));
    }
 
    // Events wait here until the search has passed them, earliest on top.
    const auto laterFirst = [&eventLess](const Event& left, const Event& right) {
        return eventLess(right, left);
    };
    std::vector<Event> queued;
    const auto queue = [&](Event event) {
        queued.push_back(std::move(event));
        std::push_heap(queued.begin(), queued.end(), laterFirst);
    };
 
    // Queues a vertex the piece passes through, unless an earlier piece of the
    // same chain already reported it.
    const auto offerVertex = [&](VertexId v) {
        if (seenVertices != nullptr && (*seenVertices)[v.index()]) {
            return;
        }
        const PointType& point = points_[v.index()];
        if (piece.intersects(point)) {
            queue(Event{false, parameterOf(point), v});
        }
    };
 
    // The walk can reach one edge from several nodes. The offers it has already
    // made live in a table sized to the offers themselves rather than to the
    // arrangement, so a query never pays for the edges it does not look at.
    static constexpr std::uint32_t noEdge = ~std::uint32_t{};
    std::vector<std::uint32_t> offered(16, noEdge);
    std::size_t offeredCount = 0;
    // Open addressing with linear probing; the table always keeps a free slot,
    // so the scan terminates.
    const auto slotOf = [](const std::vector<std::uint32_t>& table, std::uint32_t edge) {
        std::size_t slot = (edge * 2654435761U) & (table.size() - 1);
        while (table[slot] != noEdge && table[slot] != edge) {
            slot = (slot + 1) & (table.size() - 1);
        }
        return slot;
    };
    const auto offeredBefore = [&](std::uint32_t edge) {
        if (2 * (offeredCount + 1) > offered.size()) {
            std::vector<std::uint32_t> grown(2 * offered.size(), noEdge);
            for (const std::uint32_t held : offered) {
                if (held != noEdge) {
                    grown[slotOf(grown, held)] = held;
                }
            }
            offered.swap(grown);
        }
        const std::size_t slot = slotOf(offered, edge);
        if (offered[slot] == edge) {
            return true;
        }
        offered[slot] = edge;
        ++offeredCount;
        return false;
    };
 
    // Queues an edge the piece crosses, given a candidate the search turned up.
    // Candidates are a superset, so the exact predicate decides here.
    const auto offerEdge = [&](HalfedgeId h) {
        const std::size_t edge = h.index() / 2;
        if (offeredBefore(static_cast<std::uint32_t>(edge))) {
            return;
        }
        if (seenEdges != nullptr && (*seenEdges)[edge]) {
            return;
        }
        const EdgeGeometry& geometry = edgeGeometry_[edge];
 
        // Where this piece meets a finite endpoint, the vertex there is the
        // whole of the contact and represents it. Another piece meeting the
        // same edge away from its endpoints still reports the edge.
        if (piece.intersects(geometry.a) ||
            (geometry.kind == EdgeKind::segment && piece.intersects(geometry.b))) {
            return;
        }
        const bool intersects = std::visit(
            [&](const auto& value) { return piece.intersects(value); }, (*this)[h]);
        if (!intersects) {
            return;
        }
        Event event = firstOnEdge(geometry);
        event.id = h;
        queue(std::move(event));
    };
 
    // Records a cell as reported and answers whether it still needs reporting,
    // which is how a chain avoids naming the same cell from two of its pieces.
    const auto markSeen = [&](const IntersectionId& id) {
        if (const auto* v = std::get_if<VertexId>(&id)) {
            if (seenVertices != nullptr) {
                if ((*seenVertices)[v->index()]) {
                    return false;
                }
                (*seenVertices)[v->index()] = true;
            }
        } else if (seenEdges != nullptr) {
            const std::size_t edge = std::get<HalfedgeId>(id).index() / 2;
            if ((*seenEdges)[edge]) {
                return false;
            }
            (*seenEdges)[edge] = true;
        }
        return true;
    };
 
    // Hands over the events the search can no longer precede: those strictly
    // ahead of the frontier it has reached, or all of them once it is over.
    // Events tying with the frontier wait, since a candidate still to come may
    // share their parameter and order before them.
    const auto release = [&](const Parameter* frontier, bool exhausted) {
        while (!queued.empty()) {
            if (!exhausted) {
                const Event& next = queued.front();
                if (frontier == nullptr ||
                    (!next.negativeInfinity && !(next.parameter < *frontier))) {
                    break;
                }
            }
            std::pop_heap(queued.begin(), queued.end(), laterFirst);
            Event event = std::move(queued.back());
            queued.pop_back();
            if (markSeen(event.id) && detail::invokeVisitor(fn, event.id)) {
                return true;
            }
        }
        return false;
    };
 
    if (pointLocation_) {
        pointLocation_->visitCandidateVertices(*this, lowX, highX, offerVertex);
        // The walk reports earliest first, so a visitor that stops on its first
        // answer stops the search with it instead of running the query out.
        if (pointLocation_->template visitCandidateHalfedges<Parameter>(
                piece, [&](HalfedgeId h, const Parameter* frontier) {
                    if (h.valid()) {
                        offerEdge(h);
                    }
                    return release(frontier, /*exhausted=*/false);
                })) {
            return true;
        }
    } else {
        for (std::size_t i = 0; i < points_.size(); ++i) {
            offerVertex(VertexId(static_cast<std::uint32_t>(i)));
        }
        for (std::size_t edge = 0; edge < edgeGeometry_.size(); ++edge) {
            offerEdge(HalfedgeId(static_cast<std::uint32_t>(2 * edge)));
        }
    }
    return release(nullptr, /*exhausted=*/true);
}
 
template <class PointType, class TLabel>
template <PointConcept Q, class Fn>
bool Arrangement<PointType, TLabel>::visitPointIntersecting(
    const Q& point, Fn& fn, std::vector<bool>* seenVertices,
    std::vector<bool>* seenEdges) const {
    for (std::size_t i = 0; i < points_.size(); ++i) {
        const VertexId v = pointLocation_ ? pointLocation_->indexedVertex(i)
                                          : VertexId(static_cast<std::uint32_t>(i));
        if ((seenVertices == nullptr || !(*seenVertices)[v.index()]) &&
            points_[v.index()].intersects(point)) {
            if (seenVertices != nullptr) {
                (*seenVertices)[v.index()] = true;
            }
            return detail::invokeVisitor(fn, IntersectionId(v));
        }
    }
 
    for (std::size_t i = 0; i < edgeGeometry_.size(); ++i) {
        const HalfedgeId h = pointLocation_
                                 ? pointLocation_->indexedHalfedge(i)
                                 : HalfedgeId(static_cast<std::uint32_t>(2 * i));
        const std::size_t edge = h.index() / 2;
        if (seenEdges != nullptr && (*seenEdges)[edge]) {
            continue;
        }
        const EdgeGeometry& geometry = edgeGeometry_[edge];
        if (point.intersects(geometry.a) ||
            (geometry.kind == EdgeKind::segment && point.intersects(geometry.b))) {
            continue;
        }
        if (!std::visit([&](const auto& value) { return point.intersects(value); },
                        (*this)[h])) {
            continue;
        }
        if (seenEdges != nullptr) {
            (*seenEdges)[edge] = true;
        }
        if (detail::invokeVisitor(fn, IntersectionId(h))) {
            return true;
        }
    }
    return false;
}
 
template <class PointType, class TLabel>
template <class Q, class Fn>
    requires (OrientedSegmentConcept<Q> || OrientedLineConcept<Q> ||
              RayConcept<Q> || MonotoneChainConcept<Q> || PolylineConcept<Q>)
bool Arrangement<PointType, TLabel>::visitIntersecting(const Q& r, Fn fn) const {
    if constexpr (MonotoneChainConcept<Q> || PolylineConcept<Q>) {
        if (r.empty()) {
            return false;
        }
        std::vector<bool> seenVertices(points_.size(), false);
        std::vector<bool> seenEdges(edgeGeometry_.size(), false);
        if (r.size() == 1) {
            return visitPointIntersecting(r[0], fn, &seenVertices, &seenEdges);
        }
        for (const auto& edge : r.orientedEdgesView()) {
            if (visitStraightIntersecting(edge, fn, &seenVertices, &seenEdges)) {
                return true;
            }
        }
        return false;
    } else {
        return visitStraightIntersecting(r, fn, nullptr, nullptr);
    }
}
 
template <class PointType, class TLabel>
void Arrangement<PointType, TLabel>::buildPointLocation() {
    // Stable across runs: callers wanting another insertion order can use the
    // generator-taking overload.
    std::mt19937 generator(0x50474cU);  // "PGL"
    buildPointLocation(generator);
}
 
template <class PointType, class TLabel>
template <class UniformRandomBitGenerator>
void Arrangement<PointType, TLabel>::buildPointLocation(
    UniformRandomBitGenerator&& generator) {
    if (!pointLocation_) {
        pointLocation_ = std::make_shared<const TrapezoidPointLocation>(
            *this, std::forward<UniformRandomBitGenerator>(generator));
    }
}
 
template <class PointType, class TLabel>
typename Arrangement<PointType, TLabel>::FaceId
Arrangement<PointType, TLabel>::locateFace(const PointType& point) const {
    return pointLocation_ ? pointLocation_->locateFace(*this, point) : locateFaceLinear(point);
}
 
template <class PointType, class TLabel>
typename Arrangement<PointType, TLabel>::CellId
Arrangement<PointType, TLabel>::locateCell(const PointType& point) const {
    return pointLocation_ ? pointLocation_->locateCell(*this, point) : locateCellLinear(point);
}
 
template <TriangleConcept TriangleType, SegmentConcept SegmentType>
Arrangement<typename Triangulation<TriangleType, SegmentType>::PointType,
            typename Triangulation<TriangleType, SegmentType>::TriId>
Triangulation<TriangleType, SegmentType>::asArrangement() const {
    using Result = Arrangement<PointType, TriId>;
 
    // Do not pass SegmentType directly: Result uses TriId for both edge and
    // face labels, whereas SegmentType may carry an unrelated edge label.
    std::vector<Segment<PointType>> meshEdges;
    meshEdges.reserve(numEdges());
    visitEdges([&](const SegmentType& edge) { meshEdges.emplace_back(edge[0], edge[1]); });
 
    Result result(meshEdges);
    for (std::size_t i = 0; i < result.faceCount(); ++i) {
        const typename Result::FaceId face(static_cast<std::uint32_t>(i));
        if (result.isUnbounded(face)) {
            continue;
        }
        const TriId triangle = locateId(result.witness(face));
        if (triangle.valid()) {
            result.label(face) = triangle;
        }
    }
    return result;
}
 
template <TriangleConcept TriangleType, SegmentConcept SegmentType>
template <class ResultNumber>
Arrangement<Point<ResultNumber>, typename Triangulation<TriangleType, SegmentType>::PointType>
Triangulation<TriangleType, SegmentType>::voronoiDiagram() const {
    assert(!empty() && "Triangulation::voronoiDiagram requires a nonempty triangulation");
 
    using ResultPoint = Point<ResultNumber>;
    using Diagram = Arrangement<ResultPoint, PointType>;
 
    // One exact circumcenter per current real triangle. This deliberately uses
    // the internal convex-hull triangulation, including any triangles hidden by
    // a carved domain: the method's precondition says that this whole current
    // connectivity is the Delaunay triangulation of the stored sites.
    std::vector<ResultPoint> centers(static_cast<std::size_t>(firstGhost_));
    for (TriIndex t = 0; t < firstGhost_; ++t) {
        centers[static_cast<std::size_t>(t)] = ResultPoint(
            triangleValue(t).circumcircle().template center<ResultNumber>());
    }
 
    std::vector<Shape<ResultPoint>> dualEdges;
    dualEdges.reserve(segToEdge_.size());
    for (TriIndex t = 0; t < firstGhost_; ++t) {
        const Tri& triangle = triangles_[static_cast<std::size_t>(t)];
        for (int side = 0; side < 3; ++side) {
            const TriIndex neighbor = triangle.nbr[static_cast<std::size_t>(side)];
            if (!isGhost(neighbor)) {
                // Each interior primal edge is visited from both incident
                // triangles. Emit its dual only from the lower triangle id.
                if (neighbor < t) {
                    continue;
                }
                const ResultPoint& a = centers[static_cast<std::size_t>(t)];
                const ResultPoint& b = centers[static_cast<std::size_t>(neighbor)];
                if (a != b) {
                    dualEdges.emplace_back(Segment<ResultPoint>(a, b));
                }
                continue;
            }
 
            // A real triangle is counterclockwise, so its directed side
            // a -> b has the hull interior on its left. Rotating b-a clockwise
            // therefore points out of the hull and orients the dual ray.
            const PointType& a = vertices_[static_cast<std::size_t>(
                triangle.v[static_cast<std::size_t>((side + 1) % 3)])];
            const PointType& b = vertices_[static_cast<std::size_t>(
                triangle.v[static_cast<std::size_t>((side + 2) % 3)])];
            const ResultNumber dx = detail::asNumber<ResultNumber>(b.x()) -
                                    detail::asNumber<ResultNumber>(a.x());
            const ResultNumber dy = detail::asNumber<ResultNumber>(b.y()) -
                                    detail::asNumber<ResultNumber>(a.y());
            const ResultPoint& center = centers[static_cast<std::size_t>(t)];
            dualEdges.emplace_back(
                Ray<ResultPoint>(center, center + ResultPoint(dy, -dx)));
        }
    }
 
    Diagram diagram(dualEdges);
 
    // Site points are strictly inside their own cells. Build the logarithmic
    // point-location index only for this attribution pass, then release it so
    // the returned Arrangement follows the usual opt-in indexing contract.
    diagram.buildPointLocation();
    for (VertexIndex vertex = 1; vertex < static_cast<VertexIndex>(vertices_.size()); ++vertex) {
        const auto face = diagram.locateFace(ResultPoint(vertices_[static_cast<std::size_t>(vertex)]));
        diagram.label(face) = vertices_[static_cast<std::size_t>(vertex)];
    }
    diagram.clearPointLocation();
    return diagram;
}
 
// No deduction guide, deliberately: the vertex type cannot be read off the
// input, because the input's own type is usually the wrong answer. Two integral
// segments cross at a rational point, so a guide reading the type off the input
// would silently pick a vertex type that cannot hold the vertices it is about to
// compute. What `Arrangement(shapes)` deduces instead is the default vertex
// type, which is exact whatever the input is; a caller wanting another one, or
// wanting the edges to carry the input's labels, names it.
 
}  // namespace pgl