booleans.hpp

Go to the documentation of this file.

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

 
#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <cstddef>
#include <limits>
#include <map>
#include <numeric>
#include <ranges>
#include <type_traits>
#include <variant>
#include <vector>
 
namespace pgl {
 
namespace detail {

template <class ExactPoint>
void dropCollinearRingVertices(std::vector<ExactPoint>& ring) {
    bool changed = true;
    while (changed && ring.size() > 3) {
        changed = false;
        for (std::size_t i = 0; i < ring.size() && ring.size() > 3;) {
            const ExactPoint& previous = ring[(i + ring.size() - 1) % ring.size()];
            const ExactPoint& next = ring[(i + 1) % ring.size()];
            if (collinear(previous, ring[i], next)) {
                ring.erase(ring.begin() + static_cast<std::ptrdiff_t>(i));
                changed = true;
            } else {
                ++i;
            }
        }
    }
}
 
// ringOrientation and splitWalkIntoRings, which the extraction below also
// uses, live beside the arrangement in algorithm/arrangement.hpp: they are
// what turns any cell complex's boundary walks into rings.

template <class ResultPoint, class ExactPoint>
PolygonSet<ResultPoint> regularizedCellsFromKeep(
    const Arrangement<ExactPoint>& arrangement, const std::vector<char>& keep) {
    using HalfedgeId = typename Arrangement<ExactPoint>::HalfedgeId;
    using ExactNumber = typename ExactPoint::NumberType;
    using ExactPolygon = Polygon<ExactPoint>;
    using ResultPolygon = Polygon<ResultPoint>;
    constexpr std::size_t none = std::numeric_limits<std::size_t>::max();
 
    std::vector<PolygonWithHoles<ResultPoint>> result;
    const auto isKept = [&](HalfedgeId h) { return keep[arrangement.face(h).index()] != 0; };
 
    // Kept faces sharing an edge are one piece of the result. Faces meeting at a
    // single vertex are not: the result pinches shut there, and the two sides
    // have to come back as two regions, since neither a polygon nor a region may
    // have a self-touching outer ring.
    std::vector<std::size_t> parent(arrangement.faceCount());
    for (std::size_t i = 0; i < parent.size(); ++i) {
        parent[i] = i;
    }
    const auto findRoot = [&parent](std::size_t x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    };
    std::vector<HalfedgeId> boundary;
    for (std::uint32_t i = 0; i < arrangement.halfedgeCount(); ++i) {
        const HalfedgeId h(i);
        if (!isKept(h)) {
            continue;
        }
        if (isKept(arrangement.twin(h))) {
            parent[findRoot(arrangement.face(h).index())] =
                findRoot(arrangement.face(arrangement.twin(h)).index());
        } else {
            boundary.push_back(h);
        }
    }
 
    // Walking the boundary: from a halfedge with kept material on its left and
    // none across it, follow the face's own cycle, stepping over any edge whose
    // far side is kept too — that is the rotation around the shared vertex the
    // DCEL stores, and it stops at the far side of the same fan of kept faces.
    // An edge kept on both sides is interior to the result and is skipped, which
    // is what regularizes a slit away.
    const auto nextBoundary = [&](HalfedgeId h) {
        HalfedgeId ahead = arrangement.next(h);
        while (isKept(arrangement.twin(ahead))) {
            ahead = arrangement.next(arrangement.twin(ahead));
        }
        return ahead;
    };
 
    std::map<std::size_t, std::vector<std::vector<ExactPoint>>> ringsOfPiece;
    std::vector<char> walked(arrangement.halfedgeCount(), 0);
    for (const HalfedgeId start : boundary) {
        if (walked[start.index()] != 0) {
            continue;
        }
        std::vector<ExactPoint> walk;
        HalfedgeId h = start;
        do {
            walked[h.index()] = 1;
            walk.push_back(arrangement[arrangement.source(h)]);
            h = nextBoundary(h);
        } while (h != start);
        splitWalkIntoRings(walk, ringsOfPiece[findRoot(arrangement.face(start).index())]);
    }
 
    // Converting an already canonical ring into the caller's coordinates keeps
    // it canonical whenever the conversion is the identity, and only then: an
    // integral result type truncates, which can reorder the vertices or flip the
    // ring, so that case is renormalized as usual.
    const auto convert = [](const std::vector<ExactPoint>& ring) {
        std::vector<ResultPoint> converted;
        converted.reserve(ring.size());
        for (const ExactPoint& vertex : ring) {
            converted.emplace_back(vertex);
        }
        constexpr bool exact = std::is_same_v<ResultPoint, ExactPoint>;
        return ResultPolygon(std::move(converted), /*trusted=*/exact);
    };
 
    for (auto& entry : ringsOfPiece) {
        std::vector<ExactPolygon> outers;
        std::vector<ExactPolygon> holes;
        for (std::vector<ExactPoint>& ring : entry.second) {
            dropCollinearRingVertices(ring);
            const int orientation = ringOrientation(ring);
            if (orientation == 0) {
                continue;  // a ring bounding no area is not part of any region
            }
            // The walks come out of a planar subdivision, so every ring here is
            // simple and its orientation is one predicate away. Putting it into
            // canonical form by hand then keeps @ref Polygon from measuring it
            // all over again, which over rationals costs more than everything
            // else in this function put together.
            if (orientation < 0) {
                std::reverse(ring.begin(), ring.end());
            }
            std::rotate(ring.begin(), std::min_element(ring.begin(), ring.end()), ring.end());
            (orientation > 0 ? outers : holes).emplace_back(std::move(ring), /*trusted=*/true);
        }
        if (outers.empty()) {
            continue;
        }
 
        // An edge-connected piece has a single outer ring unless it pinches shut
        // at a vertex, where the walk above cuts it into one ring per side. Each
        // hole then goes to the smallest outer ring holding it, decided by a
        // witness strictly inside the hole (the rings meet at most along their
        // boundaries, so the closed containment is unambiguous).
        std::vector<std::vector<ResultPolygon>> holesOfOuter(outers.size());
        for (const ExactPolygon& hole : holes) {
            std::size_t owner = 0;
            if (outers.size() > 1) {
                const ExactPoint witness = hole.template pointInside<ExactNumber>();
                std::size_t best = none;
                for (std::size_t i = 0; i < outers.size(); ++i) {
                    if (outers[i].contains(witness) &&
                        (best == none || outers[i].twiceArea() < outers[best].twiceArea())) {
                        best = i;
                    }
                }
                owner = best == none ? 0 : best;
            }
            holesOfOuter[owner].push_back(convert(hole.vertices()));
        }
        // The same bargain `convert` strikes for a ring, struck once more for
        // the region: its holes are canonical polygons already, and not one of
        // them bounds no area — a ring whose orientation vanished was dropped
        // above — so all the region's own normalization would find to do is
        // sort them, after measuring every one of them exactly to learn what is
        // already known. Over rationals that measurement is the single most
        // expensive thing this function does. A result type that is not the
        // exact one may collapse a ring on the way, so it still normalizes.
        constexpr bool exact = std::is_same_v<ResultPoint, ExactPoint>;
        for (std::size_t i = 0; i < outers.size(); ++i) {
            if constexpr (exact) {
                std::sort(holesOfOuter[i].begin(), holesOfOuter[i].end());
            }
            result.emplace_back(convert(outers[i].vertices()), std::move(holesOfOuter[i]),
                                /*trusted=*/exact);
        }
    }
    // The pieces have pairwise disjoint interiors and share no stretch of edge —
    // an edge kept on both sides is interior to the result and never reaches the
    // boundary walk — so they are a canonical PolygonSet once sorted, and the
    // set adopts them without measuring a single area again.
    std::sort(result.begin(), result.end());
    return PolygonSet<ResultPoint>(std::move(result), /*trusted=*/true);
}

template <class ExactPoint>
Arrangement<ExactPoint> framedArrangement(const std::vector<Segment<ExactPoint>>& cuts) {
    using ExactNumber = typename ExactPoint::NumberType;
    using ExactSegment = Segment<ExactPoint>;
 
    assert(!cuts.empty());
    ExactNumber loX = cuts.front().min().x();
    ExactNumber loY = cuts.front().min().y();
    ExactNumber hiX = loX;
    ExactNumber hiY = loY;
    for (const ExactSegment& cut : cuts) {
        for (const ExactPoint& point : {cut.min(), cut.max()}) {
            loX = std::min(loX, point.x());
            loY = std::min(loY, point.y());
            hiX = std::max(hiX, point.x());
            hiY = std::max(hiY, point.y());
        }
    }
    const ExactNumber margin(1);
    const std::array<ExactPoint, 4> corners{
        ExactPoint(loX - margin, loY - margin), ExactPoint(hiX + margin, loY - margin),
        ExactPoint(hiX + margin, hiY + margin), ExactPoint(loX - margin, hiY + margin)};
    std::vector<ExactSegment> segments = cuts;
    for (std::size_t i = 0; i < corners.size(); ++i) {
        segments.emplace_back(corners[i], corners[(i + 1) % corners.size()]);
    }
    return Arrangement<ExactPoint>(segments);
}

template <class ResultPoint, class ExactPoint, class KeepWitness>
PolygonSet<ResultPoint> regularizedCells(
    const std::vector<Segment<ExactPoint>>& cuts, KeepWitness keepWitness) {
    using ExactNumber = typename ExactPoint::NumberType;
 
    if (cuts.empty()) {
        return {};  // no operand has an edge, so none has area
    }
 
    const Arrangement<ExactPoint> arrangement = framedArrangement(cuts);
    using FaceId = typename Arrangement<ExactPoint>::FaceId;
    std::vector<char> keep(arrangement.faceCount(), 0);
    for (std::uint32_t i = 0; i < arrangement.faceCount(); ++i) {
        const FaceId f(i);
        if (!arrangement.isUnbounded(f)) {
            keep[f.index()] = static_cast<char>(
                keepWitness(arrangement.template witness<ExactNumber>(f)));
        }
    }
    return regularizedCellsFromKeep<ResultPoint>(arrangement, keep);
}

template <class ResultPoint, class ShapeA, class ShapeB, class KeepCell>
PolygonSet<ResultPoint> regularizedBoolean(const ShapeA& a, const ShapeB& b, KeepCell keepCell) {
    using ExactNumber = Exact1DNumber<typename ShapeA::NumberType, typename ShapeB::NumberType>;
    using ExactPoint = Point<ExactNumber>;
 
    const bool aHasArea = !a.isDegenerate();
    const bool bHasArea = !b.isDegenerate();
 
    std::vector<Segment<ExactPoint>> cuts;
    appendCutSegments<ExactPoint>(a, cuts);
    appendCutSegments<ExactPoint>(b, cuts);
    return regularizedCells<ResultPoint>(
        cuts, [&a, &b, aHasArea, bHasArea, &keepCell](const ExactPoint& witness) {
            return keepCell(aHasArea && a.contains(witness), bHasArea && b.contains(witness));
        });
}

template <class ResultPoint, class ShapeType>
PolygonSet<ResultPoint> regularizedUnionByCoverage(
    const std::vector<ShapeType>& distinct) {
    using ShapeNumber = typename ShapeType::NumberType;
    using ExactNumber = Exact1DNumber<ShapeNumber, ShapeNumber>;
    using ExactPoint = Point<ExactNumber>;
 
    // Start outside every piece in the unbounded face, then propagate those
    // parity bits across the arrangement's face adjacency graph.
    const Arrangement<ExactPoint> arrangement(distinct, detail::simpleBoundaries);
    using HalfedgeId = typename Arrangement<ExactPoint>::HalfedgeId;
    using FaceId = typename Arrangement<ExactPoint>::FaceId;
    const std::size_t faceCount = arrangement.faceCount();
    const std::size_t pieceCount = distinct.size();
    constexpr std::size_t wordBits = 64;
    const std::size_t words = (pieceCount + wordBits - 1) / wordBits;
 
    // The halfedges bounding each face, as one pair of arrays rather than a
    // vector per face: the face count runs to hundreds of thousands here, and
    // that many separate allocations costs more than the traversal.
    std::vector<std::uint32_t> faceEdgeBegin(faceCount + 1, 0);
    for (std::uint32_t i = 0; i < arrangement.halfedgeCount(); ++i) {
        const HalfedgeId h(i);
        ++faceEdgeBegin[arrangement.face(h).index() + 1];
    }
    for (std::size_t i = 0; i < faceCount; ++i) {
        faceEdgeBegin[i + 1] += faceEdgeBegin[i];
    }
    std::vector<std::uint32_t> faceEdge(arrangement.halfedgeCount(), 0);
    {
        std::vector<std::uint32_t> cursor(faceEdgeBegin.begin(), faceEdgeBegin.end() - 1);
        for (std::uint32_t i = 0; i < arrangement.halfedgeCount(); ++i) {
            const HalfedgeId h(i);
            faceEdge[cursor[arrangement.face(h).index()]++] = h.index();
        }
    }
 
    // One shared membership word set, not one per face. Walking a spanning
    // tree of the face adjacency graph depth first, a face's membership is
    // its parent's with the crossed edge's origins toggled, so the descent
    // toggles them and the return toggles them back. The flip is its own
    // inverse in the coverage count too — a bit that was clear counts one
    // more piece, a bit that was set one fewer — so both restore exactly.
    // Storing membership per face instead would be faces x pieces bits,
    // nearly half a gigabyte on the largest shape-pair cell.
    std::vector<std::uint64_t> membership(words, 0);
    std::size_t covered = 0;
    const auto crossEdge = [&](std::uint32_t halfedge) {
        for (const std::uint32_t origin : arrangement.originsOf(HalfedgeId(halfedge))) {
            const std::size_t word = origin / wordBits;
            const std::uint64_t mask = std::uint64_t{1} << (origin % wordBits);
            if ((membership[word] & mask) == 0) {
                ++covered;
            } else {
                --covered;
            }
            membership[word] ^= mask;
        }
    };
 
    struct Frame {
        std::uint32_t face;
        std::uint32_t cursor;   // next index into faceEdge
        std::uint32_t entered;  // halfedge crossed to get here, or noEdge at the root
    };
    constexpr std::uint32_t noEdge = ~std::uint32_t{0};
 
    std::vector<std::size_t> coverage(faceCount, 0);
    std::vector<char> seen(faceCount, 0);
    std::vector<Frame> stack;
    // Face 0 is the unbounded one, which lies outside every piece.
    stack.push_back(Frame{0, faceEdgeBegin[0], noEdge});
    seen[0] = 1;
 
    while (!stack.empty()) {
        Frame& top = stack.back();
        if (top.cursor == faceEdgeBegin[top.face + 1]) {
            if (top.entered != noEdge) {
                crossEdge(top.entered);  // undo, restoring the parent's state
            }
            stack.pop_back();
            continue;
        }
        const std::uint32_t h = faceEdge[top.cursor++];
        const std::uint32_t next = arrangement.face(arrangement.twin(HalfedgeId(h))).index();
        if (seen[next] != 0) {
            continue;
        }
        seen[next] = 1;
        crossEdge(h);
        coverage[next] = covered;
        // `top` may dangle after this, so nothing above may be used again.
        stack.push_back(Frame{next, faceEdgeBegin[next], h});
    }
 
    assert(covered == 0);  // every descent undone
    assert(std::ranges::all_of(seen, [](char value) { return value != 0; }));
    std::vector<char> keep(faceCount, 0);
    for (std::uint32_t i = 0; i < arrangement.faceCount(); ++i) {
        const FaceId f(i);
        keep[f.index()] = static_cast<char>(coverage[f.index()] != 0);
    }
    return regularizedCellsFromKeep<ResultPoint>(arrangement, keep);
}
 
// The boundary of a union is a subset of the operands' boundaries. For a set
// of lattice triangles, find that subset directly: another convex triangle
// covers one interval of an edge, obtained by clipping the edge's supporting
// line against its three half-planes. Subtracting those intervals leaves only
// the exposed pieces. All the rejection and clipping arithmetic stays in
// int128; rationals are constructed only for actual interval ends on the output
// boundary. The final, usually small arrangement turns those pieces into the
// canonical PolygonSet and handles touching/collinear degeneracies centrally.
//
// The subtraction is what keeps the scan from being quadratic in earnest. An
// edge carries what is still uncovered of it, as a short sorted list of
// disjoint intervals, and a triangle that covers nothing left of it costs one
// clip and nothing else; the moment the list empties the edge is done and the
// remaining triangles are never looked at. Only an edge with a piece on the
// output boundary is ever tested against all of them, so the work is
// n x (output complexity) rather than n² whenever the union covers itself —
// which is the case a union of many overlapping pieces actually is. Testing the
// largest triangles first is the same bet: a big one covers a long stretch of
// whatever it meets, so the list empties in fewer tests.
template <class ResultPoint, class TriangleType>
std::optional<PolygonSet<ResultPoint>> regularizedUnionOfIntegralTriangles(
    const std::vector<TriangleType>& triangles) {
    using ShapeNumber = typename TriangleType::NumberType;
    using ExactNumber = Exact1DNumber<ShapeNumber, ShapeNumber>;
    using ExactPoint = Point<ExactNumber>;
    using ExactSegment = Segment<ExactPoint>;
    using Wide = int128;
 
    struct IPoint {
        std::int64_t x;
        std::int64_t y;
    };
    struct ITriangle {
        std::array<IPoint, 3> vertex;
        std::int64_t minX, minY, maxX, maxY;
        Wide twiceArea;
    };
    constexpr std::int64_t safeCoordinate = 1000000000;
    const auto narrow = [](const ShapeNumber& value) -> std::optional<std::int64_t> {
        if constexpr (is_Rational_v<ShapeNumber>) {
            if (!value.isInteger()) {
                return std::nullopt;
            }
            using Integer = rational_int_t<ShapeNumber>;
            const Integer integer = static_cast<Integer>(value);
            if (!detail::representableAs<std::int64_t>(integer)) {
                return std::nullopt;
            }
            return detail::narrowTo<std::int64_t>(integer);
        } else if constexpr (detail::extended_integral<ShapeNumber> ||
                             std::same_as<ShapeNumber, BigInt>) {
            if (!detail::representableAs<std::int64_t>(value)) {
                return std::nullopt;
            }
            return detail::narrowTo<std::int64_t>(value);
        } else {
            return std::nullopt;
        }
    };
 
    const auto cross = [](Wide ax, Wide ay, Wide bx, Wide by) {
        return ax * by - ay * bx;
    };
 
    std::vector<ITriangle> integral;
    integral.reserve(triangles.size());
    for (const TriangleType& triangle : triangles) {
        ITriangle converted{};
        const auto vertices = triangle.vertices();
        for (std::size_t i = 0; i < 3; ++i) {
            const auto x = narrow(vertices[i].x());
            const auto y = narrow(vertices[i].y());
            if (!x || !y) {
                return std::nullopt;
            }
            if (*x < -safeCoordinate || safeCoordinate < *x ||
                *y < -safeCoordinate || safeCoordinate < *y) {
                return std::nullopt;
            }
            converted.vertex[i] = {*x, *y};
        }
        converted.minX = converted.maxX = converted.vertex[0].x;
        converted.minY = converted.maxY = converted.vertex[0].y;
        for (const IPoint& point : converted.vertex) {
            converted.minX = std::min(converted.minX, point.x);
            converted.minY = std::min(converted.minY, point.y);
            converted.maxX = std::max(converted.maxX, point.x);
            converted.maxY = std::max(converted.maxY, point.y);
        }
        const IPoint& a = converted.vertex[0];
        const IPoint& b = converted.vertex[1];
        const IPoint& c = converted.vertex[2];
        converted.twiceArea =
            cross(Wide(b.x) - a.x, Wide(b.y) - a.y, Wide(c.x) - a.x, Wide(c.y) - a.y);
        if (converted.twiceArea < 0) {
            converted.twiceArea = -converted.twiceArea;
        }
        integral.push_back(converted);
    }
 
    // The order the coverers are tested in, largest first.
    std::vector<std::uint32_t> order(integral.size());
    std::iota(order.begin(), order.end(), std::uint32_t{0});
    std::sort(order.begin(), order.end(), [&integral](std::uint32_t left, std::uint32_t right) {
        return integral[left].twiceArea > integral[right].twiceArea;
    });
 
    struct Fraction {
        Wide numerator;
        Wide denominator;
    };
    const auto fraction = [](Wide numerator, Wide denominator) {
        if (denominator < 0) {
            numerator = -numerator;
            denominator = -denominator;
        }
        return Fraction{numerator, denominator};
    };
    const auto less = [](const Fraction& left, const Fraction& right) {
        return left.numerator * right.denominator <
               right.numerator * left.denominator;
    };
    const auto maximum = [&less](const Fraction& left, const Fraction& right) {
        return less(left, right) ? right : left;
    };
    const auto minimum = [&less](const Fraction& left, const Fraction& right) {
        return less(left, right) ? left : right;
    };
    const auto exact = [](const Fraction& value) {
        return ExactNumber(BigInt(value.numerator), BigInt(value.denominator));
    };
    constexpr Fraction zero{0, 1};
    constexpr Fraction one{1, 1};
    std::vector<ExactSegment> exposed;
    // What is still uncovered of the edge being scanned: sorted, disjoint, and
    // usually one interval or none.
    std::vector<std::pair<Fraction, Fraction>> uncovered;
 
    for (std::size_t owner = 0; owner < integral.size(); ++owner) {
        const ITriangle& mine = integral[owner];
        for (std::size_t edge = 0; edge < 3; ++edge) {
            const IPoint from = mine.vertex[edge];
            const IPoint to = mine.vertex[(edge + 1) % 3];
            const Wide dx = Wide(to.x) - from.x;
            const Wide dy = Wide(to.y) - from.y;
            const std::int64_t edgeMinX = std::min(from.x, to.x);
            const std::int64_t edgeMinY = std::min(from.y, to.y);
            const std::int64_t edgeMaxX = std::max(from.x, to.x);
            const std::int64_t edgeMaxY = std::max(from.y, to.y);
            uncovered.assign(1, {zero, one});
 
            for (const std::uint32_t other : order) {
                const ITriangle& theirs = integral[other];
                if (other == owner || edgeMaxX < theirs.minX || theirs.maxX < edgeMinX ||
                    edgeMaxY < theirs.minY || theirs.maxY < edgeMinY) {
                    continue;
                }
                Fraction low = zero;
                Fraction high = one;
                bool feasible = true;
                bool coincident = false;
                bool interiorOnRight = false;
                for (std::size_t wall = 0; wall < 3; ++wall) {
                    const IPoint a = theirs.vertex[wall];
                    const IPoint b = theirs.vertex[(wall + 1) % 3];
                    const Wide wx = Wide(b.x) - a.x;
                    const Wide wy = Wide(b.y) - a.y;
                    const Wide c = cross(wx, wy, Wide(from.x) - a.x,
                                         Wide(from.y) - a.y);
                    const Wide slope = cross(wx, wy, dx, dy);
                    if (slope == 0) {
                        if (c < 0) {
                            feasible = false;
                            break;
                        }
                        if (c == 0) {
                            coincident = true;
                            interiorOnRight = wx * dx + wy * dy < 0;
                        }
                        continue;
                    }
                    const Fraction at = fraction(-c, slope);
                    if (slope > 0) {
                        low = maximum(low, at);
                    } else {
                        high = minimum(high, at);
                    }
                    if (!less(low, high)) {
                        feasible = false;
                        break;
                    }
                }
                if (!feasible || !less(low, high) || (coincident && !interiorOnRight)) {
                    continue;
                }
 
                // Subtract [low, high]. Only the intervals it meets are
                // touched: those from the first one ending past `low` to the
                // last one starting before `high`. They collapse to at most a
                // head remnant and a tail remnant, so the list grows by one
                // only where a single interval is split in two.
                std::size_t first = 0;
                while (first < uncovered.size() && !less(low, uncovered[first].second)) {
                    ++first;
                }
                std::size_t last = first;
                while (last < uncovered.size() && less(uncovered[last].first, high)) {
                    ++last;
                }
                if (first == last) {
                    continue;
                }
                const Fraction head = uncovered[first].first;
                const Fraction tail = uncovered[last - 1].second;
                const bool keepHead = less(head, low);
                const bool keepTail = less(high, tail);
                if (keepHead && keepTail && last - first == 1) {
                    uncovered[first] = {head, low};
                    uncovered.insert(uncovered.begin() + static_cast<std::ptrdiff_t>(first) + 1,
                                     {high, tail});
                } else {
                    std::size_t write = first;
                    if (keepHead) {
                        uncovered[write++] = {head, low};
                    }
                    if (keepTail) {
                        uncovered[write++] = {high, tail};
                    }
                    uncovered.erase(uncovered.begin() + static_cast<std::ptrdiff_t>(write),
                                    uncovered.begin() + static_cast<std::ptrdiff_t>(last));
                }
                if (uncovered.empty()) {
                    break;
                }
            }
 
            const auto pointAt = [&](const Fraction& parameter) {
                const ExactNumber t = exact(parameter);
                return ExactPoint(ExactNumber(from.x) + ExactNumber(dx) * t,
                                  ExactNumber(from.y) + ExactNumber(dy) * t);
            };
            for (const auto& piece : uncovered) {
                exposed.emplace_back(pointAt(piece.first), pointAt(piece.second));
            }
        }
    }
 
    return regularizedCells<ResultPoint>(exposed, [&triangles](const ExactPoint& witness) {
        return std::ranges::any_of(
            triangles, [&witness](const TriangleType& triangle) { return triangle.contains(witness); });
    });
}
 
}  // namespace detail

template <class ResultPoint, class ShapeRange>
PolygonSet<ResultPoint> regularizedUnionOf(const ShapeRange& shapes,
                                           bool simpleBoundaries = false) {
    using ShapeType = std::ranges::range_value_t<ShapeRange>;
 
    // What survives to be unioned. A shape with no interior contributes nothing
    // to `closure(union of interiors)`, so dropping it cannot change the result
    // on either path below — and the parity argument on the covered one *needs*
    // it gone: such a boundary is traversed twice, but the two traversals
    // coincide and the arrangement merges them into a single edge carrying that
    // origin once, so crossing it would toggle an odd number of times and report
    // the far side as inside. A repeat contributes nothing either, and would be
    // tested again for every face it covers. Both go before the sort, which then
    // has less to order.
    std::vector<ShapeType> distinct(std::ranges::begin(shapes), std::ranges::end(shapes));
    std::erase_if(distinct, [](const ShapeType& shape) { return shape.isDegenerate(); });
    std::sort(distinct.begin(), distinct.end());
    distinct.erase(std::unique(distinct.begin(), distinct.end()), distinct.end());
 
    if (distinct.empty()) {
        return {};
    }
 
    if constexpr (detail::is_polygon_set_v<ShapeType>) {
        // A set is already a union of regions with disjoint interiors, and it is
        // those regions the paths below classify — a set is not itself something
        // an @ref Arrangement can be built from. Separating them also weakens
        // what `simpleBoundaries` has to promise: a stretch of boundary two
        // components share is two origins once they are apart, which is the two
        // crossings parity needs to see.
        std::vector<typename ShapeType::ComponentType> components;
        for (const ShapeType& set : distinct) {
            components.insert(components.end(), set.components().begin(), set.components().end());
        }
        return regularizedUnionOf<ResultPoint>(components, simpleBoundaries);
    } else if constexpr (detail::is_convex_v<ShapeType> || detail::is_triangle_v<ShapeType> ||
                         detail::is_rectangle_v<ShapeType>) {
        if constexpr (detail::is_triangle_v<ShapeType> &&
                      !std::floating_point<typename ShapeType::NumberType>) {
            if (distinct.size() >= 16) {
                if (auto result =
                        detail::regularizedUnionOfIntegralTriangles<ResultPoint>(distinct)) {
                    return std::move(*result);
                }
            }
        }
        return detail::regularizedUnionByCoverage<ResultPoint>(distinct);
    } else {
        using ShapeNumber = typename ShapeType::NumberType;
        using ExactPoint = Point<detail::Exact1DNumber<ShapeNumber, ShapeNumber>>;
 
        if (simpleBoundaries) {
            return detail::regularizedUnionByCoverage<ResultPoint>(distinct);
        }
        std::vector<Segment<ExactPoint>> cuts;
        for (const ShapeType& shape : distinct) {
            detail::appendCutSegments<ExactPoint>(shape, cuts);
        }
        return detail::regularizedCells<ResultPoint>(cuts, [&distinct](const ExactPoint& witness) {
            return std::ranges::any_of(
                distinct, [&witness](const ShapeType& shape) { return shape.contains(witness); });
        });
    }
}
 
namespace detail {

template <class OtherShape>
constexpr decltype(auto) booleanOperand(const OtherShape& other) {
    if constexpr (is_rectangle_v<OtherShape> || is_triangle_v<OtherShape> ||
                  is_convex_v<OtherShape>) {
        return other.asPolygon();
    } else {
        return (other);
    }
}

template <class ResultPoint, class ShapeA, class ShapeB>
PolygonSet<ResultPoint> regularizedDifference(const ShapeA& a, const ShapeB& b) {
    return regularizedBoolean<ResultPoint>(a, b, [](bool inA, bool inB) { return inA && !inB; });
}

template <class ResultPoint, class ShapeA, class ShapeB>
PolygonSet<ResultPoint> regularizedUnion(const ShapeA& a, const ShapeB& b) {
    return regularizedBoolean<ResultPoint>(a, b, [](bool inA, bool inB) { return inA || inB; });
}

template <class ResultPoint, class ShapeA, class ShapeB>
PolygonSet<ResultPoint> regularizedIntersection(const ShapeA& a, const ShapeB& b) {
    if (!a.bbox().interiorsIntersect(b.bbox())) {
        return {};
    }
    return regularizedBoolean<ResultPoint>(a, b, [](bool inA, bool inB) { return inA && inB; });
}

template <class ResultPoint, class ShapeA, class ShapeB>
std::vector<std::variant<ResultPoint, Polyline<ResultPoint>, PolygonWithHoles<ResultPoint>>>
literalIntersection(const ShapeA& a, const ShapeB& b) {
    using ExactNumber = Exact1DNumber<typename ShapeA::NumberType, typename ShapeB::NumberType>;
    using ExactPoint = Point<ExactNumber>;
    using ResultPolyline = Polyline<ResultPoint>;
    using ResultRegion = PolygonWithHoles<ResultPoint>;
    using Piece = std::variant<ResultPoint, ResultPolyline, ResultRegion>;
 
    std::vector<Piece> pieces;
    // Boxes that miss each other settle the whole answer, this time including
    // the boundaries: the boxes are closed, so a shared point of the operands
    // would be a shared point of them.
    if (a.empty() || b.empty() || !a.bbox().intersects(b.bbox())) {
        return pieces;
    }
 
    std::vector<Segment<ExactPoint>> cuts;
    appendCutSegments<ExactPoint>(a, cuts);
    const std::size_t cutsOfA = cuts.size();
    appendCutSegments<ExactPoint>(b, cuts);
    // An operand whose every edge has zero length has collapsed onto a single
    // point — its own bounding box — which the arrangement below would never see
    // as a vertex, there being no edge to carry it.
    if (cutsOfA == 0) {
        const ExactPoint vertex(a.bbox().min());
        if (b.contains(vertex)) {
            pieces.emplace_back(ResultPoint(vertex));
        }
        return pieces;
    }
    if (cuts.size() == cutsOfA) {
        const ExactPoint vertex(b.bbox().min());
        if (a.contains(vertex)) {
            pieces.emplace_back(ResultPoint(vertex));
        }
        return pieces;
    }
 
    const Arrangement<ExactPoint> arrangement = framedArrangement(cuts);
    using FaceId = typename Arrangement<ExactPoint>::FaceId;
    using HalfedgeId = typename Arrangement<ExactPoint>::HalfedgeId;
    using VertexId = typename Arrangement<ExactPoint>::VertexId;
 
    std::vector<char> keep(arrangement.faceCount(), 0);
    for (std::uint32_t i = 0; i < arrangement.faceCount(); ++i) {
        const FaceId f(i);
        if (!arrangement.isUnbounded(f)) {
            const ExactPoint witness = arrangement.template witness<ExactNumber>(f);
            keep[f.index()] = static_cast<char>(a.contains(witness) && b.contains(witness));
        }
    }
 
    // The strands, as an undirected graph on the arrangement's own vertices.
    std::vector<std::uint32_t> strandDegree(arrangement.vertexCount(), 0);
    std::vector<std::vector<std::uint32_t>> adjacency(arrangement.vertexCount());
    for (std::uint32_t i = 0; i < arrangement.edgeCount(); ++i) {
        const HalfedgeId h(2 * i);
        if (keep[arrangement.face(h).index()] != 0 ||
            keep[arrangement.face(arrangement.twin(h)).index()] != 0) {
            continue;
        }
        const ExactPoint witness = arrangement.template witness<ExactNumber>(h);
        if (!a.contains(witness) || !b.contains(witness)) {
            continue;
        }
        const std::uint32_t source = arrangement.source(h).index();
        const std::uint32_t target = arrangement.target(h).index();
        adjacency[source].push_back(target);
        adjacency[target].push_back(source);
        ++strandDegree[source];
        ++strandDegree[target];
    }
 
    // The isolated contact points. A vertex beside a kept face is in that
    // region piece already, and the rotational fan the DCEL carries answers that
    // without building a vector per vertex.
    for (std::uint32_t i = 0; i < arrangement.vertexCount(); ++i) {
        if (strandDegree[i] != 0) {
            continue;
        }
        const VertexId v(i);
        const HalfedgeId start = arrangement.outgoing(v);
        assert(start.valid());  // every vertex here is an endpoint of some cut
        bool besideKept = false;
        HalfedgeId h = start;
        do {
            besideKept = keep[arrangement.face(h).index()] != 0;
            h = arrangement.next(arrangement.twin(h));
        } while (h != start && !besideKept);
        if (besideKept) {
            continue;
        }
        const ExactPoint& position = arrangement[v];
        if (a.contains(position) && b.contains(position)) {
            pieces.emplace_back(ResultPoint(position));
        }
    }
 
    // Peel the open strands from their loose ends, then take whatever is left —
    // every vertex of which now carries at least two strand edges — as closed
    // walks. Peeling cannot strand a new loose end: a walk runs on through every
    // vertex it leaves with one edge, and stops only at a junction or at nothing.
    const auto takeEdge = [&adjacency](std::uint32_t from, std::uint32_t to) {
        auto& out = adjacency[from];
        out.erase(std::find(out.begin(), out.end(), to));
        auto& back = adjacency[to];
        back.erase(std::find(back.begin(), back.end(), from));
    };
    const auto walkFrom = [&](std::uint32_t start, bool stopAtJunctions) {
        std::vector<std::uint32_t> walk{start};
        std::uint32_t current = start;
        while (!adjacency[current].empty() &&
               (!stopAtJunctions || adjacency[current].size() == 1)) {
            const std::uint32_t next = adjacency[current].front();
            takeEdge(current, next);
            walk.push_back(next);
            current = next;
        }
        return walk;
    };
    std::vector<std::vector<std::uint32_t>> strands;
    for (std::uint32_t i = 0; i < adjacency.size(); ++i) {
        if (adjacency[i].size() == 1) {
            strands.push_back(walkFrom(i, /*stopAtJunctions=*/true));
        }
    }
    for (std::uint32_t i = 0; i < adjacency.size(); ++i) {
        while (!adjacency[i].empty()) {
            strands.push_back(walkFrom(i, /*stopAtJunctions=*/false));
        }
    }
 
    // A vertex the arrangement put in the middle of a straight stretch says
    // nothing about the strand, exactly as in a ring of the region pieces. Only
    // one of strand degree two may go: a junction is shared with another strand.
    const auto removable = [&](std::uint32_t before, std::uint32_t at, std::uint32_t after) {
        return strandDegree[at] == 2 &&
               collinear(arrangement[VertexId(before)], arrangement[VertexId(at)],
                         arrangement[VertexId(after)]);
    };
    for (std::vector<std::uint32_t>& walk : strands) {
        if (walk.size() > 2 && walk.front() == walk.back()) {
            // A closed strand has no distinguished first vertex, so put one that
            // survives the simplification there before simplifying.
            walk.pop_back();
            std::size_t corner = 0;
            while (corner < walk.size() &&
                   removable(walk[(corner + walk.size() - 1) % walk.size()], walk[corner],
                             walk[(corner + 1) % walk.size()])) {
                ++corner;
            }
            std::rotate(walk.begin(),
                        walk.begin() + static_cast<std::ptrdiff_t>(corner % walk.size()),
                        walk.end());
            walk.push_back(walk.front());
        }
        std::vector<ResultPoint> vertices;
        vertices.reserve(walk.size());
        vertices.emplace_back(arrangement[VertexId(walk.front())]);
        for (std::size_t i = 1; i + 1 < walk.size(); ++i) {
            if (!removable(walk[i - 1], walk[i], walk[i + 1])) {
                vertices.emplace_back(arrangement[VertexId(walk[i])]);
            }
        }
        vertices.emplace_back(arrangement[VertexId(walk.back())]);
        pieces.emplace_back(ResultPolyline(std::move(vertices)));
    }
 
    for (const ResultRegion& region : regularizedCellsFromKeep<ResultPoint>(arrangement, keep)) {
        pieces.emplace_back(region);
    }
    return pieces;
}

template <class ResultPoint, class ShapeA, class ShapeB>
PolygonSet<ResultPoint> regularizedSymmetricDifference(const ShapeA& a, const ShapeB& b) {
    return regularizedBoolean<ResultPoint>(a, b, [](bool inA, bool inB) { return inA != inB; });
}
 
}  // namespace detail
 
// Out-of-line: the boolean operations are declared in shape/polygon.hpp and
// shape/polygonwithholes.hpp, which precede this header in the layering, but
// they can only be defined once Triangulation is visible.
 
template <class PointType_, class TLabel>
template <class ResultNumber>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularized() const {
    using ResultPoint = Point<ResultNumber, typename PointType_::LabelType>;
    // Nothing with area, so closure(A°) is empty and has no pieces at all.
    if (isDegenerate()) {
        return {};
    }
    // Already the closure of its own interior: rebuilding it through the
    // arrangement could only cost time and shed collinear vertices.
    if (isRegular()) {
        return PolygonSet<ResultPoint>(PolygonWithHoles<ResultPoint>(*this));
    }
    // One arrangement over this one boundary. Every cell of it is inside the
    // region or outside it, and a slit — having no cell of its own — survives
    // in neither.
    const std::array<PolygonWithHoles, 1> self{*this};
    return regularizedUnionOf<ResultPoint>(self);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonSet<PointType_, TLabel>::regularized() const {
    using ResultPoint = Point<ResultNumber, typename PointType_::LabelType>;
    // The components meet at finitely many points at most, so no slit runs
    // between two of them and each regularizes on its own. The pieces of
    // different components keep the disjoint interiors their components had.
    std::vector<PolygonWithHoles<ResultPoint>> pieces;
    pieces.reserve(components_.size());
    for (const auto& component : components_) {
        for (const auto& piece : component.template regularized<ResultNumber>()) {
            pieces.push_back(piece);
        }
    }
    return PolygonSet<ResultPoint>(std::move(pieces));
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonConcept OtherPolygon>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::difference(const OtherPolygon& other) const {
    return detail::regularizedDifference<Point<ResultNumber, typename PointType_::LabelType>>(*this,
                                                                                              other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::difference(const OtherConvex& other) const {
    return this->template difference<ResultNumber>(other.asPolygon());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::difference(const OtherTriangle& other) const {
    return this->template difference<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::difference(const OtherRectangle& other) const {
    return this->template difference<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonWithHolesConcept OtherRegion>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::difference(const OtherRegion& other) const {
    return detail::regularizedDifference<Point<ResultNumber, typename PointType_::LabelType>>(*this,
                                                                                              other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonSetConcept OtherSet>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::difference(const OtherSet& other) const {
    return detail::regularizedDifference<Point<ResultNumber, typename PointType_::LabelType>>(*this,
                                                                                              other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonConcept OtherPolygon>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::difference(const OtherPolygon& other) const {
    return detail::regularizedDifference<Point<ResultNumber, typename PointType_::LabelType>>(*this,
                                                                                              other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::difference(const OtherConvex& other) const {
    return this->template difference<ResultNumber>(other.asPolygon());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::difference(const OtherTriangle& other) const {
    return this->template difference<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::difference(const OtherRectangle& other) const {
    return this->template difference<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonWithHolesConcept OtherRegion>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::difference(const OtherRegion& other) const {
    return detail::regularizedDifference<Point<ResultNumber, typename PointType_::LabelType>>(*this,
                                                                                              other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonSetConcept OtherSet>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::difference(const OtherSet& other) const {
    return detail::regularizedDifference<Point<ResultNumber, typename PointType_::LabelType>>(*this,
                                                                                              other);
}
 
// The three bounded convex regions take a difference against every region
// through their polygon spelling. Unlike the symmetric three, a difference has
// no higher-ranked operand to forward to — `A ∖ B` is not `B ∖ A` — so each of
// them states it against all six, and `asPolygon` is what hands the pair to the
// engine. The conversion costs nothing, the vertices already being in canonical
// polygon order, and it is the same one `detail::booleanOperand` makes for a
// set's operands.
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonalRegionConcept OtherRegion>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Rectangle<PointType_, TLabel>::difference(const OtherRegion& other) const {
    return asPolygon().template difference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonalRegionConcept OtherRegion>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Triangle<PointType_, TLabel>::difference(const OtherRegion& other) const {
    return asPolygon().template difference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonalRegionConcept OtherRegion>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Convex<PointType_, TLabel>::difference(const OtherRegion& other) const {
    return asPolygon().template difference<ResultNumber>(other);
}
 
// ---------------------------------------------------------------------------
// The unbounded subtrahends.
//
// A union or a symmetric difference with an unbounded operand is unbounded, and
// no set of regions can hold it. A *difference* against one is not: `A ∖ B` is
// contained in `A`, so it is bounded whenever the receiver is, whatever `B`
// covers — which is why these overloads exist where the symmetric ones cannot.
// The latitude is one-sided, and visibly so: it is the receiver that has to be
// bounded, so `halfplane.difference(polygon)` is nowhere to be found.
//
// Being bounded is also what makes them computable. Only the part of `B` near
// `A` can matter, so `B` is clipped to a box strictly containing `A` — which
// leaves `A ∖ B` untouched, `A` being inside the box — and what comes back is a
// convex polygon the engine already takes. The two shortcuts below are the cases
// where the clip leaves nothing with area: a subtrahend with empty interior
// removes nothing, since `A° ∖ B` is dense in `A°` when `B` is nowhere dense, so
// the answer is `closure(A°)`, which is exactly @ref PolygonWithHoles::regularized.
//
// Everything else forwards: the convex three through their polygon spelling as
// they do for a bounded operand, and a polygon through the region spelling of
// itself, so the clip is written once.
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::difference(const OtherIntersection& other) const {
    using ExactNumber = detail::region_exact_number_t<typename OtherIntersection::NumberType>;
    // Nothing with area to keep, and no bounding rectangle to clip against.
    if (isDegenerate()) {
        return {};
    }
    if (other.isDegenerate()) {
        return this->template regularized<ResultNumber>();
    }
    const auto clipped = detail::regionClippedToBox(other, bbox());
    if (clipped.isDegenerate()) {
        return this->template regularized<ResultNumber>();
    }
    return this->template difference<ResultNumber>(clipped.template asConvex<ExactNumber>());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::difference(const OtherHalfplane& other) const {
    return this->template difference<ResultNumber>(other.asHalfplaneIntersection());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonSet<PointType_, TLabel>::difference(const OtherIntersection& other) const {
    using ExactNumber = detail::region_exact_number_t<typename OtherIntersection::NumberType>;
    if (isDegenerate()) {
        return {};
    }
    if (other.isDegenerate()) {
        return this->template regularized<ResultNumber>();
    }
    const auto clipped = detail::regionClippedToBox(other, bbox());
    if (clipped.isDegenerate()) {
        return this->template regularized<ResultNumber>();
    }
    return this->template difference<ResultNumber>(clipped.template asConvex<ExactNumber>());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonSet<PointType_, TLabel>::difference(const OtherHalfplane& other) const {
    return this->template difference<ResultNumber>(other.asHalfplaneIntersection());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::difference(const OtherIntersection& other) const {
    return asPolygonWithHoles().template difference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::difference(const OtherHalfplane& other) const {
    return asPolygonWithHoles().template difference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Rectangle<PointType_, TLabel>::difference(const OtherIntersection& other) const {
    return asPolygon().template difference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Rectangle<PointType_, TLabel>::difference(const OtherHalfplane& other) const {
    return asPolygon().template difference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Triangle<PointType_, TLabel>::difference(const OtherIntersection& other) const {
    return asPolygon().template difference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Triangle<PointType_, TLabel>::difference(const OtherHalfplane& other) const {
    return asPolygon().template difference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Convex<PointType_, TLabel>::difference(const OtherIntersection& other) const {
    return asPolygon().template difference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Convex<PointType_, TLabel>::difference(const OtherHalfplane& other) const {
    return asPolygon().template difference<ResultNumber>(other);
}
 
// The three bounded convex regions unite among themselves through the polygon
// engine, each pair on the higher-ranked of its two operands. The union of two
// convex shapes is not convex in general, so there is nothing a convex operand
// could contribute that its outline does not; going through `asPolygon` is the
// same conversion `detail::booleanOperand` makes for a set's operands, and it
// costs nothing, the vertices already being in canonical polygon order.
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Rectangle<PointType_, TLabel>::regularizedUnion(const OtherRectangle& other) const {
    return asPolygon().template regularizedUnion<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Triangle<PointType_, TLabel>::regularizedUnion(const OtherTriangle& other) const {
    return asPolygon().template regularizedUnion<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Triangle<PointType_, TLabel>::regularizedUnion(const OtherRectangle& other) const {
    return asPolygon().template regularizedUnion<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Convex<PointType_, TLabel>::regularizedUnion(const OtherConvex& other) const {
    return asPolygon().template regularizedUnion<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Convex<PointType_, TLabel>::regularizedUnion(const OtherTriangle& other) const {
    return asPolygon().template regularizedUnion<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Convex<PointType_, TLabel>::regularizedUnion(const OtherRectangle& other) const {
    return asPolygon().template regularizedUnion<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonConcept OtherPolygon>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::regularizedUnion(const OtherPolygon& other) const {
    using ExactNumber = detail::Exact1DNumber<NumberType, typename OtherPolygon::NumberType>;
    using ExactPoint = Point<ExactNumber>;
    const std::array<Polygon<ExactPoint>, 2> operands{Polygon<ExactPoint>(*this),
                                                      Polygon<ExactPoint>(other)};
    return regularizedUnionOf<Point<ResultNumber, typename PointType_::LabelType>>(
        operands, true);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::regularizedUnion(const OtherConvex& other) const {
    return this->template regularizedUnion<ResultNumber>(other.asPolygon());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::regularizedUnion(const OtherTriangle& other) const {
    return this->template regularizedUnion<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::regularizedUnion(const OtherRectangle& other) const {
    return this->template regularizedUnion<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonWithHolesConcept OtherRegion>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::regularizedUnion(const OtherRegion& other) const {
    return detail::regularizedUnion<Point<ResultNumber, typename PointType_::LabelType>>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonConcept OtherPolygon>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedUnion(const OtherPolygon& other) const {
    return detail::regularizedUnion<Point<ResultNumber, typename PointType_::LabelType>>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedUnion(const OtherConvex& other) const {
    return this->template regularizedUnion<ResultNumber>(other.asPolygon());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedUnion(const OtherTriangle& other) const {
    return this->template regularizedUnion<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedUnion(const OtherRectangle& other) const {
    return this->template regularizedUnion<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonWithHolesConcept OtherRegion>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedUnion(const OtherRegion& other) const {
    return detail::regularizedUnion<Point<ResultNumber, typename PointType_::LabelType>>(*this, other);
}
 
// A symmetric difference is symmetric, so the three convex regions own exactly
// the same pairs among themselves that they own for the union — each on the
// higher-ranked of its two operands — and reach everything above them through
// the same rank forwarder. The polygon spelling is the engine's operand here too.
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Rectangle<PointType_, TLabel>::symmetricDifference(const OtherRectangle& other) const {
    return asPolygon().template symmetricDifference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Triangle<PointType_, TLabel>::symmetricDifference(const OtherTriangle& other) const {
    return asPolygon().template symmetricDifference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Triangle<PointType_, TLabel>::symmetricDifference(const OtherRectangle& other) const {
    return asPolygon().template symmetricDifference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Convex<PointType_, TLabel>::symmetricDifference(const OtherConvex& other) const {
    return asPolygon().template symmetricDifference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Convex<PointType_, TLabel>::symmetricDifference(const OtherTriangle& other) const {
    return asPolygon().template symmetricDifference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Convex<PointType_, TLabel>::symmetricDifference(const OtherRectangle& other) const {
    return asPolygon().template symmetricDifference<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonConcept OtherPolygon>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::symmetricDifference(const OtherPolygon& other) const {
    return detail::regularizedSymmetricDifference<Point<ResultNumber, typename PointType_::LabelType>>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::symmetricDifference(const OtherConvex& other) const {
    return this->template symmetricDifference<ResultNumber>(other.asPolygon());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::symmetricDifference(const OtherTriangle& other) const {
    return this->template symmetricDifference<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::symmetricDifference(const OtherRectangle& other) const {
    return this->template symmetricDifference<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonWithHolesConcept OtherRegion>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
Polygon<PointType_, TLabel>::symmetricDifference(const OtherRegion& other) const {
    return detail::regularizedSymmetricDifference<Point<ResultNumber, typename PointType_::LabelType>>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonConcept OtherPolygon>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::symmetricDifference(const OtherPolygon& other) const {
    return detail::regularizedSymmetricDifference<Point<ResultNumber, typename PointType_::LabelType>>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::symmetricDifference(const OtherConvex& other) const {
    return this->template symmetricDifference<ResultNumber>(other.asPolygon());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::symmetricDifference(const OtherTriangle& other) const {
    return this->template symmetricDifference<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::symmetricDifference(const OtherRectangle& other) const {
    return this->template symmetricDifference<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonWithHolesConcept OtherRegion>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::symmetricDifference(const OtherRegion& other) const {
    return detail::regularizedSymmetricDifference<Point<ResultNumber, typename PointType_::LabelType>>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonConcept OtherPolygon>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedIntersection(const OtherPolygon& other) const {
    return detail::regularizedIntersection<Point<ResultNumber, typename PointType_::LabelType>>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedIntersection(const OtherConvex& other) const {
    return this->template regularizedIntersection<ResultNumber>(other.asPolygon());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedIntersection(const OtherTriangle& other) const {
    return this->template regularizedIntersection<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedIntersection(const OtherRectangle& other) const {
    return this->template regularizedIntersection<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonWithHolesConcept OtherRegion>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedIntersection(const OtherRegion& other) const {
    return detail::regularizedIntersection<Point<ResultNumber, typename PointType_::LabelType>>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedIntersection(const OtherIntersection& other) const {
    using ExactNumber = detail::region_exact_number_t<typename OtherIntersection::NumberType>;
    // Neither a region without area nor a half-plane intersection without
    // interior can contribute to closure(A° ∩ B°).
    if (isDegenerate() || other.isDegenerate()) {
        return {};
    }
    // The clip only has to preserve A ∩ B, and A lies strictly inside the box.
    const auto clipped = detail::regionClippedToBox(other, bbox());
    if (clipped.isDegenerate()) {
        return {};
    }
    return this->template regularizedIntersection<ResultNumber>(clipped.template asConvex<ExactNumber>());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonWithHoles<PointType_, TLabel>::regularizedIntersection(const OtherHalfplane& other) const {
    return this->template regularizedIntersection<ResultNumber>(other.asHalfplaneIntersection());
}
 
 
// The literal intersection: the one operation here that keeps what
// regularization drops. It takes the same operand grid as
// regularizedIntersection and makes the same reductions onto the engine's own
// shapes — a bounded convex operand goes in as its polygon, an unbounded one is
// clipped against this region first — and differs only in the engine it calls.
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonConcept OtherPolygon>
std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType_::LabelType>>,
                         PolygonWithHoles<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherPolygon& other) const {
    return detail::literalIntersection<Point<ResultNumber, typename PointType_::LabelType>>(
        *this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType_::LabelType>>,
                         PolygonWithHoles<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherConvex& other) const {
    return this->template intersection<ResultNumber>(other.asPolygon());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType_::LabelType>>,
                         PolygonWithHoles<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherTriangle& other) const {
    return this->template intersection<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType_::LabelType>>,
                         PolygonWithHoles<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherRectangle& other) const {
    return this->template intersection<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonWithHolesConcept OtherRegion>
std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType_::LabelType>>,
                         PolygonWithHoles<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherRegion& other) const {
    return detail::literalIntersection<Point<ResultNumber, typename PointType_::LabelType>>(
        *this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType_::LabelType>>,
                         PolygonWithHoles<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherIntersection& other) const {
    using ExactNumber = detail::region_exact_number_t<typename OtherIntersection::NumberType>;
    if (empty() || other.empty()) {
        return {};
    }
    // The clip only has to preserve A ∩ B, and A lies strictly inside the box.
    // A clip that comes back without interior is kept and handed to the engine
    // as the segment or point it collapsed to: unlike the regularized answer,
    // a contact of that dimension is a piece of this one.
    const auto clipped = detail::regionClippedToBox(other, bbox());
    return detail::literalIntersection<Point<ResultNumber, typename PointType_::LabelType>>(
        *this, clipped.template asConvex<ExactNumber>());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType_::LabelType>>,
                         PolygonWithHoles<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonWithHoles<PointType_, TLabel>::intersection(const OtherHalfplane& other) const {
    return this->template intersection<ResultNumber>(other.asHalfplaneIntersection());
}
 
 
// ---------------------------------------------------------------------------
// The set of regions, where the four operations close.
//
// One definition per operation over every operand the engine takes, another set
// included. There is nothing per-operand about them: the engine asks a shape for
// its cut segments and for one containment test per cell, which a set answers
// like anything else. Note in particular that a set operand goes in whole rather
// than one component at a time — folding would build one arrangement per step.
 
template <class PointType_, class TLabel>
template <class ResultNumber, detail::SetBooleanOperandConcept OtherShape>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonSet<PointType_, TLabel>::difference(const OtherShape& other) const {
    return detail::regularizedDifference<Point<ResultNumber, typename PointType_::LabelType>>(
        *this, detail::booleanOperand(other));
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, detail::SetBooleanOperandConcept OtherShape>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonSet<PointType_, TLabel>::regularizedUnion(const OtherShape& other) const {
    return detail::regularizedUnion<Point<ResultNumber, typename PointType_::LabelType>>(
        *this, detail::booleanOperand(other));
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, detail::SetBooleanOperandConcept OtherShape>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonSet<PointType_, TLabel>::regularizedIntersection(const OtherShape& other) const {
    return detail::regularizedIntersection<Point<ResultNumber, typename PointType_::LabelType>>(
        *this, detail::booleanOperand(other));
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, detail::SetBooleanOperandConcept OtherShape>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonSet<PointType_, TLabel>::symmetricDifference(const OtherShape& other) const {
    return detail::regularizedSymmetricDifference<
        Point<ResultNumber, typename PointType_::LabelType>>(*this, detail::booleanOperand(other));
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonSet<PointType_, TLabel>::regularizedIntersection(const OtherIntersection& other) const {
    using ExactNumber = detail::region_exact_number_t<typename OtherIntersection::NumberType>;
    // Neither a set without area nor a half-plane intersection without interior
    // can contribute to closure(A° ∩ B°).
    if (isDegenerate() || other.isDegenerate()) {
        return {};
    }
    // The clip only has to preserve A ∩ B, and A lies strictly inside the box.
    const auto clipped = detail::regionClippedToBox(other, bbox());
    if (clipped.isDegenerate()) {
        return {};
    }
    return this->template regularizedIntersection<ResultNumber>(clipped.template asConvex<ExactNumber>());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
PolygonSet<PointType_, TLabel>::regularizedIntersection(const OtherHalfplane& other) const {
    return this->template regularizedIntersection<ResultNumber>(other.asHalfplaneIntersection());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, detail::SetBooleanOperandConcept OtherShape>
std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType_::LabelType>>,
                         PolygonWithHoles<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonSet<PointType_, TLabel>::intersection(const OtherShape& other) const {
    return detail::literalIntersection<Point<ResultNumber, typename PointType_::LabelType>>(
        *this, detail::booleanOperand(other));
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType_::LabelType>>,
                         PolygonWithHoles<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonSet<PointType_, TLabel>::intersection(const OtherIntersection& other) const {
    using ExactNumber = detail::region_exact_number_t<typename OtherIntersection::NumberType>;
    if (empty() || other.empty()) {
        return {};
    }
    // The clip only has to preserve A ∩ B, and A lies strictly inside the box.
    // A clip without interior still carries pieces of this answer; see
    // PolygonWithHoles::intersection(const OtherIntersection&) const.
    const auto clipped = detail::regionClippedToBox(other, bbox());
    return detail::literalIntersection<Point<ResultNumber, typename PointType_::LabelType>>(
        *this, clipped.template asConvex<ExactNumber>());
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
std::vector<std::variant<Point<ResultNumber, typename PointType_::LabelType>,
                         Polyline<Point<ResultNumber, typename PointType_::LabelType>>,
                         PolygonWithHoles<Point<ResultNumber, typename PointType_::LabelType>>>>
PolygonSet<PointType_, TLabel>::intersection(const OtherHalfplane& other) const {
    return this->template intersection<ResultNumber>(other.asHalfplaneIntersection());
}
 
}  // namespace pgl