triangulation.hpp

Go to the documentation of this file.

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

 
#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <deque>
#include <limits>
#include <map>
#include <memory>
#include <optional>
#include <queue>
#include <random>
#include <set>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <variant>
#include <vector>
 
namespace pgl {
 
namespace detail {
// The entity's LabelType, or NoLabel if it has none. Lets the triangulation
// pick up a triangle label the moment pgl::Triangle exposes one, without
// breaking while it still has none.
template <class T, class = void>
struct optional_label {
    using type = NoLabel;
};
template <class T>
struct optional_label<T, std::void_t<typename T::LabelType>> {
    using type = typename T::LabelType;
};
template <class T>
using optional_label_t = typename optional_label<T>::type;
 
// A directed or undirected segment. The triangulation's traversal queries
// accept either, tracing from endpoint [0] to endpoint [1] (for an oriented
// segment that is source -> target; for a plain segment, its sorted endpoints).
template <class T>
concept SegmentOrOriented = SegmentConcept<T> || OrientedSegmentConcept<T>;
 
// A doubly-infinite straight query: a line or an oriented line. It is traced in
// order by the same directed walk as a segment — the only differences are that
// it enters the hull at a ghost found by a directional descent (rather than at a
// finite source) and never stops at a finite target, so it runs until it leaves
// the hull. An oriented line keeps its own direction; a plain line's order
// follows its (arbitrary) defining-point direction.
template <class T>
concept LineOrOriented = LineConcept<T> || OrientedLineConcept<T>;
 
// Queries traced in order by the directed walk of visitTrianglesIntersecting:
// segments, oriented segments, lines, oriented lines, and rays. A ray keeps a
// finite source (entered like a segment, or at a directional ghost when the
// source is outside the hull) but is unbounded forwards, so like a line it runs
// until it leaves the hull instead of stopping at a target.
template <class T>
concept DirectedTraversal = SegmentOrOriented<T> || LineOrOriented<T> || RayConcept<T>;
 
// A connected convex query shape accepted by the region-traversal overloads of
// visitTrianglesIntersecting / trianglesIntersecting. Segments, oriented
// segments, (oriented) lines, and rays are deliberately excluded — they have
// their own directed walk. The set is convex so its intersection with the
// (convex) triangulated hull stays connected, which is what lets a single seed
// flood-fill find every triangle it meets. (Non-convex polygons are future work.)
template <class T>
concept TriangulationRegionQuery =
    PointConcept<T> || TriangleConcept<T> || RectangleConcept<T> || ConvexConcept<T> ||
    DiskConcept<T> || HalfplaneConcept<T>;
 
// A polygonal chain query: a polyline or a monotone chain. It is neither straight
// (so no single directed walk covers it) nor convex (its intersection with the
// hull may fall apart into several arcs, so no single flood fill covers it
// either), and it gets its own traversal: the directed segment walk run over each
// edge in turn, in chain order.
template <class T>
concept ChainTraversal = PolylineConcept<T> || MonotoneChainConcept<T>;
 
// A query shape bounded by a closed polygonal chain: a triangle, rectangle,
// convex polygon, or (possibly non-convex) polygon. Containment of
// such a shape in the triangulated domain reduces to containment of its boundary
// edges: the domain is a simply connected region (a simple polygon, or the convex
// hull), so its complement is connected — a point of the shape outside the domain
// could be joined to infinity without ever meeting the domain, and that path has
// to cross the shape's boundary, which is inside it. Contradiction.
template <class T>
concept PolygonalRegion =
    TriangleConcept<T> || RectangleConcept<T> || ConvexConcept<T> || PolygonConcept<T>;
 
// Any shape the triangulation's intersection queries accept: a directed-traversal
// shape (traced in order by the directed walk), a chain (traced edge by edge, in
// order), or a region-query shape (grown by the flood fill, unspecified order).
// The derived wrappers — trianglesIntersecting, the interior-intersecting
// variants, the edge variants — are all built on visitTrianglesIntersecting plus
// a per-shape predicate, so they take this whole set; only the three
// visitTrianglesIntersecting overloads carry distinct walks.
template <class T>
concept TriangulationQuery =
    DirectedTraversal<T> || ChainTraversal<T> || TriangulationRegionQuery<T>;
 
// Narrow access point for the paper-specific graph used by
// Polygon::convexCovering(). Defined after Triangulation is complete.
struct ConvexCoverBuilder;
}  // namespace detail

template <TriangleConcept TriangleType_,
          SegmentConcept SegmentType_ = typename TriangleType_::template BoundaryType<false>>
struct Triangulation {
  private:
    // Tags that keep the two handle families distinct types, so a vertex handle
    // can never be passed where a triangle one is meant.
    struct VertexTag;
    struct TriTag;
 
  public:
    using TriangleType = TriangleType_;
    using SegmentType = SegmentType_;
    using PointType = typename TriangleType::PointType;
    using NumberType = typename PointType::NumberType;
    using SegmentLabel = typename SegmentType::LabelType;
    using TriangleLabel = detail::optional_label_t<TriangleType>;
    using VertexId = detail::Handle<VertexTag>;
    using TriId = detail::Handle<TriTag>;

    Triangulation() = default;

    template <class TriangleRange>
        requires TriangleConcept<typename TriangleRange::value_type>
    explicit Triangulation(const TriangleRange& tris) {
        std::unordered_map<PointType, VertexIndex> vid;
        const auto idOfPoint = makeVertexInterner(vid);
        // First pass: intern every vertex so vertices_ is complete, then keep it
        // in Hilbert order for cache locality (see the point-set constructor) and
        // rebuild the point->id map against the reordered vertices.
        for (const auto& tr : tris) {
            idOfPoint(tr[0]);
            idOfPoint(tr[1]);
            idOfPoint(tr[2]);
        }
        hilbertSort(vertices_);
        syncVertexApproximations();
        vid.clear();
        for (VertexIndex i = 0; i < static_cast<VertexIndex>(vertices_.size()); ++i) {
            vid.emplace(vertices_[i], i);
        }
        // Second pass: resolve each triangle's vertices against the reordered ids.
        std::vector<std::array<VertexIndex, 3>> triples;
        std::vector<TriangleLabel> triLabels;
        for (const auto& tr : tris) {
            triples.push_back({vid.at(tr[0]), vid.at(tr[1]), vid.at(tr[2])});
            triLabels.push_back(detail::copyLabel<TriangleLabel>(tr));
        }
        buildFromTriples(triples, triLabels);
    }

    template <class SegmentRange>
        requires SegmentConcept<typename SegmentRange::value_type>
    explicit Triangulation(const SegmentRange& segs) {
        // Intern endpoints so vertices_ is complete, then keep it in Hilbert
        // order for cache locality (see the point-set constructor) and rebuild
        // the point->id map against the reordered vertices.
        std::unordered_map<PointType, VertexIndex> vid;
        const auto idOfPoint = makeVertexInterner(vid);
        for (const auto& s : segs) {
            idOfPoint(s[0]);
            idOfPoint(s[1]);
        }
        hilbertSort(vertices_);
        syncVertexApproximations();
        vid.clear();
        for (VertexIndex i = 0; i < static_cast<VertexIndex>(vertices_.size()); ++i) {
            vid.emplace(vertices_[i], i);
        }
        // Collect the edges as vertex-id pairs against the reordered ids.
        std::vector<std::pair<VertexIndex, VertexIndex>> elist;
        for (const auto& s : segs) {
            VertexIndex a = vid.at(s[0]);
            VertexIndex b = vid.at(s[1]);
            if (a != b) {
                elist.emplace_back(a, b);
            }
        }
        const VertexIndex n = static_cast<VertexIndex>(vertices_.size());
 
        // Undirected adjacency, deduplicated.
        std::vector<std::vector<VertexIndex>> adj(static_cast<std::size_t>(n));
        std::set<std::pair<VertexIndex, VertexIndex>> seen;
        for (auto [a, b] : elist) {
            auto k = a < b ? std::pair{a, b} : std::pair{b, a};
            if (seen.insert(k).second) {
                adj[a].push_back(b);
                adj[b].push_back(a);
            }
        }
 
        // Sort each vertex's neighbors counterclockwise; record their position.
        std::vector<std::unordered_map<VertexIndex, int>> posIn(static_cast<std::size_t>(n));
        for (VertexIndex v = 0; v < n; ++v) {
            std::vector<PointType> nbr;
            nbr.reserve(adj[v].size());
            for (VertexIndex w : adj[v]) {
                nbr.push_back(vertices_[w]);
            }
            sortAround(nbr, vertices_[v]);
            adj[v].clear();
            for (const auto& p : nbr) {
                adj[v].push_back(vid.at(p));
            }
            for (int i = 0; i < static_cast<int>(adj[v].size()); ++i) {
                posIn[v][adj[v][i]] = i;
            }
        }
 
        // next(u->v): the half-edge leaving v that keeps the face on its left,
        // i.e. v -> (neighbor of v immediately clockwise from u). Tracing it
        // walks each face; bounded faces come out CCW, the outer face CW.
        const auto nextHE = [&](VertexIndex u, VertexIndex v) -> std::pair<VertexIndex, VertexIndex> {
            const auto& a = adj[v];
            const int deg = static_cast<int>(a.size());
            const int pu = posIn[v].at(u);
            return {v, a[(pu + deg - 1) % deg]};
        };
 
        std::vector<std::array<VertexIndex, 3>> triples;
        std::set<std::pair<VertexIndex, VertexIndex>> visited;
        for (VertexIndex v = 0; v < n; ++v) {
            for (VertexIndex w : adj[v]) {
                std::pair<VertexIndex, VertexIndex> h{v, w};
                if (visited.count(h)) {
                    continue;
                }
                std::vector<VertexIndex> cycle;
                while (visited.insert(h).second) {
                    cycle.push_back(h.first);
                    h = nextHE(h.first, h.second);
                }
                if (cycle.size() == 3 &&
                    orientationSign(vertices_[cycle[0]], vertices_[cycle[1]], vertices_[cycle[2]]) > 0) {
                    triples.push_back({cycle[0], cycle[1], cycle[2]});
                }
            }
        }
        buildFromTriples(triples, std::vector<TriangleLabel>(triples.size()));
 
        // Carry each input segment's label onto its edge record.
        if constexpr (detail::has_label_v<SegmentLabel>) {
            for (const auto& s : segs) {
                auto it = segmentMap().find(s);
                if (it != segmentMap().end()) {
                    it->second.segLabel = detail::copyLabel<SegmentLabel>(s);
                }
            }
        }
    }

    template <class PointRange>
        requires PointConcept<typename PointRange::value_type>
    explicit Triangulation(const PointRange& pts) {
        std::unordered_map<PointType, VertexIndex> vid;
        const auto idOfPoint = makeVertexInterner(vid);
        for (const auto& p : pts) {
            idOfPoint(PointType(p));
        }
        // Store the vertices in Hilbert-curve order: spatially close points then
        // sit close together both in vertices_ and — because triangles are
        // created in insertion order — in triangles_. That keeps the incremental
        // build's point-location walks short (each is seeded from the previous
        // insertion) and improves cache locality for later query walks too.
        // Vertex order is purely internal, so this is transparent downstream.
        hilbertSort(vertices_);
        syncVertexApproximations();
        auto triples = delaunayTriples(vertices_, vertexApproximations_);
        buildFromTriples(triples, std::vector<TriangleLabel>(triples.size()));
    }

    template <class PointRange, class SegmentRange>
        requires PointConcept<typename PointRange::value_type> &&
                 SegmentConcept<typename SegmentRange::value_type>
    Triangulation(const PointRange& pts, const SegmentRange& segments) {
        std::unordered_map<PointType, VertexIndex> vid;
        const auto idOfPoint = makeVertexInterner(vid);
        for (const auto& p : pts) {
            idOfPoint(PointType(p));
        }
        for (const auto& s : segments) {
            idOfPoint(PointType(s[0]));
            idOfPoint(PointType(s[1]));
        }
        // Keep vertices_ in Hilbert order (see the point-set constructor).
        hilbertSort(vertices_);
        syncVertexApproximations();
        auto triples = delaunayTriples(vertices_, vertexApproximations_);
        buildFromTriples(triples, std::vector<TriangleLabel>(triples.size()));
 
        // Resolve the constraint endpoints against the final ids — the
        // interner's are stale after the reorder and the ghost prepend (see
        // constructConstrained) — then force each segment in as a constrained
        // edge and restore the constrained Delaunay property.
        vid.clear();
        for (VertexIndex i = 1; i < static_cast<VertexIndex>(vertices_.size()); ++i) {
            vid.emplace(vertices_[i], i);
        }
        for (const auto& s : segments) {
            const VertexIndex a = vid.at(PointType(s[0]));
            const VertexIndex b = vid.at(PointType(s[1]));
            if (a != b) {
                insertConstraint(a, b);
            }
        }
        restoreConstrainedDelaunay();
 
        // Carry each constraint segment's label onto its edge record; the
        // edges exist (constrained edges are never flipped away).
        if constexpr (detail::has_label_v<SegmentLabel>) {
            for (const auto& s : segments) {
                auto it = segmentMap().find(SegmentType(PointType(s[0]), PointType(s[1])));
                if (it != segmentMap().end()) {
                    it->second.segLabel = detail::copyLabel<SegmentLabel>(s);
                }
            }
        }
    }

    explicit Triangulation(const Polygon<PointType>& poly) {
        constructConstrained({poly}, std::array<PointType, 0>{}, std::array<SegmentType, 0>{});
    }

    template <class PointRange>
        requires PointConcept<typename PointRange::value_type>
    Triangulation(const Polygon<PointType>& poly, const PointRange& points) {
        constructConstrained({poly}, points, std::array<SegmentType, 0>{});
    }

    template <class SegmentRange>
        requires SegmentConcept<typename SegmentRange::value_type>
    Triangulation(const Polygon<PointType>& poly, const SegmentRange& segments) {
        constructConstrained({poly}, std::array<PointType, 0>{}, segments);
    }

    template <class PointRange, class SegmentRange>
        requires PointConcept<typename PointRange::value_type> &&
                 SegmentConcept<typename SegmentRange::value_type>
    Triangulation(const Polygon<PointType>& poly, const PointRange& points,
                  const SegmentRange& segments) {
        constructConstrained({poly}, points, segments);
    }

    explicit Triangulation(const PolygonWithHoles<PointType>& region) {
        constructConstrained({region.outer()}, std::array<PointType, 0>{},
                             std::array<SegmentType, 0>{}, region.holes());
    }

    template <class PointRange>
        requires PointConcept<typename PointRange::value_type>
    Triangulation(const PolygonWithHoles<PointType>& region, const PointRange& points) {
        constructConstrained({region.outer()}, points, std::array<SegmentType, 0>{}, region.holes());
    }

    template <class SegmentRange>
        requires SegmentConcept<typename SegmentRange::value_type>
    Triangulation(const PolygonWithHoles<PointType>& region, const SegmentRange& segments) {
        constructConstrained({region.outer()}, std::array<PointType, 0>{}, segments, region.holes());
    }

    template <class PointRange, class SegmentRange>
        requires PointConcept<typename PointRange::value_type> &&
                 SegmentConcept<typename SegmentRange::value_type>
    Triangulation(const PolygonWithHoles<PointType>& region, const PointRange& points,
                  const SegmentRange& segments) {
        constructConstrained({region.outer()}, points, segments, region.holes());
    }

    explicit Triangulation(const PolygonSet<PointType>& set) {
        constructConstrained(setOuters(set), std::array<PointType, 0>{},
                             std::array<SegmentType, 0>{}, setHoles(set));
    }

    template <class PointRange>
        requires PointConcept<typename PointRange::value_type>
    Triangulation(const PolygonSet<PointType>& set, const PointRange& points) {
        constructConstrained(setOuters(set), points, std::array<SegmentType, 0>{}, setHoles(set));
    }

    template <class SegmentRange>
        requires SegmentConcept<typename SegmentRange::value_type>
    Triangulation(const PolygonSet<PointType>& set, const SegmentRange& segments) {
        constructConstrained(setOuters(set), std::array<PointType, 0>{}, segments, setHoles(set));
    }

    template <class PointRange, class SegmentRange>
        requires PointConcept<typename PointRange::value_type> &&
                 SegmentConcept<typename SegmentRange::value_type>
    Triangulation(const PolygonSet<PointType>& set, const PointRange& points,
                  const SegmentRange& segments) {
        constructConstrained(setOuters(set), points, segments, setHoles(set));
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] Arrangement<Point<ResultNumber>, PointType> voronoiDiagram() const;
 
    // ---- sizes -----------------------------------------------------------

    [[nodiscard]] std::size_t numVertices() const {
        // A nonempty vertices_ always holds the ghost vertex at index GHOST.
        return vertices_.empty() ? 0 : vertices_.size() - 1;
    }

    [[nodiscard]] std::size_t numTriangles() const { return domainTriangleCount_; }

    [[nodiscard]] std::size_t numEdges() const {
        // Counted off the triangles rather than off the segment map, which
        // holds exactly these edges but which a triangulation asked only for a
        // count should not have to fill. An edge between two real triangles is
        // reached from both and counted from the lower-numbered one; every
        // other edge has one real side and is reached once.
        std::size_t count = 0;
        for (TriIndex t = 0; t < firstGhost_; ++t) {
            for (std::int8_t s = 0; s < 3; ++s) {
                const TriIndex other = triangles_[t].nbr[s];
                if (other != NO_TRI && other < firstGhost_ && other < t) {
                    continue;
                }
                if (edgeInDomain(Edge{t, s})) {
                    ++count;
                }
            }
        }
        return count;
    }

    [[nodiscard]] bool empty() const { return domainTriangleCount_ == 0; }
 
    // ---- membership ------------------------------------------------------

    [[nodiscard]] bool has(const TriangleType& t) const { return inDomain(idOf(t)); }

    [[nodiscard]] bool has(const SegmentType& s) const {
        auto se = segmentMap().find(s);
        return se != segmentMap().end() && edgeInDomain(se->second);
    }
 
    // ---- navigation ------------------------------------------------------

    [[nodiscard]] std::optional<TriangleType> otherTriangle(const TriangleType& t,
                                                       const SegmentType& shared) const {
        auto se = segmentMap().find(shared);
        if (se == segmentMap().end()) {
            return std::nullopt;
        }
        const TriIndex given = idOf(t);
        if (given == NO_TRI) {
            return std::nullopt;
        }
        const Edge e = se->second;
        const TriIndex i1 = e.tri;
        const TriIndex i2 = mirror(e).tri;
        TriIndex other = (given == i1) ? i2 : (given == i2 ? i1 : NO_TRI);
        if (!inDomain(other)) {
            return std::nullopt;  // shared not on t, or boundary
        }
        return triangleValue(other);
    }

    [[nodiscard]] std::vector<TriangleType> edgeAdjacentTriangles(const TriangleType& t) const {
        return trianglesOf(edgeAdjacentTriangles(triHandle(idOf(t))));
    }

    [[nodiscard]] std::vector<TriangleType> vertexAdjacentTriangles(const TriangleType& t) const {
        return trianglesOf(vertexAdjacentTriangles(triHandle(idOf(t))));
    }

    [[nodiscard]] std::vector<TriangleType> incidentTriangles(const SegmentType& s) const {
        std::vector<TriangleType> out;
        auto se = segmentMap().find(s);
        if (se == segmentMap().end()) {
            return out;
        }
        const Edge e = se->second;
        if (inDomain(e.tri)) {
            out.push_back(triangleValue(e.tri));
        }
        const TriIndex other = mirror(e).tri;
        if (inDomain(other)) {
            out.push_back(triangleValue(other));
        }
        return out;
    }

    [[nodiscard]] std::vector<TriangleType> incidentTriangles(const PointType& p) const {
        return trianglesOf(incidentTriangles(getId(p)));
    }

    template <class Fn>
    bool visitTriangles(Fn fn) const {
        for (TriIndex t = 0; t < firstGhost_; ++t) {
            if (!inDomain(t)) {
                continue;
            }
            if (reportTriangle(fn, t)) {
                return true;
            }
        }
        return false;
    }

    template <class Fn>
    bool visitEdges(Fn fn) const {
        for (const auto& [seg, e] : segmentMap()) {
            (void)seg;
            if (!edgeInDomain(e)) {
                continue;
            }
            if (detail::invokeVisitor(fn, edgeSegment(e))) {
                return true;
            }
        }
        return false;
    }

    [[nodiscard]] std::vector<TriangleType> triangles() const {
        std::vector<TriangleType> out;
        out.reserve(numTriangles());
        visitTriangles([&](const TriangleType& t) { out.push_back(t); });
        std::sort(out.begin(), out.end());
        return out;
    }

    [[nodiscard]] std::vector<SegmentType> edges() const {
        std::vector<SegmentType> out;
        out.reserve(numEdges());
        visitEdges([&](const SegmentType& s) { out.push_back(s); });
        std::sort(out.begin(), out.end());
        return out;
    }

    [[nodiscard]] Graph<PointType> asGraph() const {
        Graph<PointType> result;
        for (VertexIndex v = GHOST + 1; v < static_cast<VertexIndex>(vertices_.size()); ++v) {
            result.addVertex(vertices_[v]);
        }
        visitEdges([&](const SegmentType& s) { result.addEdge(s[0], s[1]); });
        return result;
    }

    [[nodiscard]] Arrangement<PointType, TriId> asArrangement() const;

    void buildPointLocation();

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

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

    [[nodiscard]] bool hasCurrentPointLocation() const noexcept {
        return pointLocation_ && pointLocationRevision_ == revision_;
    }
 
    // ---- low-level navigation --------------------------------------------

    [[nodiscard]] std::vector<TriId> triangleIds() const {
        std::vector<TriId> out;
        out.reserve(numTriangles());
        visitTriangles([&](TriId t) { out.push_back(t); });
        return out;
    }

    [[nodiscard]] std::vector<VertexId> vertexIds() const {
        std::vector<VertexId> out;
        out.reserve(numVertices());
        for (VertexIndex v = GHOST + 1; v < static_cast<VertexIndex>(vertices_.size()); ++v) {
            out.push_back(vertexHandle(v));
        }
        return out;
    }

    [[nodiscard]] std::size_t triangleIndexBound() const {
        return static_cast<std::size_t>(firstGhost_);
    }

    [[nodiscard]] std::size_t vertexIndexBound() const { return vertices_.size(); }

    [[nodiscard]] TriId locateId(const PointType& p) const {
        // The hierarchy lands on the triangle itself rather than beside it, so
        // where it does the walk has nothing left to do and is not run: what it
        // would return is a triangle whose closure holds p, which this is.
        bool answered = false;
        const TriIndex seed = pointLocationSeed(p, answered);
        if (answered) {
            hint_ = seed;
            return triHandle(inDomain(seed) ? seed : NO_TRI);
        }
        const TriIndex id = locateIndex(p, seed);
        return triHandle(inDomain(id) ? id : NO_TRI);
    }

    template <PointConcept QueryPoint>
    [[nodiscard]] TriId locateId(const QueryPoint& p) const {
        if constexpr (std::same_as<typename QueryPoint::NumberType, NumberType> &&
                      std::constructible_from<PointType, const QueryPoint&>) {
            return locateId(PointType(p));
        }
        const TriIndex id = locateIndex(p);
        return triHandle(inDomain(id) ? id : NO_TRI);
    }

    [[nodiscard]] TriangleType getShape(TriId t) const {
        assert(has(t) && "getShape(): the handle is not a triangle of the triangulation");
        return triangleValue(indexOf(t));
    }

    [[nodiscard]] const PointType& getShape(VertexId v) const {
        assert(has(v) && "getShape(): the handle is not a vertex of the triangulation");
        return vertices_[static_cast<std::size_t>(indexOf(v))];
    }

    [[nodiscard]] TriangleType operator[](TriId t) const { return getShape(t); }

    [[nodiscard]] const PointType& operator[](VertexId v) const { return getShape(v); }

    [[nodiscard]] TriId getId(const TriangleType& t) const {
        const TriIndex id = idOf(t);
        return triHandle(inDomain(id) ? id : NO_TRI);
    }

    [[nodiscard]] VertexId getId(const PointType& p) const {
        return vertexHandle(vertexIndexAt(p));
    }

    [[nodiscard]] bool has(TriId t) const { return inDomain(indexOf(t)); }

    [[nodiscard]] bool has(VertexId v) const { return realVertex(indexOf(v)); }

    [[nodiscard]] std::array<VertexId, 3> vertices(TriId t) const {
        assert(has(t) && "vertices(): the handle is not a triangle of the triangulation");
        const TriIndex id = indexOf(t);
        const auto& v = triangles_[static_cast<std::size_t>(id)].v;
        const int first = firstVertex(id);
        return {vertexHandle(v[first]), vertexHandle(v[(first + 1) % 3]),
                vertexHandle(v[(first + 2) % 3])};
    }

    [[nodiscard]] std::optional<TriId> otherTriangle(TriId t, int side) const {
        assert(has(t) && "otherTriangle(): the handle is not a triangle of the triangulation");
        assert(side >= 0 && side < 3 && "otherTriangle(): side is not one of 0, 1, 2");
        const TriIndex id = indexOf(t);
        const TriIndex nb = triangles_[static_cast<std::size_t>(id)].nbr[internalSide(id, side)];
        return inDomain(nb) ? std::optional<TriId>(triHandle(nb)) : std::nullopt;
    }

    [[nodiscard]] bool isConstrained(TriId t, int side) const {
        assert(has(t) && "isConstrained(): the handle is not a triangle of the triangulation");
        assert(side >= 0 && side < 3 && "isConstrained(): side is not one of 0, 1, 2");
        const TriIndex id = indexOf(t);
        return bit(triangles_[static_cast<std::size_t>(id)].constrainedMask,
                   internalSide(id, side));
    }

    void setConstrained(TriId t, int side, bool value = true) {
        assert(has(t) && "setConstrained(): the handle is not a triangle of the triangulation");
        assert(side >= 0 && side < 3 && "setConstrained(): side is not one of 0, 1, 2");
        const TriIndex id = indexOf(t);
        const Edge e{id, static_cast<std::int8_t>(internalSide(id, side))};
        setBit(triangles_[static_cast<std::size_t>(id)].constrainedMask, e.side, value);
        const Edge m = mirror(e);
        if (m.tri != NO_TRI) {
            setBit(triangles_[static_cast<std::size_t>(m.tri)].constrainedMask, m.side, value);
        }
    }

    [[nodiscard]] std::optional<TriId> otherTriangle(TriId t, VertexId a, VertexId b) const {
        const TriIndex id = indexOf(t);
        if (!inDomain(id)) {
            return std::nullopt;
        }
        const VertexIndex first = indexOf(a);
        const VertexIndex second = indexOf(b);
        const auto& v = triangles_[static_cast<std::size_t>(id)].v;
        for (int s = 0; s < 3; ++s) {
            const VertexIndex left = v[(s + 1) % 3];
            const VertexIndex right = v[(s + 2) % 3];
            if ((left == first && right == second) || (left == second && right == first)) {
                const TriIndex nb = triangles_[static_cast<std::size_t>(id)].nbr[s];
                return inDomain(nb) ? std::optional<TriId>(triHandle(nb)) : std::nullopt;
            }
        }
        return std::nullopt;  // a b is not an edge of t
    }

    [[nodiscard]] std::vector<TriId> edgeAdjacentTriangles(TriId t) const {
        std::vector<TriId> out;
        const TriIndex id = indexOf(t);
        if (!realTriangle(id)) {
            return out;
        }
        for (int s = 0; s < 3; ++s) {
            const TriIndex nb = triangles_[static_cast<std::size_t>(id)].nbr[s];
            if (inDomain(nb)) {
                out.push_back(triHandle(nb));
            }
        }
        return out;
    }

    [[nodiscard]] std::vector<TriId> vertexAdjacentTriangles(TriId t) const {
        std::vector<TriId> out;
        const TriIndex self = indexOf(t);
        if (!realTriangle(self)) {
            return out;
        }
        // "Already listed" is tracked by the per-triangle walkMark bit (as in
        // visitTrianglesIntersecting): marking `self` first skips t itself, and
        // the guard clears every set bit on the way out so the const query
        // leaves the triangulation pristine.
        std::vector<TriIndex> marked{self};
        triangles_[self].walkMark = 1;
        struct MarkClearer {
            const std::vector<Tri>& tris;
            const std::vector<TriIndex>& marked;
            ~MarkClearer() { for (TriIndex m : marked) tris[m].walkMark = 0; }
        } markClearer{triangles_, marked};
        for (const VertexIndex w : triangles_[self].v) {
            visitVertexFan(self, w, [&](TriIndex cur) {
                if (inDomain(cur) && !triangles_[cur].walkMark) {
                    triangles_[cur].walkMark = 1;
                    marked.push_back(cur);
                    out.push_back(triHandle(cur));
                }
            });
        }
        return out;
    }

    [[nodiscard]] std::vector<TriId> incidentTriangles(VertexId v) const {
        std::vector<TriId> out;
        const VertexIndex w = indexOf(v);
        if (!realVertex(w)) {
            return out;
        }
        const TriIndex start = fanSeedOf(w);
        if (start == NO_TRI) {
            return out;
        }
        visitVertexFan(start, w, [&](TriIndex cur) {
            if (inDomain(cur)) {
                out.push_back(triHandle(cur));
            }
        });
        return out;
    }

    template <class L = TriangleLabel>
        requires(detail::has_label_v<L>)
    [[nodiscard]] L& label(TriId t) {
        assert(has(t) && "label(): the handle is not a triangle of the triangulation");
        return triangles_[static_cast<std::size_t>(indexOf(t))].triLabel;
    }

    template <class L = TriangleLabel>
        requires(detail::has_label_v<L>)
    [[nodiscard]] const L& label(TriId t) const {
        assert(has(t) && "label(): the handle is not a triangle of the triangulation");
        return triangles_[static_cast<std::size_t>(indexOf(t))].triLabel;
    }
 
    // ---- visibility ------------------------------------------------------

    [[nodiscard]] Graph<PointType> visibilityGraph() const;

    [[nodiscard]] Graph<PointType> clearVisibilityGraph() const;

    [[nodiscard]] Graph<PointType> reducedVisibilityGraph() const;

    [[nodiscard]] std::vector<PointType> visibleVertices(const PointType& query) const;

    [[nodiscard]] std::vector<PointType> clearlyVisibleVertices(const PointType& query) const;

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] Polygon<Point<ResultNumber>> regularizedVisiblePolygon(
        const PointType& query) const;
 
  private:
    friend struct detail::ConvexCoverBuilder;
 
    // Full visibility of two mesh triangles needs only the new boundary edges
    // introduced by their convex hull. A hull edge whose endpoints both belong
    // to one input triangle already lies in that triangle and therefore in the
    // domain; only the (at most two) bridges between the triangles need a mesh
    // walk. The generic contains(Convex) path would locate both endpoints and
    // walk every hull edge, including all of those known-contained chains.
    [[nodiscard]] bool trianglesFullyVisible(
        std::int32_t a, std::int32_t b, const std::vector<std::uint8_t>& visible) const {
        const TriangleType first = triangleValue(a);
        const TriangleType second = triangleValue(b);
        std::array<PointType, 6> points{
            first.a(), first.b(), first.c(), second.a(), second.b(), second.c()
        };
        std::sort(points.begin(), points.end());
        const auto uniqueEnd = std::unique(points.begin(), points.end());
        const std::size_t pointCount =
            static_cast<std::size_t>(std::distance(points.begin(), uniqueEnd));
 
        // Andrew's monotone chain on six points, kept entirely on the stack.
        // Twelve slots cover the lower and upper chains before their repeated
        // endpoint is discarded.
        std::array<PointType, 12> hull;
        std::size_t hullSize = 0;
        for (std::size_t i = 0; i < pointCount; ++i) {
            while (hullSize >= 2 &&
                   orientationSign(hull[hullSize - 2], hull[hullSize - 1], points[i]) <= 0) {
                --hullSize;
            }
            hull[hullSize++] = points[i];
        }
        const std::size_t lowerSize = hullSize;
        for (std::size_t i = pointCount - 1; i-- != 0;) {
            while (hullSize > lowerSize &&
                   orientationSign(hull[hullSize - 2], hull[hullSize - 1], points[i]) <= 0) {
                --hullSize;
            }
            hull[hullSize++] = points[i];
        }
        assert(hullSize >= 2);
        --hullSize;  // the first vertex was repeated at the end
 
        const auto isVertexOf = [&](TriIndex triangle, const PointType& point) {
            const auto& ids = triangles_[triangle].v;
            return vertices_[ids[0]] == point || vertices_[ids[1]] == point ||
                   vertices_[ids[2]] == point;
        };
 
        // The candidate triangle is already visible if each of its sides that
        // became internal to the joint hull borders a visible triangle. Their
        // union fills the candidate-facing neighborhood of every such side.
        // Conversely, an internalized side on the domain boundary proves that
        // the hull escapes immediately. These are the paper implementation's
        // inexpensive positive and negative tests.
        bool boundedByVisibleTriangles = true;
        for (std::int8_t side = 0; side < 3; ++side) {
            const PointType u = vertices_[triangles_[b].v[(side + 1) % 3]];
            const PointType v = vertices_[triangles_[b].v[(side + 2) % 3]];
            bool onHullBoundary = false;
            for (std::size_t i = 0; i < hullSize; ++i) {
                const Segment<PointType> hullSide(hull[i], hull[(i + 1) % hullSize]);
                if (hullSide.contains(u) && hullSide.contains(v)) {
                    onHullBoundary = true;
                    break;
                }
            }
            if (onHullBoundary) {
                continue;
            }
            const TriIndex across = triangles_[b].nbr[side];
            if (!inDomain(across)) {
                return false;
            }
            if (!visible[static_cast<std::size_t>(across)]) {
                boundedByVisibleTriangles = false;
            }
        }
        if (boundedByVisibleTriangles && holeWitnesses_.empty()) {
            return true;
        }
 
        for (std::size_t i = 0; i < hullSize; ++i) {
            const PointType u = hull[i];
            const PointType v = hull[(i + 1) % hullSize];
            if ((isVertexOf(a, u) && isVertexOf(a, v)) ||
                (isVertexOf(b, u) && isVertexOf(b, v))) {
                continue;
            }
            if (!segmentInteriorContained(Segment<PointType>(u, v))) {
                return false;
            }
        }
 
        // This helper is currently used by Polygon::convexCovering(), whose
        // domain has no holes. Retain the general triangulation semantics here:
        // a contained boundary must not have enclosed a hole.
        if (!holeWitnesses_.empty()) {
            const std::vector<PointType> hullVertices(hull.begin(), hull.begin() + hullSize);
            const Convex<PointType> convexHull(hullVertices, true);
            for (const TriangleType& witness : holeWitnesses_) {
                if (convexHull.contains(witness)) {
                    return false;
                }
            }
        }
        return true;
    }

    [[nodiscard]] Graph<TriangleType> convexCoverVisibilityGraph() const {
        Graph<TriangleType> result;
        std::vector<TriIndex> sources;
        sources.reserve(domainTriangleCount_);
        for (TriIndex source = 0; source < firstGhost_; ++source) {
            if (inDomain(source)) {
                sources.push_back(source);
                result.addVertex(triangleValue(source));
            }
        }
 
        // A depth-first source order tends to make a previously processed
        // source a neighbor of the next one. Its already known symmetric edges
        // can then seed more of the next BFS without another geometric test.
        std::vector<TriIndex> sourceOrder;
        sourceOrder.reserve(sources.size());
        std::vector<bool> ordered(static_cast<std::size_t>(firstGhost_), false);
        std::vector<TriIndex> stack;
        for (const TriIndex seed : sources) {
            if (ordered[static_cast<std::size_t>(seed)]) {
                continue;
            }
            stack.push_back(seed);
            while (!stack.empty()) {
                const TriIndex current = stack.back();
                stack.pop_back();
                if (ordered[static_cast<std::size_t>(current)]) {
                    continue;
                }
                ordered[static_cast<std::size_t>(current)] = true;
                sourceOrder.push_back(current);
                for (const TriIndex neighbor : triangles_[current].nbr) {
                    if (inDomain(neighbor) &&
                        !ordered[static_cast<std::size_t>(neighbor)]) {
                        stack.push_back(neighbor);
                    }
                }
            }
        }
 
        std::vector<std::uint8_t> visited(static_cast<std::size_t>(firstGhost_));
        std::vector<std::uint8_t> visible(static_cast<std::size_t>(firstGhost_));
        std::queue<TriIndex> queue;
        for (const TriIndex source : sourceOrder) {
            std::fill(visited.begin(), visited.end(), false);
            std::fill(visible.begin(), visible.end(), false);
            visited[static_cast<std::size_t>(source)] = true;
            visible[static_cast<std::size_t>(source)] = true;
            queue.push(source);
            const TriangleType sourceTriangle = triangleValue(source);
 
            while (!queue.empty()) {
                const TriIndex current = queue.front();
                queue.pop();
                for (const TriIndex neighbor : triangles_[current].nbr) {
                    if (!inDomain(neighbor) ||
                        visited[static_cast<std::size_t>(neighbor)]) {
                        continue;
                    }
                    visited[static_cast<std::size_t>(neighbor)] = true;
                    const TriangleType neighborTriangle = triangleValue(neighbor);
 
                    bool fullyVisible = result.containsEdge(sourceTriangle, neighborTriangle);
                    if (!fullyVisible) {
                        fullyVisible = trianglesFullyVisible(source, neighbor, visible);
                    }
                    if (fullyVisible) {
                        result.addEdge(sourceTriangle, neighborTriangle);
                        visible[static_cast<std::size_t>(neighbor)] = true;
                        queue.push(neighbor);
                    }
                }
            }
        }
 
        return result;
    }
 
  public:
    // ---- convex partition ------------------------------------------------

    [[nodiscard]] std::vector<Convex<PointType>> convexPartition() const {
        std::vector<Convex<PointType>> pieces;
        if (domainTriangleCount_ == 0) {
            return pieces;
        }
 
        // Pieces are grown by union-find over the triangles, and each piece
        // carries its boundary as a circular list of half-edges. A half-edge is
        // one (triangle, side) pair — `3 * tri + side` — which is exactly the
        // directed edge `v[(side+1)%3] -> v[(side+2)%3]`, CCW around the triangle
        // because the vertices are. A triangle starts as its own piece, so its
        // three sides start as its ring, in that same CCW order.
        std::vector<TriIndex> parent(static_cast<std::size_t>(firstGhost_));
        for (TriIndex t = 0; t < firstGhost_; ++t) {
            parent[static_cast<std::size_t>(t)] = t;
        }
        const auto findRoot = [&parent](TriIndex x) {
            while (parent[static_cast<std::size_t>(x)] != x) {
                parent[static_cast<std::size_t>(x)] =
                    parent[static_cast<std::size_t>(parent[static_cast<std::size_t>(x)])];
                x = parent[static_cast<std::size_t>(x)];
            }
            return x;
        };
 
        const std::size_t slots = static_cast<std::size_t>(firstGhost_) * 3;
        std::vector<std::int32_t> succ(slots, -1);
        std::vector<std::int32_t> pred(slots, -1);
        std::vector<std::int32_t> diagonals;
        for (TriIndex t = 0; t < firstGhost_; ++t) {
            if (!inDomain(t)) {
                continue;
            }
            for (int s = 0; s < 3; ++s) {
                const std::int32_t h = 3 * t + s;
                succ[static_cast<std::size_t>(h)] = 3 * t + (s + 1) % 3;
                pred[static_cast<std::size_t>(3 * t + (s + 1) % 3)] = h;
                // Each interior unconstrained edge is a candidate once, from the
                // lower-numbered of the two triangles that meet along it.
                const TriIndex n = triangles_[static_cast<std::size_t>(t)].nbr[s];
                if (n > t && inDomain(n) &&
                    !bit(triangles_[static_cast<std::size_t>(t)].constrainedMask, s)) {
                    diagonals.push_back(h);
                }
            }
        }
 
        const auto originOf = [this](std::int32_t h) {
            return triangles_[static_cast<std::size_t>(h / 3)].v[static_cast<std::size_t>((h % 3 + 1) % 3)];
        };
        const auto targetOf = [this](std::int32_t h) {
            return triangles_[static_cast<std::size_t>(h / 3)].v[static_cast<std::size_t>((h % 3 + 2) % 3)];
        };
        const auto neighborOf = [this](std::int32_t h) {
            return triangles_[static_cast<std::size_t>(h / 3)].nbr[static_cast<std::size_t>(h % 3)];
        };
        // How many edges the piece holding `start` shares with the piece rooted
        // at `otherRoot`. Two convex sets with disjoint interiors are separated by
        // a line, so they can only share more than one edge if those edges are
        // collinear — which takes a slit or collinear vertices, and is rare
        // enough to pay a ring walk for. Splicing across one of two shared edges
        // would leave a boundary pinched shut at the other, which is not a ring.
        const auto sharedEdges = [&](std::int32_t start, TriIndex otherRoot) {
            int count = 0;
            std::int32_t h = start;
            do {
                const TriIndex n = neighborOf(h);
                if (n != NO_TRI && inDomain(n) && findRoot(n) == otherRoot) {
                    ++count;
                }
                h = succ[static_cast<std::size_t>(h)];
            } while (h != start);
            return count;
        };
 
        for (const std::int32_t h : diagonals) {
            const TriIndex t = h / 3;
            const TriIndex n = neighborOf(h);
            const TriIndex rootA = findRoot(t);
            const TriIndex rootB = findRoot(n);
            if (rootA == rootB) {
                continue;  // an earlier merge already swallowed this diagonal
            }
            const std::int32_t twin = 3 * n + findSide(n, t);
            // The ring runs ... -> p -> u -> v -> q -> ... on this side and
            // ... -> r -> v -> u -> w -> ... on the other, so deleting the two
            // half-edges u->v and v->u leaves u followed by w and v followed by q.
            const VertexIndex u = originOf(h);
            const VertexIndex v = targetOf(h);
            const VertexIndex p = originOf(pred[static_cast<std::size_t>(h)]);
            const VertexIndex q = targetOf(succ[static_cast<std::size_t>(h)]);
            const VertexIndex r = originOf(pred[static_cast<std::size_t>(twin)]);
            const VertexIndex w = targetOf(succ[static_cast<std::size_t>(twin)]);
            const auto convexAt = [this](VertexIndex before, VertexIndex at, VertexIndex after) {
                return orientationSign(vertices_[static_cast<std::size_t>(before)],
                                       vertices_[static_cast<std::size_t>(at)],
                                       vertices_[static_cast<std::size_t>(after)]) >= 0;
            };
            // Collinear counts as convex: the merged piece is still convex and
            // carries one fewer corner, which is the whole point.
            if (!convexAt(p, u, w) || !convexAt(r, v, q)) {
                continue;
            }
            if (sharedEdges(h, rootB) != 1) {
                continue;
            }
            const std::int32_t beforeU = pred[static_cast<std::size_t>(h)];
            const std::int32_t afterU = succ[static_cast<std::size_t>(twin)];
            const std::int32_t beforeV = pred[static_cast<std::size_t>(twin)];
            const std::int32_t afterV = succ[static_cast<std::size_t>(h)];
            succ[static_cast<std::size_t>(beforeU)] = afterU;
            pred[static_cast<std::size_t>(afterU)] = beforeU;
            succ[static_cast<std::size_t>(beforeV)] = afterV;
            pred[static_cast<std::size_t>(afterV)] = beforeV;
            succ[static_cast<std::size_t>(h)] = -1;
            pred[static_cast<std::size_t>(h)] = -1;
            succ[static_cast<std::size_t>(twin)] = -1;
            pred[static_cast<std::size_t>(twin)] = -1;
            parent[static_cast<std::size_t>(rootA)] = rootB;
        }
 
        // Read each surviving ring out once, in the canonical form a trusted
        // `Convex` asks for: CCW, no vertex in the middle of a straight stretch,
        // starting at the lexicographically smallest. The rings are convex, so a
        // vertex is interior to a straight stretch exactly when its own two
        // neighbours in the ring are collinear with it, and dropping every such
        // vertex in one pass is correct however many run together.
        std::vector<char> read(slots, 0);
        std::vector<PointType> ring;
        std::vector<PointType> kept;
        for (std::int32_t start = 0; start < static_cast<std::int32_t>(slots); ++start) {
            if (succ[static_cast<std::size_t>(start)] < 0 || read[static_cast<std::size_t>(start)]) {
                continue;
            }
            ring.clear();
            std::int32_t h = start;
            do {
                read[static_cast<std::size_t>(h)] = 1;
                ring.push_back(vertices_[static_cast<std::size_t>(originOf(h))]);
                h = succ[static_cast<std::size_t>(h)];
            } while (h != start);
 
            kept.clear();
            for (std::size_t i = 0; i < ring.size(); ++i) {
                const PointType& before = ring[(i + ring.size() - 1) % ring.size()];
                const PointType& after = ring[(i + 1) % ring.size()];
                if (orientationSign(before, ring[i], after) != 0) {
                    kept.push_back(ring[i]);
                }
            }
            if (kept.size() < 3) {
                continue;  // no area, so nothing of the domain to carry
            }
            std::rotate(kept.begin(), std::min_element(kept.begin(), kept.end()), kept.end());
            pieces.push_back(Convex<PointType>(kept, /*trusted=*/true));
        }
        std::sort(pieces.begin(), pieces.end());
        return pieces;
    }

    [[nodiscard]] std::vector<Convex<PointType>> convexCovering() const {
        std::vector<Convex<PointType>> result;
        if (domainTriangleCount_ == 0) {
            return result;
        }
 
        const std::size_t triangleSlots = static_cast<std::size_t>(firstGhost_);
        const std::size_t halfEdgeSlots = triangleSlots * 3;
        const auto originOf = [this](std::int32_t h) {
            return triangles_[static_cast<std::size_t>(h / 3)]
                .v[static_cast<std::size_t>((h % 3 + 1) % 3)];
        };
        const auto targetOf = [this](std::int32_t h) {
            return triangles_[static_cast<std::size_t>(h / 3)]
                .v[static_cast<std::size_t>((h % 3 + 2) % 3)];
        };
        const auto convexAt = [this](VertexIndex before, VertexIndex at, VertexIndex after) {
            return orientationSign(vertices_[static_cast<std::size_t>(before)],
                                   vertices_[static_cast<std::size_t>(at)],
                                   vertices_[static_cast<std::size_t>(after)]) >= 0;
        };
 
        struct Candidate {
            std::vector<TriIndex> triangles;
            Convex<PointType> convex;
        };
        std::vector<Candidate> candidates;
        candidates.reserve(domainTriangleCount_);
 
        for (TriIndex seed = 0; seed < firstGhost_; ++seed) {
            if (!inDomain(seed)) {
                continue;
            }
 
            // Only the growing candidate has a live ring. Newly fused triangles
            // receive their original three-edge ring immediately before it is
            // spliced into that candidate.
            std::vector<std::int32_t> succ(halfEdgeSlots, -1);
            std::vector<std::int32_t> pred(halfEdgeSlots, -1);
            std::vector<char> fused(triangleSlots, 0);
            std::vector<TriIndex> fusedTriangles{seed};
            fused[static_cast<std::size_t>(seed)] = 1;
 
            const auto initializeRing = [&](TriIndex t) {
                for (int s = 0; s < 3; ++s) {
                    const std::int32_t h = 3 * t + s;
                    const std::int32_t next = 3 * t + (s + 1) % 3;
                    succ[static_cast<std::size_t>(h)] = next;
                    pred[static_cast<std::size_t>(next)] = h;
                }
            };
            initializeRing(seed);
 
            std::deque<std::int32_t> frontier;
            for (int s = 0; s < 3; ++s) {
                frontier.push_back(3 * seed + s);
            }
 
            while (!frontier.empty()) {
                const std::int32_t h = frontier.front();
                frontier.pop_front();
                if (succ[static_cast<std::size_t>(h)] < 0) {
                    continue;  // removed by an earlier fusion
                }
 
                const TriIndex t = h / 3;
                const int side = h % 3;
                const TriIndex neighbor = triangles_[static_cast<std::size_t>(t)].nbr[side];
                if (!inDomain(neighbor) || fused[static_cast<std::size_t>(neighbor)] ||
                    bit(triangles_[static_cast<std::size_t>(t)].constrainedMask, side)) {
                    continue;
                }
 
                // Splicing one pair of twins is valid only while this triangle
                // meets the candidate along exactly that one edge. Once it meets
                // along two edges, that count can only grow, so it need not be
                // reconsidered.
                int sharedEdges = 0;
                for (const TriIndex adjacent :
                     triangles_[static_cast<std::size_t>(neighbor)].nbr) {
                    if (inDomain(adjacent) && fused[static_cast<std::size_t>(adjacent)]) {
                        ++sharedEdges;
                    }
                }
                if (sharedEdges != 1) {
                    continue;
                }
 
                const std::int32_t twin = 3 * neighbor + findSide(neighbor, t);
                initializeRing(neighbor);
 
                const VertexIndex u = originOf(h);
                const VertexIndex v = targetOf(h);
                const VertexIndex p = originOf(pred[static_cast<std::size_t>(h)]);
                const VertexIndex q = targetOf(succ[static_cast<std::size_t>(h)]);
                const VertexIndex r = originOf(pred[static_cast<std::size_t>(twin)]);
                const VertexIndex w = targetOf(succ[static_cast<std::size_t>(twin)]);
                if (!convexAt(p, u, w) || !convexAt(r, v, q)) {
                    // Initializing a rejected triangle must not leave a live ring.
                    for (int s = 0; s < 3; ++s) {
                        const std::size_t rejected =
                            static_cast<std::size_t>(3 * neighbor + s);
                        succ[rejected] = -1;
                        pred[rejected] = -1;
                    }
                    continue;
                }
 
                const std::int32_t beforeU = pred[static_cast<std::size_t>(h)];
                const std::int32_t afterU = succ[static_cast<std::size_t>(twin)];
                const std::int32_t beforeV = pred[static_cast<std::size_t>(twin)];
                const std::int32_t afterV = succ[static_cast<std::size_t>(h)];
                succ[static_cast<std::size_t>(beforeU)] = afterU;
                pred[static_cast<std::size_t>(afterU)] = beforeU;
                succ[static_cast<std::size_t>(beforeV)] = afterV;
                pred[static_cast<std::size_t>(afterV)] = beforeV;
                succ[static_cast<std::size_t>(h)] = -1;
                pred[static_cast<std::size_t>(h)] = -1;
                succ[static_cast<std::size_t>(twin)] = -1;
                pred[static_cast<std::size_t>(twin)] = -1;
                fused[static_cast<std::size_t>(neighbor)] = 1;
                fusedTriangles.push_back(neighbor);
 
                // These are precisely the live edges whose endpoint angles are
                // new, including both newly exposed sides of `neighbor`.
                frontier.push_back(beforeU);
                frontier.push_back(afterU);
                frontier.push_back(beforeV);
                frontier.push_back(afterV);
            }
 
            std::int32_t start = -1;
            for (std::int32_t h = 0; h < static_cast<std::int32_t>(halfEdgeSlots); ++h) {
                if (succ[static_cast<std::size_t>(h)] >= 0) {
                    start = h;
                    break;
                }
            }
            assert(start >= 0);
 
            std::vector<PointType> ring;
            std::int32_t h = start;
            do {
                ring.push_back(vertices_[static_cast<std::size_t>(originOf(h))]);
                h = succ[static_cast<std::size_t>(h)];
            } while (h != start);
 
            std::vector<PointType> kept;
            kept.reserve(ring.size());
            for (std::size_t i = 0; i < ring.size(); ++i) {
                const PointType& before = ring[(i + ring.size() - 1) % ring.size()];
                const PointType& after = ring[(i + 1) % ring.size()];
                if (orientationSign(before, ring[i], after) != 0) {
                    kept.push_back(ring[i]);
                }
            }
            assert(kept.size() >= 3);
            std::rotate(kept.begin(), std::min_element(kept.begin(), kept.end()), kept.end());
            candidates.push_back(
                Candidate{std::move(fusedTriangles), Convex<PointType>(kept, /*trusted=*/true)});
        }
 
        // Maintain each candidate's uncovered gain through the inverse
        // triangle-to-candidate incidence lists. Every incidence is processed at
        // most once, when its triangle first becomes covered.
        std::vector<std::vector<std::size_t>> coverers(triangleSlots);
        std::vector<std::size_t> gain(candidates.size());
        for (std::size_t c = 0; c < candidates.size(); ++c) {
            gain[c] = candidates[c].triangles.size();
            for (const TriIndex t : candidates[c].triangles) {
                coverers[static_cast<std::size_t>(t)].push_back(c);
            }
        }
 
        std::vector<char> covered(triangleSlots, 0);
        std::vector<char> selected(candidates.size(), 0);
        std::vector<std::size_t> selectionOrder;
        std::size_t remaining = domainTriangleCount_;
        while (remaining != 0) {
            std::size_t best = candidates.size();
            std::size_t bestGain = 0;
            for (std::size_t c = 0; c < candidates.size(); ++c) {
                if (!selected[c] && gain[c] > bestGain) {
                    best = c;
                    bestGain = gain[c];
                }
            }
            assert(best != candidates.size() && bestGain != 0);
            selected[best] = 1;
            selectionOrder.push_back(best);
            for (const TriIndex t : candidates[best].triangles) {
                const std::size_t ti = static_cast<std::size_t>(t);
                if (covered[ti]) {
                    continue;
                }
                covered[ti] = 1;
                --remaining;
                for (const std::size_t c : coverers[ti]) {
                    assert(gain[c] != 0);
                    --gain[c];
                }
            }
        }
 
        std::vector<std::size_t> coverageCount(triangleSlots, 0);
        for (const std::size_t c : selectionOrder) {
            for (const TriIndex t : candidates[c].triangles) {
                ++coverageCount[static_cast<std::size_t>(t)];
            }
        }
        for (auto it = selectionOrder.rbegin(); it != selectionOrder.rend(); ++it) {
            const std::size_t c = *it;
            const bool redundant = std::ranges::all_of(
                candidates[c].triangles, [&](TriIndex t) {
                    return coverageCount[static_cast<std::size_t>(t)] > 1;
                });
            if (!redundant) {
                continue;
            }
            selected[c] = 0;
            for (const TriIndex t : candidates[c].triangles) {
                --coverageCount[static_cast<std::size_t>(t)];
            }
        }
 
        result.reserve(selectionOrder.size());
        for (const std::size_t c : selectionOrder) {
            if (selected[c]) {
                result.push_back(std::move(candidates[c].convex));
            }
        }
        std::sort(result.begin(), result.end());
        return result;
    }
 
    // ---- traversal along a segment or (oriented) line --------------------

    template <detail::DirectedTraversal OS, class Fn>
    bool visitTrianglesIntersecting(const OS& s, Fn f) const {
        return visitTriangleIdsIntersecting<false>(
            s, [&](TriIndex t) { return reportTriangle(f, t); });
    }
 
  private:
    template <bool AllTriangles, detail::DirectedTraversal OS, class Fn>
    bool visitTriangleIdsIntersecting(const OS& s, Fn f) const {
        if (firstGhost_ == 0) return false;
        const auto a = s[0];
        const auto b = s[1];
        // The query's two defining points are read by every sign this walk
        // takes — a dozen of them, once per triangle of the path and once per
        // edge of each — so they are converted here rather than at each.
        const auto fa = filteredPoint(a);
        const auto fb = filteredPoint(b);
        // How the directed query extends past its defining points `a`, `b` (which
        // always give its supporting line and forward direction a->b):
        //  - unboundedBack: no finite source (a line / oriented line). Entered at
        //    the ghost where a->b crosses into the hull, found by a directional
        //    descent, rather than by locating a source.
        //  - unboundedFront: no finite target (a line / oriented line / ray). The
        //    walk never stops at `b`; it runs until it leaves the hull, and picks
        //    each forward exit edge by direction rather than by segment bounds.
        constexpr bool unboundedBack = detail::LineOrOriented<OS>;
        constexpr bool unboundedFront = detail::LineOrOriented<OS> || RayConcept<OS>;
 
        // "Already emitted" is tracked by a per-triangle mark bit rather than a
        // hash set: O(1) test, no per-lookup allocation. `marked` records which
        // bits were set; the guard clears every one of them on the way out —
        // normal return, early exit, or an exception from f — so the const walk
        // leaves the triangulation pristine. NOTE: because the mark lives on the
        // triangulation, the walk is not reentrant: f must not start another
        // walk on the same Triangulation.
        std::vector<TriIndex> marked;
        struct MarkClearer {
            const std::vector<Tri>& tris;
            const std::vector<TriIndex>& marked;
            ~MarkClearer() { for (TriIndex t : marked) tris[t].walkMark = 0; }
        } markClearer{triangles_, marked};
        bool stop = false;
        // The triangle the walk came from, so it never steps straight back out
        // through the edge it just crossed. Set by the entry helpers below too.
        TriIndex prev = NO_TRI;
        // Emit a triangle to f, once. Ghost triangles are walked through for
        // navigation but never reported here, and out-of-domain (polygon
        // hull-fill) ones only under AllTriangles, so the reported set is exactly
        // the (in-domain, or all real) triangles s meets.
        const auto emit = [&](TriIndex t) {
            const bool report = AllTriangles ? (t != NO_TRI && !isGhost(t)) : inDomain(t);
            if (!stop && report && !triangles_[t].walkMark) {
                triangles_[t].walkMark = 1;
                marked.push_back(t);
                if (f(t)) stop = true;
            }
        };
        // The next triangle when rotating around vertex w, having arrived from
        // `from` (NO_TRI to start). Orientation-independent (ghost triangles are
        // not CCW-normalized): leave through the w-incident edge we didn't enter.
        const auto rotateAround = [&](TriIndex cur, VertexIndex w, TriIndex from) -> TriIndex {
            const int lw = localIndex(cur, w);
            int s1 = -1, s2 = -1;
            for (int s = 0; s < 3; ++s) {
                if (s != lw) {
                    if (s1 < 0) s1 = s; else s2 = s;
                }
            }
            return triangles_[cur].nbr[triangles_[cur].nbr[s1] == from ? s2 : s1];
        };
        // Every real triangle of vertex w's fan (rotate around w through ghosts).
        const auto emitFan = [&](VertexIndex w, TriIndex startTri) {
            TriIndex cur = startTri, from = NO_TRI;
            std::size_t g = 0, lim = triangles_.size() + 1;
            do {
                emit(cur);
                const TriIndex next = rotateAround(cur, w, from);
                from = cur;
                cur = next;
            } while (cur != startTri && cur != NO_TRI && !stop && ++g < lim);
        };
        // Continuing forward from p — a point on the closure of t — does the query
        // enter t's interior? 1 = enters, 0 = no, -1 = runs collinearly along an
        // edge of t. Forward is the *direction* a->b, not "towards the point b":
        // b is the far end only for a segment, while a line or a ray runs past its
        // second defining point, and beyond it "towards b" points backwards.
        // cross(w - u, b - a) is that direction's side of the edge (u,w), and the
        // triangle being CCW its interior is the positive side.
        const auto rayEnters = [&](TriIndex t, const auto& p) -> int {
            const auto& v = triangles_[t].v;
            const auto fp = filteredPoint(p);
            for (int k = 0; k < 3; ++k) {
                const auto& u = vertices_[v[(k + 1) % 3]];
                const auto& w = vertices_[v[(k + 2) % 3]];
                if (detail::orientationSignOf(filteredVertex(v[(k + 1) % 3]),
                                              filteredVertex(v[(k + 2) % 3]), fp).value() == 0) {
                    const auto forward =
                        orientationDeterminant(u, w, b) - orientationDeterminant(u, w, a);
                    if (forward < 0) return 0;
                    if (forward == 0) return -1;
                }
            }
            return 1;
        };
        // The fan triangle around w whose interior the ray w->b enters (or NO_TRI).
        const auto forwardAround = [&](VertexIndex w, TriIndex startTri) -> TriIndex {
            TriIndex cur = startTri, from = NO_TRI;
            std::size_t g = 0, lim = triangles_.size() + 1;
            do {
                if (!isGhost(cur) && rayEnters(cur, vertices_[w]) == 1) return cur;
                const TriIndex next = rotateAround(cur, w, from);
                from = cur;
                cur = next;
            } while (cur != startTri && cur != NO_TRI && ++g < lim);
            return NO_TRI;
        };
        // Three-way order of two points along the directed line a->b: the sign
        // of q - p projected on that direction. Every question the trace asks
        // about progress along the segment is a comparison of two such
        // projections, and their difference is a single dot product -- so the
        // signed progress of each point never has to be formed. Comparing
        // against the target b is how the walk stops at the far end.
        const auto alongOrder = [&](const auto& p, const auto& q) {
            return dotSign(p, q, a, b);
        };
        const auto vertexOrder = [&](VertexIndex u, VertexIndex w) {
            return alongOrder(vertices_[u], vertices_[w]);
        };
        // From a vertex w on the segment, emit its fan and follow the segment: if
        // it enters a triangle interior return that triangle (resume the trace);
        // if it continues collinearly along an edge, hop to the next on-segment
        // vertex and repeat; return NO_TRI when the segment leaves/ends here.
        const auto advanceFromVertex = [&](VertexIndex w, TriIndex anchor) -> TriIndex {
            std::size_t hops = 0, lim = vertices_.size() + 1;
            while (!stop && ++hops < lim) {
                emitFan(w, anchor);
                if (stop) return NO_TRI;
                if constexpr (!unboundedFront) {
                    if (alongOrder(b, vertices_[w]) >= 0) return NO_TRI;  // reached the target end
                }
                const TriIndex interior = forwardAround(w, anchor);
                if (interior != NO_TRI) return interior;
                // Look for a forward neighbor x with edge {w,x} collinear with a->b,
                // not past the target.
                VertexIndex nextV = NO_TRI;
                TriIndex nextAnchor = NO_TRI;
                TriIndex cur = anchor, from = NO_TRI;
                std::size_t g = 0, glim = triangles_.size() + 1;
                do {
                    if (!isGhost(cur)) {
                        for (VertexIndex y : triangles_[cur].v) {
                            if (y != w &&
                                detail::orientationSignOf(fa, fb, filteredVertex(y)).value() == 0 &&
                                vertexOrder(w, y) > 0 &&
                                (unboundedFront || alongOrder(b, vertices_[y]) <= 0)) {
                                nextV = y;
                                nextAnchor = cur;
                                break;
                            }
                        }
                    }
                    if (nextV != NO_TRI) break;
                    const TriIndex next = rotateAround(cur, w, from);
                    from = cur;
                    cur = next;
                } while (cur != anchor && cur != NO_TRI && ++g < glim);
                if (nextV == NO_TRI) return NO_TRI;  // segment leaves or ends at w
                w = nextV;
                anchor = nextAnchor;
            }
            return NO_TRI;
        };
 
        // Whether the directed line a->b crosses ghost g's hull edge from outside
        // into the hull (rather than exiting it). The interior lies on the side of
        // the hull edge holding the inside triangle's apex; the line enters iff its
        // forward direction a->b points toward that same side. (Line mode only.)
        [[maybe_unused]] const auto lineEnters = [&](TriIndex g) -> bool {
            const VertexIndex va = triangles_[g].v[0];
            const VertexIndex vb = triangles_[g].v[1];
            const auto& iv = triangles_[triangles_[g].nbr[2]].v;  // inside (real) triangle
            VertexIndex apex = iv[0];
            for (int i = 1; i < 3; ++i) {
                if (iv[i] != va && iv[i] != vb) apex = iv[i];
            }
            const auto insideSide =
                orientationDeterminant(vertices_[va], vertices_[vb], vertices_[apex]);
            // (vb - va) x (b - a): the side of the hull edge that a->b points toward.
            const auto dirCross = orientationDeterminant(vertices_[va], vertices_[vb], b) -
                                  orientationDeterminant(vertices_[va], vertices_[vb], a);
            return dirCross != 0 && (dirCross > 0) == (insideSide > 0);
        };
 
        // The ghost whose hull edge the directed query crosses to ENTER the hull,
        // or NO_TRI if it never reaches the hull at all. Used by every query whose
        // start lies outside: the direction a->b is what picks the entry out of the
        // (at most two) hull edges the supporting line meets, so a segment starting
        // outside is entered exactly like a line.
        //
        // Ghost triangles share a single placeholder apex (they are topological,
        // not Euclidean), so navigation uses only the real hull endpoints: each step
        // moves toward whichever endpoint of the current hull edge is closer to the
        // supporting line a->b (smaller |orientation determinant|), so the descent
        // walks O(arc) ghosts toward a crossing rather than scanning all of them.
        //
        // That descent lands on *a* crossing, which need not be the entry: it may be
        // where the query leaves the hull, and a crossing where the line merely
        // touches a vertex the two hull edges share is no entry at all. So the ring
        // is then followed until a hull edge the line genuinely enters through.
        // A line running *along* a hull edge enters through no edge, yet still meets
        // the triangles lining it, so that collinear edge is kept as a fallback.
        // @p g0 seeds the descent: the ghost the query's own source located, when it
        // has one, so the walk starts near its crossing and stays O(arc).
        [[maybe_unused]] const auto lineEntry = [&](TriIndex g0) -> TriIndex {
            TriIndex g = g0, crossing = NO_TRI;
            std::size_t guard = 0, lim = triangles_.size() + 1;
            while (g != NO_TRI && isGhost(g) && ++guard < lim) {
                const VertexIndex va = triangles_[g].v[0];
                const VertexIndex vb = triangles_[g].v[1];
                const auto da = orientationDeterminant(a, b, vertices_[va]);
                const auto db = orientationDeterminant(a, b, vertices_[vb]);
                if (da == 0 || db == 0 || (da > 0) != (db > 0)) { crossing = g; break; }
                const auto absA = da < 0 ? -da : da;
                const auto absB = db < 0 ? -db : db;
                g = triangles_[g].nbr[absA < absB ? 1 : 0];
            }
            if (crossing == NO_TRI) return NO_TRI;  // the supporting line misses the hull
            TriIndex h = crossing, from = NO_TRI, collinear = NO_TRI;
            guard = 0;
            while (h != NO_TRI && isGhost(h) && ++guard < lim) {
                const VertexIndex va = triangles_[h].v[0];
                const VertexIndex vb = triangles_[h].v[1];
                const auto da = orientationDeterminant(a, b, vertices_[va]);
                const auto db = orientationDeterminant(a, b, vertices_[vb]);
                const bool straddle = da == 0 || db == 0 || (da > 0) != (db > 0);
                if (da == 0 && db == 0) {
                    collinear = h;  // the query runs along this hull edge
                } else if (straddle && lineEnters(h)) {
                    return h;
                }
                const TriIndex next = triangles_[h].nbr[0] == from ? triangles_[h].nbr[1]
                                                                : triangles_[h].nbr[0];
                from = h;
                h = next;
                if (h == crossing) break;  // came full circle
            }
            return collinear;
        };
 
        // Begin the trace at the hull edge of ghost `g`, which the query crosses to
        // enter the triangulated region. Crossing through an *endpoint* of that edge
        // — a mesh vertex — is not the same as crossing through its relative
        // interior: there the query meets that vertex's whole fan, and may only
        // graze the hull there or run on along its boundary. So it enters exactly
        // like any other on-vertex advance (fan, then follow), which is also what
        // makes the entry independent of which of the vertex's two ghosts was found.
        // Otherwise it simply enters the triangle just inside the edge.
        [[maybe_unused]] const auto enterThroughGhost = [&](TriIndex g) -> TriIndex {
            const VertexIndex va = triangles_[g].v[0];
            const VertexIndex vb = triangles_[g].v[1];
            const bool onA = detail::orientationSignOf(fa, fb, filteredVertex(va)).value() == 0;
            const bool onB = detail::orientationSignOf(fa, fb, filteredVertex(vb)).value() == 0;
            if (onA || onB) {
                VertexIndex w = onA ? va : vb;
                if (onA && onB) {
                    // Along the hull edge: start at the endpoint met first, never at
                    // one behind a finite source.
                    const VertexIndex first = vertexOrder(va, vb) >= 0 ? va : vb;
                    const VertexIndex second = (first == va) ? vb : va;
                    w = (unboundedBack || alongOrder(a, vertices_[first]) >= 0) ? first : second;
                }
                prev = NO_TRI;
                return advanceFromVertex(w, g);
            }
            prev = g;  // entry boundary edge; don't step back out through it
            const TriIndex inside = triangles_[g].nbr[2];
            return (inside != NO_TRI && !isGhost(inside)) ? inside : NO_TRI;
        };
 
        // Pick the triangle to begin tracing the component reached at `start`,
        // which contains point `p`; visit p's contact triangles along the way.
        [[maybe_unused]] const auto enterAt = [&](const auto& p, TriIndex start) -> TriIndex {
            const auto& v = triangles_[start].v;
            const auto fp = filteredPoint(p);
            int zeros = 0, z0 = -1, z1 = -1;
            for (int k = 0; k < 3; ++k) {
                if (detail::orientationSignOf(filteredVertex(v[(k + 1) % 3]),
                                              filteredVertex(v[(k + 2) % 3]), fp).value() == 0) {
                    ++zeros;
                    if (z0 < 0) z0 = k; else z1 = k;
                }
            }
            if (zeros >= 2) {  // p is a vertex: visit its fan, then continue
                return advanceFromVertex(v[3 - z0 - z1], start);
            }
            if (zeros == 1) {  // p on an edge: visit both sides, continue into the forward one
                const TriIndex other = triangles_[start].nbr[z0];
                emit(start);
                emit(other);
                if (rayEnters(start, p) == 1) return start;
                if (other != NO_TRI && !isGhost(other) && rayEnters(other, p) == 1) return other;
                const VertexIndex e1 = v[(z0 + 1) % 3];
                const VertexIndex e2 = v[(z0 + 2) % 3];
                // Continue only if the segment runs ALONG the shared edge (hop to
                // its forward endpoint). Otherwise the segment crosses the edge
                // into the other side already emitted — if that side is outside
                // the hull, nothing more lies on this side.
                if (detail::orientationSignOf(filteredVertex(e1), filteredVertex(e2), fb)
                        .value() != 0) {
                    return NO_TRI;
                }
                const VertexIndex fwd = vertexOrder(e2, e1) > 0 ? e1 : e2;
                if constexpr (!unboundedFront) {
                    // The target lies strictly inside this edge: s stops before the
                    // forward endpoint, so its fan must not be visited. The two
                    // triangles sharing the edge are all s meets, and both are out.
                    if (alongOrder(b, vertices_[fwd]) > 0) {
                        return NO_TRI;
                    }
                }
                return advanceFromVertex(fwd, start);
            }
            emit(start);  // p strictly inside
            return start;
        };
 
        // Emit exactly the triangle(s) whose closure contains the target b — the
        // one triangle it is interior to, both triangles sharing the edge it lies
        // on, or the whole fan if b is a mesh vertex — then the walk ends. (Unlike
        // enterAt, this never continues, so it cannot over-emit a vertex fan when
        // b merely lies on an edge.)
        [[maybe_unused]] const auto emitTargetContacts = [&](TriIndex t) {
            const auto& v = triangles_[t].v;
            int zeros = 0, z0 = -1, z1 = -1;
            for (int k = 0; k < 3; ++k) {
                if (detail::orientationSignOf(filteredVertex(v[(k + 1) % 3]),
                                              filteredVertex(v[(k + 2) % 3]), fb).value() == 0) {
                    ++zeros;
                    if (z0 < 0) z0 = k; else z1 = k;
                }
            }
            if (zeros >= 2) {  // b is a mesh vertex: every incident triangle touches it
                emitFan(v[3 - z0 - z1], t);
                return;
            }
            emit(t);
            if (zeros == 1) {  // b on an edge: the triangle on the far side also touches it
                emit(triangles_[t].nbr[z0]);
            }
        };
 
        TriIndex t;
        if constexpr (unboundedBack) {
            // A line has no finite source: enter the hull at the ghost where the
            // directed line a->b crosses in, then trace forward like a segment. With
            // no source to locate, the descent starts at an arbitrary ghost.
            const TriIndex g = lineEntry(firstGhost_);
            if (g == NO_TRI) {
                return false;  // the line misses the triangulated region
            }
            t = enterThroughGhost(g);
        } else {
            t = locateIndex(a);
            if (t == NO_TRI) {
                return false;  // empty triangulation
            }
            if (isGhost(t)) {
                // Source outside the hull: the hull boundary is crossed twice, so it
                // is the direction a->b that says which crossing is the way in — for
                // a segment and a ray no less than for a line — and the query must
                // then actually reach that edge to meet the region at all. The ghost
                // the source landed in seeds the descent, keeping it O(arc).
                const TriIndex g = lineEntry(t);
                if (g == NO_TRI || !s.intersects(edgeSegment(Edge{g, 2}))) {
                    return false;  // s never meets the triangulated region
                }
                t = enterThroughGhost(g);
            } else {
                t = enterAt(a, t);  // source inside the hull: start there
            }
        }
 
        const std::size_t cap = triangles_.size() * 4 + 16;
        std::size_t guard = 0;
        while (!stop && t != NO_TRI) {
            if (++guard > cap) break;  // safety net
            emit(t);
            if (stop) break;
            if constexpr (!unboundedFront) {
                if (pointInClosure(b, t)) {  // target reached: visit its contact triangles
                    emitTargetContacts(t);
                    break;
                }
            }
            const auto& v = triangles_[t].v;
            // Forward exit edge: the supporting line ab crosses it, it isn't where
            // we came from, and we leave forward through it — for a segment that is
            // "the segment ab straddles the edge"; for an unbounded-forward query
            // (ray / line) it is "a->b points out of the triangle across the edge"
            // (no finite endpoint to bound the crossing).
            int exitK = -1;
            for (int k = 0; k < 3; ++k) {
                if (triangles_[t].nbr[k] == prev) continue;
                const auto& u = vertices_[v[(k + 1) % 3]];
                const auto& w = vertices_[v[(k + 2) % 3]];
                const auto fu = filteredVertex(v[(k + 1) % 3]);
                const auto fw = filteredVertex(v[(k + 2) % 3]);
                const auto du = detail::orientationSignOf(fa, fb, fu).value();
                const auto dw = detail::orientationSignOf(fa, fb, fw).value();
                const bool straddleLine = (du > 0 && dw < 0) || (du < 0 && dw > 0);
                if (!straddleLine) continue;
                bool forward;
                if constexpr (unboundedFront) {
                    // a->b leaves t outward across (u,w): its direction points to the
                    // side of (u,w) away from the apex v[k] (the triangle interior).
                    const auto apexSide =
                        detail::orientationSignOf(fu, fw, filteredVertex(v[k])).value();
                    const auto dirCross = orientationDeterminant(u, w, b) -
                                          orientationDeterminant(u, w, a);
                    forward = dirCross != 0 && (dirCross > 0) != (apexSide > 0);
                } else {
                    const auto ea = detail::orientationSignOf(fu, fw, fa).value();
                    const auto eb = detail::orientationSignOf(fu, fw, fb).value();
                    forward = (ea > 0 && eb < 0) || (ea < 0 && eb > 0);
                }
                if (forward) { exitK = k; break; }
            }
            if (exitK != -1) {
                prev = t;
                t = triangles_[t].nbr[exitK];
                if (isGhost(t)) {
                    // Left the triangulated region; by convexity s cannot re-enter.
                    // Under AllTriangles report the ghost we stepped into, so the
                    // caller can tell this exit from a mere touch of the boundary.
                    if constexpr (AllTriangles) {
                        (void)f(t);
                    }
                    break;
                }
                continue;
            }
            // Degenerate: s leaves t through a vertex or runs along an edge.
            int onCount = 0, on0 = -1, on1 = -1;
            for (int m = 0; m < 3; ++m) {
                if (detail::orientationSignOf(fa, fb, filteredVertex(v[m])).value() == 0) {
                    ++onCount;
                    if (on0 < 0) on0 = m; else on1 = m;
                }
            }
            if (onCount >= 1) {
                // s leaves t through the on-line vertex farthest along a->b; from
                // there, advance (fan + collinear hops) until it re-enters an interior.
                const int m = (onCount >= 2 && vertexOrder(v[on0], v[on1]) > 0) ? on1 : on0;
                t = advanceFromVertex(v[m], t);
                prev = NO_TRI;
                continue;
            }
            break;  // no progress (should not happen for a valid triangulation)
        }
        return stop;
    }
 
  public:
    template <detail::TriangulationQuery OS>
    [[nodiscard]] std::vector<TriangleType> trianglesIntersecting(const OS& s) const {
        std::vector<TriangleType> out;
        visitTrianglesIntersecting(s, [&](const TriangleType& t) { out.push_back(t); });
        return out;
    }

    template <detail::TriangulationQuery OS, class Fn>
    bool visitTrianglesInteriorIntersecting(const OS& s, Fn f) const {
        return visitTrianglesIntersecting(s, [&](TriId id) -> bool {
            const TriangleType t = triangleValue(indexOf(id));
            if (!t.interiorsIntersect(s)) {
                return false;  // only boundary contact: skip, but keep walking
            }
            if constexpr (std::is_invocable_v<Fn&, const TriangleType&>) {
                return detail::invokeVisitor(f, t);
            } else {
                return detail::invokeVisitor(f, id);
            }
        });
    }

    template <detail::TriangulationQuery OS>
    [[nodiscard]] std::vector<TriangleType> trianglesInteriorIntersecting(const OS& s) const {
        std::vector<TriangleType> out;
        visitTrianglesInteriorIntersecting(s, [&](const TriangleType& t) { out.push_back(t); });
        return out;
    }

    template <detail::TriangulationQuery OS, class Fn>
    bool visitEdgesIntersecting(const OS& s, Fn f) const {
        std::unordered_set<SegmentType> seen;
        bool stop = false;
        visitTrianglesIntersecting(s, [&](const TriangleType& t) -> bool {
            for (const auto& edge : t.edges()) {
                if (!s.intersects(edge)) {
                    continue;
                }
                // Canonical edge (the constructor re-sorts the endpoints), used
                // for reporting and for deduplicating across the two incident
                // triangles. Only when edges carry a label do we look it up in
                // segToEdge_ to recover the stored label; otherwise the segment
                // built from the endpoints already is the value to report.
                SegmentType seg(edge[0], edge[1]);
                if constexpr (detail::has_label_v<SegmentLabel>) {
                    auto se = segmentMap().find(seg);
                    if (se == segmentMap().end()) {
                        continue;  // not a stored edge (should not happen)
                    }
                    seg = edgeSegment(se->second);
                }
                if (seen.insert(seg).second && detail::invokeVisitor(f, seg)) {
                    stop = true;
                    return true;
                }
            }
            return false;
        });
        return stop;
    }

    template <detail::TriangulationQuery OS, class Fn>
    bool visitEdgesInteriorIntersecting(const OS& s, Fn f) const {
        return visitEdgesIntersecting(s, [&](const SegmentType& e) -> bool {
            if (s.interiorsIntersect(e)) {
                return detail::invokeVisitor(f, e);
            }
            return false;  // endpoint-only contact: skip, but keep walking
        });
    }

    template <detail::TriangulationQuery OS>
    [[nodiscard]] std::vector<SegmentType> edgesIntersecting(const OS& s) const {
        std::vector<SegmentType> out;
        visitEdgesIntersecting(s, [&](const SegmentType& e) { out.push_back(e); });
        return out;
    }

    template <detail::TriangulationQuery OS>
    [[nodiscard]] std::vector<SegmentType> edgesInteriorIntersecting(const OS& s) const {
        std::vector<SegmentType> out;
        visitEdgesInteriorIntersecting(s, [&](const SegmentType& e) { out.push_back(e); });
        return out;
    }
 
    // ---- traversal over a region (a connected convex query shape) --------

    template <detail::TriangulationRegionQuery Q, class Fn>
    bool visitTrianglesIntersecting(const Q& shape, Fn f) const {
        if (firstGhost_ == 0) return false;  // no real triangles
        const std::vector<TriIndex> seeds = seedTrianglesIntersecting(shape);
        if (seeds.empty()) return false;
        return floodTriangleIdsIntersecting<false>(
            shape, seeds, [&](TriIndex t) { return reportTriangle(f, t); });
    }
 
    // The materialized form — trianglesIntersecting(shape) — and the
    // interior-intersecting and edge variants are the shared wrappers above,
    // constrained on detail::TriangulationQuery, so they accept a region shape
    // here just as they accept a segment.
 
    // ---- traversal along a chain (a polyline or a monotone chain) --------

    template <detail::ChainTraversal C, class Fn>
    bool visitTrianglesIntersecting(const C& c, Fn f) const {
        if (firstGhost_ == 0) return false;  // no real triangles
        if (c.size() < 2) {
            return c.empty() ? false : visitTrianglesIntersecting(c[0], f);
        }
        // Triangles already reported, so an edge re-entering a triangle an
        // earlier edge met does not report it again (the per-edge walks each
        // deduplicate only within themselves).
        std::unordered_set<TriId> seen;
        bool stop = false;
        for (const auto& e : c.orientedEdgesView()) {
            visitTrianglesIntersecting(e, [&](TriId t) -> bool {
                if (!seen.insert(t).second) {
                    return false;  // already reported by an earlier edge
                }
                stop = reportTriangle(f, indexOf(t));
                return stop;
            });
            if (stop) {
                break;
            }
        }
        return stop;
    }
 
    // As for a region shape, the materialized form — trianglesIntersecting(c) —
    // and the interior-intersecting and edge variants are the shared wrappers
    // above, constrained on detail::TriangulationQuery.
 
    // ---- point location --------------------------------------------------

    template <PointConcept QueryPoint>
    [[nodiscard]] std::optional<TriangleType> locate(const QueryPoint& p) const {
        const TriId id = locateId(p);
        if (!id.valid()) {
            return std::nullopt;  // outside the region, or in a hull-fill triangle
        }
        return getShape(id);
    }
 
    // ---- predicates against the domain -----------------------------------

    template <PointConcept QueryPoint>
    [[nodiscard]] bool contains(const QueryPoint& shape) const {
        // The domain is closed, so containing a point and meeting it are the same
        // question. Note that locate() would not answer it: for a point on the
        // domain boundary the visibility walk may land on the out-of-domain side,
        // and the point is contained all the same.
        return intersects(shape);
    }

    template <detail::SegmentOrOriented S>
    [[nodiscard]] bool contains(const S& shape) const {
        if (!contains(shape[0]) || !contains(shape[1])) {
            return false;
        }
        return segmentInteriorContained(shape);
    }
 
  private:
    // Segment containment once its endpoints are known to lie in the closed
    // domain. Kept separate because mesh vertices satisfy that precondition and
    // the convex-cover visibility hot path must not locate them repeatedly.
    template <detail::SegmentOrOriented S>
    [[nodiscard]] bool segmentInteriorContained(const S& shape) const {
        // Both endpoints are in the domain, so the segment can only escape between
        // them, and the walk sees every triangle it crosses on the way. It escapes
        // iff it (a) enters the interior of an out-of-domain hull-fill triangle,
        // (b) runs along an edge with no in-domain side — a diagonal of the fill
        // region, which no interior is entered by — or (c) leaves the triangulated
        // region altogether, which the walk reports as a step into a ghost. Merely
        // touching a fill triangle (along a boundary edge, or at a boundary vertex
        // whose fan the walk sweeps) is not an escape: the boundary belongs to the
        // domain.
        return !visitTriangleIdsIntersecting<true>(shape, [&](TriIndex t) {
            if (isGhost(t)) {
                return true;  // (c)
            }
            if (inDomain(t)) {
                return false;
            }
            if (triangleValue(t).interiorsIntersect(shape)) {
                return true;  // (a)
            }
            for (std::int8_t side = 0; side < 3; ++side) {
                const Edge e{t, side};
                if (!edgeInDomain(e) && shape.interiorsIntersect(edgeSegment(e))) {
                    return true;  // (b)
                }
            }
            return false;
        });
    }
 
  public:

    template <detail::ChainTraversal C>
    [[nodiscard]] bool contains(const C& shape) const {
        if (shape.empty()) {
            return true;
        }
        if (shape.size() < 2) {
            return contains(shape[0]);
        }
        return containsBoundary(shape.orientedEdgesView(), shape);
    }

    template <detail::PolygonalRegion Q>
    [[nodiscard]] bool contains(const Q& shape) const {
        // A closed boundary inside the domain encloses nothing outside it — see
        // detail::PolygonalRegion — so the edges decide, as long as the domain is
        // simply connected.
        if (!containsBoundary(shape.edges(), shape)) {
            return false;
        }
        // A region domain is not simply connected, and the one thing its holes
        // let a contained boundary enclose is a hole. The shape's boundary is
        // inside the domain, hence clear of every hole interior, so each hole is
        // wholly inside the shape or wholly outside it and one witness triangle
        // per hole settles it.
        for (const TriangleType& witness : holeWitnesses_) {
            if (shape.contains(witness)) {
                return false;
            }
        }
        return true;
    }

    template <PolygonWithHolesConcept Q>
    [[nodiscard]] bool contains(const Q& shape) const {
        if (shape.outer().empty()) {
            return true;
        }
        bool traced = false;
        for (const auto& edge : shape.edges()) {
            if (edge[0] == edge[1]) {
                continue;
            }
            traced = true;
            if (!contains(edge)) {
                return false;
            }
        }
        if (!traced && !contains(shape.outer()[0])) {
            return false;
        }
        for (const TriangleType& witness : holeWitnesses_) {
            if (shape.contains(witness)) {
                return false;
            }
        }
        return true;
    }

    template <PolygonSetConcept Q>
    [[nodiscard]] bool contains(const Q& shape) const {
        for (const auto& component : shape.components()) {
            if (!contains(component)) {
                return false;
            }
        }
        return true;
    }

    template <HalfplaneIntersectionConcept Q>
    [[nodiscard]] bool contains(const Q& shape) const {
        if (shape.empty()) {
            return true;
        }
        if (!shape.isBounded()) {
            return false;
        }
        return contains(shape.template asConvex<division_result_t<typename Q::NumberType>>());
    }

    template <DiskConcept D>
    [[nodiscard]] bool contains(const D& shape) const {
        // A disk has no edges to trace, so it is grown by the region flood fill
        // instead. The domain is closed, so the disk lies in it iff its *interior*
        // does, and the open disk escapes iff it reaches into the interior of an
        // out-of-domain hull-fill triangle or crosses out through a boundary edge
        // of the triangulated region (one whose far side is a ghost). Neither
        // fires on a disk that only touches the domain boundary — an open disk on
        // one side of a line never meets that line — so containment also demands
        // that the open disk actually reach inside: a disk resting against the
        // boundary from the outside, meeting the domain in a single point, escapes
        // through no triangle and no edge, and is still not contained.
        const std::vector<TriIndex> seeds = seedTrianglesIntersecting(shape);
        if (seeds.empty()) {
            return false;  // the disk misses the triangulated region entirely
        }
        bool reachesInside = false;
        const bool escapes = floodTriangleIdsIntersecting<true>(shape, seeds, [&](TriIndex t) {
            const TriangleType triangle = triangleValue(t);
            if (triangle.interiorsIntersect(shape)) {
                if (!inDomain(t)) {
                    return true;  // reaches into the hull-fill region
                }
                reachesInside = true;
            }
            for (std::int8_t side = 0; side < 3; ++side) {
                if (isGhost(triangles_[t].nbr[side]) &&
                    shape.interiorsIntersect(edgeSegment(Edge{t, side}))) {
                    return true;  // crosses out through the boundary of the region
                }
            }
            return false;
        });
        // reachesInside is complete only when the fill ran to the end, which is
        // exactly when nothing escaped.
        return !escapes && reachesInside;
    }

    template <class U>
        requires detail::LineOrOriented<U> || RayConcept<U> || HalfplaneConcept<U>
    [[nodiscard]] bool contains(const U&) const {
        return false;
    }

    template <EmptyShapeConcept E>
    [[nodiscard]] bool contains(const E&) const {
        return true;
    }

    [[nodiscard]] bool contains(const Shape<PointType>& shape) const {
        return std::visit([this](const auto& s) { return this->contains(s); }, shape.variant());
    }

    template <detail::TriangulationQuery Q>
    [[nodiscard]] bool intersects(const Q& shape) const {
        return visitTrianglesIntersecting(shape, [](const TriangleType&) { return true; });
    }

    template <PolygonConcept Q>
    [[nodiscard]] bool intersects(const Q& shape) const {
        // A polygon is not a traversal query (it need not be convex), so it is
        // taken edge by edge. If its boundary misses the domain entirely then
        // each connected piece of the domain is either inside the polygon or
        // disjoint from it, and one vertex of each settles which.
        for (const auto& e : shape.edges()) {
            if (e[0] != e[1] && intersects(e)) {
                return true;
            }
        }
        return anyDomainComponent(
            [&](TriIndex t) { return shape.contains(vertices_[triangles_[t].v[0]]); });
    }

    template <PolygonWithHolesConcept Q>
    [[nodiscard]] bool intersects(const Q& shape) const {
        if (shape.outer().empty()) {
            return false;
        }
        for (const auto& edge : shape.edges()) {
            if (edge[0] != edge[1] && intersects(edge)) {
                return true;
            }
        }
        return anyDomainComponent(
            [&](TriIndex t) { return shape.contains(vertices_[triangles_[t].v[0]]); });
    }

    template <PolygonSetConcept Q>
    [[nodiscard]] bool intersects(const Q& shape) const {
        for (const auto& component : shape.components()) {
            if (intersects(component)) {
                return true;
            }
        }
        return false;
    }

    template <HalfplaneIntersectionConcept Q>
    [[nodiscard]] bool intersects(const Q& shape) const {
        for (TriIndex t = 0; t < static_cast<TriIndex>(triangles_.size()); ++t) {
            if (inDomain(t) && shape.intersects(triangleValue(t))) {
                return true;
            }
        }
        return false;
    }

    template <EmptyShapeConcept E>
    [[nodiscard]] bool intersects(const E&) const {
        return false;
    }

    [[nodiscard]] bool intersects(const Shape<PointType>& shape) const {
        return std::visit([this](const auto& s) { return this->intersects(s); }, shape.variant());
    }

    template <detail::TriangulationQuery Q>
    [[nodiscard]] bool interiorContains(const Q& shape) const {
        return contains(shape) && !meetsDomainBoundary(shape);
    }

    template <PolygonConcept Q>
    [[nodiscard]] bool interiorContains(const Q& shape) const {
        // As for contains(), the edges decide: a polygon lying in the domain can
        // only touch ∂D along its own boundary. An interior point of the polygon
        // touching ∂D would have points outside the domain arbitrarily close to it,
        // yet a whole neighbourhood of it lies in the polygon, hence in the domain.
        return containsBoundary<true>(shape.edges(), shape);
    }

    template <PolygonWithHolesConcept Q>
    [[nodiscard]] bool interiorContains(const Q& shape) const {
        if (shape.outer().empty()) {
            return true;
        }
        bool traced = false;
        for (const auto& edge : shape.edges()) {
            if (edge[0] == edge[1]) {
                continue;
            }
            traced = true;
            if (!interiorContains(edge)) {
                return false;
            }
        }
        if (!traced && !interiorContains(shape.outer()[0])) {
            return false;
        }
        for (const TriangleType& witness : holeWitnesses_) {
            if (shape.contains(witness)) {
                return false;
            }
        }
        return true;
    }

    template <PolygonSetConcept Q>
    [[nodiscard]] bool interiorContains(const Q& shape) const {
        for (const auto& component : shape.components()) {
            if (!interiorContains(component)) {
                return false;
            }
        }
        return true;
    }

    template <HalfplaneIntersectionConcept Q>
    [[nodiscard]] bool interiorContains(const Q& shape) const {
        if (shape.empty()) {
            return true;
        }
        if (!shape.isBounded()) {
            return false;
        }
        if (!contains(shape)) {
            return false;
        }
        return interiorContains(shape.template asConvex<division_result_t<typename Q::NumberType>>());
    }

    template <EmptyShapeConcept E>
    [[nodiscard]] bool interiorContains(const E&) const {
        return true;
    }

    [[nodiscard]] bool interiorContains(const Shape<PointType>& shape) const {
        return std::visit([this](const auto& s) { return this->interiorContains(s); },
                          shape.variant());
    }

    template <PointConcept QueryPoint>
    [[nodiscard]] bool interiorsIntersect(const QueryPoint& shape) const {
        return interiorContains(shape);  // a point's interior is the point itself
    }

    template <detail::TriangulationQuery Q>
        requires(!PointConcept<Q>)
    [[nodiscard]] bool interiorsIntersect(const Q& shape) const {
        // The interior of the domain is the open triangles plus the relative
        // interiors of the interior edges (those with an in-domain triangle on
        // both sides) plus the interior vertices — and the vertices are redundant
        // here: a shape of positive dimension whose interior reaches one of them
        // also reaches, arbitrarily close by, an open triangle or an interior
        // edge. (For a point it would not, which is why a point is its own case.)
        return visitTrianglesIntersecting(shape, [&](const TriangleType& t) {
            if (t.interiorsIntersect(shape)) {
                return true;
            }
            return anyEdgeOf(t, [&](const SegmentType& e, Edge handle) {
                return interiorEdge(handle) && e.interiorsIntersect(shape);
            });
        });
    }

    template <PolygonConcept Q>
    [[nodiscard]] bool interiorsIntersect(const Q& shape) const {
        // Edge by edge, as for the other polygon overloads. If no edge reaches the
        // domain's interior, then the polygon's boundary misses it entirely, and
        // the interiors meet iff some connected piece of the domain lies inside
        // the polygon — one triangle of each piece settles it.
        for (const auto& e : shape.edges()) {
            if (e[0] != e[1] && interiorsIntersect(e)) {
                return true;
            }
        }
        return anyDomainComponent(
            [&](TriIndex t) { return shape.interiorsIntersect(triangleValue(t)); });
    }

    template <PolygonWithHolesConcept Q>
    [[nodiscard]] bool interiorsIntersect(const Q& shape) const {
        if (shape.outer().empty()) {
            return false;
        }
        for (const auto& edge : shape.edges()) {
            if (edge[0] != edge[1] && interiorsIntersect(edge)) {
                return true;
            }
        }
        return anyDomainComponent(
            [&](TriIndex t) { return shape.interiorsIntersect(triangleValue(t)); });
    }

    template <PolygonSetConcept Q>
    [[nodiscard]] bool interiorsIntersect(const Q& shape) const {
        for (const auto& component : shape.components()) {
            if (interiorsIntersect(component)) {
                return true;
            }
        }
        return false;
    }

    template <HalfplaneIntersectionConcept Q>
    [[nodiscard]] bool interiorsIntersect(const Q& shape) const {
        for (TriIndex t = 0; t < static_cast<TriIndex>(triangles_.size()); ++t) {
            if (inDomain(t) && triangleValue(t).interiorsIntersect(shape)) {
                return true;
            }
        }
        return false;
    }

    template <EmptyShapeConcept E>
    [[nodiscard]] bool interiorsIntersect(const E&) const {
        return false;
    }

    [[nodiscard]] bool interiorsIntersect(const Shape<PointType>& shape) const {
        return std::visit([this](const auto& s) { return this->interiorsIntersect(s); },
                          shape.variant());
    }
 
    // ---- constrained edges ----------------------------------------------

    [[nodiscard]] bool isConstrained(const SegmentType& s) const {
        auto se = segmentMap().find(s);
        return se != segmentMap().end() && bit(triangles_[se->second.tri].constrainedMask, se->second.side);
    }

    void setConstrained(const SegmentType& s, bool value = true) {
        auto se = segmentMap().find(s);
        if (se == segmentMap().end()) {
            return;
        }
        const Edge e = se->second;
        setBit(triangles_[e.tri].constrainedMask, e.side, value);
        const Edge m = mirror(e);
        if (m.tri != NO_TRI) {
            setBit(triangles_[m.tri].constrainedMask, m.side, value);
        }
    }
 
    // ---- labels ----------------------------------------------------------

    template <class L = TriangleLabel>
        requires(detail::has_label_v<L>)
    [[nodiscard]] L& label(const TriangleType& t) {
        const TriIndex id = idOf(t);
        assert(inDomain(id) && "label(): triangle is not part of the triangulation");
        return triangles_[id].triLabel;
    }

    template <class L = TriangleLabel>
        requires(detail::has_label_v<L>)
    [[nodiscard]] const L& label(const TriangleType& t) const {
        const TriIndex id = idOf(t);
        assert(inDomain(id) && "label(): triangle is not part of the triangulation");
        return triangles_[id].triLabel;
    }

    template <class L = SegmentLabel>
        requires(detail::has_label_v<L>)
    [[nodiscard]] L& label(const SegmentType& s) {
        auto se = segmentMap().find(s);
        assert(se != segmentMap().end() && "label(): segment is not an edge of the triangulation");
        return se->second.segLabel;
    }

    template <class L = SegmentLabel>
        requires(detail::has_label_v<L>)
    [[nodiscard]] const L& label(const SegmentType& s) const {
        auto se = segmentMap().find(s);
        assert(se != segmentMap().end() && "label(): segment is not an edge of the triangulation");
        return se->second.segLabel;
    }
 
    // ---- mutation --------------------------------------------------------

    [[nodiscard]] bool flippable(const SegmentType& s) const {
        auto se = segmentMap().find(s);
        return se != segmentMap().end() && flippableEdge(se->second);
    }

    std::optional<SegmentType> flip(const SegmentType& s) {
        auto se = segmentMap().find(s);
        if (se == segmentMap().end()) {
            return std::nullopt;
        }
        const Edge e = se->second;
        if (!flippableEdge(e)) {
            return std::nullopt;
        }
        const TriIndex t = e.tri;
        const TriIndex t2 = mirror(e).tri;
        flipEdge(e);
        // The shared edge is gone; re-register the six sides of the two
        // rewritten triangles (the four surrounding edges get fresh handles and
        // the new diagonal is added).
        segmentMap().erase(s);
        registerSides(t);
        registerSides(t2);
        ++revision_;
        return edgeSegment(Edge{t, 1});  // new diagonal (side opposite a in t)
    }

    template <class EdgeRange>
        requires SegmentConcept<typename EdgeRange::value_type>
    [[nodiscard]] bool flippable(const EdgeRange& edges) const {
        std::unordered_set<TriIndex> claimed;  // triangles already covered by some quad
        for (const auto& s : edges) {
            const auto se = segmentMap().find(SegmentType(s[0], s[1]));
            if (se == segmentMap().end() || !flippableEdge(se->second)) {
                return false;
            }
            const Edge e = se->second;
            if (!claimed.insert(e.tri).second || !claimed.insert(mirror(e).tri).second) {
                return false;  // a quad shares a triangle with an earlier one
            }
        }
        return true;
    }

    template <class EdgeRange>
        requires SegmentConcept<typename EdgeRange::value_type>
    std::optional<std::vector<SegmentType>> flip(const EdgeRange& edges) {
        if (!flippable(edges)) {
            return std::nullopt;
        }
        std::vector<SegmentType> diagonals;
        for (const auto& s : edges) {
            diagonals.push_back(*flip(SegmentType(s[0], s[1])));  // disjoint: each succeeds
        }
        return diagonals;
    }
 
    // ---- point insertion -------------------------------------------------

    bool insert(const PointType& p) {
        const bool inserted = insertVertexImpl(p).has_value();
        if (inserted) {
            ++revision_;
        }
        return inserted;
    }

    bool insertDelaunay(const PointType& p) {
        const auto inserted = insertVertexImpl(p);
        if (!inserted) {
            return false;
        }
        const auto [vp, start] = *inserted;
        // The suspect edges are the new vertex's link: the side opposite vp in
        // every real triangle of its fan.
        std::vector<SegmentType> suspect;
        visitVertexFan(start, vp, [&](TriIndex cur) {
            if (!isGhost(cur)) {
                suspect.push_back(
                    edgeSegment(Edge{cur, static_cast<std::int8_t>(localIndex(cur, vp))}));
            }
        });
        legalize(suspect);
        ++revision_;
        return true;
    }
 
    // ---- validation ------------------------------------------------------

    [[nodiscard]] bool checkInvariants() const {
        for (TriIndex t = 0; t < static_cast<TriIndex>(triangles_.size()); ++t) {
            const auto& T = triangles_[t];
            if (!isGhost(t) &&
                !(orientationSign(vertices_[T.v[0]], vertices_[T.v[1]], vertices_[T.v[2]]) > 0)) {
                return false;
            }
            for (int i = 0; i < 3; ++i) {
                const TriIndex t2 = T.nbr[i];
                if (t2 == NO_TRI) {
                    continue;
                }
                int j = -1;
                for (int s = 0; s < 3; ++s) {
                    if (triangles_[t2].nbr[s] == t) {
                        j = s;
                    }
                }
                if (j < 0) {
                    return false;  // link not mutual
                }
                const VertexIndex a = T.v[(i + 1) % 3];
                const VertexIndex b = T.v[(i + 2) % 3];
                const VertexIndex c = triangles_[t2].v[(j + 1) % 3];
                const VertexIndex d = triangles_[t2].v[(j + 2) % 3];
                if (!((a == c && b == d) || (a == d && b == c))) {
                    return false;  // links disagree on the shared edge
                }
            }
        }
        return true;
    }

    friend Canvas& operator<<(Canvas& canvas, const Triangulation& triangulation) {
        for (const auto &t : triangulation.triangles()) {
            canvas << t;
        }
        return canvas;
    }
 
  private:
    // ---- internal vocabulary (never exposed) -----------------------------
 
    using VertexIndex = std::int32_t;  // index into vertices_
    using TriIndex = std::int32_t;     // index into triangles_
    static constexpr TriIndex NO_TRI = -1;
    static constexpr VertexIndex NO_VERTEX = -1;
 
    // A directed reference to one side of one triangle: `side` selects the edge
    // opposite vertex `side`, i.e. { v[(side+1)%3], v[(side+2)%3] }. As the value
    // stored in segToEdge_ it also carries that edge's label.
    struct Edge {
        TriIndex tri = NO_TRI;
        std::int8_t side = 0;
        [[no_unique_address]] SegmentLabel segLabel{};
    };
 
    // One triangle record: three CCW vertices, three neighbors (neighbor i is
    // across the edge opposite v[i]), a constrained-edge mask, an out-of-domain
    // flag (set for hull-fill triangles outside a polygon), a transient
    // walk-visited mark (used by visitTrianglesIntersecting in lieu of a hash
    // set; always cleared again before that const method returns), and the
    // label. `walkMark` is mutable because the walk is a const query.
    struct Tri {
        std::array<VertexIndex, 3> v{NO_TRI, NO_TRI, NO_TRI};
        std::array<TriIndex, 3> nbr{NO_TRI, NO_TRI, NO_TRI};
        std::uint8_t constrainedMask = 0;
        std::uint8_t outOfDomain = 0;
        mutable std::uint8_t walkMark = 0;
        [[no_unique_address]] TriangleLabel triLabel{};
    };
    static_assert(sizeof(Tri) == 28 || detail::has_label_v<TriangleLabel>,
                  "Tri should stay 28 bytes when unlabeled");
 
    // The coordinate type the vertex approximations are kept for: the in-circle
    // test promotes furthest of the predicates run here, so it is the one whose
    // gate decides whether keeping them is worth it at all. A predicate that
    // promotes less reads its own gate and simply ignores them.
    using VertexCoordinate = detail::incircle_coordinate_t<typename PointType::NumberType>;
 
    // Rebuilds every approximation, after vertices_ is reordered or refilled.
    void syncVertexApproximations() {
        if constexpr (detail::filtersSign<VertexCoordinate>) {
            vertexApproximations_.clear();
            vertexApproximations_.reserve(vertices_.size());
            for (const auto& vertex : vertices_) {
                vertexApproximations_.push_back(detail::approximatePoint(vertex));
            }
        }
    }
 
    // Appends a vertex and its approximation together.
    VertexIndex appendVertex(const PointType& p) {
        const VertexIndex id = static_cast<VertexIndex>(vertices_.size());
        vertices_.push_back(p);
        if constexpr (detail::filtersSign<VertexCoordinate>) {
            vertexApproximations_.push_back(detail::approximatePoint(p));
        }
        return id;
    }
 
    // A vertex paired with the approximation kept for it, for the sign
    // predicates in detail:: that take filtered operands.
    [[nodiscard]] auto filteredVertex(VertexIndex v) const {
        const auto index = static_cast<std::size_t>(v);
        return detail::filtered<VertexCoordinate>(vertices_[index], vertexApproximations_, index);
    }
 
    // A point that is not a vertex — a query, or one about to be inserted —
    // converted once for however many signs read it. Its own coordinate type
    // joins the gate: a rational query point against `int` vertices promotes
    // the predicates that read it into the range where filtering pays, even
    // though the vertices alone would not.
    template <class QueryPoint>
    [[nodiscard]] static auto filteredPoint(const QueryPoint& p) {
        using QueryCoordinate =
            detail::incircle_coordinate_t<typename PointType::NumberType,
                                          typename QueryPoint::NumberType>;
        return detail::filtered<QueryCoordinate>(p);
    }
 
    std::vector<PointType> vertices_;  // ghost vertex at index 0 (when nonempty),
                                       // then the real vertices; inserts append
    // Approximations of vertices_, parallel to it, or empty where the filter
    // would not pay for itself. Every sign predicate here reads its operands
    // out of vertices_ and the same vertices come up again and again — the walk
    // revisits them, the flood-fill revisits them, each legalization revisits
    // them — so converting a coordinate to double once per vertex rather than
    // once per predicate is most of what the filter costs. Kept in step with
    // vertices_ by syncVertexApproximations and appendVertex.
    std::vector<detail::ApproximatePoint> vertexApproximations_;
    std::vector<Tri> triangles_;       // real triangles [0,firstGhost_), then ghost triangles
    // Outside edge -> internal handle, reached through @ref segmentMap. Mutable
    // and paired with `mapStale_` because a bulk build leaves it unfilled: see
    // @ref materializeSegmentMap.
    mutable std::unordered_map<SegmentType, Edge> segToEdge_;
    mutable bool mapStale_ = false;
    // One triangle inside each hole of a region domain, empty otherwise. A
    // closed boundary drawn inside such a domain can enclose a hole, which the
    // boundary alone never reveals, so the region queries consult these.
    std::vector<TriangleType> holeWitnesses_;
    // One real triangle incident to each vertex, indexed by VertexIndex — the seed
    // a fan rotation needs to reach a vertex's neighbourhood without scanning
    // every triangle. Only a *hint*: it is refreshed by registerSides and
    // verified by incidentTriangleOf, so an entry left stale by an edit that
    // does not re-register costs a fallback, never an answer.
    std::vector<TriIndex> vertexTri_;
    static constexpr VertexIndex GHOST = 0;  // index of the ghost vertex
    TriIndex firstGhost_ = 0;             // first ghost triangle index
    std::size_t domainTriangleCount_ = 0;  // in-domain real triangles (<= firstGhost_)
    mutable TriIndex hint_ = NO_TRI;      // last located triangle (walk seed)
    mutable std::mt19937 rng_;         // drives the stochastic walk in locateIndex
    // The Kirkpatrick hierarchy: a stack of triangulations of the same box, the
    // finest being the mesh (plus the ring filling the box around it) and each
    // coarser one obtained by removing an independent set of vertices and
    // retriangulating the holes. A query descends it one cell per level.
    //
    // Every cell above level 0 is a triangle of the retriangulation of one
    // removed vertex's star, so the cells below it are consecutive triangles of
    // that star's fan, and which of them holds a point already known to be in
    // the parent is decided by where the point falls in the fan — one sign per
    // fan edge crossed, and none at all where the parent has a single child.
    // That is what a cell stores: the removed vertex, and the run of fan edges
    // and cells between them. It needs no triangle of its own; only the top
    // level, whose cells a query is not yet known to be inside, keeps one.
    struct Kirkpatrick {
        // Set in runCount when the run closes the whole fan, which happens for
        // the one triangle of a retriangulation that holds the removed vertex:
        // the run then spans a full turn instead of an arc, so a query cannot
        // assume it starts inside it. On a cell with no children — where there
        // is no run to mark — the same bit says the cell is one of those
        // filling the box rather than a mesh triangle, so what it carries only
        // starts the walk instead of answering it.
        static constexpr std::uint32_t FULL_TURN = 0x80000000u;
        struct Cell {
            std::uint32_t apex = 0;      // the removed vertex the fan turns around
            std::uint32_t runBegin = 0;  // where the run starts, or a leaf's seed
            std::uint32_t runCount = 0;  // children, 0 for a leaf, plus FULL_TURN
        };
        std::vector<Cell> cells;
        // The fan runs, one per cell with children: the vertices bounding the
        // fan triangles, interleaved with the cells between them, as
        // `link, cell, link, cell, ..., link` — one more vertex than cells.
        std::vector<std::uint32_t> run;
        std::vector<std::uint32_t> roots;                    // the top level
        std::vector<std::array<VertexIndex, 3>> rootShape;   // and its triangles
        // The four box corners, addressed by vertex indices vertices_.size() + i.
        std::vector<PointType> extra;
        std::vector<detail::ApproximatePoint> extraApprox;
    };
    // Immutable once built, so copies share it rather than redrawing it.
    std::shared_ptr<const Kirkpatrick> pointLocation_;
 
    // Whether the descent's signs are exact, and its answer therefore final. An
    // integer coordinate's are — the promotion holds every product the
    // orientation forms — and so are a coordinate type carrying its own exact
    // arithmetic. Everything else, a fixed-width rational (whose products
    // overflow) as much as a floating-point coordinate, can put the descent in
    // a triangle beside the right one, and there the walk finishes the query
    // exactly as it does for a mesh the hierarchy has fallen behind.
    static constexpr bool exactDescent =
        detail::extended_integral<NumberType> || detail::arbitraryPrecision<NumberType>;
 
    // The largest star a level takes off. A Delaunay vertex has six triangles
    // around it on average and all but a handful have ten or fewer, so a bound
    // here only turns away the rare crowded vertex — which its own neighbours
    // would usually have blocked anyway. Query time measures the same over
    // 6..20 at 100,000 vertices, and from 12 up the hierarchy comes out
    // identical; ten sits in the middle of that plateau.
    static constexpr std::size_t pointLocationMaxStar = 10;
 
    // Where the hierarchy stops growing. Its top level is scanned triangle by
    // triangle, so it has to stay small; anything from 2 to 16 measures the
    // same, the levels that would replace the scan costing about what the scan
    // does, and above that the scan starts to show (64 triangles cost a tenth
    // of the query).
    static constexpr std::size_t pointLocationTopSize = 8;
 
    // Counts the structural edits — the ones that move vertices or connectivity
    // on from what the index was drawn against. Nothing invalidates the index,
    // which stays correct across any edit; this only tells buildPointLocation
    // whether redrawing would find anything new, and answers
    // hasCurrentPointLocation. Constraint flags do not count: they change no
    // geometry a walk can see.
    std::size_t revision_ = 0;
    std::size_t pointLocationRevision_ = 0;
 
    // Where the query lands in the hierarchy: the mesh triangle holding it, and
    // whether that is the answer or only where the walk should start. It is the
    // answer when the descent settled the query — it reached the triangle
    // itself, or left the box altogether — and was exact and drawn against the
    // mesh as it now stands. Otherwise the triangle is still one of this mesh,
    // and still beside the query, so the walk goes on from there; and where
    // there is neither index nor triangle, the walk starts where it would have
    // without one.
    [[nodiscard]] TriIndex pointLocationSeed(const PointType& p, bool& answered) const {
        answered = false;
        if (!pointLocation_) {
            return NO_TRI;
        }
        bool settled = false;
        const TriIndex seed = kirkpatrickSeed(p, settled);
        answered = settled && exactDescent && pointLocationRevision_ == revision_;
        return realTriangle(seed) ? seed : NO_TRI;
    }
 
    // A hierarchy vertex, filtered: the mesh's own vertices keep the stored
    // approximation, the four box corners the hierarchy's own.
    [[nodiscard]] auto kirkpatrickVertex(const Kirkpatrick& kp, VertexIndex v) const {
        const std::size_t i = static_cast<std::size_t>(v);
        const std::size_t meshVertices = vertices_.size();
        return i < meshVertices
                   ? filteredVertex(v)
                   : detail::filtered<VertexCoordinate>(kp.extra[i - meshVertices],
                                                        kp.extraApprox, i - meshVertices);
    }
 
    // Where the descent lands: one cell per level, ending on the level-0 cell
    // whose closure holds the query. @p settled says the descent decided the
    // query rather than only narrowing it — it ended on a mesh triangle, or off
    // the box, which the hull is strictly inside of — as against ending on a
    // triangle filling the box, whose closure holds a point of the hull's
    // boundary too, and whose answer is therefore the neighbouring mesh
    // triangle it seeds the walk with.
    template <class QueryPoint>
    [[nodiscard]] TriIndex kirkpatrickSeed(const QueryPoint& p, bool& settled) const {
        const Kirkpatrick& kp = *pointLocation_;
        settled = false;
        // The descent tests p against several fan edges per level, so p is
        // converted once here rather than once per sign.
        const auto q = filteredPoint(p);
        // The fan scan turns about one vertex, so that one is converted once per
        // cell and only the edge's far end is fetched per sign.
        const auto leftOfFrom = [&](const auto& a, VertexIndex b) {
            return !(detail::orientationSignOf(a, kirkpatrickVertex(kp, b), q).value() < 0);
        };
        const auto leftOf = [&](VertexIndex a, VertexIndex b) {
            return leftOfFrom(kirkpatrickVertex(kp, a), b);
        };
 
        // The top level is the only one a query has to be placed in rather than
        // handed down into, so it is the only one whose triangles are tested.
        std::uint32_t cell = ~0u;
        for (std::size_t i = 0; i < kp.roots.size(); ++i) {
            const auto& corners = kp.rootShape[i];
            if (leftOf(corners[0], corners[1]) && leftOf(corners[1], corners[2]) &&
                leftOf(corners[2], corners[0])) {
                cell = kp.roots[i];
                break;
            }
        }
        if (cell == ~0u) {
            settled = true;
            return NO_TRI;  // outside the box, so outside the mesh
        }
 
        for (;;) {
            const auto& current = kp.cells[cell];
            const std::uint32_t count = current.runCount & ~Kirkpatrick::FULL_TURN;
            if (count == 0) {
                settled = current.runCount == 0;
                return static_cast<TriIndex>(current.runBegin);
            }
            const std::uint32_t* const fan = kp.run.data() + current.runBegin;
            cell = fan[2 * count - 1];  // the last child, where the scan runs out
            if (count > 1) {
                const auto apex =
                    kirkpatrickVertex(kp, static_cast<VertexIndex>(current.apex));
                // p is already known to be past the run's first edge — the
                // parent starts there — unless the run is a whole turn, which
                // starts nowhere in particular.
                bool inside = (current.runCount & Kirkpatrick::FULL_TURN) == 0 ||
                              leftOfFrom(apex, static_cast<VertexIndex>(fan[0]));
                for (std::uint32_t m = 1; m < count; ++m) {
                    const bool beyond = leftOfFrom(apex, static_cast<VertexIndex>(fan[2 * m]));
                    if (inside && !beyond) {
                        cell = fan[2 * m - 1];
                        break;
                    }
                    inside = beyond;
                }
            }
        }
    }
 
    // ---- small helpers ---------------------------------------------------
 
    // ---- handles <-> internal indices ------------------------------------
 
    // The public handle for an internal index, and back. An invalid handle
    // holds ~0u, which is exactly NO_TRI / NO_VERTEX read as unsigned, so both
    // directions are the plain cast and the invalid handle needs no special
    // case. A handle a caller made up out of range survives the round trip
    // unchanged and is then rejected by realTriangle / realVertex / inDomain.
    static constexpr TriId triHandle(TriIndex t) {
        return TriId(static_cast<std::uint32_t>(t));
    }
    static constexpr VertexId vertexHandle(VertexIndex v) {
        return VertexId(static_cast<std::uint32_t>(v));
    }
    static constexpr TriIndex indexOf(TriId t) { return static_cast<TriIndex>(t.index()); }
    static constexpr VertexIndex indexOf(VertexId v) {
        return static_cast<VertexIndex>(v.index());
    }
 
    // Position within the stored triple of the vertex the public vertex order
    // starts at. The stored order is counterclockwise but starts anywhere;
    // rotating it onto the lexicographically smallest vertex is exactly what the
    // Triangle constructor does to the same three points, so the public order is
    // the one triangleValue() reports.
    [[nodiscard]] int firstVertex(TriIndex t) const {
        const auto& v = triangles_[static_cast<std::size_t>(t)].v;
        int first = 0;
        for (int i = 1; i < 3; ++i) {
            if (vertices_[static_cast<std::size_t>(v[i])] <
                vertices_[static_cast<std::size_t>(v[first])]) {
                first = i;
            }
        }
        return first;
    }
 
    // The stored side index of public side `side`: the public sides follow the
    // public vertex order (side i joins public vertices i and i+1, as
    // Triangle::edges does), while a stored side is the one *opposite* its
    // vertex, which is the third one.
    [[nodiscard]] int internalSide(TriIndex t, int side) const {
        return (firstVertex(t) + side + 2) % 3;
    }
 
    // True if the index refers to a stored real (non-ghost) triangle, whether
    // in the domain or a hull-fill one outside it.
    [[nodiscard]] bool realTriangle(TriIndex t) const { return t >= 0 && t < firstGhost_; }
 
    // True if the index refers to a stored vertex other than the ghost.
    [[nodiscard]] bool realVertex(VertexIndex v) const {
        return v > GHOST && static_cast<std::size_t>(v) < vertices_.size();
    }
 
    // Hands one triangle to a user visitor: a callable that accepts a Triangle
    // gets the value, one that accepts only a TriId gets the handle (a generic
    // callable accepts both and gets the value). Every triangle visitor of the
    // public interface reports through here, so the rule is the same for all.
    template <class Fn>
    bool reportTriangle(Fn& f, TriIndex t) const {
        if constexpr (std::is_invocable_v<Fn&, const TriangleType&>) {
            return detail::invokeVisitor(f, triangleValue(t));
        } else {
            return detail::invokeVisitor(f, triHandle(t));
        }
    }
 
    // The triangle values of a list of handles, in the same order.
    [[nodiscard]] std::vector<TriangleType> trianglesOf(const std::vector<TriId>& ids) const {
        std::vector<TriangleType> out;
        out.reserve(ids.size());
        for (const TriId id : ids) {
            out.push_back(triangleValue(indexOf(id)));
        }
        return out;
    }
 
    // The vertex p is stored as, or NO_VERTEX when p is not a vertex of the
    // triangulation (or lies outside the hull): locate p, then match it against
    // the vertices of the triangle found.
    [[nodiscard]] VertexIndex vertexIndexAt(const PointType& p) const {
        const TriIndex start = locateIndex(p);
        if (start == NO_TRI || isGhost(start)) {
            return NO_VERTEX;  // empty triangulation, or p outside the hull
        }
        const auto& sv = triangles_[start].v;
        for (int i = 0; i < 3; ++i) {
            if (vertices_[sv[i]] == p) {
                return sv[i];
            }
        }
        return NO_VERTEX;
    }
 
    // A real triangle of vertex w's fan, or NO_TRI when w carries none. The
    // per-vertex hint answers in O(1) when it still holds; a hint left stale by
    // an edit that did not re-register falls back on locating w's position.
    [[nodiscard]] TriIndex fanSeedOf(VertexIndex w) const {
        const TriIndex hinted = incidentTriangleOf(w);
        if (hinted != NO_TRI) {
            return hinted;
        }
        const TriIndex start = locateIndex(vertices_[static_cast<std::size_t>(w)]);
        if (start == NO_TRI) {
            return NO_TRI;
        }
        const auto& v = triangles_[start].v;
        return (v[0] == w || v[1] == w || v[2] == w) ? start : NO_TRI;
    }
 
    // Reads bit `i` of a 3-bit edge mask (e.g. constrainedMask).
    static constexpr bool bit(std::uint8_t mask, int i) { return (mask >> i) & 1; }
    // Sets or clears bit `i` of a 3-bit edge mask.
    static constexpr void setBit(std::uint8_t& mask, int i, bool value) {
        mask = static_cast<std::uint8_t>(value ? (mask | (1u << i)) : (mask & ~(1u << i)));
    }
    // Packs three per-edge flags into a 3-bit mask (bit i is edge i).
    static constexpr std::uint8_t mask(bool b0, bool b1, bool b2) {
        return static_cast<std::uint8_t>((b0 ? 1 : 0) | (b1 ? 2 : 0) | (b2 ? 4 : 0));
    }
 
    // True if t is one of the boundary-closing ghost triangles (stored last).
    [[nodiscard]] bool isGhost(TriIndex t) const { return t >= firstGhost_; }
 
    // A real triangle that is part of the visible triangulation: not a ghost and
    // not a hull-fill triangle outside a polygon's boundary. The public view
    // (sizes, navigation, locate) speaks only of in-domain triangles, while the
    // internal walks still traverse every real triangle, including fill ones.
    [[nodiscard]] bool inDomain(TriIndex t) const {
        return t != NO_TRI && t < firstGhost_ && !triangles_[t].outOfDomain;
    }
 
    // ---- visibility internals (defined in implementation/visibilitygraph.hpp) ---
 
    // True if side `s` of in-domain triangle `t` stops sight: a constrained edge
    // is an opaque wall, and so is the boundary of the domain.
    [[nodiscard]] bool blocksVisibility(TriIndex t, int s) const {
        return bit(triangles_[t].constrainedMask, s) || !inDomain(triangles_[t].nbr[s]);
    }
 
    // One step of the expansion: cross `side` of `tri` carrying the open cone
    // whose clockwise bound is the ray origin->right and whose counterclockwise
    // bound is origin->left. Both bounds are vertices, so every test the
    // traversal makes is one orientation predicate on stored points and no
    // coordinate is ever constructed. The crossed edge runs from its clockwise
    // end `tri.v[(side+1)%3]` to its counterclockwise end `tri.v[(side+2)%3]`.
    struct VisibilityCone {
        TriIndex tri;
        std::int8_t side;
        VertexIndex right;
        VertexIndex left;
    };
 
    // Where an expansion starts from a query point: the cones covering the
    // directions that immediately enter the domain, sorted counterclockwise by
    // their clockwise bound and grouped into contiguous arcs — one arc per lobe
    // the visible region reaches the query point along. `direct` holds the
    // vertices visible without any expansion at all, each flagged when sight to
    // it only grazes (runs along the boundary or a wall) rather than passing
    // through the interior.
    struct VisibilitySeeds {
        std::vector<VisibilityCone> cones;
        std::vector<std::pair<VertexIndex, bool>> direct;
        std::vector<std::size_t> arcs;  // index into `cones` where each arc opens
        bool located = false;           // the query point lies in the domain
        bool fullTurn = false;          // one arc, covering every direction
    };
 
    [[nodiscard]] VisibilitySeeds visibilitySeeds(const PointType& query) const;
 
    // An in-domain triangle whose closure holds `query`, starting from the
    // triangle `start` that locateIndex stopped at. That walk may stop on a ghost or
    // a hull-fill triangle when the query sits on the domain boundary — which
    // every polygon vertex does — so a boundary query needs this to find the
    // triangle it belongs to. NO_TRI when the query really is outside.
    [[nodiscard]] TriIndex inDomainTriangleAt(const PointType& query, TriIndex start) const;
 
    // Drains the open-cone expansion of `start`, calling onVertex(VertexIndex) for
    // each clearly visible vertex met and onBlocked(VisibilityCone) for each
    // blocking edge a cone runs into. Leaves come out counterclockwise, which is
    // what lets regularizedVisiblePolygon assemble its ring without sorting.
    // `scratch` is the traversal stack, passed in so repeated calls reuse it.
    template <class OnVertex, class OnBlocked>
    void expandVisibility(const PointType& origin, VisibilityCone start,
                          std::vector<VisibilityCone>& scratch,
                          OnVertex onVertex, OnBlocked onBlocked) const;
 
    // The vertices clearly visible from `query`, by id. `grazing` selects whether
    // the seeds reached along the boundary come too, which is what separates
    // visibleVertices from clearlyVisibleVertices.
    [[nodiscard]] std::vector<VertexIndex> visibleIds(const PointType& query,
                                                   const VisibilitySeeds& seeds,
                                                   bool grazing) const;
 
    // Clear-visibility adjacency indexed by vertex id (slot GHOST stays empty),
    // built by one triangular expansion per vertex.
    [[nodiscard]] std::vector<std::vector<VertexIndex>> clearVisibleAdjacency() const;
 
    // clearVisibleAdjacency plus the mesh's own blocking edges — together the
    // pairs that see each other with no vertex in between — closed along
    // collinear chains, which is the full visibility relation.
    [[nodiscard]] std::vector<std::vector<VertexIndex>> visibleAdjacency() const;
 
    // The other endpoint of every wall incident to each vertex, indexed by
    // vertex id. Drives the tangency test of reducedVisibilityGraph.
    [[nodiscard]] std::vector<std::vector<VertexIndex>> wallNeighbors() const;
 
    // Whether some line through vertex `m` stays in the domain on both sides of
    // it, so that a visibility segment can pass straight through. Only a
    // strictly convex corner of the domain fails. Deciding whether the collinear
    // closure need look at `m` at all, it is allowed to answer optimistically —
    // never the other way round.
    [[nodiscard]] bool passesThrough(VertexIndex m) const;
 
    // The next mesh vertex met by the ray leaving `current` in the direction
    // `current - tail`, when the segment reaching it stays in the domain and
    // crosses no wall; GHOST when the ray leaves the domain first. `tail` is a
    // point rather than a vertex so a chain can start at a query point.
    [[nodiscard]] VertexIndex nextVertexAlongRay(const PointType& tail, VertexIndex current) const;
 
    // Side of triangle x whose neighbor is `target` (x shares <=1 edge with it).
    [[nodiscard]] std::int8_t findSide(TriIndex x, TriIndex target) const {
        const auto& n = triangles_[x].nbr;
        return static_cast<std::int8_t>(n[0] == target ? 0 : (n[1] == target ? 1 : 2));
    }
 
    // An edge belongs to the visible triangulation if at least one of its two
    // incident triangles is in-domain.
    [[nodiscard]] bool edgeInDomain(Edge e) const {
        return inDomain(e.tri) || inDomain(mirror(e).tri);
    }
 
    // An edge of the domain's interior: an in-domain triangle on both sides, so a
    // neighbourhood of its relative interior lies in the domain.
    [[nodiscard]] bool interiorEdge(Edge e) const {
        return inDomain(e.tri) && inDomain(mirror(e).tri);
    }
 
    // An edge of the domain's boundary (∂D): in-domain on exactly one side, the
    // other being a ghost or a hull-fill triangle.
    [[nodiscard]] bool boundaryEdge(Edge e) const { return edgeInDomain(e) && !interiorEdge(e); }
 
    // Calls f(edge, handle) on the three edges of mesh triangle `t` until one
    // returns true; returns whether one did.
    template <class Fn>
    [[nodiscard]] bool anyEdgeOf(const TriangleType& t, Fn f) const {
        for (const auto& edge : t.edges()) {
            const SegmentType seg(edge[0], edge[1]);
            const auto se = segmentMap().find(seg);
            if (se != segmentMap().end() && f(seg, se->second)) {
                return true;
            }
        }
        return false;
    }
 
    // True if `shape` touches the boundary of the domain. Sound for a shape known
    // to lie in the domain, which is where interiorContains() uses it: every
    // boundary edge such a shape meets belongs to an in-domain triangle it meets,
    // and the traversal reports all of those.
    template <detail::TriangulationQuery Q>
    [[nodiscard]] bool meetsDomainBoundary(const Q& shape) const {
        return visitTrianglesIntersecting(shape, [&](const TriangleType& t) {
            return anyEdgeOf(t, [&](const SegmentType& e, Edge handle) {
                return boundaryEdge(handle) && e.intersects(shape);
            });
        });
    }
 
    // contains() — or, with Interior, interiorContains() — for a shape given by
    // its boundary edges: a chain, or a shape bounded by a closed polygonal chain.
    // Every edge of it must be contained. Zero-length edges are skipped, so a flat
    // rectangle (or a chain with a repeated vertex) is traced along the edges it
    // does have; one whose edges are all degenerate has collapsed to a single
    // point, which then decides.
    template <bool Interior = false, class EdgeRange, class Q>
    [[nodiscard]] bool containsBoundary(const EdgeRange& edges, const Q& shape) const {
        const auto held = [&](const auto& piece) {
            if constexpr (Interior) {
                return interiorContains(piece);
            } else {
                return contains(piece);
            }
        };
        bool traced = false;
        for (const auto& e : edges) {
            if (e[0] == e[1]) {
                continue;
            }
            traced = true;
            if (!held(e)) {
                return false;
            }
        }
        return traced || held(shape[0]);
    }
 
    // Rewrites a ring's vertex loop so that consecutive entries span an
    // unobstructed edge: rings are allowed to touch, and where they do a vertex
    // of one lands in the relative interior of an edge of the other, which
    // insertConstraint cannot force through. Every such vertex is spliced into
    // the loop, in order along the edge carrying it.
    //
    // Only vertices of the *other* rings can obstruct an edge — a ring is simple,
    // so a vertex of its own inside one of its edges would be a self-crossing —
    // but scanning all of @p candidates is simpler and costs the same test.
    std::vector<VertexIndex> expandRing(const std::vector<VertexIndex>& ring,
                                     const std::vector<VertexIndex>& candidates) const {
        std::vector<VertexIndex> expanded;
        expanded.reserve(ring.size());
        std::vector<VertexIndex> onEdge;
        for (std::size_t i = 0; i < ring.size(); ++i) {
            const VertexIndex a = ring[i];
            const VertexIndex b = ring[(i + 1) % ring.size()];
            expanded.push_back(a);
            const SegmentType edge(vertices_[a], vertices_[b]);
            onEdge.clear();
            for (const VertexIndex v : candidates) {
                if (v != a && v != b && edge.contains(vertices_[v])) {
                    onEdge.push_back(v);
                }
            }
            if (onEdge.empty()) {
                continue;
            }
            // The obstructing vertices are collinear with the edge, so ordering
            // them lexicographically orders them along it, forward or backward.
            std::sort(onEdge.begin(), onEdge.end(),
                      [this](VertexIndex p, VertexIndex q) { return vertices_[p] < vertices_[q]; });
            if (vertices_[b] < vertices_[a]) {
                std::reverse(onEdge.begin(), onEdge.end());
            }
            expanded.insert(expanded.end(), onEdge.begin(), onEdge.end());
        }
        return expanded;
    }
 
    // The triangle lying to the left of the directed edge a -> b, or NO_TRI when
    // that edge is not in the triangulation. Triangle vertices are stored
    // counterclockwise and side s spans v[s+1] -> v[s+2], so the triangle
    // carrying the edge in that direction is the one on its left.
    [[nodiscard]] TriIndex triangleLeftOf(VertexIndex a, VertexIndex b) const {
        const auto it = segmentMap().find(SegmentType(vertices_[a], vertices_[b]));
        if (it == segmentMap().end()) {
            return NO_TRI;
        }
        for (const Edge& e : {it->second, mirror(it->second)}) {
            if (e.tri == NO_TRI) {
                continue;
            }
            const auto& v = triangles_[e.tri].v;
            if (v[(e.side + 1) % 3] == a && v[(e.side + 2) % 3] == b) {
                return e.tri;
            }
        }
        return NO_TRI;
    }
 
    // The same edge seen from the triangle on its other side (an invalid Edge if
    // there is none). mirror(e).tri is the neighbor across e.
    [[nodiscard]] Edge mirror(Edge e) const {
        const TriIndex t2 = triangles_[e.tri].nbr[e.side];
        if (t2 == NO_TRI) {
            return Edge{NO_TRI, 0};
        }
        return Edge{t2, findSide(t2, e.tri)};
    }
 
    // Materializes triangle t as a public Triangle value (with its label, if any).
    [[nodiscard]] TriangleType triangleValue(TriIndex t) const {
        const auto& T = triangles_[t];
        TriangleType tri(vertices_[T.v[0]], vertices_[T.v[1]], vertices_[T.v[2]]);
        // Activates automatically once pgl::Triangle gains a mutable label():
        // until then TriangleLabel is NoLabel and this branch is discarded.
        if constexpr (detail::has_label_v<TriangleLabel>) {
            tri.label() = T.triLabel;
        }
        return tri;
    }
 
    // Materializes the edge as a Segment, attaching the label stored on @p e.
    // (Transient handles from mirror()/registerSides carry a default label, used
    // only for key construction where labels are ignored.)
    [[nodiscard]] SegmentType edgeSegment(Edge e) const {
        const auto& T = triangles_[e.tri];
        SegmentType s(vertices_[T.v[(e.side + 1) % 3]], vertices_[T.v[(e.side + 2) % 3]]);
        if constexpr (detail::has_label_v<SegmentLabel>) {
            s.label() = e.segLabel;
        }
        return s;
    }
 
    // The vertex of t not lying on edge s (its apex relative to s).
    [[nodiscard]] static PointType apexOf(const TriangleType& t, const SegmentType& s) {
        for (const auto& p : t.vertices()) {
            if (p != s[0] && p != s[1]) {
                return p;
            }
        }
        return t[0];  // unreachable when s is an edge of t
    }
 
    // Resolves a public Triangle to its internal id, or NO_TRI if absent: look
    // up one of its edges, then pick the incident side whose apex matches t.
    [[nodiscard]] TriIndex idOf(const TriangleType& t) const {
        const auto edges = t.edges();
        auto se = segmentMap().find(edges[0]);
        if (se == segmentMap().end()) {
            return NO_TRI;
        }
        const PointType apex = apexOf(t, edges[0]);
        const Edge e = se->second;
        if (vertices_[triangles_[e.tri].v[e.side]] == apex) {
            return e.tri;
        }
        const Edge m = mirror(e);
        if (m.tri != NO_TRI && !isGhost(m.tri) && vertices_[triangles_[m.tri].v[m.side]] == apex) {
            return m.tri;
        }
        return NO_TRI;
    }
 
    // True if p lies in the closed triangle t (interior or boundary). Assumes t
    // is CCW; @p p may use a different point type.
    template <class QueryPoint>
    [[nodiscard]] bool pointInClosure(const QueryPoint& p, TriIndex t) const {
        const auto& v = triangles_[t].v;
        const auto fq = filteredPoint(p);
        const auto f0 = filteredVertex(v[0]);
        const auto f1 = filteredVertex(v[1]);
        const auto f2 = filteredVertex(v[2]);
        return detail::orientationSignOf(f0, f1, fq).value() >= 0 &&
               detail::orientationSignOf(f1, f2, fq).value() >= 0 &&
               detail::orientationSignOf(f2, f0, fq).value() >= 0;
    }
 
    // Position (0,1,2) of vertex w within triangle t.
    [[nodiscard]] int localIndex(TriIndex t, VertexIndex w) const {
        const auto& v = triangles_[t].v;
        return v[0] == w ? 0 : (v[1] == w ? 1 : 2);
    }
 
    // The neighbour of `cur` reached by rotating around vertex `w`, having
    // arrived from `from` (NO_TRI to start): leave through the w-incident edge we
    // did not enter. Orientation-independent (ghosts are not CCW-normalized), so
    // it steps through the ghost triangles of w's fan as well. Mirrors the lambda
    // in the segment walk; shared by the region flood fill.
    [[nodiscard]] TriIndex rotateAroundVertex(TriIndex cur, VertexIndex w, TriIndex from) const {
        const int lw = localIndex(cur, w);
        int s1 = -1, s2 = -1;
        for (int s = 0; s < 3; ++s) {
            if (s != lw) {
                if (s1 < 0) s1 = s; else s2 = s;
            }
        }
        return triangles_[cur].nbr[triangles_[cur].nbr[s1] == from ? s2 : s1];
    }
 
    // Calls fn(TriIndex) on every triangle of vertex w's fan, in rotational order
    // starting from `start` (which must have w as a vertex). The rotation steps
    // through ghost and fill triangles too — that is what closes the ring — so
    // fn also sees those and must filter with inDomain if it wants only visible
    // triangles.
    template <class Fn>
    void visitVertexFan(TriIndex start, VertexIndex w, Fn fn) const {
        TriIndex cur = start, from = NO_TRI;
        std::size_t g = 0, lim = triangles_.size() + 1;
        do {
            fn(cur);
            const TriIndex next = rotateAroundVertex(cur, w, from);
            from = cur;
            cur = next;
        } while (cur != start && cur != NO_TRI && ++g < lim);
    }
 
    // The first in-domain real triangle, or NO_TRI if the domain is empty.
    [[nodiscard]] TriIndex firstInDomainTriangle() const {
        for (TriIndex t = 0; t < firstGhost_; ++t) {
            if (!triangles_[t].outOfDomain) return t;
        }
        return NO_TRI;
    }
 
    // Calls f(TriIndex) once per connected component of the domain, stopping at the
    // first true. This is what a query falls back on when the shape's boundary
    // misses the domain entirely: every component then lies wholly inside the
    // shape or wholly outside it, and one triangle of each settles which.
    //
    // Only a region's domain can come apart — a hole spanning it leaves two
    // slabs, and two holes meeting at a point pinch it likewise — so a domain
    // without holes keeps answering from its first triangle at O(1). Components
    // are the edge-connected classes of in-domain triangles: sharing an edge
    // puts the open edge in the domain's interior and joins them, sharing only a
    // vertex does not.
    template <class Fn>
    bool anyDomainComponent(Fn&& f) const {
        const TriIndex first = firstInDomainTriangle();
        if (first == NO_TRI) {
            return false;
        }
        if (holeWitnesses_.empty()) {
            return f(first);
        }
        std::vector<char> seen(triangles_.size(), 0);
        std::vector<TriIndex> stack;
        for (TriIndex t = 0; t < firstGhost_; ++t) {
            if (!inDomain(t) || seen[t]) {
                continue;
            }
            if (f(t)) {
                return true;
            }
            seen[t] = 1;
            stack.push_back(t);
            while (!stack.empty()) {
                const TriIndex cur = stack.back();
                stack.pop_back();
                for (int s = 0; s < 3; ++s) {
                    const TriIndex nb = triangles_[cur].nbr[s];
                    if (nb != NO_TRI && !seen[nb] && inDomain(nb)) {
                        seen[nb] = 1;
                        stack.push_back(nb);
                    }
                }
            }
        }
        return false;
    }
 
    // Finds one or more triangles meeting `shape` to seed the region flood fill,
    // by navigation rather than a full scan. Returns an empty vector when `shape`
    // does not meet the triangulated region at all. (Helper for the
    // detail::TriangulationRegionQuery overload of visitTrianglesIntersecting.)
    template <class Q>
    [[nodiscard]] std::vector<TriIndex> seedTrianglesIntersecting(const Q& shape) const {
        std::vector<TriIndex> seeds;
        if constexpr (PointConcept<Q>) {
            // A point: the triangle locate lands in already meets it (its closure
            // contains the point). Outside the hull, locate returns a ghost.
            const TriIndex t = locateIndex(shape);
            if (t != NO_TRI && !isGhost(t)) {
                seeds.push_back(t);
            }
        } else if constexpr (HalfplaneConcept<Q>) {
            // Unbounded shape: it meets the triangulated region iff it meets some
            // convex-hull edge. Scan the ghost ring for the first such edge; the
            // real triangle just inside it is a seed. (O(hull).) Segments, lines,
            // oriented lines, and rays are not region queries — they have their
            // own ordered walk.
            for (TriIndex g = firstGhost_; g < static_cast<TriIndex>(triangles_.size()); ++g) {
                if (edgeSegment(Edge{g, 2}).intersects(shape)) {
                    seeds.push_back(triangles_[g].nbr[2]);  // the real triangle inside
                    break;
                }
            }
        } else if constexpr (DiskConcept<Q>) {
            // A disk is bounded but has no straight edges to trace. One of its
            // defining boundary points, shape[0], lies in the (closed) disk, so
            // the triangle locate lands in is a seed when that point is inside
            // the hull. Otherwise the disk reaches in from outside the hull, so —
            // as for an unbounded shape — a hull edge it crosses gives the seed.
            const TriIndex t = locateIndex(shape[0]);
            if (t != NO_TRI && !isGhost(t)) {
                seeds.push_back(t);
            } else {
                for (TriIndex g = firstGhost_; g < static_cast<TriIndex>(triangles_.size()); ++g) {
                    if (edgeSegment(Edge{g, 2}).intersects(shape)) {
                        seeds.push_back(triangles_[g].nbr[2]);
                        break;
                    }
                }
            }
        } else {
            // Bounded shape with straight edges (triangle, rectangle, convex
            // polygon): any triangle its boundary meets is a seed. Trace each edge
            // with the segment walk and take the first triangle it reports.
            for (const auto& e : shape.edges()) {
                TriIndex found = NO_TRI;
                visitTrianglesIntersecting(e, [&](const TriangleType& t) {
                    found = idOf(t);
                    return true;  // the first triangle suffices
                });
                if (found != NO_TRI) {
                    seeds.push_back(found);
                    break;
                }
            }
            // No edge met the domain: each connected piece of it is either
            // inside `shape` or disjoint from it, and one vertex of each decides.
            // Every piece inside gets a seed — a hole can split the domain, and
            // the flood does not cross from one piece to another.
            if (seeds.empty() && numVertices() > 0) {
                anyDomainComponent([&](TriIndex t) {
                    if (shape.contains(vertices_[triangles_[t].v[0]])) {
                        seeds.push_back(t);
                    }
                    return false;  // every component, not just the first
                });
            }
        }
        return seeds;
    }
 
    // Grows the set of triangles meeting `shape` outward from `seeds` by a flood
    // fill over edge/vertex adjacency, emitting the in-domain ones to `f` (every
    // real one under AllTriangles, as in the directed walk — what contains() needs
    // to see the fill region a shape reaches into). `f` takes a TriIndex and returns
    // bool. Every real (non-ghost) triangle meeting `shape` and reachable from a
    // seed through other meeting triangles is visited; ghosts are walls. A
    // per-triangle mark bit stands in for a visited set and is cleared on every
    // exit path (normal, early-stop, or an exception from f), so the const query
    // leaves the mesh pristine. Like the segment walk, it is not reentrant: f must
    // not start another walk on the same Triangulation. Returns whether f stopped
    // early.
    template <bool AllTriangles, class Q, class Fn>
    bool floodTriangleIdsIntersecting(const Q& shape, const std::vector<TriIndex>& seeds, Fn f) const {
        std::vector<TriIndex> marked;
        struct MarkClearer {
            const std::vector<Tri>& tris;
            const std::vector<TriIndex>& marked;
            ~MarkClearer() { for (TriIndex t : marked) tris[t].walkMark = 0; }
        } markClearer{triangles_, marked};
 
        std::vector<TriIndex> stack;
        // Marks t seen; if it is a real triangle meeting `shape`, queues it for
        // expansion. Rejected triangles are marked too, so each is tested once.
        const auto consider = [&](TriIndex t) {
            if (t == NO_TRI || isGhost(t) || triangles_[t].walkMark) {
                return;
            }
            triangles_[t].walkMark = 1;
            marked.push_back(t);
            if (triangleValue(t).intersects(shape)) {
                stack.push_back(t);
            }
        };
        for (const TriIndex s : seeds) {
            consider(s);
        }
 
        bool stop = false;
        while (!stop && !stack.empty()) {
            const TriIndex t = stack.back();
            stack.pop_back();
            if ((AllTriangles || inDomain(t)) && f(t)) {
                stop = true;
                break;
            }
            // Expand through every triangle sharing a vertex with t (which covers
            // its edge neighbours too): the meeting set is connected under this
            // adjacency, so this reaches all of it without scanning the mesh.
            for (const VertexIndex w : triangles_[t].v) {
                TriIndex cur = t, from = NO_TRI;
                std::size_t g = 0, lim = triangles_.size() + 1;
                do {
                    consider(cur);
                    const TriIndex next = rotateAroundVertex(cur, w, from);
                    from = cur;
                    cur = next;
                } while (cur != t && cur != NO_TRI && ++g < lim);
            }
        }
        return stop;
    }
 
    // Inserts t's three edges into the segment-to-edge map, keyed by Segment,
    // and records t as the fan seed of each of its vertices. The two go
    // together: every edit that rewrites a triangle's vertices re-registers it,
    // so hooking the seed here is what keeps it fresh across a flip — the flip
    // re-registers both rewritten triangles, and a vertex of the quad ends up in
    // at least one of them.
    void registerSides(TriIndex t) {
        for (std::int8_t s = 0; s < 3; ++s) {
            segmentMap()[edgeSegment(Edge{t, s})] = Edge{t, s};
        }
        noteVertexIncidence(t);
    }
 
    // Records t as the fan seed of each of its real vertices. Ghost triangles
    // are skipped: a fan rotation may pass through them but never starts there.
    void noteVertexIncidence(TriIndex t) {
        if (t < 0 || t >= firstGhost_) {
            return;
        }
        if (vertexTri_.size() < vertices_.size()) {
            vertexTri_.resize(vertices_.size(), NO_TRI);
        }
        for (const VertexIndex w : triangles_[t].v) {
            if (w != GHOST) {
                vertexTri_[static_cast<std::size_t>(w)] = t;
            }
        }
    }
 
    // A real triangle having w as a vertex, or NO_TRI when none is recorded.
    // The record is verified against the triangle's current vertices, so it is
    // authoritative when it answers: a triangle that still lists w really is in
    // w's fan.
    [[nodiscard]] TriIndex incidentTriangleOf(VertexIndex w) const {
        if (static_cast<std::size_t>(w) >= vertexTri_.size()) {
            return NO_TRI;
        }
        const TriIndex t = vertexTri_[static_cast<std::size_t>(w)];
        if (t == NO_TRI || t >= firstGhost_) {
            return NO_TRI;
        }
        const auto& v = triangles_[t].v;
        return (v[0] == w || v[1] == w || v[2] == w) ? t : NO_TRI;
    }
 
    // Fills the segment-to-edge map from the mesh, over all real triangles.
    //
    // A bulk build leaves the map unfilled and comes here on the first lookup
    // instead. Registering an edge costs a hash insert keyed by a segment --
    // two exact points copied and hashed -- three per triangle, and a
    // triangulation that is only ever asked about triangles, points or
    // neighbours never reads one of them. What the build cannot defer is the
    // per-triangle vertex-incidence hint, which the walks do read; that is why
    // the two are registered apart.
    //
    // Only the edge handles are recovered here. An edge's segment label is set
    // on the record afterwards, by whoever had the labelled segment, and every
    // site that does so goes through the accessor first.
    void materializeSegmentMap() const {
        segToEdge_.clear();
        // Three sides a triangle, each shared by at most two of them: the map
        // ends up between one and a half and three entries per triangle, and
        // sizing it once beats rehashing it as it fills.
        segToEdge_.reserve(static_cast<std::size_t>(std::max(firstGhost_, TriIndex(0))) * 2);
        for (TriIndex t = 0; t < firstGhost_; ++t) {
            for (std::int8_t s = 0; s < 3; ++s) {
                segToEdge_[edgeSegment(Edge{t, s})] = Edge{t, s};
            }
        }
        mapStale_ = false;
    }
 
    // The segment-to-edge map, filled first if a bulk build left it deferred.
    [[nodiscard]] std::unordered_map<SegmentType, Edge>& segmentMap() const {
        if (mapStale_) {
            materializeSegmentMap();
        }
        return segToEdge_;
    }
 
    // Returns a function that maps a point to its vertex id, appending it to
    // vertices_ (and to `vid`) the first time it is seen.
    auto makeVertexInterner(std::unordered_map<PointType, VertexIndex>& vid) {
        return [this, &vid](const PointType& p) -> VertexIndex {
            auto it = vid.find(p);
            if (it != vid.end()) {
                return it->second;
            }
            VertexIndex id = static_cast<VertexIndex>(vertices_.size());
            vertices_.push_back(p);
            vid.emplace(p, id);
            return id;
        };
    }
 
    // Exact Delaunay triangle triples (CCW, vertex indices into `pts`) via
    // incremental Bowyer–Watson insertion. A single symbolic "vertex at infinity"
    // (INF) closes the convex hull instead of a containing super-triangle: every
    // hull edge a->b carries a ghost triangle {b,a,INF} that tiles the exterior,
    // so the working triangulation is a closed surface in which every triangle has
    // three neighbours. No far-away coordinates are ever constructed, so every
    // in-circle / orientation test runs on the real point coordinates at the
    // library's normal exactness — no BigInt, and no magnitude assumption (which
    // is why a super-triangle is unsound for rational points: its required scale
    // is unbounded).
    //
    // Each point is located by a visibility walk over the neighbour links (the
    // same walk locateIndex() runs on the finished structure) and inserted by
    // carving out the triangles whose open circumdisk contains it — found by a
    // local flood-fill from the located triangle, not a global scan — then
    // re-fanning the star-shaped cavity to the new vertex. With the random
    // insertion order the walk is short (seeded from the previously inserted
    // triangle), so the build is ~O(n^1.5) here rather than the O(n^2) of testing
    // every triangle against every point.
    static std::vector<std::array<VertexIndex, 3>>
    delaunayTriples(const std::vector<PointType>& pts,
                    const std::vector<detail::ApproximatePoint>& approximations) {
        const VertexIndex n = static_cast<VertexIndex>(pts.size());
        std::vector<std::array<VertexIndex, 3>> out;
        if (n < 3) {
            return out;
        }
        const VertexIndex INF = n;  // the symbolic vertex at infinity
 
        // Every predicate below reads its operands out of `pts`, and the walk
        // and the flood-fill revisit the same vertices over and over, so the
        // approximations the caller keeps are what stop each of those reads
        // from converting a coordinate to double again.
        const auto fp = [&](VertexIndex v) {
            const auto index = static_cast<std::size_t>(v);
            return detail::filtered<VertexCoordinate>(pts[index], approximations, index);
        };
 
        // Local closed triangulation: CCW vertices (ghosts contain INF) and three
        // neighbours each (nbr[i] is across the edge opposite v[i]). Killed
        // triangles keep `dead` set and recycle their slots through freeList; no
        // live triangle ever references a dead slot, so `dead` doubles as the
        // per-insertion "already in the cavity" mark during the flood-fill.
        struct LTri {
            std::array<VertexIndex, 3> v{};
            std::array<int, 3> nbr{-1, -1, -1};
            bool dead = false;
        };
        std::vector<LTri> tri;
        std::vector<int> freeList;
 
        auto newTri = [&](VertexIndex a, VertexIndex b, VertexIndex d) -> int {
            int id;
            if (!freeList.empty()) {
                id = freeList.back();
                freeList.pop_back();
                tri[id] = LTri{{a, b, d}, {-1, -1, -1}, false};
            } else {
                id = static_cast<int>(tri.size());
                tri.push_back(LTri{{a, b, d}, {-1, -1, -1}, false});
            }
            return id;
        };
 
        auto isGhost = [&](int t) {
            const auto& q = tri[t].v;
            return q[0] == INF || q[1] == INF || q[2] == INF;
        };
 
        // Inside-circumdisk test, finite or ghost. The circle through two finite
        // points and INF degenerates to the line through them, so for a ghost
        // {.,.,INF} "inside the open disk" reduces to "left of the finite directed
        // edge" (the edge opposite INF, in CCW order so its left is the hull
        // exterior).
        //
        // The collinear boundary case (orientation 0) is decided so the re-fan
        // never builds a zero-area triangle: p exactly on the open hull-edge
        // segment counts as inside, pulling the ghost into the cavity so the edge
        // splits in two instead of spanning a degenerate {edge, p}. p on the
        // edge's line but beyond an endpoint stays outside — the hull just extends
        // straight along the line, no degenerate triangle either way.
        auto inDisk = [&](int t, VertexIndex p) -> bool {
            const auto& q = tri[t].v;
            const int inf = q[0] == INF ? 0 : (q[1] == INF ? 1 : (q[2] == INF ? 2 : -1));
            if (inf < 0) {
                return detail::inCircleSignOf(fp(q[0]), fp(q[1]),
                                              fp(q[2]), fp(p)) ==
                       std::partial_ordering::greater;
            }
            const VertexIndex u = q[(inf + 1) % 3];
            const VertexIndex w = q[(inf + 2) % 3];
            const auto side =
                detail::orientationSignOf(fp(u), fp(w), fp(p)).value();
            if (side > 0) {
                return true;
            }
            if (side == 0) {  // collinear: inside iff strictly between u and w
                // Betweenness only reads the signs of the two dot products, and
                // dotSign takes them in the promoted coordinate type — the bare
                // dot product multiplies in the coordinate type, where a wrap
                // would silently corrupt the triangulation.
                return dotSign(pts[p] - pts[u], pts[w] - pts[u]) > 0 &&
                       dotSign(pts[p] - pts[w], pts[u] - pts[w]) > 0;
            }
            return false;
        };
 
        // Seed with the first non-collinear triple, oriented CCW, plus the three
        // ghost triangles covering its hull edges. Points 2..c-1 (if any) are
        // collinear with 0 and 1 and get inserted in the main loop like any other.
        VertexIndex c = 2;
        while (c < n &&
               detail::orientationSignOf(fp(0), fp(1), fp(c)).value() == 0) {
            ++c;
        }
        if (c == n) {
            return out;  // all points collinear: the Delaunay triangulation is empty
        }
        const std::array<VertexIndex, 3> seed =
            detail::orientationSignOf(fp(0), fp(1), fp(c)).value() > 0
                ? std::array<VertexIndex, 3>{0, 1, c}
                : std::array<VertexIndex, 3>{1, 0, c};
        const int seedTris[4] = {
            newTri(seed[0], seed[1], seed[2]),
            newTri(seed[1], seed[0], INF),  // ghost outside edge seed0->seed1
            newTri(seed[2], seed[1], INF),  // ghost outside edge seed1->seed2
            newTri(seed[0], seed[2], INF),  // ghost outside edge seed2->seed0
        };
        // Link the seed: match each directed edge (a,b) to the neighbour carrying
        // its reverse (b,a). The four triangles tile the sphere, so every side
        // finds a partner.
        for (int t : seedTris) {
            for (int s = 0; s < 3; ++s) {
                const VertexIndex a = tri[t].v[(s + 1) % 3];
                const VertexIndex b = tri[t].v[(s + 2) % 3];
                for (int u : seedTris) {
                    for (int q = 0; q < 3; ++q) {
                        if (tri[u].v[(q + 1) % 3] == b && tri[u].v[(q + 2) % 3] == a) {
                            tri[t].nbr[s] = u;
                        }
                    }
                }
            }
        }
 
        // Visibility walk: from triangle t, step across whichever edge p lies
        // strictly to the right of (outside), until p is inside the current
        // triangle (no such edge) or a ghost is reached (p outside the hull). A
        // randomised start edge guarantees termination; the step cap is a safety
        // net only.
        std::uint64_t rngState = 0x9e3779b97f4a7c15ULL;
        auto walk = [&](VertexIndex p, int t) -> int {
            int from = -1;
            const int64_t cap = int64_t(3) * static_cast<int64_t>(tri.size()) + 16;
            for (int64_t step = 0; step < cap; ++step) {
                if (isGhost(t)) {
                    return t;
                }
                rngState = rngState * 6364136223846793005ULL + 1442695040888963407ULL;
                const int begin = static_cast<int>((rngState >> 33) % 3);
                int next = -1;
                for (int k = 0; k < 3; ++k) {
                    const int s = (begin + k) % 3;
                    if (tri[t].nbr[s] == from) {
                        continue;
                    }
                    const VertexIndex ea = tri[t].v[(s + 1) % 3];
                    const VertexIndex eb = tri[t].v[(s + 2) % 3];
                    if (detail::orientationSignOf(fp(ea), fp(eb),
                                                  fp(p)).value() < 0) {
                        next = tri[t].nbr[s];
                        break;
                    }
                }
                if (next < 0) {
                    return t;
                }
                from = t;
                t = next;
            }
            return t;
        };
 
        std::vector<int> cavity;
        // A boundary edge of the cavity: its directed edge (a,b), the surviving
        // triangle behind it, and that triangle's side facing the cavity.
        struct Bnd {
            VertexIndex a, b;
            int surv, survSide;
        };
        std::vector<Bnd> boundary;
        // Spoke edges {vertex,p} awaiting their partner among the new triangles.
        struct Spoke {
            VertexIndex vertex;
            int tri, side;
        };
        std::vector<Spoke> spokes;
 
        int hint = seedTris[0];
        for (VertexIndex i = 0; i < n; ++i) {
            if (i == seed[0] || i == seed[1] || i == seed[2]) {
                continue;
            }
            const int start = walk(i, hint);
            if (!inDisk(start, i)) {
                continue;  // i lies in no open circumdisk (collinear/duplicate): skip
            }
 
            // Flood-fill the cavity: every triangle whose open disk contains i,
            // reachable from `start` through neighbours. Mark them dead as we go
            // (so they double as the visited set); record each edge where the
            // flood meets a surviving triangle.
            cavity.clear();
            boundary.clear();
            tri[start].dead = true;
            cavity.push_back(start);
            for (std::size_t qi = 0; qi < cavity.size(); ++qi) {
                const int t = cavity[qi];
                for (int s = 0; s < 3; ++s) {
                    const int nb = tri[t].nbr[s];
                    if (tri[nb].dead) {
                        continue;  // already carved into the cavity this insertion
                    }
                    if (inDisk(nb, i)) {
                        tri[nb].dead = true;
                        cavity.push_back(nb);
                    } else {
                        const int survSide =
                            tri[nb].nbr[0] == t ? 0 : (tri[nb].nbr[1] == t ? 1 : 2);
                        boundary.push_back(
                            {tri[t].v[(s + 1) % 3], tri[t].v[(s + 2) % 3], nb, survSide});
                    }
                }
            }
            for (int t : cavity) {
                freeList.push_back(t);
            }
 
            // Re-fan: one new triangle {a,b,i} per boundary edge, linked across its
            // base {a,b} to the surviving triangle and across its two spokes
            // {b,i}/{a,i} to the adjacent new triangles, paired by shared vertex.
            spokes.clear();
            int newReal = -1;
            for (const Bnd& e : boundary) {
                const int nt = newTri(e.a, e.b, i);
                tri[nt].nbr[2] = e.surv;  // side 2 (opposite i) is the base edge {a,b}
                tri[e.surv].nbr[e.survSide] = nt;
                if (newReal < 0 && !isGhost(nt)) {
                    newReal = nt;
                }
                // side 0 (opposite a) is edge {b,i}; side 1 (opposite b) is {a,i}.
                for (const auto& [vertex, side] :
                     {std::pair<VertexIndex, int>{e.b, 0}, std::pair<VertexIndex, int>{e.a, 1}}) {
                    bool paired = false;
                    for (std::size_t k = 0; k < spokes.size(); ++k) {
                        if (spokes[k].vertex == vertex) {
                            tri[nt].nbr[side] = spokes[k].tri;
                            tri[spokes[k].tri].nbr[spokes[k].side] = nt;
                            spokes[k] = spokes.back();
                            spokes.pop_back();
                            paired = true;
                            break;
                        }
                    }
                    if (!paired) {
                        spokes.push_back({vertex, nt, side});
                    }
                }
            }
            if (newReal >= 0) {
                hint = newReal;
            }
        }
 
        for (int t = 0; t < static_cast<int>(tri.size()); ++t) {
            if (tri[t].dead) {
                continue;
            }
            const auto& q = tri[t].v;
            if (q[0] != INF && q[1] != INF && q[2] != INF) {
                out.push_back({q[0], q[1], q[2]});
            }
        }
        return out;
    }
 
    // ---- constrained Delaunay (polygon constructor) ----------------------
 
    // True if the open segments AB and CD cross properly (no shared endpoint,
    // no collinear contact). Exact via orientationSign.
    [[nodiscard]] bool properCross(const PointType& a, const PointType& b,
                                   const PointType& c, const PointType& d) const {
        const auto s1 = orientationSign(a, b, c);
        const auto s2 = orientationSign(a, b, d);
        const auto s3 = orientationSign(c, d, a);
        const auto s4 = orientationSign(c, d, b);
        if (s1 == 0 || s2 == 0 || s3 == 0 || s4 == 0) {
            return false;
        }
        return ((s1 > 0) != (s2 > 0)) && ((s3 > 0) != (s4 > 0));
    }
 
    // The current internal handle of edge {p,q}, or an invalid edge if absent.
    [[nodiscard]] Edge edgeHandle(VertexIndex p, VertexIndex q) const {
        auto se = segmentMap().find(SegmentType(vertices_[p], vertices_[q]));
        return se == segmentMap().end() ? Edge{NO_TRI, 0} : se->second;
    }
 
    [[nodiscard]] bool edgeExists(VertexIndex p, VertexIndex q) const {
        return segmentMap().contains(SegmentType(vertices_[p], vertices_[q]));
    }
 
    // The interior edges that the open segment va->vb crosses, as vertex pairs,
    // in order from va to vb. Empty when {va,vb} is already an edge. Assumes no
    // vertex lies in the interior of the segment (true for simple-polygon edges).
    [[nodiscard]] std::vector<std::pair<VertexIndex, VertexIndex>>
    collectCrossings(VertexIndex va, VertexIndex vb) const {
        const PointType& A = vertices_[va];
        const PointType& B = vertices_[vb];
        std::vector<std::pair<VertexIndex, VertexIndex>> out;
 
        // Find the triangle incident to va that the segment first enters. Only
        // va's own fan can hold it, so the search rotates around va rather than
        // scanning the mesh: the scan is what made inserting k constraints cost
        // O(k * triangles), which dominated every constrained triangulation
        // here. Under this method's precondition the answer is unique — the
        // link of va is star-shaped about it, so a segment leaving va properly
        // crosses exactly one link edge — and the visit order therefore does not
        // matter.
        TriIndex t = NO_TRI;
        std::pair<VertexIndex, VertexIndex> entry{NO_TRI, NO_TRI};
        const auto enterFrom = [&](TriIndex k) {
            if (t != NO_TRI || k == NO_TRI || isGhost(k)) {
                return;  // already found, or a ghost triangle closing the fan
            }
            const auto& v = triangles_[k].v;
            const int i = localIndex(k, va);
            const VertexIndex p = v[(i + 1) % 3];
            const VertexIndex q = v[(i + 2) % 3];
            if (properCross(A, B, vertices_[p], vertices_[q])) {
                t = triangles_[k].nbr[i];
                entry = {p, q};
                out.push_back({p, q});
            }
        };
        const TriIndex seed = incidentTriangleOf(va);
        if (seed != NO_TRI) {
            visitVertexFan(seed, va, enterFrom);
        } else {
            // No seed recorded (a mesh built by a path that never registered
            // one): fall back to the scan, which needs no incidence at all.
            for (TriIndex k = 0; k < firstGhost_ && t == NO_TRI; ++k) {
                const auto& v = triangles_[k].v;
                if (v[0] == va || v[1] == va || v[2] == va) {
                    enterFrom(k);
                }
            }
        }
        if (t == NO_TRI) {
            return out;  // {va,vb} is already an edge (or va not incident anywhere)
        }
 
        // Walk across the crossed edges until a triangle holds vb.
        std::size_t guard = 0, cap = triangles_.size() + 1;
        while (++guard < cap) {
            const auto& v = triangles_[t].v;
            if (v[0] == vb || v[1] == vb || v[2] == vb) {
                break;
            }
            int exitK = -1;
            for (int k = 0; k < 3; ++k) {
                const VertexIndex p = v[(k + 1) % 3];
                const VertexIndex q = v[(k + 2) % 3];
                const bool sameEntry = (p == entry.first && q == entry.second) ||
                                       (p == entry.second && q == entry.first);
                if (sameEntry) {
                    continue;
                }
                if (properCross(A, B, vertices_[p], vertices_[q])) {
                    exitK = k;
                    entry = {p, q};
                    out.push_back({p, q});
                    break;
                }
            }
            if (exitK < 0) {
                break;  // should not happen for a valid triangulation
            }
            t = triangles_[t].nbr[exitK];
        }
        return out;
    }
 
    // Inserts the segment {va,vb} as an edge by flipping the edges it crosses
    // (Sloan's flip algorithm), then flags it constrained on both sides.
    //
    // A FIFO queue of crossing edges is processed front-to-back: a convex
    // crossing edge is flipped, and if its new diagonal still crosses the
    // constraint it goes to the back; a non-convex one is deferred to the back
    // until its quad becomes convex. This ordering guarantees progress (a naive
    // "flip the first convex one" can oscillate, repeatedly flipping a diagonal
    // back and forth).
    void insertConstraint(VertexIndex va, VertexIndex vb) {
        if (va == vb || edgeExists(va, vb)) {
            setConstrained(SegmentType(vertices_[va], vertices_[vb]), true);
            return;
        }
        const PointType& A = vertices_[va];
        const PointType& B = vertices_[vb];
        std::deque<std::pair<VertexIndex, VertexIndex>> queue;
        for (const auto& pq : collectCrossings(va, vb)) {
            queue.push_back(pq);
        }
        std::size_t guard = 0, cap = (queue.size() + 1) * (triangles_.size() + 1) * 4 + 64;
        while (!queue.empty() && ++guard < cap) {
            auto [p, q] = queue.front();
            queue.pop_front();
            const Edge e = edgeHandle(p, q);
            if (e.tri == NO_TRI || !properCross(A, B, vertices_[p], vertices_[q])) {
                continue;  // edge gone, or no longer crosses the constraint
            }
            if (!flippableEdge(e)) {
                queue.push_back({p, q});  // quad not yet convex; revisit later
                continue;
            }
            const VertexIndex r = triangles_[e.tri].v[e.side];        // apex on one side
            const Edge m = mirror(e);
            const VertexIndex l = triangles_[m.tri].v[m.side];        // apex on the other
            flip(SegmentType(vertices_[p], vertices_[q]));
            if (properCross(A, B, vertices_[r], vertices_[l])) {
                queue.push_back({r, l});  // new diagonal still crosses; re-queue
            }
        }
        setConstrained(SegmentType(vertices_[va], vertices_[vb]), true);
    }
 
    // Restores the constrained Delaunay property: flip every non-constrained
    // interior edge whose opposite apex lies inside the incident circumcircle,
    // until none remain. Constrained edges are never flipped.
    //
    // A triangulation is Delaunay exactly when every edge is locally Delaunay,
    // and an edge's local test reads only its two incident triangles — so
    // @ref legalize's work queue reaches the same fixpoint as a sweep, for the
    // cost of the flips instead of the cost of the mesh per flip. Seeding it
    // with every edge is what makes it a *global* restore: both callers reach
    // here having forced constraints into a mesh whose edges have not been
    // tested since, so none can be assumed legal.
    void restoreConstrainedDelaunay() {
        std::vector<SegmentType> suspect;
        suspect.reserve(segmentMap().size());
        for (const auto& [seg, handle] : segmentMap()) {
            (void)handle;
            suspect.push_back(seg);
        }
        legalize(suspect);
    }
 
    // Flood-fills from the ghost (outside) across non-constrained edges, marking
    // every real triangle reachable without crossing the polygon boundary as
    // out-of-domain (the hull-fill triangles between polygon and convex hull).
    //
    // A region with holes seeds the same flood once more per hole, from a
    // triangle inside it: hole interiors are fenced off from the outside by the
    // outer boundary, so nothing else would reach them.
    void markOutOfDomain(const std::vector<TriIndex>& holeSeeds = {}) {
        std::vector<char> seen(triangles_.size(), 0);
        std::vector<TriIndex> stack;
        for (TriIndex g = firstGhost_; g < static_cast<TriIndex>(triangles_.size()); ++g) {
            seen[g] = 1;
            stack.push_back(g);
        }
        std::size_t marked = 0;
        for (const TriIndex seed : holeSeeds) {
            if (seed == NO_TRI || seen[seed]) {
                continue;  // two rings bounding the same triangle: seeded already
            }
            seen[seed] = 1;
            triangles_[seed].outOfDomain = 1;
            ++marked;
            stack.push_back(seed);
        }
        while (!stack.empty()) {
            const TriIndex t = stack.back();
            stack.pop_back();
            for (int s = 0; s < 3; ++s) {
                if (bit(triangles_[t].constrainedMask, s)) {
                    continue;  // the polygon boundary fences off the interior
                }
                const TriIndex nb = triangles_[t].nbr[s];
                if (nb == NO_TRI || seen[nb]) {
                    continue;
                }
                seen[nb] = 1;
                if (!isGhost(nb)) {
                    triangles_[nb].outOfDomain = 1;
                    ++marked;
                }
                stack.push_back(nb);
            }
        }
        domainTriangleCount_ = static_cast<std::size_t>(firstGhost_) - marked;
    }
 
    // The outer rings and, respectively, all the hole rings of a set's
    // components, flattened for @ref constructConstrained — which asks for the
    // two ring kinds separately because only holes need a flood seed of their
    // own.
    static std::vector<Polygon<PointType>> setOuters(const PolygonSet<PointType>& set) {
        std::vector<Polygon<PointType>> outers;
        outers.reserve(set.componentCount());
        for (const auto& component : set) {
            outers.push_back(component.outer());
        }
        return outers;
    }
 
    static std::vector<Polygon<PointType>> setHoles(const PolygonSet<PointType>& set) {
        std::vector<Polygon<PointType>> holes;
        holes.reserve(set.holeCount());
        for (const auto& component : set) {
            for (const auto& hole : component.holes()) {
                holes.push_back(hole);
            }
        }
        return holes;
    }
 
    // Shared implementation of the polygon, region and region-set constructors.
    // Triangulates the union of the ring vertices, the extra interior points,
    // and the constraint segment endpoints; constrains every ring and every
    // interior segment; then marks the exterior (between the rings and their
    // convex hull) out of domain, along with each hole. Interior constraints
    // never reach the exterior flood — the boundary fences it off — so they stay
    // in-domain. @p extraPoints and @p constraintSegments are assumed to lie in
    // the domain (not checked).
    //
    // @p outers holds one outer ring per piece. Several of them need nothing
    // extra: the gap between two pieces is reached by the same exterior flood,
    // which the pieces' own rings fence it out of, and a piece stranded inside
    // another's hole is fenced out of that hole's flood the same way.
    template <class PointRange, class SegmentRange>
    void constructConstrained(const std::vector<Polygon<PointType>>& outers,
                              const PointRange& extraPoints,
                              const SegmentRange& constraintSegments,
                              const std::vector<Polygon<PointType>>& holes = {}) {
        std::unordered_map<PointType, VertexIndex> vid;
        const auto idOfPoint = makeVertexInterner(vid);
        for (const auto& outer : outers) {
            for (std::size_t i = 0; i < outer.size(); ++i) {
                idOfPoint(outer[i]);
            }
        }
        for (const auto& hole : holes) {
            for (std::size_t i = 0; i < hole.size(); ++i) {
                idOfPoint(hole[i]);
            }
        }
        for (const auto& p : extraPoints) {
            idOfPoint(PointType(p));
        }
        for (const auto& s : constraintSegments) {
            idOfPoint(PointType(s[0]));
            idOfPoint(PointType(s[1]));
        }
        // Keep vertices_ in Hilbert order (see the point-set constructor).
        hilbertSort(vertices_);
        syncVertexApproximations();
        auto triples = delaunayTriples(vertices_, vertexApproximations_);
        buildFromTriples(triples, std::vector<TriangleLabel>(triples.size()));
 
        // The interner's ids went stale twice — the Hilbert reorder and the
        // ghost vertex buildFromTriples prepended — so resolve the boundary
        // loop and constraint endpoints against the final ids (reals start
        // at 1; index 0 is the ghost, whose placeholder coordinates must not
        // shadow a real vertex in the map).
        vid.clear();
        for (VertexIndex i = 1; i < static_cast<VertexIndex>(vertices_.size()); ++i) {
            vid.emplace(vertices_[i], i);
        }
        std::vector<std::vector<VertexIndex>> outerLoops;
        outerLoops.reserve(outers.size());
        for (const auto& outer : outers) {
            std::vector<VertexIndex> loop;
            loop.reserve(outer.size());
            for (std::size_t i = 0; i < outer.size(); ++i) {
                loop.push_back(vid.at(outer[i]));
            }
            outerLoops.push_back(std::move(loop));
        }
 
        std::vector<std::vector<VertexIndex>> holeLoops;
        holeLoops.reserve(holes.size());
        for (const auto& hole : holes) {
            std::vector<VertexIndex> ring;
            ring.reserve(hole.size());
            for (std::size_t i = 0; i < hole.size(); ++i) {
                ring.push_back(vid.at(hole[i]));
            }
            holeLoops.push_back(std::move(ring));
        }
 
        // Where rings touch, one ring's vertex can sit inside another ring's
        // edge; splice those in so every constrained edge is unobstructed. Only
        // several rings can produce them, so a lone polygon skips the scan.
        if (!holes.empty() || outers.size() > 1) {
            std::vector<VertexIndex> ringVertices;
            for (const auto& ring : outerLoops) {
                ringVertices.insert(ringVertices.end(), ring.begin(), ring.end());
            }
            for (const auto& ring : holeLoops) {
                ringVertices.insert(ringVertices.end(), ring.begin(), ring.end());
            }
            for (auto& ring : outerLoops) {
                ring = expandRing(ring, ringVertices);
            }
            for (auto& ring : holeLoops) {
                ring = expandRing(ring, ringVertices);
            }
        }
 
        // Constrain every outer ring, every hole ring, and every interior
        // segment, restore the constrained Delaunay property, then carve away
        // the exterior and the hole interiors.
        for (const auto& ring : outerLoops) {
            for (std::size_t i = 0; i < ring.size(); ++i) {
                insertConstraint(ring[i], ring[(i + 1) % ring.size()]);
            }
        }
        for (const auto& ring : holeLoops) {
            for (std::size_t i = 0; i < ring.size(); ++i) {
                insertConstraint(ring[i], ring[(i + 1) % ring.size()]);
            }
        }
        for (const auto& s : constraintSegments) {
            const VertexIndex a = vid.at(PointType(s[0]));
            const VertexIndex b = vid.at(PointType(s[1]));
            if (a != b) {
                insertConstraint(a, b);
            }
        }
        restoreConstrainedDelaunay();
 
        // One seed per hole, taken from a directed ring edge: a hole ring is
        // counterclockwise, so the triangle on the left of any of its edges lies
        // inside it. Recorded as a triangle value too — that is what tells a
        // region query later whether it has swallowed a hole, and unlike a
        // triangle id it survives the edits that follow.
        std::vector<TriIndex> holeSeeds;
        holeSeeds.reserve(holeLoops.size());
        for (const auto& ring : holeLoops) {
            const TriIndex seed = triangleLeftOf(ring[0], ring[1 % ring.size()]);
            if (seed != NO_TRI && !isGhost(seed)) {
                holeSeeds.push_back(seed);
                holeWitnesses_.push_back(triangleValue(seed));
            }
        }
        markOutOfDomain(holeSeeds);
 
        // Carry each constraint segment's label onto its edge record. Constrained
        // edges are never flipped, so they are still in segToEdge_ (keyed by
        // coordinates; the label is ignored for the lookup). Mirrors the
        // segment-range constructor.
        if constexpr (detail::has_label_v<SegmentLabel>) {
            for (const auto& s : constraintSegments) {
                auto it = segmentMap().find(SegmentType(PointType(s[0]), PointType(s[1])));
                if (it != segmentMap().end()) {
                    it->second.segLabel = detail::copyLabel<SegmentLabel>(s);
                }
            }
        }
    }
 
    // Prepends the ghost vertex at index 0, materializes the triangle records
    // (normalized to CCW) from vertex-index triples — 0-based into the
    // pre-ghost vertices_, so shifted by one here — and their parallel labels,
    // then links adjacency and the edge map. With the ghost first, the real
    // vertices are the contiguous range [1, vertices_.size()) and insertions
    // can append real vertices without disturbing it.
    void buildFromTriples(std::vector<std::array<VertexIndex, 3>>& triples,
                          const std::vector<TriangleLabel>& triLabels) {
        vertices_.insert(vertices_.begin(), PointType{});  // ghost (GHOST); coordinates unused
        syncVertexApproximations();
 
        for (std::size_t k = 0; k < triples.size(); ++k) {
            VertexIndex x = triples[k][0] + 1, y = triples[k][1] + 1, z = triples[k][2] + 1;
            if (orientationSign(vertices_[x], vertices_[y], vertices_[z]) < 0) {
                std::swap(y, z);
            }
            assert(orientationSign(vertices_[x], vertices_[y], vertices_[z]) > 0 &&
                   "Triangulation: degenerate triangle");
            triangles_.push_back(Tri{{x, y, z}, {NO_TRI, NO_TRI, NO_TRI}, 0, 0, 0, triLabels[k]});
        }
        firstGhost_ = static_cast<TriIndex>(triangles_.size());
        domainTriangleCount_ = static_cast<std::size_t>(firstGhost_);
        buildAdjacency();
        for (TriIndex t = 0; t < firstGhost_; ++t) {
            noteVertexIncidence(t);
        }
        mapStale_ = true;  // materialized on the first lookup that needs it
        assert(checkInvariants());
    }
 
    // ---- internal predicates / mutation ----------------------------------
 
    // True if edge e can be flipped: unconstrained, interior (both sides real),
    // and the two incident triangles form a strictly convex quadrilateral.
    [[nodiscard]] bool flippableEdge(Edge e) const {
        const TriIndex t = e.tri;
        if (t == NO_TRI || bit(triangles_[t].constrainedMask, e.side)) {
            return false;
        }
        const TriIndex t2 = triangles_[t].nbr[e.side];
        if (t2 == NO_TRI || isGhost(t) || isGhost(t2)) {
            return false;
        }
        const Edge m = mirror(e);
        const VertexIndex c = triangles_[t].v[e.side];
        const VertexIndex a = triangles_[t].v[(e.side + 1) % 3];
        const VertexIndex b = triangles_[t].v[(e.side + 2) % 3];
        const VertexIndex d = triangles_[t2].v[m.side];
        const auto fc = filteredVertex(c);
        const auto fd = filteredVertex(d);
        const auto oa = detail::orientationSignOf(fc, fd, filteredVertex(a)).value();
        const auto ob = detail::orientationSignOf(fc, fd, filteredVertex(b)).value();
        return (oa > 0 && ob < 0) || (oa < 0 && ob > 0);  // strictly convex quad
    }
 
    // Replaces edge e by the opposite diagonal of its quadrilateral, rewriting
    // the two incident triangle records and relinking the four outer neighbors.
    // Surrounding constrained flags are carried over; the new diagonal is
    // unconstrained. Does not touch segToEdge_ (callers re-register). Returns
    // false if e is not flippable.
    bool flipEdge(Edge e) {
        if (!flippableEdge(e)) {
            return false;
        }
        const TriIndex t = e.tri;
        const int i = e.side;
        const Edge m = mirror(e);
        const TriIndex t2 = m.tri;
        const int j = m.side;
 
        const VertexIndex c = triangles_[t].v[i];
        const VertexIndex a = triangles_[t].v[(i + 1) % 3];
        const VertexIndex b = triangles_[t].v[(i + 2) % 3];
        const VertexIndex d = triangles_[t2].v[j];
 
        const TriIndex nCA = triangles_[t].nbr[(i + 2) % 3];
        const TriIndex nBC = triangles_[t].nbr[(i + 1) % 3];
        const TriIndex nDB = triangles_[t2].nbr[(j + 2) % 3];
        const TriIndex nAD = triangles_[t2].nbr[(j + 1) % 3];
        const bool cCA = bit(triangles_[t].constrainedMask, (i + 2) % 3);
        const bool cBC = bit(triangles_[t].constrainedMask, (i + 1) % 3);
        const bool cDB = bit(triangles_[t2].constrainedMask, (j + 2) % 3);
        const bool cAD = bit(triangles_[t2].constrainedMask, (j + 1) % 3);
 
        const int sAD = (nAD != NO_TRI) ? findSide(nAD, t2) : -1;
        const int sBC = (nBC != NO_TRI) ? findSide(nBC, t) : -1;
 
        triangles_[t].v = {c, a, d};
        triangles_[t].nbr = {nAD, t2, nCA};
        triangles_[t].constrainedMask = mask(cAD, false, cCA);
        triangles_[t].triLabel = TriangleLabel{};  // flipped-in triangle has no source label
 
        triangles_[t2].v = {c, d, b};
        triangles_[t2].nbr = {nDB, nBC, t};
        triangles_[t2].constrainedMask = mask(cDB, cBC, false);
        triangles_[t2].triLabel = TriangleLabel{};
 
        if (sAD >= 0) {
            triangles_[nAD].nbr[sAD] = t;
        }
        if (sBC >= 0) {
            triangles_[nBC].nbr[sBC] = t2;
        }
        // assert(checkInvariants());  // O(n): uncomment when debugging
        return true;
    }
 
    // ---- point insertion internals ----------------------------------------
 
    // Re-registers the three sides of real triangle t in segToEdge_, keeping
    // any label already stored under a persisting edge key. (registerSides
    // would reset those labels — acceptable for flip, whose contract says so,
    // but an insertion only re-points surviving edges and must not wipe them.)
    void reRegisterSides(TriIndex t) {
        for (std::int8_t s = 0; s < 3; ++s) {
            const SegmentType key = edgeSegment(Edge{t, s});
            auto it = segmentMap().find(key);
            if (it == segmentMap().end()) {
                segmentMap().emplace(key, Edge{t, s});
            } else {
                it->second.tri = t;
                it->second.side = s;
            }
        }
    }
 
    // Ensures room for `extra` more elements without reallocation, growing
    // geometrically when needed. (A bare reserve(size + k) reallocates to
    // exactly that capacity, which would make repeated insertions quadratic.)
    template <class T>
    static void reserveExtra(std::vector<T>& v, std::size_t extra) {
        if (v.capacity() < v.size() + extra) {
            v.reserve(std::max(v.size() + extra, v.capacity() * 2));
        }
    }
 
    // Frees k slots directly above the real block — for the real triangles an
    // insertion creates — by relocating the k lowest ghost triangles to the
    // end of triangles_ and bumping firstGhost_. Only neighbor links reference
    // ghost triangles (segToEdge_ and hint_ never do), so relocation rewires
    // those and nothing else. Requires k ghosts (any nonempty triangulation
    // has >= 3) and spare capacity (callers reserve). Returns the first freed
    // slot; the freed slots hold stale copies the caller must overwrite.
    TriIndex makeRoom(int k) {
        const TriIndex base = firstGhost_;
        const TriIndex oldSize = static_cast<TriIndex>(triangles_.size());
        assert(oldSize - base >= k);
        for (int j = 0; j < k; ++j) {
            triangles_.push_back(triangles_[base + j]);
        }
        for (int j = 0; j < k; ++j) {
            const TriIndex moved = oldSize + j;
            for (int s = 0; s < 3; ++s) {
                TriIndex& nb = triangles_[moved].nbr[s];
                assert(nb != NO_TRI);  // ghosts always have three neighbors
                if (nb >= base && nb < base + k) {
                    nb = nb - base + oldSize;  // the neighbor was relocated too
                }
                for (int q = 0; q < 3; ++q) {
                    if (triangles_[nb].nbr[q] == base + j) {
                        triangles_[nb].nbr[q] = moved;
                    }
                }
            }
        }
        firstGhost_ += k;
        return base;
    }
 
    // Structural vertex insertion shared by insert and insertDelaunay: locates
    // @p p and subdivides the triangle (1->3) or edge (2->4; 2->3 on the hull,
    // where the ghost across splits too) containing it, or grows the hull when
    // p is outside. Returns the new vertex and one real triangle incident to
    // it, or nullopt — with the triangulation unchanged — when p is already a
    // vertex or the triangulation is empty. Splits inherit their parents'
    // out-of-domain flags, so an insertion inside a polygon's domain keeps the
    // carved-away region carved (a point outside the closed polygon is a
    // precondition violation; see insert).
    std::optional<std::pair<VertexIndex, TriIndex>> insertVertexImpl(const PointType& p) {
        const TriIndex t0 = locateIndex(p);
        if (t0 == NO_TRI) {
            return std::nullopt;  // empty triangulation
        }
        if (isGhost(t0)) {
            // p is strictly outside the hull (a point on the hull boundary
            // lands in a real triangle's closure), so it cannot be a duplicate.
            return growHull(t0, p);
        }
        const auto& tv = triangles_[t0].v;
        for (int k = 0; k < 3; ++k) {
            if (vertices_[tv[k]] == p) {
                return std::nullopt;  // p is already a vertex
            }
        }
        // p is in t0's closure (locateIndex stopped here) and is not a vertex, so
        // it lies strictly inside either the triangle or exactly one edge.
        int onSide = -1;
        const auto fp = filteredPoint(p);
        for (int k = 0; k < 3; ++k) {
            if (detail::orientationSignOf(filteredVertex(tv[(k + 1) % 3]),
                                          filteredVertex(tv[(k + 2) % 3]), fp).value() == 0) {
                onSide = k;
            }
        }
        if (onSide < 0) {
            return splitTriangle(t0, p);
        }
        return splitEdge(t0, onSide, p);
    }
 
    // 1->3 subdivision: replaces the in-domain triangle t, whose interior
    // strictly contains p, by the fan from its three edges to the new vertex.
    std::pair<VertexIndex, TriIndex> splitTriangle(TriIndex t, const PointType& p) {
        reserveExtra(vertices_, 1);
        reserveExtra(triangles_, 2);
        const VertexIndex vp = appendVertex(p);
        const TriIndex n1 = makeRoom(2);
        const TriIndex n2 = n1 + 1;
 
        // Copy the record after makeRoom, whose relocation already fixed the
        // neighbor links of every triangle adjacent to a moved ghost.
        const Tri old = triangles_[t];
        const VertexIndex a = old.v[0];
        const VertexIndex b = old.v[1];
        const VertexIndex c = old.v[2];
        // Children (CCW because p is strictly interior); each keeps one parent
        // edge — with that edge's constrained flag — as its side 2, opposite
        // vp, and inherits the parent's out-of-domain flag.
        triangles_[t] = Tri{{a, b, vp},
                            {n1, n2, old.nbr[2]},
                            mask(false, false, bit(old.constrainedMask, 2)),
                            old.outOfDomain, 0, TriangleLabel{}};
        triangles_[n1] = Tri{{b, c, vp},
                             {n2, t, old.nbr[0]},
                             mask(false, false, bit(old.constrainedMask, 0)),
                             old.outOfDomain, 0, TriangleLabel{}};
        triangles_[n2] = Tri{{c, a, vp},
                             {t, n1, old.nbr[1]},
                             mask(false, false, bit(old.constrainedMask, 1)),
                             old.outOfDomain, 0, TriangleLabel{}};
        triangles_[old.nbr[0]].nbr[findSide(old.nbr[0], t)] = n1;
        triangles_[old.nbr[1]].nbr[findSide(old.nbr[1], t)] = n2;
        domainTriangleCount_ += old.outOfDomain ? 0 : 2;
        for (const TriIndex x : {t, n1, n2}) {
            reRegisterSides(x);
        }
        hint_ = t;
        // assert(checkInvariants() && checkEdgeMap());  // O(n): uncomment when debugging
        return {vp, t};
    }
 
    // 2->4 (interior edge) or 2->3 (hull edge, where the ghost across splits
    // too) subdivision: p lies strictly inside the edge of triangle t opposite
    // t.v[s]. The edge splits into two collinear halves, which inherit its
    // constrained flag and label; each incident triangle splits in two, its
    // children inheriting its out-of-domain flag.
    std::pair<VertexIndex, TriIndex> splitEdge(TriIndex t, int s, const PointType& p) {
        const bool ghostSide = isGhost(triangles_[t].nbr[s]);
        reserveExtra(vertices_, 1);
        reserveExtra(triangles_, 2);
 
        const bool cUW = bit(triangles_[t].constrainedMask, s);
        const VertexIndex vp = appendVertex(p);
        const TriIndex n1 = makeRoom(ghostSide ? 1 : 2);
 
        // Read the records after makeRoom: it relocates the lowest ghosts, so a
        // ghost id read before it is stale (and now names a freed real slot),
        // while its relocation already fixed the neighbor links of every
        // triangle adjacent to a moved ghost.
        const TriIndex across = triangles_[t].nbr[s];
        const Tri oldT = triangles_[t];
        const VertexIndex apex = oldT.v[s];
        const VertexIndex u = oldT.v[(s + 1) % 3];
        const VertexIndex w = oldT.v[(s + 2) % 3];
        const TriIndex nWApex = oldT.nbr[(s + 1) % 3];  // neighbor across {w, apex}
        const TriIndex nApexU = oldT.nbr[(s + 2) % 3];  // neighbor across {apex, u}
        const bool cWApex = bit(oldT.constrainedMask, (s + 1) % 3);
        const bool cApexU = bit(oldT.constrainedMask, (s + 2) % 3);
 
        // Both halves inherit the split edge's label; the old key is dropped
        // now and the halves are registered (then labeled) below.
        SegmentLabel halfLabel{};
        {
            const auto it = segmentMap().find(SegmentType(vertices_[u], vertices_[w]));
            assert(it != segmentMap().end());
            halfLabel = it->second.segLabel;
            segmentMap().erase(it);
        }
 
        if (!ghostSide) {
            // Interior edge: t = (apex,u,w) and t2 = (apex2,w,u) become the
            // four triangles fanning around vp.
            const TriIndex n2 = n1 + 1;
            const TriIndex t2 = across;
            const int j = findSide(t2, t);
            const Tri oldT2 = triangles_[t2];
            const VertexIndex apex2 = oldT2.v[j];
            assert(oldT2.v[(j + 1) % 3] == w && oldT2.v[(j + 2) % 3] == u);
            const TriIndex nApex2W = oldT2.nbr[(j + 2) % 3];  // across {apex2, w}
            const TriIndex nUApex2 = oldT2.nbr[(j + 1) % 3];  // across {u, apex2}
            const bool cApex2W = bit(oldT2.constrainedMask, (j + 2) % 3);
            const bool cUApex2 = bit(oldT2.constrainedMask, (j + 1) % 3);
 
            triangles_[t] = Tri{{apex, u, vp},
                                {n2, n1, nApexU},
                                mask(cUW, false, cApexU),
                                oldT.outOfDomain, 0, TriangleLabel{}};
            triangles_[n1] = Tri{{apex, vp, w},
                                 {t2, nWApex, t},
                                 mask(cUW, cWApex, false),
                                 oldT.outOfDomain, 0, TriangleLabel{}};
            triangles_[t2] = Tri{{apex2, w, vp},
                                 {n1, n2, nApex2W},
                                 mask(cUW, false, cApex2W),
                                 oldT2.outOfDomain, 0, TriangleLabel{}};
            triangles_[n2] = Tri{{apex2, vp, u},
                                 {t, nUApex2, t2},
                                 mask(cUW, cUApex2, false),
                                 oldT2.outOfDomain, 0, TriangleLabel{}};
            triangles_[nWApex].nbr[findSide(nWApex, t)] = n1;
            triangles_[nUApex2].nbr[findSide(nUApex2, t2)] = n2;
            domainTriangleCount_ += (oldT.outOfDomain ? 0 : 1) + (oldT2.outOfDomain ? 0 : 1);
            for (const TriIndex x : {t, n1, t2, n2}) {
                reRegisterSides(x);
            }
        } else {
            // Hull edge: t splits in two and so does the ghost across, keeping
            // the ghost convention v = {real, real, GHOST} with nbr[2] real.
            const TriIndex g = across;
            const Tri oldG = triangles_[g];
            assert(oldG.v[0] == u && oldG.v[1] == w && oldG.v[2] == GHOST && oldG.nbr[2] == t);
            const TriIndex gw = oldG.nbr[0];  // ghost-ring neighbor across {w, ghost}
            const TriIndex gu = oldG.nbr[1];  // ghost-ring neighbor across {u, ghost}
 
            const TriIndex g2 = static_cast<TriIndex>(triangles_.size());
            triangles_.push_back(Tri{{vp, w, GHOST},
                                     {gw, g, n1},
                                     mask(false, false, cUW),
                                     0, 0, TriangleLabel{}});
            triangles_[t] = Tri{{apex, u, vp},
                                {g, n1, nApexU},
                                mask(cUW, false, cApexU),
                                oldT.outOfDomain, 0, TriangleLabel{}};
            triangles_[n1] = Tri{{apex, vp, w},
                                 {g2, nWApex, t},
                                 mask(cUW, cWApex, false),
                                 oldT.outOfDomain, 0, TriangleLabel{}};
            triangles_[g] = Tri{{u, vp, GHOST},
                                {g2, gu, t},
                                mask(false, false, cUW),
                                0, 0, TriangleLabel{}};
            triangles_[nWApex].nbr[findSide(nWApex, t)] = n1;
            triangles_[gw].nbr[findSide(gw, g)] = g2;
            domainTriangleCount_ += oldT.outOfDomain ? 0 : 1;
            for (const TriIndex x : {t, n1}) {
                reRegisterSides(x);
            }
        }
        if constexpr (detail::has_label_v<SegmentLabel>) {
            segmentMap().at(SegmentType(vertices_[u], vertices_[vp])).segLabel = halfLabel;
            segmentMap().at(SegmentType(vertices_[vp], vertices_[w])).segLabel = halfLabel;
        }
        hint_ = t;
        // assert(checkInvariants() && checkEdgeMap());  // O(n): uncomment when debugging
        return std::pair{vp, t};
    }
 
    // Hull growth: p lies strictly outside the convex hull, and the locate
    // walk stopped at ghost g0, whose base hull edge p strictly sees. Joins p
    // to the maximal contiguous chain of hull edges visible from it — the
    // fan is the only triangulation of the pocket between hull and point, as
    // any other diagonal would cut into the old hull — replacing the chain's
    // m ghosts by m real triangles, closed by two new ghosts along the new
    // hull edges {u_0, p} and {p, u_m}. Visibility is strict, so hull edges
    // collinear with p are never in the chain and no degenerate triangle can
    // arise. The new triangles are in-domain; a base hull edge that is
    // constrained stays constrained and simply becomes interior. (For a
    // polygon triangulation an outside point is a precondition violation —
    // see insert — so no domain fencing happens here.)
    std::pair<VertexIndex, TriIndex> growHull(TriIndex g0, const PointType& p) {
        const auto fp = filteredPoint(p);
        const auto visible = [&](TriIndex g) {
            const auto& gv = triangles_[g].v;
            return detail::orientationSignOf(filteredVertex(gv[0]), filteredVertex(gv[1]), fp)
                       .value() < 0;
        };
        assert(isGhost(g0) && visible(g0));
 
        // Rewind to the first visible ghost, then record the chain through its
        // inner anchors (r, s): ghost ids go stale across makeRoom, the inner
        // real triangles do not. For a ghost {a, b, GHOST}, nbr[1] is the ring
        // predecessor (sharing a) and nbr[0] the successor (sharing b).
        std::size_t guard = 0;
        const std::size_t ringCap = triangles_.size() + 1;
        TriIndex gStart = g0;
        while (visible(triangles_[gStart].nbr[1]) && ++guard < ringCap) {
            gStart = triangles_[gStart].nbr[1];
        }
        std::vector<TriIndex> inner;            // r_i: the real triangle inside base i
        std::vector<std::int8_t> innerSide;  // its side facing that base
        std::vector<VertexIndex> u;             // hull chain u_0 -> ... -> u_m
        TriIndex g = gStart;
        guard = 0;
        while (visible(g) && ++guard < ringCap) {
            const TriIndex r = triangles_[g].nbr[2];
            inner.push_back(r);
            innerSide.push_back(findSide(r, g));
            u.push_back(triangles_[g].v[0]);
            g = triangles_[g].nbr[0];
        }
        u.push_back(triangles_[g].v[0]);  // == v[1] of the last chain ghost
        const int m = static_cast<int>(inner.size());
        assert(m >= 1);
 
        reserveExtra(vertices_, 1);
        reserveExtra(triangles_, m == 1 ? 2 : static_cast<std::size_t>(m));
        const VertexIndex vp = appendVertex(p);
        const TriIndex nr = makeRoom(m);  // slots for the m new real triangles
 
        // Re-resolve the (possibly relocated) chain ghosts and the ring ends.
        std::vector<TriIndex> dead(static_cast<std::size_t>(m));
        for (int i = 0; i < m; ++i) {
            dead[i] = triangles_[inner[i]].nbr[innerSide[i]];
        }
        const TriIndex prevG = triangles_[dead.front()].nbr[1];
        const TriIndex nextG = triangles_[dead.back()].nbr[0];
 
        // The two new ghosts reuse dead chain slots (plus a push_back when
        // only one ghost died); leftover dead slots are compacted away below,
        // so the vector grows by exactly two triangles in every case.
        const TriIndex ga = dead[0];
        TriIndex gb;
        if (m == 1) {
            gb = static_cast<TriIndex>(triangles_.size());
            triangles_.push_back(Tri{});  // written below
        } else {
            gb = dead[1];
        }
 
        // The m pocket triangles {u_{i+1}, u_i, vp} fanning vp: side 2 is the
        // base hull edge (keeping its constrained flag), sides 0/1 the spokes.
        for (int i = 0; i < m; ++i) {
            const bool cBase = bit(triangles_[inner[i]].constrainedMask, innerSide[i]);
            triangles_[nr + i] = Tri{{u[i + 1], u[i], vp},
                                     {i > 0 ? nr + i - 1 : ga,
                                      i < m - 1 ? nr + i + 1 : gb, inner[i]},
                                     mask(false, false, cBase),
                                     0, 0, TriangleLabel{}};
            triangles_[inner[i]].nbr[innerSide[i]] = nr + i;
        }
        triangles_[ga] = Tri{{u.front(), vp, GHOST}, {gb, prevG, nr}, 0, 0, 0, TriangleLabel{}};
        triangles_[gb] =
            Tri{{vp, u.back(), GHOST}, {nextG, ga, nr + m - 1}, 0, 0, 0, TriangleLabel{}};
        triangles_[prevG].nbr[0] = ga;
        triangles_[nextG].nbr[1] = gb;
 
        if (m >= 3) {
            // m ghosts died but only two slots were reused: fill the remaining
            // holes with live ghosts taken from the end of triangles_ (the
            // same link rewiring as makeRoom, in the other direction), then
            // drop the all-dead tail.
            std::vector<TriIndex> holes(dead.begin() + 2, dead.end());
            std::sort(holes.begin(), holes.end());
            const auto isHole = [&](TriIndex t) {
                return std::binary_search(holes.begin(), holes.end(), t);
            };
            TriIndex last = static_cast<TriIndex>(triangles_.size()) - 1;
            for (std::size_t h = 0; h < holes.size() && holes[h] < last;) {
                if (isHole(last)) {
                    --last;  // already dead: it will be truncated
                    continue;
                }
                const TriIndex hole = holes[h];
                triangles_[hole] = triangles_[last];
                for (int s = 0; s < 3; ++s) {
                    const TriIndex nb = triangles_[hole].nbr[s];
                    for (int q = 0; q < 3; ++q) {
                        if (triangles_[nb].nbr[q] == last) {
                            triangles_[nb].nbr[q] = hole;
                        }
                    }
                }
                ++h;
                --last;
            }
            triangles_.resize(triangles_.size() - static_cast<std::size_t>(m - 2));
        }
 
        for (int i = 0; i < m; ++i) {
            reRegisterSides(nr + i);
        }
        domainTriangleCount_ += static_cast<std::size_t>(m);
        hint_ = nr;
        // assert(checkInvariants() && checkEdgeMap());  // O(n): uncomment when debugging
        return {vp, nr};
    }
 
    // Restores the local Delaunay property after an insertion by Lawson flips:
    // pops suspect edges, flips any that are flippable (so never a constrained
    // edge) with an apex strictly inside the opposite circumcircle, and
    // re-queues the rewritten triangles' sides. Terminates by Lawson's
    // argument; the guard is a safety net only.
    void legalize(std::vector<SegmentType>& suspect) {
        std::size_t guard = 0;
        const std::size_t cap = triangles_.size() * triangles_.size() + 64;
        while (!suspect.empty() && ++guard < cap) {
            const SegmentType s = suspect.back();
            suspect.pop_back();
            const auto se = segmentMap().find(s);
            if (se == segmentMap().end() || !flippableEdge(se->second)) {
                continue;  // edge gone, constrained, or quad not convex
            }
            const Edge e = se->second;
            const Edge m = mirror(e);
            const auto& tv = triangles_[e.tri].v;
            const VertexIndex d = triangles_[m.tri].v[m.side];
            if (detail::inCircleSignOf(filteredVertex(tv[0]), filteredVertex(tv[1]),
                                       filteredVertex(tv[2]), filteredVertex(d)) !=
                std::partial_ordering::greater) {
                continue;  // locally Delaunay (the test is symmetric across e)
            }
            const TriIndex t = e.tri;
            const TriIndex t2 = m.tri;
            flip(s);
            // The rewritten triangles' sides — the four quad edges plus the new
            // diagonal, which the test above now accepts — are suspect again.
            for (const TriIndex x : {t, t2}) {
                for (std::int8_t q = 0; q < 3; ++q) {
                    suspect.push_back(edgeSegment(Edge{x, q}));
                }
            }
        }
        assert(suspect.empty() && "Triangulation: legalization did not terminate");
    }
 
    // Debug validation of segToEdge_: every real triangle side is registered
    // under its edge key, and every entry references a real triangle that
    // still owns the keyed edge (no stale handles or leftover keys). Kept
    // separate from checkInvariants because flipEdge checks invariants at a
    // point where the map is deliberately stale (its callers re-register).
    [[nodiscard]] bool checkEdgeMap() const {
        for (TriIndex t = 0; t < firstGhost_; ++t) {
            for (std::int8_t s = 0; s < 3; ++s) {
                if (!segmentMap().contains(edgeSegment(Edge{t, s}))) {
                    return false;
                }
            }
        }
        for (const auto& [seg, e] : segmentMap()) {
            if (e.tri < 0 || e.tri >= firstGhost_) {
                return false;
            }
            const auto& tv = triangles_[e.tri].v;
            if (SegmentType(vertices_[tv[(e.side + 1) % 3]],
                            vertices_[tv[(e.side + 2) % 3]]) != seg) {
                return false;
            }
        }
        return true;
    }
 
    // Locates the triangle containing p by a stochastic visibility walk (random
    // start edge guarantees termination). Returns a ghost triangle if p is
    // outside the triangulated region, NO_TRI only if empty. Starts at @p start
    // when that is a real triangle, else seeds from hint_; updates hint_ either
    // way. @p p may use a different point type.
    template <class QueryPoint>
    [[nodiscard]] TriIndex locateIndex(const QueryPoint& p, TriIndex start = NO_TRI) const {
        if (triangles_.empty()) {
            return NO_TRI;
        }
        TriIndex t = realTriangle(start) ? start
                                         : ((hint_ != NO_TRI && !isGhost(hint_)) ? hint_ : 0);
        TriIndex from = NO_TRI;
        const std::size_t cap = triangles_.size() * 3 + 16;
        // The walk tests p against three edges of every triangle it steps
        // through, so p is converted once here rather than once per test.
        const auto fq = filteredPoint(p);
        for (std::size_t step = 0; step < cap; ++step) {
            if (isGhost(t)) {
                hint_ = NO_TRI;
                return t;  // outside the triangulated region
            }
            const auto& T = triangles_[t];
            const int begin = static_cast<int>(rng_() % 3);  // random start: provably terminates
            TriIndex next = NO_TRI;
            for (int k = 0; k < 3; ++k) {
                const int s = (begin + k) % 3;
                if (T.nbr[s] == from) {
                    continue;
                }
                const VertexIndex ea = T.v[(s + 1) % 3];
                const VertexIndex eb = T.v[(s + 2) % 3];
                if (detail::orientationSignOf(filteredVertex(ea), filteredVertex(eb), fq)
                        .value() < 0) {
                    next = T.nbr[s];
                    break;
                }
            }
            if (next == NO_TRI) {
                hint_ = t;
                return t;
            }
            from = t;
            t = next;
        }
        return t;
    }
 
    // Links neighbor pointers between real triangles sharing an edge, then closes
    // every still-unmatched (boundary) edge with a ghost triangle to the ghost
    // vertex and links those ghosts into a ring, so every edge has two sides.
    void buildAdjacency() {
        // Undirected edge key (endpoints sorted) for matching the two sides.
        const auto key = [](VertexIndex u, VertexIndex w) {
            return u < w ? std::pair<VertexIndex, VertexIndex>{u, w} : std::pair<VertexIndex, VertexIndex>{w, u};
        };
 
        // Every side of every triangle, sorted by its edge key, so the two
        // sides of an interior edge land next to each other and a boundary
        // edge's single side lands alone. Sorting an array beats matching them
        // through a search tree, which allocates a node per edge and revisits
        // it to erase.
        struct Side {
            std::pair<VertexIndex, VertexIndex> edge;
            TriIndex tri;
            int side;
        };
        std::vector<Side> sides;
        sides.reserve(static_cast<std::size_t>(std::max(firstGhost_, TriIndex(0))) * 3);
        for (TriIndex t = 0; t < firstGhost_; ++t) {
            for (int i = 0; i < 3; ++i) {
                sides.push_back(Side{key(triangles_[t].v[(i + 1) % 3],
                                         triangles_[t].v[(i + 2) % 3]),
                                     t, i});
            }
        }
        std::sort(sides.begin(), sides.end(),
                  [](const Side& a, const Side& b) { return a.edge < b.edge; });
 
        // The still-unmatched sides, which are the boundary ones.
        std::vector<Side> edges;
        for (std::size_t i = 0; i < sides.size();) {
            if (i + 1 < sides.size() && sides[i].edge == sides[i + 1].edge) {
                const Side& one = sides[i];
                const Side& other = sides[i + 1];
                triangles_[one.tri].nbr[one.side] = other.tri;
                triangles_[other.tri].nbr[other.side] = one.tri;
                i += 2;
            } else {
                edges.push_back(sides[i]);
                ++i;
            }
        }
 
        std::map<VertexIndex, std::pair<TriIndex, int>> ghostEdges;
        for (const Side& unmatched : edges) {
            const TriIndex t = unmatched.tri;
            const int i = unmatched.side;
            const VertexIndex a = triangles_[t].v[(i + 1) % 3];
            const VertexIndex b = triangles_[t].v[(i + 2) % 3];
            const TriIndex g = static_cast<TriIndex>(triangles_.size());
            triangles_.push_back(Tri{{a, b, GHOST}, {NO_TRI, NO_TRI, NO_TRI}, 0});
            triangles_[g].nbr[2] = t;  // side 2 (opposite ghost) is the shared edge {a,b}
            triangles_[t].nbr[i] = g;
            for (auto [realVertex, side] :
                 {std::pair<VertexIndex, int>{b, 0}, std::pair<VertexIndex, int>{a, 1}}) {
                auto it = ghostEdges.find(realVertex);
                if (it == ghostEdges.end()) {
                    ghostEdges.emplace(realVertex, std::pair<TriIndex, int>{g, side});
                } else {
                    auto [g2, s2] = it->second;
                    triangles_[g].nbr[side] = g2;
                    triangles_[g2].nbr[s2] = g;
                    ghostEdges.erase(it);
                }
            }
        }
        assert(ghostEdges.empty() && "Triangulation: open boundary (input is not a triangulation)");
    }
};
 
// The Kirkpatrick hierarchy. Level 0 is the mesh, extended by a ring of
// triangles filling an enclosing box, so that every vertex but the four box
// corners is interior and can be removed. Each further level removes an
// independent set of low-degree vertices and retriangulates the hole each star
// leaves; the new triangles record the ones they cover, which is what a query
// descends. The four corners survive every level, so the top is the box itself
// as a couple of triangles.
template <TriangleConcept TriangleType, SegmentConcept SegmentType>
void Triangulation<TriangleType, SegmentType>::buildPointLocation() {
    if (hasCurrentPointLocation()) {
        return;  // already drawn against this mesh; redrawing would find nothing
    }
    pointLocation_.reset();
    const std::size_t meshCells = static_cast<std::size_t>(firstGhost_);
    if (meshCells == 0 || triangles_.size() <= meshCells) {
        return;  // no triangle to index, or no ghost ring to read the hull off
    }
    auto kp = std::make_shared<Kirkpatrick>();
    using Cell = typename Kirkpatrick::Cell;
    // The triangle of every cell. A finished hierarchy keeps only the top
    // level's, a query below it being handed down rather than placed, but the
    // construction reads all of them.
    std::vector<std::array<VertexIndex, 3>> shape;
 
    // ---- level 0, first half: the mesh triangles, at their own indices ----
    kp->cells.resize(meshCells);
    shape.resize(meshCells);
    for (std::size_t t = 0; t < meshCells; ++t) {
        shape[t] = triangles_[t].v;
        kp->cells[t].runBegin = static_cast<std::uint32_t>(t);
    }
 
    // ---- the hull, read off the ghost ring -------------------------------
    // ring[i] -> ring[i + 1] is a boundary edge with the mesh on its left, and
    // ringTri[i] is the mesh triangle there.
    std::vector<VertexIndex> ring;
    std::vector<TriIndex> ringTri;
    {
        const TriIndex first = firstGhost_;
        TriIndex g = first;
        do {
            const Tri& ghost = triangles_[static_cast<std::size_t>(g)];
            ring.push_back(ghost.v[0]);
            ringTri.push_back(ghost.nbr[2]);
            g = ghost.nbr[0];  // the ghost across {ghost.v[1], GHOST}
        } while (g != first && ring.size() <= triangles_.size());
        if (g != first || ring.size() < 3) {
            return;
        }
    }
    const std::uint32_t hullSize = static_cast<std::uint32_t>(ring.size());
 
    // ---- the enclosing box -----------------------------------------------
    NumberType xlo = vertices_[static_cast<std::size_t>(ring[0])].x();
    NumberType xhi = xlo;
    NumberType ylo = vertices_[static_cast<std::size_t>(ring[0])].y();
    NumberType yhi = ylo;
    for (const VertexIndex v : ring) {
        const PointType& p = vertices_[static_cast<std::size_t>(v)];
        if (p.x() < xlo) xlo = p.x();
        if (xhi < p.x()) xhi = p.x();
        if (p.y() < ylo) ylo = p.y();
        if (yhi < p.y()) yhi = p.y();
    }
    if constexpr (std::numeric_limits<NumberType>::is_specialized &&
                  std::numeric_limits<NumberType>::is_integer &&
                  std::numeric_limits<NumberType>::is_bounded) {
        // Stepping the box out would wrap: leave the mesh unindexed rather than
        // fold the outside onto the inside.
        if (xlo == std::numeric_limits<NumberType>::lowest() ||
            ylo == std::numeric_limits<NumberType>::lowest() ||
            xhi == std::numeric_limits<NumberType>::max() ||
            yhi == std::numeric_limits<NumberType>::max()) {
            return;
        }
    }
    const NumberType unit(1);
    const NumberType bxlo = xlo - unit;
    const NumberType bylo = ylo - unit;
    const NumberType bxhi = xhi + unit;
    const NumberType byhi = yhi + unit;
    if (!(bxlo < xlo) || !(bylo < ylo) || !(xhi < bxhi) || !(yhi < byhi)) {
        return;  // coordinates too large to step away from (inexact types)
    }
    kp->extra = {PointType(bxlo, bylo), PointType(bxhi, bylo), PointType(bxhi, byhi),
                 PointType(bxlo, byhi)};
    if constexpr (detail::filtersSign<VertexCoordinate>) {
        kp->extraApprox.reserve(kp->extra.size());
        for (const PointType& corner : kp->extra) {
            kp->extraApprox.push_back(detail::approximatePoint(corner));
        }
    }
    const auto cornerVertex = [&](std::size_t j) {
        return static_cast<VertexIndex>(vertices_.size() + j);
    };
 
    // ---- the geometry the construction runs on ---------------------------
    const auto orient = [&](VertexIndex a, VertexIndex b, VertexIndex c) {
        return detail::orientationSignOf(kirkpatrickVertex(*kp, a), kirkpatrickVertex(*kp, b),
                                         kirkpatrickVertex(*kp, c))
            .value();
    };
    const auto holds = [&](VertexIndex a, VertexIndex b, VertexIndex c, VertexIndex p) {
        return !(orient(a, b, p) < 0) && !(orient(b, c, p) < 0) && !(orient(c, a, p) < 0);
    };
    const auto addLeaf = [&](VertexIndex a, VertexIndex b, VertexIndex c, TriIndex seed) {
        Cell cell;
        cell.runBegin = static_cast<std::uint32_t>(seed);
        cell.runCount = Kirkpatrick::FULL_TURN;  // fills the box; seeds, never answers
        kp->cells.push_back(cell);
        shape.push_back({a, b, c});
    };
 
    // Ear clipping of a simple counterclockwise polygon, reporting each ear as
    // the triple of *positions* it cut, in clipping order. Only a strictly
    // convex corner is cut, so every corner of the polygon survives in some
    // triangle: a collinear one dropped instead would leave the vertex on an
    // edge of the retriangulation without being a corner of it, which no later
    // level could then remove. @p mountain says the polygon is a monotone
    // mountain over its first edge: its two base corners are never cut, and
    // every convex corner is then an ear, which is what lets the pass skip
    // testing a corner against the rest of the polygon.
    std::vector<std::array<std::uint32_t, 3>> ears;
    std::vector<std::uint32_t> prevAt;
    std::vector<std::uint32_t> nextAt;
    std::vector<std::int8_t> turnAt;
    const auto earClip = [&](const std::vector<VertexIndex>& poly, bool mountain) {
        ears.clear();
        const std::uint32_t m = static_cast<std::uint32_t>(poly.size());
        if (m < 3) {
            return false;
        }
        prevAt.resize(m);
        nextAt.resize(m);
        turnAt.resize(m);
        for (std::uint32_t i = 0; i < m; ++i) {
            prevAt[i] = (i + m - 1) % m;
            nextAt[i] = (i + 1) % m;
        }
        const auto turn = [&](std::uint32_t i) -> std::int8_t {
            const auto side = orient(poly[prevAt[i]], poly[i], poly[nextAt[i]]);
            return side > 0 ? std::int8_t{1} : (side < 0 ? std::int8_t{-1} : std::int8_t{0});
        };
        std::uint32_t blockingCount = 0;  // corners that can sit inside an ear
        for (std::uint32_t i = 0; i < m; ++i) {
            turnAt[i] = turn(i);
            blockingCount += turnAt[i] <= 0 ? 1u : 0u;
        }
        std::uint32_t remaining = m;
        std::uint32_t cursor = mountain ? 2 : 0;
        std::uint32_t skipped = 0;
        while (remaining > 3) {
            if (skipped > remaining) {
                return false;  // a whole turn with no ear: not a simple polygon
            }
            const std::uint32_t i = cursor;
            cursor = nextAt[cursor];
            if (turnAt[i] <= 0 || (mountain && i < 2)) {
                ++skipped;
                continue;
            }
            const std::uint32_t p = prevAt[i];
            const std::uint32_t n = nextAt[i];
            if (!mountain && blockingCount != 0) {
                bool clean = true;
                for (std::uint32_t k = nextAt[n]; k != p; k = nextAt[k]) {
                    if (turnAt[k] <= 0 && holds(poly[p], poly[i], poly[n], poly[k])) {
                        clean = false;
                        break;
                    }
                }
                if (!clean) {
                    ++skipped;
                    continue;
                }
            }
            ears.push_back({p, i, n});
            nextAt[p] = n;
            prevAt[n] = p;
            --remaining;
            for (const std::uint32_t z : {p, n}) {
                const std::int8_t was = turnAt[z];
                turnAt[z] = turn(z);
                blockingCount += (turnAt[z] <= 0 ? 1u : 0u) - (was <= 0 ? 1u : 0u);
            }
            cursor = p;
            skipped = 0;
        }
        const std::uint32_t a = cursor;
        const std::uint32_t b = nextAt[a];
        ears.push_back({a, b, nextAt[b]});
        return true;
    };
 
    // ---- level 0, second half: the ring filling the box ------------------
    // The hull vertex extreme in a diagonal direction anchors a spoke to the box
    // corner facing it: the supporting line separates the two, so the spoke runs
    // outside the hull. The four spokes cut the ring into monotone mountains —
    // a box side as base, a monotone chain above it — which ear clipping
    // triangulates in one pass.
    std::array<std::uint32_t, 4> anchor{0, 0, 0, 0};
    {
        using Wide = detail::promoted_number_t<NumberType>;
        const auto wx = [&](std::uint32_t i) {
            return detail::asNumber<Wide>(vertices_[static_cast<std::size_t>(ring[i])].x());
        };
        const auto wy = [&](std::uint32_t i) {
            return detail::asNumber<Wide>(vertices_[static_cast<std::size_t>(ring[i])].y());
        };
        Wide leastSum = wx(0) + wy(0);
        Wide mostSum = leastSum;
        for (std::uint32_t i = 1; i < hullSize; ++i) {
            const Wide sum = wx(i) + wy(i);
            if (sum < leastSum) {
                leastSum = sum;
                anchor[0] = i;
            }
            if (mostSum < sum) {
                mostSum = sum;
                anchor[2] = i;
            }
            // The extremes of x - y, compared without forming a difference.
            if (wx(anchor[1]) + wy(i) < wx(i) + wy(anchor[1])) {
                anchor[1] = i;
            }
            if (wx(i) + wy(anchor[3]) < wx(anchor[3]) + wy(i)) {
                anchor[3] = i;
            }
        }
    }
    {
        std::vector<VertexIndex> poly;
        std::vector<TriIndex> edgeSeed;
        for (std::size_t j = 0; j < 4; ++j) {
            poly.clear();
            edgeSeed.clear();
            poly.push_back(cornerVertex(j));
            edgeSeed.push_back(NO_TRI);
            poly.push_back(cornerVertex((j + 1) % 4));
            edgeSeed.push_back(NO_TRI);
            // The hull walked backwards, from the next corner's anchor to this
            // one's: the filling runs counterclockwise where the hull runs
            // clockwise, since it lies on the hull's other side.
            for (std::uint32_t k = anchor[(j + 1) % 4];; k = (k + hullSize - 1) % hullSize) {
                poly.push_back(ring[k]);
                const bool last = k == anchor[j];
                edgeSeed.push_back(last ? NO_TRI : ringTri[(k + hullSize - 1) % hullSize]);
                if (last) {
                    break;
                }
            }
            if (!earClip(poly, /*mountain=*/true)) {
                return;
            }
            for (const auto& ear : ears) {
                // The seed a query landing in this triangle walks from: the mesh
                // triangle across whichever of its edges lies on the hull, which
                // the diagonal it leaves behind then carries to its neighbors.
                TriIndex seed = edgeSeed[ear[0]];
                if (seed == NO_TRI) seed = edgeSeed[ear[1]];
                if (seed == NO_TRI) seed = edgeSeed[ear[2]];
                addLeaf(poly[ear[0]], poly[ear[1]], poly[ear[2]], seed);
                edgeSeed[ear[0]] = seed;
            }
        }
    }
 
    // ---- the levels above ------------------------------------------------
    std::vector<std::uint32_t> active(kp->cells.size());
    for (std::uint32_t i = 0; i < static_cast<std::uint32_t>(active.size()); ++i) {
        active[i] = i;
    }
 
    std::vector<std::int32_t> slotOf(vertices_.size() + kp->extra.size(), -1);
    std::vector<VertexIndex> used;
    std::vector<std::uint32_t> starBegin;
    std::vector<std::uint32_t> starFill;
    std::vector<std::uint32_t> starCells;
    std::vector<VertexIndex> starFrom;
    std::vector<VertexIndex> starTo;
    std::vector<std::uint8_t> blocked;
    std::vector<std::uint32_t> order;
    std::vector<std::uint8_t> replaced;
    std::vector<std::uint32_t> fan;
    std::vector<VertexIndex> link;
    std::vector<std::uint32_t> fresh;
    std::vector<std::uint32_t> kept;
    std::size_t degreeLimit = pointLocationMaxStar;
 
    while (active.size() > pointLocationTopSize) {
        // The vertices this level still has, and the cells around each of them.
        used.clear();
        for (const std::uint32_t c : active) {
            for (const VertexIndex v : shape[c]) {
                if (slotOf[static_cast<std::size_t>(v)] < 0) {
                    slotOf[static_cast<std::size_t>(v)] = static_cast<std::int32_t>(used.size());
                    used.push_back(v);
                }
            }
        }
        const std::uint32_t vertexCount = static_cast<std::uint32_t>(used.size());
        starBegin.assign(vertexCount + 1, 0);
        for (const std::uint32_t c : active) {
            for (const VertexIndex v : shape[c]) {
                ++starBegin[static_cast<std::size_t>(slotOf[static_cast<std::size_t>(v)]) + 1];
            }
        }
        for (std::uint32_t s = 0; s < vertexCount; ++s) {
            starBegin[s + 1] += starBegin[s];
        }
        starFill = starBegin;
        starCells.resize(starBegin[vertexCount]);
        for (const std::uint32_t c : active) {
            for (const VertexIndex v : shape[c]) {
                starCells[starFill[static_cast<std::size_t>(
                    slotOf[static_cast<std::size_t>(v)])]++] = c;
            }
        }
 
        // Candidates, smallest star first: a small star makes few triangles and
        // blocks few other candidates, so taking those first leaves more room.
        // This, rather than either constant above, is what decides how far a
        // level gets — taking the candidates in storage order instead costs 16%
        // more cells, two more levels and 28% of the query time at 100,000
        // vertices, a maximal independent set being that much smaller when a
        // crowded vertex blocks its neighbours first.
        order.clear();
        {
            std::vector<std::uint32_t> tally(degreeLimit + 2, 0);
            for (std::uint32_t s = 0; s < vertexCount; ++s) {
                const std::size_t degree = starBegin[s + 1] - starBegin[s];
                if (degree >= 3 && degree <= degreeLimit) {
                    ++tally[degree + 1];
                }
            }
            for (std::size_t d = 1; d < tally.size(); ++d) {
                tally[d] += tally[d - 1];
            }
            order.resize(tally.back());
            for (std::uint32_t s = 0; s < vertexCount; ++s) {
                const std::size_t degree = starBegin[s + 1] - starBegin[s];
                if (degree >= 3 && degree <= degreeLimit) {
                    order[tally[degree]++] = s;
                }
            }
        }
 
        blocked.assign(vertexCount, 0);
        replaced.assign(kp->cells.size(), 0);
        fresh.clear();
        for (const std::uint32_t slot : order) {
            if (blocked[slot]) {
                continue;
            }
            const VertexIndex v = used[slot];
            const std::uint32_t begin = starBegin[slot];
            const std::uint32_t degree = starBegin[slot + 1] - begin;
 
            // Walk the star into a fan: cell m of it is (v, link[m], link[m+1]).
            // A vertex whose cells do not close into a single turn — a box
            // corner, or one a neighboring cell only touches — is not one this
            // can remove, and the walk finding no successor is how it says so.
            starFrom.resize(degree);
            starTo.resize(degree);
            bool closes = true;
            for (std::uint32_t k = 0; k < degree && closes; ++k) {
                const auto& corners = shape[starCells[begin + k]];
                std::uint32_t position = 0;
                while (position < 3 && corners[position] != v) {
                    ++position;
                }
                closes = position < 3;
                if (closes) {
                    starFrom[k] = corners[(position + 1) % 3];
                    starTo[k] = corners[(position + 2) % 3];
                }
            }
            fan.clear();
            link.clear();
            for (std::uint32_t step = 0, at = 0; closes && step < degree; ++step) {
                fan.push_back(starCells[begin + at]);
                link.push_back(starFrom[at]);
                const VertexIndex after = starTo[at];
                if (step + 1 == degree) {
                    closes = after == link[0];
                    break;
                }
                std::uint32_t following = degree;
                for (std::uint32_t k = 0; k < degree; ++k) {
                    if (starFrom[k] == after) {
                        following = k;
                        break;
                    }
                }
                closes = following < degree;
                at = closes ? following : 0;
            }
            if (!closes || link.size() != degree || !earClip(link, /*mountain=*/false)) {
                blocked[slot] = 1;  // not removable; its neighbors stay free
                continue;
            }
 
            for (const auto& ear : ears) {
                // The fan cells this triangle covers. Consecutive link vertices
                // turn around v, so a triangle that does not hold v spans one
                // arc of them — the two gaps other than the one wider than half
                // a turn — and it covers exactly the cells in that arc.
                std::uint32_t from = 0;
                std::uint32_t count = degree;
                int outward = -1;
                bool decided = true;
                for (int g = 0; g < 3 && decided; ++g) {
                    const auto side =
                        orient(v, link[ear[static_cast<std::size_t>(g)]],
                               link[ear[static_cast<std::size_t>((g + 1) % 3)]]);
                    if (side == 0) {
                        decided = false;
                    } else if (side < 0) {
                        decided = outward < 0;
                        outward = g;
                    }
                }
                if (decided && outward >= 0) {
                    from = ear[static_cast<std::size_t>((outward + 1) % 3)];
                    const std::uint32_t stop = ear[static_cast<std::size_t>(outward)];
                    count = (stop + degree - from) % degree;
                    if (count == 0) {
                        count = degree;
                    }
                }
                Cell cell;
                cell.apex = static_cast<std::uint32_t>(v);
                cell.runBegin = static_cast<std::uint32_t>(kp->run.size());
                cell.runCount = count | (count == degree ? Kirkpatrick::FULL_TURN : 0u);
                for (std::uint32_t k = 0; k < count; ++k) {
                    kp->run.push_back(static_cast<std::uint32_t>(link[(from + k) % degree]));
                    kp->run.push_back(fan[(from + k) % degree]);
                }
                kp->run.push_back(static_cast<std::uint32_t>(link[(from + count) % degree]));
                fresh.push_back(static_cast<std::uint32_t>(kp->cells.size()));
                kp->cells.push_back(cell);
                shape.push_back({link[ear[0]], link[ear[1]], link[ear[2]]});
            }
 
            blocked[slot] = 1;
            for (const VertexIndex u : link) {
                blocked[static_cast<std::size_t>(slotOf[static_cast<std::size_t>(u)])] = 1;
            }
            for (const std::uint32_t cell : fan) {
                replaced[cell] = 1;
            }
        }
 
        for (const VertexIndex v : used) {
            slotOf[static_cast<std::size_t>(v)] = -1;
        }
 
        if (fresh.empty()) {
            // Nothing came off at this degree: let the next pass reach further,
            // and give up once even a generous star finds nothing.
            if (degreeLimit >= 64) {
                break;
            }
            degreeLimit += 4;
            continue;
        }
        kept.clear();
        kept.reserve(active.size());
        for (const std::uint32_t c : active) {
            if (!replaced[c]) {
                kept.push_back(c);
            }
        }
        kept.insert(kept.end(), fresh.begin(), fresh.end());
        active.swap(kept);
        degreeLimit = pointLocationMaxStar;
    }
 
    // ---- lay the cells out along the descent -----------------------------
    // A query reads one cell per level, and the cells it reads are scattered
    // over the order they were made in — which follows the removal order,
    // itself sorted by star size. Renumbering them depth first from the top
    // instead puts a cell beside the child it hands the query to, so a descent
    // walks memory forwards rather than jumping over the whole hierarchy.
    {
        // The top level keeps its triangles, being the one a query has to be
        // placed in; the rest go now, before the relaid arrays double what the
        // hierarchy holds.
        kp->roots.reserve(active.size());
        kp->rootShape.reserve(active.size());
        for (const std::uint32_t c : active) {
            kp->rootShape.push_back(shape[c]);
        }
        shape.clear();
        shape.shrink_to_fit();
 
        const std::uint32_t cellCount = static_cast<std::uint32_t>(kp->cells.size());
        std::vector<std::uint32_t> relabel(cellCount, ~0u);
        std::vector<Cell> laid;
        laid.reserve(cellCount);
        std::vector<std::uint32_t> laidRun;
        laidRun.reserve(kp->run.size());
        std::vector<std::uint32_t> stack(active.rbegin(), active.rend());
        while (!stack.empty()) {
            const std::uint32_t c = stack.back();
            stack.pop_back();
            if (relabel[c] != ~0u) {
                continue;
            }
            relabel[c] = static_cast<std::uint32_t>(laid.size());
            laid.push_back(kp->cells[c]);
            const Cell& cell = kp->cells[c];
            for (std::uint32_t k = cell.runCount & ~Kirkpatrick::FULL_TURN; k-- > 0;) {
                stack.push_back(kp->run[cell.runBegin + 2 * k + 1]);
            }
        }
        for (Cell& cell : laid) {
            const std::uint32_t count = cell.runCount & ~Kirkpatrick::FULL_TURN;
            if (count == 0) {
                continue;  // runBegin is this leaf's seed, not a run
            }
            const std::uint32_t from = cell.runBegin;
            cell.runBegin = static_cast<std::uint32_t>(laidRun.size());
            for (std::uint32_t k = 0; k < count; ++k) {
                laidRun.push_back(kp->run[from + 2 * k]);
                laidRun.push_back(relabel[kp->run[from + 2 * k + 1]]);
            }
            laidRun.push_back(kp->run[from + 2 * count]);
        }
        for (const std::uint32_t c : active) {
            kp->roots.push_back(relabel[c]);
        }
        kp->cells.swap(laid);
        kp->run.swap(laidRun);
    }
    pointLocation_ = std::move(kp);
    pointLocationRevision_ = revision_;
}
 
namespace detail {

struct ConvexCoverBuilder {
    template <TriangleConcept TriangleType, SegmentConcept SegmentType>
    [[nodiscard]] static Graph<TriangleType> visibilityGraph(
        const Triangulation<TriangleType, SegmentType>& triangulation) {
        return triangulation.convexCoverVisibilityGraph();
    }
};
 
}  // namespace detail
 
// Deduction guides: the stored triangle type is not deducible from the
// container-templated constructors on their own. From a triangle container it is
// the element type (keeping its labels); from a segment container it is
// `Triangle<PointType>`. The `requires` clauses keep the two guides disjoint.
template <class TriangleRange>
    requires TriangleConcept<typename TriangleRange::value_type>
Triangulation(const TriangleRange&)
    -> Triangulation<typename TriangleRange::value_type>;
 
template <class SegmentRange>
    requires SegmentConcept<typename SegmentRange::value_type>
Triangulation(const SegmentRange&)
    -> Triangulation<Triangle<typename SegmentRange::value_type::PointType>,
                     typename SegmentRange::value_type>;
 
template <class PointRange>
    requires PointConcept<typename PointRange::value_type>
Triangulation(const PointRange&)
    -> Triangulation<Triangle<Point<typename PointRange::value_type::NumberType>>>;
 
// Point set with constraint segments (conforming constrained Delaunay): the
// point type comes from the point range, the stored edge type takes the
// segments' label so it survives. Disjoint from the polygon guides below —
// Polygon is not a point range.
template <class PointRange, class SegmentRange>
    requires PointConcept<typename PointRange::value_type> &&
             SegmentConcept<typename SegmentRange::value_type>
Triangulation(const PointRange&, const SegmentRange&)
    -> Triangulation<Triangle<Point<typename PointRange::value_type::NumberType>>,
                     Segment<Point<typename PointRange::value_type::NumberType>,
                             typename SegmentRange::value_type::LabelType>>;
 
template <class PointType>
Triangulation(const Polygon<PointType>&) -> Triangulation<Triangle<PointType>>;
 
// Polygon with extra interior points and/or constraint segments: the point and
// triangle types come from the polygon. When constraint segments are present the
// stored edge type takes their label (over the polygon's own point type), so the
// segments' labels survive; otherwise the edge type defaults like the
// polygon-only guide. The `requires` clauses disambiguate the two two-argument
// forms (a points range vs a segment range), mirroring the disjoint guides above.
template <class PolyPoint, class PointRange>
    requires PointConcept<typename PointRange::value_type>
Triangulation(const Polygon<PolyPoint>&, const PointRange&)
    -> Triangulation<Triangle<PolyPoint>>;
 
template <class PolyPoint, class SegmentRange>
    requires SegmentConcept<typename SegmentRange::value_type>
Triangulation(const Polygon<PolyPoint>&, const SegmentRange&)
    -> Triangulation<Triangle<PolyPoint>,
                     Segment<PolyPoint, typename SegmentRange::value_type::LabelType>>;
 
template <class PolyPoint, class PointRange, class SegmentRange>
    requires PointConcept<typename PointRange::value_type> &&
             SegmentConcept<typename SegmentRange::value_type>
Triangulation(const Polygon<PolyPoint>&, const PointRange&, const SegmentRange&)
    -> Triangulation<Triangle<PolyPoint>,
                     Segment<PolyPoint, typename SegmentRange::value_type::LabelType>>;
 
// Region guides, mirroring the polygon ones above.
template <class RegionPoint>
Triangulation(const PolygonWithHoles<RegionPoint>&) -> Triangulation<Triangle<RegionPoint>>;
 
template <class RegionPoint, class PointRange>
    requires PointConcept<typename PointRange::value_type>
Triangulation(const PolygonWithHoles<RegionPoint>&, const PointRange&)
    -> Triangulation<Triangle<RegionPoint>>;
 
template <class RegionPoint, class SegmentRange>
    requires SegmentConcept<typename SegmentRange::value_type>
Triangulation(const PolygonWithHoles<RegionPoint>&, const SegmentRange&)
    -> Triangulation<Triangle<RegionPoint>,
                     Segment<RegionPoint, typename SegmentRange::value_type::LabelType>>;
 
template <class RegionPoint, class PointRange, class SegmentRange>
    requires PointConcept<typename PointRange::value_type> &&
             SegmentConcept<typename SegmentRange::value_type>
Triangulation(const PolygonWithHoles<RegionPoint>&, const PointRange&, const SegmentRange&)
    -> Triangulation<Triangle<RegionPoint>,
                     Segment<RegionPoint, typename SegmentRange::value_type::LabelType>>;
 
// Region-set guides, mirroring the region ones above.
template <class SetPoint>
Triangulation(const PolygonSet<SetPoint>&) -> Triangulation<Triangle<SetPoint>>;
 
template <class SetPoint, class PointRange>
    requires PointConcept<typename PointRange::value_type>
Triangulation(const PolygonSet<SetPoint>&, const PointRange&)
    -> Triangulation<Triangle<SetPoint>>;
 
template <class SetPoint, class SegmentRange>
    requires SegmentConcept<typename SegmentRange::value_type>
Triangulation(const PolygonSet<SetPoint>&, const SegmentRange&)
    -> Triangulation<Triangle<SetPoint>,
                     Segment<SetPoint, typename SegmentRange::value_type::LabelType>>;
 
template <class SetPoint, class PointRange, class SegmentRange>
    requires PointConcept<typename PointRange::value_type> &&
             SegmentConcept<typename SegmentRange::value_type>
Triangulation(const PolygonSet<SetPoint>&, const PointRange&, const SegmentRange&)
    -> Triangulation<Triangle<SetPoint>,
                     Segment<SetPoint, typename SegmentRange::value_type::LabelType>>;
 
// Out-of-line: Polygon::triangulation is declared in shape/polygon.hpp (which
// precedes this header in the layering) but can only be defined once
// Triangulation and its deduction guides are visible.
template <class PointType_, class TLabel>
auto Polygon<PointType_, TLabel>::triangulation() const {
    return Triangulation(*this);
}
 
template <class PointType_, class TLabel>
template <class SegmentRange>
auto Polygon<PointType_, TLabel>::triangulation(const SegmentRange& segments) const {
    return Triangulation(*this, segments);
}
 
template <class PointType_, class TLabel>
template <class PointRange, class SegmentRange>
auto Polygon<PointType_, TLabel>::triangulation(const PointRange& points,
                                                const SegmentRange& segments) const {
    return Triangulation(*this, points, segments);
}
 
template <class PointType_, class TLabel>
std::vector<Convex<PointType_>> Polygon<PointType_, TLabel>::convexPartition() const {
    return triangulation().convexPartition();
}
 
template <class PointType_, class TLabel>
std::vector<Convex<PointType_>> Polygon<PointType_, TLabel>::convexCovering() const {
    const auto partition = triangulation();
    const auto triangles = partition.triangles();
    const auto cliques = detail::ConvexCoverBuilder::visibilityGraph(partition).cliqueCover();
 
    std::vector<Convex<PointType_>> result;
    result.reserve(cliques.size());
    for (const auto& clique : cliques) {
        std::vector<PointType_> vertices;
        vertices.reserve(3 * clique.size());
        for (const auto& triangle : clique) {
            vertices.push_back(triangle.a());
            vertices.push_back(triangle.b());
            vertices.push_back(triangle.c());
        }
        result.emplace_back(vertices);
    }
 
    // A clique partition covers every source triangle, but one clique hull may
    // also cover triangles assigned to other cliques. Remove such redundant
    // hulls while preserving coverage of the whole triangulation.
    for (std::size_t i = result.size(); i-- > 0;) {
        bool redundant = true;
        for (const auto& triangle : triangles) {
            bool coveredElsewhere = false;
            for (std::size_t j = 0; j < result.size(); ++j) {
                if (j != i && result[j].contains(triangle)) {
                    coveredElsewhere = true;
                    break;
                }
            }
            if (!coveredElsewhere) {
                redundant = false;
                break;
            }
        }
        if (redundant) {
            result.erase(result.begin() + static_cast<std::ptrdiff_t>(i));
        }
    }
 
    std::sort(result.begin(), result.end());
    return result;
}
 
// Out-of-line for the same reason: declared in shape/polygonwithholes.hpp.
template <class PointType_, class TLabel>
auto PolygonWithHoles<PointType_, TLabel>::triangulation() const {
    return Triangulation(*this);
}
 
template <class PointType_, class TLabel>
std::vector<Convex<PointType_>> PolygonWithHoles<PointType_, TLabel>::convexPartition() const {
    return triangulation().convexPartition();
}
 
template <class PointType_, class TLabel>
std::vector<Convex<PointType_>> PolygonWithHoles<PointType_, TLabel>::convexCovering() const {
    return triangulation().convexCovering();
}
 
template <class PointType_, class TLabel>
template <class SegmentRange>
auto PolygonWithHoles<PointType_, TLabel>::triangulation(const SegmentRange& segments) const {
    return Triangulation(*this, segments);
}
 
template <class PointType_, class TLabel>
auto PolygonSet<PointType_, TLabel>::triangulation() const {
    return Triangulation(*this);
}
 
template <class PointType_, class TLabel>
template <class SegmentRange>
auto PolygonSet<PointType_, TLabel>::triangulation(const SegmentRange& segments) const {
    return Triangulation(*this, segments);
}
 
template <class PointType_, class TLabel>
std::vector<Convex<PointType_>> PolygonSet<PointType_, TLabel>::convexPartition() const {
    return triangulation().convexPartition();
}
 
template <class PointType_, class TLabel>
std::vector<Convex<PointType_>> PolygonSet<PointType_, TLabel>::convexCovering() const {
    return triangulation().convexCovering();
}
 
// A set's interior is the union of its components', so any component's own
// witness serves — no triangulation of the whole set needed.
template <class PointType_, class TLabel>
template <class ResultNumber>
Point<ResultNumber> PolygonSet<PointType_, TLabel>::pointInside() const {
    if (components_.empty()) {
        return Point<ResultNumber>();
    }
    return components_.front().template pointInside<ResultNumber>();
}
 
// pointInside also lives here rather than in measures.hpp.  It first tries the
// outer ring's cheap ear/diagonal witness. A hole can occupy that ear or
// interrupt that diagonal, however, so a witness in or on a hole falls back to
// the triangulated domain: every one of its triangles is inside the region by
// construction.
template <class PointType_, class TLabel>
template <class ResultNumber>
Point<ResultNumber> PolygonWithHoles<PointType_, TLabel>::pointInside() const {
    if (isDegenerate()) {
        // No interior to point at (UB per the library contract); fall back to a
        // representative point rather than triangulating nothing.
        return verticesCentroid<ResultNumber>();
    }
 
    const auto outerWitness = outer_.template pointInside<ResultNumber>();
    bool insideHole = false;
    for (const auto& hole : holes_) {
        if (hole.contains(outerWitness)) {
            insideHole = true;
            break;
        }
    }
    if (!insideHole) {
        return outerWitness;
    }
 
    // Any domain triangle serves, so stop at the first non-degenerate one the
    // visit meets rather than materializing and sorting the whole mesh.
    const auto mesh = triangulation();
    Point<ResultNumber> witness = verticesCentroid<ResultNumber>();
    mesh.visitTriangles([&](const auto& triangle) {
        if (triangle.isDegenerate()) {
            return false;
        }
        witness = triangle.template pointInside<ResultNumber>();
        return true;
    });
    return witness;
}
 
template <class PointType_, class TLabel>
template <class OtherShape>
bool PolygonWithHoles<PointType_, TLabel>::pointInsideInteriorContainedIn(const OtherShape& shape) const {
    const auto witness = pointInside<NumberType>();
    if (interiorContains(witness)) {
        return shape.interiorContains(witness);
    }
    // pointInside() divides by 4 (Triangle::pointInside), so integer truncation
    // can round it onto the boundary; scaling by 4 makes it exact without
    // changing containment. Scaling is a similarity, so the scaled region's
    // triangulation is the scaled triangulation and the witness scales with it.
    return (shape * 4).interiorContains((*this * 4).template pointInside<NumberType>());
}
 
}  // namespace pgl