separates.hpp

Go to the documentation of this file.

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

 
#include <algorithm>
#include <array>
#include <compare>
#include <cstddef>
#include <map>
#include <optional>
#include <stdexcept>
#include <limits>
#include <type_traits>
#include <utility>
#include <vector>
#include "shape/orientedline.hpp"
#include "shape/orientedsegment.hpp"
#include "implementation/orientation.hpp"
#include "predicates_helpers.hpp"
 
 
namespace pgl {
 
namespace detail {


template <PointConcept HitPoint>
struct LineHit {
    HitPoint u{};
    HitPoint v{};
    bool crossing = false;
};

template <PointConcept LinePoint, PointConcept APoint, PointConcept BPoint>
constexpr std::partial_ordering lineHitOrder(
    const LinePoint& from, const LinePoint& to,
    const LineHit<APoint>& a, const LineHit<BPoint>& b) {
    if (!b.crossing) {
        if (!a.crossing) {
            // Two exact points: order by projection on the line direction.
            return dotSign(from, to, b.u, a.u);
        }
        const auto reversed = lineHitOrder(from, to, b, a);
        return reversed < 0 ? std::partial_ordering::greater
             : reversed > 0 ? std::partial_ordering::less
                            : std::partial_ordering::equivalent;
    }
    // The side of edge b from which the directed line arrives at b's
    // traversal point.
    const auto arrival = crossSign(b.u, b.v, from, to);
    if (!a.crossing) {
        const auto side = orientationSign(b.u, b.v, a.u);
        if (side == 0) {
            return std::partial_ordering::equivalent;
        }
        return ((side > 0) != (arrival > 0)) ? std::partial_ordering::less
                                             : std::partial_ordering::greater;
    }
    // Substitute a's traversal point p = a.u + t (a.v - a.u), with
    // t = du / (du - dv), into the determinant against edge b and multiply
    // through by (du - dv): sign(det(b, p)) = sign(ov du - ou dv) sign(du - dv).
    const auto du = orientationDeterminant(from, to, a.u);
    const auto dv = orientationDeterminant(from, to, a.v);
    const auto ou = orientationDeterminant(b.u, b.v, a.u);
    const auto ov = orientationDeterminant(b.u, b.v, a.v);
    using Wide = promoted_number_t<std::common_type_t<decltype(du), decltype(ou)>>;
    const auto scaled = threeWay(detail::asNumber<Wide>(ov) * detail::asNumber<Wide>(du),
                                 detail::asNumber<Wide>(ou) * detail::asNumber<Wide>(dv));
    if (scaled == 0) {
        return std::partial_ordering::equivalent;
    }
    const bool side_positive = (scaled > 0) == (du > dv);
    return (side_positive != (arrival > 0)) ? std::partial_ordering::less
                                            : std::partial_ordering::greater;
}

template <PointConcept VertexPoint, PointConcept LinePoint>
constexpr bool interiorWedgeContainsDirection(
    const VertexPoint& vertex, const VertexPoint& previous, const VertexPoint& next,
    const LinePoint& from, const LinePoint& to) {
    const auto outgoing_side = crossSign(vertex, next, from, to);
    if (outgoing_side == 0 && dotSign(vertex, next, from, to) > 0) {
        return false;  // the direction continues along the outgoing edge
    }
    const auto incoming_side = crossSign(from, to, vertex, previous);
    if (incoming_side == 0 && dotSign(vertex, previous, from, to) > 0) {
        return false;  // the direction doubles back along the incoming edge
    }
    const auto wedge = orientationSign(vertex, next, previous);
    if (wedge > 0) {  // convex corner: wedge narrower than a halfplane
        return outgoing_side > 0 && incoming_side > 0;
    }
    if (wedge < 0) {  // reflex corner: wedge wider than a halfplane
        return outgoing_side > 0 || incoming_side > 0;
    }
    return outgoing_side > 0;  // straight vertex: open left halfplane
}

template <PointConcept LinePoint, PolygonConcept Polygon>
constexpr bool lineSectionSeparatesPolygon(
    const LinePoint& from, const LinePoint& to,
    const Polygon& other, bool bounded_above) {
    using PolygonPoint = typename Polygon::PointType;
    using Hit = LineHit<PolygonPoint>;
    const LineHit<LinePoint> low{from, from, false};
    const LineHit<LinePoint> high{to, to, false};
 
    const std::ptrdiff_t n = other.size();
    // Visits every point where the boundary meets the line; stops early when
    // the visitor returns true.
    const auto scan_line_hits = [&](auto&& visit) {
        auto side = orientationSign(from, to, other[0]);
        for (std::ptrdiff_t i = 0; i < n; ++i) {
            const auto next_side = orientationSign(from, to, other.get(i + 1));
            if (side == 0) {
                if (visit(Hit{other[i], other[i], false}, i)) {
                    return true;
                }
            } else if (next_side != 0 && side != next_side) {
                if (visit(Hit{other[i], other.get(i + 1), true}, i)) {
                    return true;
                }
            }
            side = next_side;
        }
        return false;
    };
 
    Hit opening{};
    bool found = false;
    scan_line_hits([&](const Hit& hit, std::ptrdiff_t i) {
        const bool opens_interior = hit.crossing
            ? crossSign(hit.u, hit.v, from, to) > 0
            : interiorWedgeContainsDirection(
                  hit.u, other.get(i - 1), other.get(i + 1), from, to);
        if (opens_interior &&
            lineHitOrder(from, to, hit, low) >= 0 &&
            (!bounded_above || lineHitOrder(from, to, hit, high) <= 0) &&
            (!found || lineHitOrder(from, to, hit, opening) < 0)) {
            opening = hit;
            found = true;
        }
        return false;
    });
    if (!found) {
        return false;
    }
 
    return scan_line_hits([&](const Hit& hit, std::ptrdiff_t) {
        return lineHitOrder(from, to, hit, opening) > 0 &&
               (!bounded_above || lineHitOrder(from, to, hit, high) <= 0);
    });
}

 
}  // namespace detail

 
template <class Number, class Label>
template<SegmentConcept OtherSegment>
constexpr bool Point<Number, Label>::separates(const OtherSegment& other) const {
    return other.interiorContains(*this);
}
 
template <class Number, class Label>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Point<Number, Label>::separates(const OtherOrientedSegment& other) const {
    return static_cast<Segment<typename OtherOrientedSegment::PointType>>(other).interiorContains(*this);
}
 
template <class Number, class Label>
template<LineConcept OtherLine>
constexpr bool Point<Number, Label>::separates(const OtherLine& other) const {
    return other.contains(*this);
}
 
template <class Number, class Label>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool Point<Number, Label>::separates(const OtherOrientedLine& other) const {
    return other.asLine().contains(*this);
}
 
template <class Number, class Label>
template<RayConcept OtherRay>
constexpr bool Point<Number, Label>::separates(const OtherRay& other) const {
    return other.interiorContains(*this);
}
 
template <class Number, class Label>
template<PointConcept OtherPoint>
constexpr bool Point<Number, Label>::separates(const OtherPoint&) const {
    return false;
}
 
template <class Number, class Label>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Point<Number, Label>::separates(const OtherHalfplane&) const {
    return false;
}
 
template <class Number, class Label>
template<DiskConcept OtherDisk>
constexpr bool Point<Number, Label>::separates(const OtherDisk&) const {
    return false;
}
 
// A point can only separate a 2D region when that region is degenerate to a
// segment; its convex hull then has size 2, whose whole segment is the
// region's boundary, so the point separates it iff that segment contains it.
template <class Number, class Label>
template<ConvexConcept OtherConvex>
constexpr bool Point<Number, Label>::separates(const OtherConvex& other) const {
    return other.size() == 2 &&
           Segment<typename OtherConvex::PointType>(other[0], other[1]).contains(*this);
}
 
template <class Number, class Label>
template<RectangleConcept OtherRectangle>
constexpr bool Point<Number, Label>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return separates(static_cast<Convex<typename OtherRectangle::PointType>>(other));
}
 
template <class Number, class Label>
template<TriangleConcept OtherTriangle>
constexpr bool Point<Number, Label>::separates(const OtherTriangle& other) const {
    return separates(static_cast<Convex<typename OtherTriangle::PointType>>(other));
}
 
template <class Number, class Label>
template<PolygonConcept OtherPolygon>
constexpr bool Point<Number, Label>::separates(const OtherPolygon& other) const {
    return separates(Convex<typename OtherPolygon::PointType>(other.vertices()));
}

 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool Segment<PointType, LabelType>::separates(const OtherPoint&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool Segment<PointType, LabelType>::separates(const OtherSegment& other) const {
    using Coordinate = detail::sign_coordinate_t<NumberType, typename OtherSegment::NumberType>;
    const auto& a = min();
    const auto& b = max();
    const auto& c = other.min();
    const auto& d = other.max();
 
    // Four orientation signs over four endpoints, each endpoint in three of
    // them; see Segment::intersects for what the filtered wrappers buy.
    const auto fa = detail::filtered<Coordinate>(a);
    const auto fb = detail::filtered<Coordinate>(b);
    const auto fc = detail::filtered<Coordinate>(c);
    const auto fd = detail::filtered<Coordinate>(d);
    const auto s1 = detail::orientationSignOf(fa, fb, fc);
    const auto s2 = detail::orientationSignOf(fa, fb, fd);
    const auto s3 = detail::orientationSignOf(fc, fd, fa);
    const auto s4 = detail::orientationSignOf(fc, fd, fb);
 
    if (detail::allDecided(s1, s2, s3, s4)) {
        return s1.value() != s2.value() && s3.value() != s4.value();
    }
 
    const int cross = boundingBoxesCross(other);
    if (cross == 0) {
        return false;
    }
    if (cross == 2) {
        return true;
    }
    const auto d1 = s1.value();
    const auto d2 = s2.value();
    if (d1 == 0 && d2 == 0) {
        return other.interiorContains(a) && other.interiorContains(b);
    }
    const bool other_endpoints_are_on_strictly_opposite_sides =
        d1 != 0 && d2 != 0 && d1 != d2;
    if (!other_endpoints_are_on_strictly_opposite_sides) {
        return false;
    }
    const auto d3 = s3.value();
    if (d3 == 0 && other.containsCollinear(a)) {
        return true;
    }
    const auto d4 = s4.value();
    if (d4 == 0 && other.containsCollinear(b)) {
        return true;
    }
    return d3 != 0 && d4 != 0 && d3 != d4;
}
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Segment<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return separates(other.asSegment());
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool Segment<PointType, LabelType>::separates(const OtherLine& other) const {
    return intersects(other);
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool Segment<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return intersects(other);
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool Segment<PointType, LabelType>::separates(const OtherRay& other) const {
    // The segment splits the ray only when it meets the ray ahead of the
    // source: a piece then survives between the source and the segment, and
    // another runs to infinity. If the source lies on the segment, the near
    // piece is empty and the ray stays connected.
    return intersects(other) && !contains(other.source());
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool Segment<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return other.interiorsIntersect(*this) && !other.interiorContains(min()) && !other.interiorContains(max());
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Segment<PointType, LabelType>::separates(const OtherHalfplane& other) const {
    return other.separates(*this);
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool Segment<PointType, LabelType>::separates(const OtherTriangle& other) const {
    if (other.isDegenerate()) {
        // A collinear triangle is the segment or point it spans: the segment
        // overload answers for the former, and a point is never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
    if (isDegenerate()) {
        return false;
    }
    if (other.interiorContains(min()) || other.interiorContains(max())) {
        return false;
    }
 
    const auto triangle_edges = other.edges();
    for (const auto& edge : triangle_edges) {
        if (parallel(edge) &&
            (edge.contains(min()) ||
             edge.contains(max()) ||
             contains(edge.min()) ||
             contains(edge.max()))) {
            return false;
        }
    }
 
    int boundary_contact_count = 0;
    if (other.boundaryContains(min())) {
        ++boundary_contact_count;
    }
    if (other.boundaryContains(max()) && max() != min()) {
        ++boundary_contact_count;
    }
 
    const auto triangle_vertices = other.vertices();
    for (const auto& vertex : triangle_vertices) {
        if (contains(vertex) && vertex != min() && vertex != max()) {
            ++boundary_contact_count;
        }
    }
 
    for (const auto& edge : triangle_edges) {
        if (interiorsIntersect(edge)) {
            ++boundary_contact_count;
        }
    }
 
    return boundary_contact_count >= 2;
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool Segment<PointType, LabelType>::separates(const OtherConvex& other) const {
    // The segment separates the polygon iff its intersection with the polygon
    // is a true chord through the interior — neither endpoint lies strictly
    // inside (otherwise the segment ends midway and leaves a slit, not a
    // split) and the segment actually crosses the interior (otherwise it
    // lies along a single boundary edge). Both checks are O(log n).
    if (other.isDegenerate()) {
        // A convex with fewer than three vertices is the segment or point it
        // spans: the segment overload answers for the former, and a point is
        // never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
    return !isDegenerate()
           && !other.interiorContains(min())
           && !other.interiorContains(max())
           && other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool Segment<PointType, LabelType>::separates(const OtherPolygon& other) const {
    if (other.isDegenerate()) {
        // A polygon that collapses to a segment is cut as that segment is; a
        // point is never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
    if (isDegenerate()) {
        return false;
    }
    // The segment separates the polygon iff it covers whole some connected
    // component of supporting-line ∩ interior; see the engine for details.
    return detail::lineSectionSeparatesPolygon(min(), max(), other, true);
}
 
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool Segment<PointType, LabelType>::separates(const OtherDisk& other) const {
    return    !other.interiorContains(min())
           && !other.interiorContains(max())
           && other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
constexpr bool Segment<PointType, LabelType>::separates(const Shape<PointType>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}

 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherSegment& other) const {
    return !other.isDegenerate() && !contains(other.min()) && !contains(other.max()) && intersects(other);
}
 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherPoint&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return separates(static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherLine& other) const {
    return !other.isDegenerate() && intersects(other);
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return separates(other.asLine());
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherRay& other) const {
    // Removing the triangle splits the ray iff the ray reaches it with its
    // source outside: the far end runs to infinity (always outside), so any
    // contact -- including a tangential touch of the boundary -- leaves a piece
    // on each side.
    return !other.isDegenerate() && !isDegenerate()
           && !contains(other.source()) && intersects(other);
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherHalfplane&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    // A non-degenerate rectangle is its convex view; use the general convex-body
    // algorithm so the result matches separates(Convex) exactly. The previous
    // bespoke edge/vertex formula diverged from it.
    return separates(other.asConvex());
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherTriangle& other) const {
    // See separates(Rectangle): defer to the general convex-body algorithm.
    return separates(other.asConvex());
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherConvex& other) const {
    if (other.isDegenerate()) {
        // A convex with fewer than three vertices is the segment or point it
        // spans: the segment overload answers for the former, and a point is
        // never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
    if (isDegenerate()) {
        // What this degenerate shape spans does the cutting: a segment answers
        // as one, and a point disconnects nothing.
        if (const auto spanned = getIfSegment()) {
            return spanned->separates(other);
        }
        return false;
    }
 
    int interior_vertex_count = 0;
    for (const auto& vertex : vertices()) {
        if (other.interiorContains(vertex)) {
            interior_vertex_count++;
            if (interior_vertex_count >= 2) {
                return false;
            }
        }
    }
    if (interior_vertex_count == 1) {
        for (const auto& edge : edges()) {
            if (edge.separates(other)) {
                return true;
            }
        }
        return false;
    }
 
    int separating_edge_count = 0;
    for (const auto& edge : edges()) {
        if (edge.separates(other)) {
            separating_edge_count++;
        }
    }
    return separating_edge_count >= 2;
}
 
// Disk: specialized rather than delegated to Convex. For a triangle, removing it
// disconnects the disk iff the boundary crosses the circle >= 4 times. With k =
// vertices strictly inside the disk and c = edges carrying a full chord, the
// crossing count is fully determined by (k, c): a triangle can never weave in and
// out four times via single crossings, so the answer is a closed-form table:
//   k == 0 -> separates iff c >= 2
//   k == 1 -> separates iff c >= 1
//   k >= 2 -> never (the disk bulges out on a single arc)
// This avoids building a Convex and short-circuits on k before any chord test.
// A degenerate disk is a point and disconnects nothing; a collinear triangle is
// the segment it spans, which answers for it.
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherDisk& other) const {
    if (other.isDegenerate()) {
        return false;
    }
    if (isDegenerate()) {
        // What this degenerate shape spans does the cutting: a segment answers
        // as one, and a point disconnects nothing.
        if (const auto spanned = getIfSegment()) {
            return spanned->separates(other);
        }
        return false;
    }
 
    int interior_vertex_count = 0;
    for (const auto& vertex : vertices()) {
        interior_vertex_count += other.interiorContains(vertex) ? 1 : 0;
    }
    if (interior_vertex_count >= 2) {
        return false;
    }
 
    const int needed = interior_vertex_count == 0 ? 2 : 1;
    int separating_edge_count = 0;
    for (const auto& edge : edges()) {
        if (edge.separates(other) && ++separating_edge_count >= needed) {
            return true;
        }
    }
    return false;
}
 
template <class PointType, class LabelType>
constexpr bool Triangle<PointType, LabelType>::separates(const Shape<PointType>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}

 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherPoint& other) const {
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherSegment& other) const {
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return this->asSegment().separates(other.asSegment());
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherLine& other) const {
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherRay& other) const {
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherTriangle& other) const {
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherHalfplane& other) const {
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherConvex& other) const {
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherDisk& other) const {
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const Shape<PointType>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}

 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool Line<PointType, LabelType>::separates(const OtherPoint&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool Line<PointType, LabelType>::separates(const OtherSegment& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
    const auto first_side = orientationSign(min(), max(), other.min());
    const auto second_side = orientationSign(min(), max(), other.max());
    return first_side != 0 && second_side != 0 && first_side != second_side;
}
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Line<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return separates(other.asSegment());
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool Line<PointType, LabelType>::separates(const OtherLine& other) const {
    return !isDegenerate() && !other.isDegenerate() &&
           intersects(other) && !collinear(other);
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool Line<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return separates(other.asLine());
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool Line<PointType, LabelType>::separates(const OtherRay& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
    const auto source_side = orientationSign(min(), max(), other.source());
    const auto target_side = orientationSign(min(), max(), other.target());
    if (source_side == 0) {
        return false;
    }
    return target_side == 0 || target_side != source_side;
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool Line<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool Line<PointType, LabelType>::separates(const OtherTriangle& other) const {
    if (isDegenerate()) {
        return false;
    }
    return other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Line<PointType, LabelType>::separates(const OtherHalfplane& other) const {
    return other.separates(*this);
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool Line<PointType, LabelType>::separates(const OtherConvex& other) const {
    if (isDegenerate()) {
        return false;
    }
    return other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool Line<PointType, LabelType>::separates(const OtherDisk& other) const {
    if (isDegenerate()) {
        return false;
    }
    return other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
constexpr bool Line<PointType, LabelType>::separates(const Shape<PointType>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}

 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherPoint& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherSegment& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherLine& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherRay& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherTriangle& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherHalfplane& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherConvex& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherDisk& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
constexpr bool OrientedLine<PointType, LabelType>::separates(const Shape<PointType>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}

 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool Ray<PointType, LabelType>::separates(const OtherPoint&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool Ray<PointType, LabelType>::separates(const OtherSegment& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
 
    const auto first_side = orientationSign(source(), target(), other.min());
    const auto second_side = orientationSign(source(), target(), other.max());
    if (first_side == 0 && second_side == 0) {
        return false;
    }
    if (other.interiorContains(source())) {
        return true;
    }
    if (first_side == 0 || second_side == 0 || first_side == second_side) {
        return false;
    }
 
    // The segment crosses the ray's supporting line at a point interior to the
    // segment; it cuts the segment in two exactly when that crossing lies on the
    // ray itself (anywhere from the source onward, not just up to the target).
    return intersects(other);
}
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Ray<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return separates(other.asSegment());
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool Ray<PointType, LabelType>::separates(const OtherLine& other) const {
    return !isDegenerate() && !other.isDegenerate() &&
           intersects(other) && !collinear(other);
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool Ray<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return separates(other.asLine());
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool Ray<PointType, LabelType>::separates(const OtherRay& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
    if (collinear(other)) {
        return false;
    }
    if (other.interiorContains(source())) {
        return true;
    }
 
    const auto other_source_side = orientationSign(source(), target(), other.source());
    const auto other_target_side = orientationSign(source(), target(), other.target());
    if (other_source_side == 0) {
        return false;
    }
    if (other_target_side != 0 && other_target_side == other_source_side) {
        return false;
    }
 
    const auto source_side = orientationSign(other.source(), other.target(), source());
    const auto target_side = orientationSign(other.source(), other.target(), target());
    return source_side != 0 && (target_side == 0 || target_side != source_side);
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Ray<PointType, LabelType>::separates(const OtherHalfplane& other) const {
    return other.separates(*this);
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool Ray<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return !other.interiorContains(source()) && other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool Ray<PointType, LabelType>::separates(const OtherTriangle& other) const {
    if (isDegenerate() || other.isDegenerate() || other.interiorContains(source())) {
        return false;
    }
 
    const auto triangle_edges = other.edges();
    for (const auto& edge : triangle_edges) {
        if (parallel(edge) && (contains(edge.min()) || contains(edge.max()))) {
            return false;
        }
    }
 
    int boundary_contact_count = 0;
    if (other.boundaryContains(source())) {
        ++boundary_contact_count;
    }
    const auto triangle_vertices = other.vertices();
    for (const auto& vertex : triangle_vertices) {
        if (contains(vertex) && vertex != source()) {
            ++boundary_contact_count;
        }
    }
 
    for (const auto& edge : triangle_edges) {
        if (interiorsIntersect(edge)) {
            ++boundary_contact_count;
        }
    }
 
    return boundary_contact_count >= 2;
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool Ray<PointType, LabelType>::separates(const OtherConvex& other) const {
    if (other.isDegenerate()) {
        // A convex with fewer than three vertices is the segment or point it
        // spans: the segment overload answers for the former, and a point is
        // never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
 
    if (other.interiorContains(source())) {
        return false;
    }
    return other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool Ray<PointType, LabelType>::separates(const OtherDisk& other) const {
    // The disk is convex, so removing the ray disconnects it exactly when the ray
    // runs clear through as a full chord: the source must not be strictly inside
    // (an interior source leaves only a slit, which stays connected) and the ray's
    // interior must reach the disk's interior, after which the half-infinite far
    // end guarantees a second boundary crossing.
    if (other.isDegenerate()) {
        return false;
    }
    if (other.interiorContains(source())) {
        return false;
    }
    return other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool Ray<PointType, LabelType>::separates(const OtherPolygon& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
    // Same chord criterion as Segment::separates(Polygon) with the far end at
    // infinity: only the source can leave a component end uncovered, so the
    // covered-component search is bounded below by the source alone.
    return detail::lineSectionSeparatesPolygon(source(), target(), other, false);
}
 
template <class PointType, class LabelType>
constexpr bool Ray<PointType, LabelType>::separates(const Shape<PointType>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}

 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherRectangle& other) const {
    // The empty set disconnects nothing and cannot be disconnected, and the
    // inverted corners of an empty rectangle can pass the tests below, so both
    // branches that reach true rule emptiness out. The check trails the
    // geometry rather than guarding the function because false is the common
    // answer, and that path then never pays for it.
    const bool splits_horizontally =
        other.min().x() < min().x() &&
        max().x() < other.max().x() &&
        !(other.min().y() < min().y()) &&
        !(max().y() < other.max().y());
    if (splits_horizontally) {
        return !empty() && !other.empty();
    }
    const bool splits_vertically =
        other.min().y() < min().y() &&
        max().y() < other.max().y() &&
        !(other.min().x() < min().x()) &&
        !(max().x() < other.max().x());
    return splits_vertically && !empty() && !other.empty();
}
 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherPoint&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherLine& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return intersects(other) && !contains(other);
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return separates(other.asLine());
}
 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherSegment& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    // Both endpoints outside the closed rectangle and the segment touching it
    // anywhere (boundary contact included) leave a piece on each side.
    return !other.isDegenerate() && !contains(other.min()) && !contains(other.max()) && intersects(other);
}
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return separates(static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherRay& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return !other.isDegenerate() &&
           !contains(other.source()) &&
           intersects(other);
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherHalfplane& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherTriangle& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    if (other.isDegenerate()) {
        // A collinear triangle is the segment or point it spans: the segment
        // overload answers for the former, and a point is never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
 
    const auto target_edges = other.edges();
    if (contains(other.a()) && separates(target_edges[1])) {
        return true;
    }
    if (contains(other.b()) && separates(target_edges[2])) {
        return true;
    }
    if (contains(other.c()) && separates(target_edges[0])) {
        return true;
    }
 
    int separated_edges = 0;
    for (const auto& edge : target_edges) {
        separated_edges += separates(edge) ? 1 : 0;
    }
    return separated_edges >= 2;
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherConvex& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return asConvex().separates(other);
}
 
// Disk: specialized rather than delegated to Convex (see Triangle's overload for
// the reasoning). A rectangle's four corners are concyclic, so a disk can never
// contain exactly the two corners of a diagonal; hence the answer is again a
// closed-form table in k = corners strictly inside and c = full-chord edges:
//   k == 0      -> separates iff c >= 2
//   k in {1, 2} -> separates iff c >= 1   (k == 2 is always two *adjacent* corners)
//   k >= 3      -> never
// Degenerate inputs are undefined here; we report false.
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherDisk& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    if (other.isDegenerate()) {
        return false;
    }
    if (isDegenerate()) {
        // What this degenerate shape spans does the cutting: a segment answers
        // as one, and a point disconnects nothing.
        if (const auto spanned = getIfSegment()) {
            return spanned->separates(other);
        }
        return false;
    }
 
    int interior_vertex_count = 0;
    for (const auto& vertex : vertices()) {
        interior_vertex_count += other.interiorContains(vertex) ? 1 : 0;
    }
    if (interior_vertex_count >= 3) {
        return false;
    }
 
    const int needed = interior_vertex_count == 0 ? 2 : 1;
    int separating_edge_count = 0;
    for (const auto& edge : edges()) {
        if (edge.separates(other) && ++separating_edge_count >= needed) {
            return true;
        }
    }
    return false;
}
 
template <class PointType, class LabelType>
constexpr bool Rectangle<PointType, LabelType>::separates(const Shape<PointType>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}

 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherPoint&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherLine& other) const {
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherSegment& other) const {
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherRay& other) const {
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherHalfplane& other) const {
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherTriangle& other) const {
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherConvex&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherDisk& other) const {
    // Removing a closed half-plane from a disk leaves the circular segment on the
    // far side of the boundary line, which is always a single connected piece, so
    // a half-plane never disconnects a (convex) disk.
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherPolygon& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
 
    if (other.size() < 3) {
        return false;
    }
 
    const ptrdiff_t m = other.size();
 
    ptrdiff_t start = 0;
    for (; start < m && !contains(other[start]) ; ++start) {
    }
 
    if (start == m) {
        return false; // No vertex in the halfplane
    }
    // Now we know that other[start] is in the halfplane
 
    int arcs = 0;
    pgl::OrientedLine<PointType> boundary = asOrientedLine();
    pgl::Line<typename OtherPolygon::PointType> leaving;
    start++; // Now other[start-1] is in the halfplane
    bool prev_in = true;
 
    for (ptrdiff_t i = start; i < start+m; ++i) {
        const bool cur_in = contains(other.get(i));
        if (prev_in && !cur_in) {
            // Just went outside
            leaving = pgl::Line<typename OtherPolygon::PointType>(other.get(i), other.get(i-1));
        }
        if (cur_in && !prev_in) {
            // Just came inside, must check order
            pgl::Line<typename OtherPolygon::PointType> entering(other.get(i), other.get(i-1));
            if (boundary.crossingOrder(entering, leaving) >= 0) {
                ++arcs;
                if (arcs >= 2) {
                    return true;
                }
            }
        }
 
        prev_in = cur_in;
    }
    return false;
 
}
 
template <class PointType, class LabelType>
constexpr bool Halfplane<PointType, LabelType>::separates(const Shape<PointType>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}
 
 
// ---------------------------------------------------------------------------
// Convex
 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool Convex<PointType, LabelType>::separates(const OtherPoint&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool Convex<PointType, LabelType>::separates(const OtherSegment& other) const {
    // A contact anywhere on the convex polygon -- boundary touch included --
    // leaves a piece of the segment on each side when both endpoints are
    // outside, so this gates on closed intersection, not interior crossing.
    return !other.isDegenerate() && !contains(other.min()) && !contains(other.max()) && intersects(other);
}
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Convex<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return separates(static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool Convex<PointType, LabelType>::separates(const OtherLine& other) const {
    return !other.isDegenerate() && intersects(other);
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool Convex<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return separates(static_cast<Line<typename OtherOrientedLine::PointType>>(other));
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool Convex<PointType, LabelType>::separates(const OtherRay& other) const {
    return !other.isDegenerate() && !isDegenerate() &&
           !contains(other.source()) && intersects(other);
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Convex<PointType, LabelType>::separates(const OtherHalfplane&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool Convex<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return separates(other.asConvex());
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool Convex<PointType, LabelType>::separates(const OtherTriangle& other) const {
    return separates(other.asConvex());
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool Convex<PointType, LabelType>::separates(const OtherConvex& other) const {
    if (other.isDegenerate()) {
        // A convex with fewer than three vertices is the segment or point it
        // spans: the segment overload answers for the former, and a point is
        // never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
    if (isDegenerate()) {
        // What this degenerate shape spans does the cutting: a segment answers
        // as one, and a point disconnects nothing.
        if (const auto spanned = getIfSegment()) {
            return spanned->separates(other);
        }
        return false;
    }
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
 
    if (other.size() <= 2*size()) {
        const ptrdiff_t m = other.size();
        int arcs = 0;
        bool prev_in = contains(other[m - 1]);
        for (ptrdiff_t i = 0; i < m; ++i) {
            const bool cur_in = contains(other[i]);
            if (!prev_in && !cur_in) {
                // Both endpoints are out, but maybe the edge went through this
                Segment<typename OtherConvex::PointType> edge(other.get(i-1),other[i]);
                if (intersects(edge)) {
                    ++arcs;
                    if (arcs >= 2) {
                        return true;
                    }
                }
            }
            else if (prev_in && !cur_in) {
                // Just went outside
                ++arcs;
                if (arcs >= 2) {
                    return true;
                }
            }
 
            prev_in = cur_in;
        }
    }
    else {
        const ptrdiff_t n = size();
        int arcs = 0;
        bool prev_in = other.interiorContains((*this)[n - 1]);
        for (ptrdiff_t i = 0; i < n; ++i) {
            const bool cur_in = other.interiorContains((*this)[i]);
            if (!prev_in && !cur_in) {
                // Both endpoints are out, but maybe the edge went through this
                Segment<typename OtherConvex::PointType> edge(get(i-1),(*this)[i]);
                if (edge.separates(other)) {
                    ++arcs;
                    if (arcs >= 2) {
                        return true;
                    }
                }
            }
            else if (prev_in && !cur_in) {
                // Just went outside
                ++arcs;
                if (arcs >= 2) {
                    return true;
                }
            }
 
            prev_in = cur_in;
        }
    }
 
    return false;
}
 
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool Convex<PointType, LabelType>::separates(const OtherDisk& other) const {
    if (other.isDegenerate()) {
        return false;
    }
    if (isDegenerate()) {
        // What this degenerate shape spans does the cutting: a segment answers
        // as one, and a point disconnects nothing.
        if (const auto spanned = getIfSegment()) {
            return spanned->separates(other);
        }
        return false;
    }
 
    // Removing the polygon disconnects the disk iff its boundary crosses the
    // circle at least four times (>= 2 in/out arcs). Each edge contributes:
    //   - one endpoint inside, one outside -> 1 crossing,
    //   - both endpoints outside and the edge carries a full chord -> 2,
    //   - otherwise 0.
    // (A segment with both endpoints outside meets the circle 0 or 2 times, so
    // an interior crossing is exactly a full chord.) Counting per edge avoids
    // the earlier scheme that ignored crossings split across single-crossing
    // edges -- which under-reported, e.g. a thin rhombus laid across the disk.
    int crossings = 0;
    for (const auto& edge : edgesView()) {
        const bool min_inside = other.interiorContains(edge.min());
        const bool max_inside = other.interiorContains(edge.max());
        if (min_inside != max_inside) {
            ++crossings;
        } else if (!min_inside && edge.separates(other)) {
            crossings += 2;
        }
        if (crossings >= 4) {
            return true;
        }
    }
 
    return false;
}

 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool Disk<PointType, LabelType>::separates(const OtherPoint&) const {
    return false;
}
 
template <class PointType, class LabelType>
template <SegmentConcept OtherSegment>
constexpr bool Disk<PointType, LabelType>::separates(const OtherSegment& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
    return !contains(other.min()) && !contains(other.max()) && intersects(other);
}
 
template <class PointType, class LabelType>
template <OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Disk<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return separates(other.asSegment());
}
 
template <class PointType, class LabelType>
template <LineConcept OtherLine>
constexpr bool Disk<PointType, LabelType>::separates(const OtherLine& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
    return intersects(other);
}
 
template <class PointType, class LabelType>
template <OrientedLineConcept OtherOrientedLine>
constexpr bool Disk<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return separates(other.asLine());
}
 
template <class PointType, class LabelType>
template <ConvexConcept OtherConvex>
constexpr bool Disk<PointType, LabelType>::separates(const OtherConvex& other) const {
    if (other.isDegenerate()) {
        // A convex with fewer than three vertices is the segment or point it
        // spans: the segment overload answers for the former, and a point is
        // never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
    if (isDegenerate()) {
        return false;
    }
 
    const ptrdiff_t m = other.size();
    int arcs = 0;
    bool prev_in = contains(other[m - 1]);
    for (ptrdiff_t i = 0; i < m; ++i) {
        const bool cur_in = contains(other[i]);
        if (!prev_in && !cur_in) {
            // Both endpoints are out, but maybe the edge went through this
            Segment<typename OtherConvex::PointType> edge(other.get(i-1),other[i]);
            if (intersects(edge)) {
                ++arcs;
                if (arcs >= 2) {
                    return true;
                }
            }
        }
        else if (prev_in && !cur_in) {
            // Just went outside
            ++arcs;
            if (arcs >= 2) {
                return true;
            }
        }
 
        prev_in = cur_in;
    }
    return false;
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool Disk<PointType, LabelType>::separates(const OtherRay& other) const {
    return !contains(other.source()) && intersects(other);
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Disk<PointType, LabelType>::separates(const OtherHalfplane&) const {
    return false;  // A disk never separates a halfplane
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool Disk<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    int count = 0;
    for (int i = 0; i < 4; i++) {
        pgl::Segment<typename OtherRectangle::PointType> edge(other.get(i),other.get(i+1));
        if (separates(edge)) {
            count++;
            typename OtherRectangle::PointType opposite1 = other.get(i+2);
            typename OtherRectangle::PointType opposite2 = other.get(i+3);
            if (contains(opposite1) || contains(opposite2)) {
                return true;
            }
        }
    }
    return count >= 2;
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool Disk<PointType, LabelType>::separates(const OtherTriangle& other) const {
    int count = 0;
    for (int i = 0; i < 3; i++) {
        pgl::Segment<typename OtherTriangle::PointType> edge(other.get(i),other.get(i+1));
        if (separates(edge)) {
            count++;
            typename OtherTriangle::PointType opposite = other.get(i+2);
            if (contains(opposite)) {
                return true;
            }
        }
    }
    return count >= 2;
}
 
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool Disk<PointType, LabelType>::separates(const OtherDisk&) const {
    return false;  // a disk never separates another disk
}
 
template <class PointType, class LabelType>
template <PointConcept OtherPoint>
constexpr bool Convex<PointType, LabelType>::separates(const Shape<OtherPoint>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}

 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherPoint&) const {
    return false;
}
 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherSegment& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
 
    std::optional<pgl::OrientedSegment<PointType>> minSeg, maxSeg;
 
    auto crossingOrder = [](const pgl::Segment<typename OtherSegment::PointType> &s, const pgl::OrientedSegment<PointType> &a, const pgl::OrientedSegment<PointType> &b) {
        bool gte = pgl::Triangle<PointType>(s[0],a[0],a[1]).interiorsIntersect(pgl::Triangle<PointType>(s[1],b[0],b[1]));
        bool lte = pgl::Triangle<PointType>(s[0],b[0],b[1]).interiorsIntersect(pgl::Triangle<PointType>(s[1],a[0],a[1]));
        if (gte && !lte) return std::partial_ordering::greater;
        if (lte && !gte) return std::partial_ordering::less;
        return std::partial_ordering::equivalent;
    };
 
    for (ptrdiff_t i = 0; i < (ptrdiff_t) size(); ++i) {
        pgl::OrientedSegment<PointType> edge(get(i), get(i+1));
        if (edge.separates(other) && !edge.collinear(other)) {
            auto h = edge.leftHalfplane();
            if (h.contains(other[0])) {
                if (!maxSeg || crossingOrder(other, *maxSeg, edge) < 0) {
                    if (other.contains(edge[0])) {
                        pgl::OrientedSegment<PointType> previous(get(i-1), get(i));
                        if (previous.collinear(other) && !h.contains(previous)) continue;
                        auto hprevious = previous.leftHalfplane();
                        if (!previous.collinear(other) && hprevious.contains(other[1]) && !hprevious.contains(edge[1])) continue;
                    }
                    else if (other.contains(edge[1])) {
                        pgl::OrientedSegment<PointType> next(get(i+1), get(i+2));
                        if (next.collinear(other) && !h.contains(next)) continue;
                        auto hnext = next.leftHalfplane();
                        if (!next.collinear(other) && hnext.contains(other[1]) && !hnext.contains(edge[0])) continue;
                    }
                    maxSeg = edge;
                    if (minSeg && maxSeg && crossingOrder(other, *minSeg, *maxSeg) <= 0) {
                        return true;
                    }
                }
            }
            else {
                if (!minSeg || crossingOrder(other, *minSeg, edge) > 0) {
                    if (other.contains(edge[0])) {
                        pgl::OrientedSegment<PointType> previous(get(i-1), get(i));
                        if (previous.collinear(other) && !h.contains(previous)) continue;
                        auto hprevious = previous.leftHalfplane();
                        if (!previous.collinear(other) && hprevious.contains(other[0]) && !hprevious.contains(edge[1])) continue;
                    }
                    else if (other.contains(edge[1])) {
                        pgl::OrientedSegment<PointType> next(get(i+1), get(i+2));
                        if (next.collinear(other) && !h.contains(next)) continue;
                        auto hnext = next.leftHalfplane();
                        if (!next.collinear(other) && hnext.contains(other[0]) && !hnext.contains(edge[0])) continue;
                    }
                    minSeg = edge;
                    if (minSeg && maxSeg && crossingOrder(other, *minSeg, *maxSeg) <= 0) {
                        return true;
                    }
                }
            }
        }
    }
    return false;
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherRay& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
 
    pgl::OrientedLine<typename OtherRay::PointType> ol(other[0],other[1]);
 
    for (ptrdiff_t i = 0; i < (ptrdiff_t) size(); ++i) {
        pgl::OrientedSegment<PointType> edge(get(i), get(i+1));
        if (edge.separates(other) && !edge.collinear(other)) {
            auto h = edge.leftHalfplane();
            if (!h.contains(other.source())) {
                if (other.contains(edge[0])) {
                    pgl::OrientedSegment<PointType> previous(get(i-1), get(i));
                    if (previous.collinear(other) && !h.contains(previous)) continue;
                    auto hprevious = previous.leftHalfplane();
                    if (!previous.collinear(other) && hprevious.contains(other[0]) && !hprevious.contains(edge[1])) continue;
                }
                else if (other.contains(edge[1])) {
                    pgl::OrientedSegment<PointType> next(get(i+1), get(i+2));
                    if (next.collinear(other) && !h.contains(next)) continue;
                    auto hnext = next.leftHalfplane();
                    if (!next.collinear(other) && hnext.contains(other[0]) && !hnext.contains(edge[0])) continue;
                }
                return true;
            }
        }
    }
    return false;
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherLine& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
    return intersects(other);
}
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return separates(other.asSegment());
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return separates(other.asLine());
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherHalfplane&) const {
    return false;  // A polygon never separates a halfplane
}
 
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherPolygon& other) const {
    return this->asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherPolygon& other) const {
    return this->asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool Line<PointType, LabelType>::separates(const OtherPolygon& other) const {
    if (isDegenerate()) {
        return false;
    }
    return other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherPolygon& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return asConvex().separates(other);
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherPolygon& other) const {
    return asConvex().separates(other);
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool Convex<PointType, LabelType>::separates(const OtherPolygon& other) const {
    if (other.isDegenerate()) {
        // A polygon that collapses to a segment is cut as that segment is; a
        // point is never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
    if (isDegenerate()) {
        // What this degenerate shape spans does the cutting: a segment answers
        // as one, and a point disconnects nothing.
        if (const auto spanned = getIfSegment()) {
            return spanned->separates(other);
        }
        return false;
    }
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
 
    // Removing the convex body C from the polygon P disconnects P iff some
    // connected component of C ∩ P touches ∂P in two or more pieces. Counting
    // boundary arcs (as the Convex overload does) is not enough here: a reflex
    // polygon can dip into C through several separate pockets while P \ C
    // stays connected. Instead, walk around ∂C: P is cut exactly when some
    // maximal arc of ∂C through the open interior of P joins two contacts
    // belonging to *different* components of ∂P ∩ C. This also catches a
    // convex body inside P that pinches ∂P at two isolated touch points.
    //
    // Events are the contacts of ∂P with ∂C, each labelled with the component
    // of ∂P ∩ C it belongs to and located on its convex edge by a line that
    // crosses the edge at the contact, so no intersection point is ever
    // constructed. Each event also records the local feature of ∂P at the
    // contact (the crossing edge, or the wedge at a polygon vertex), so that
    // whether the ∂C arc leaving a contact enters the interior of P is decided
    // by orientation signs at that contact alone. Stretches where the two
    // boundaries overlap are bounded by events of the same component and can
    // therefore never fire the test.
 
    using CommonNumber = std::common_type_t<typename PointType::NumberType,
                                            typename OtherPolygon::NumberType>;
    using CommonPoint = Point<CommonNumber>;
    using PosLine = Line<CommonPoint>;
 
    const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(other.size());
    const std::ptrdiff_t m = static_cast<std::ptrdiff_t>(size());
 
    const auto common = [](const auto& p) {
        return CommonPoint(static_cast<CommonNumber>(p.x()),
                           static_cast<CommonNumber>(p.y()));
    };
 
    // A line through `at` perpendicular to convex edge k; it crosses that
    // edge exactly at `at`, placing a contact point that is known explicitly.
    const auto perpendicularAt = [&](const CommonPoint& at, std::ptrdiff_t k) {
        const CommonPoint a = common(get(k));
        const CommonPoint b = common(get(k + 1));
        return PosLine(at, CommonPoint(at.x() - (b.y() - a.y()),
                                       at.y() + (b.x() - a.x())));
    };
 
    struct Event {
        std::ptrdiff_t cedge;  // convex edge carrying the contact
        PosLine where;         // crosses that edge at the contact point
        int component;         // component of ∂P ∩ C the contact belongs to
        CommonPoint a, b, c;   // ∂P at the contact: edge a->b, or wedge a->b->c
        bool atVertex;         // contact is the polygon vertex b
    };
    std::vector<Event> events;
 
    // The first (mode first) or last contact of polygon edge pa->pb with ∂C.
    // Contacts are ordered along pa->pb by where the convex edge's line
    // crosses it; that line is parallel to pa->pb only for a collinear
    // overlap, whose near and far ends are explicit endpoints. A contact
    // landing exactly on the run's vertex inside C is dropped when requested:
    // the vertex's own wedge events describe that contact better.
    const auto extremeContact = [&](const typename OtherPolygon::PointType& pa, const typename OtherPolygon::PointType& pb,
                                    bool first, bool skipAtRunVertex)
        -> std::optional<std::pair<std::ptrdiff_t, PosLine>> {
        const CommonPoint a = common(pa);
        const CommonPoint b = common(pb);
        const Segment<CommonPoint> s(a, b);
        const OrientedLine<CommonPoint> axis(a, b);
        std::optional<std::pair<std::ptrdiff_t, PosLine>> best;
        std::optional<PosLine> bestAlong;
        const auto consider = [&](std::ptrdiff_t k, const PosLine& where,
                                  const PosLine& along) {
            if (bestAlong) {
                const auto order = axis.crossingOrder(along, *bestAlong);
                const bool improves =
                    first ? order == std::partial_ordering::less
                          : order == std::partial_ordering::greater;
                if (!improves) {
                    return;
                }
            }
            best.emplace(k, where);
            bestAlong = along;
        };
        for (std::ptrdiff_t k = 0; k < m; ++k) {
            const Segment<CommonPoint> edge(common(get(k)), common(get(k + 1)));
            if (!edge.intersects(s)) {
                continue;
            }
            if (s.parallel(edge)) {
                const CommonPoint lo = std::max(s.min(), edge.min());
                const CommonPoint hi = std::min(s.max(), edge.max());
                const CommonPoint at = (first == (a < b)) ? lo : hi;
                const PosLine across = perpendicularAt(at, k);
                consider(k, across, across);
            } else {
                consider(k, PosLine(a, b),
                         PosLine(common(get(k)), common(get(k + 1))));
            }
        }
        if (best && skipAtRunVertex) {
            const CommonPoint inVertex = first ? b : a;
            const PosLine atInVertex(
                inVertex, CommonPoint(inVertex.x() - (b.y() - a.y()),
                                      inVertex.y() + (b.x() - a.x())));
            if (axis.crossingOrder(*bestAlong, atInVertex) ==
                std::partial_ordering::equivalent) {
                return std::nullopt;
            }
        }
        return best;
    };
 
    std::vector<char> in(static_cast<std::size_t>(n));
    std::ptrdiff_t start = -1;
    for (std::ptrdiff_t i = 0; i < n; ++i) {
        in[static_cast<std::size_t>(i)] = contains(other.get(i));
        if (!in[static_cast<std::size_t>(i)]) {
            start = i;
        }
    }
    if (start < 0) {
        return false;  // the whole polygon lies inside the convex body
    }
 
    int component = 0;
 
    // A polygon vertex inside C that lies on ∂C is a contact in its own
    // right; the pair of wedge events splits the ∂C arcs it sits between
    // (the two-touch pinch case) without affecting any other arc.
    const auto addVertexContact = [&](std::ptrdiff_t i) {
        const typename OtherPolygon::PointType vertex = other.get(i);
        if (!boundaryContains(vertex)) {
            return;
        }
        const CommonPoint v = common(vertex);
        for (std::ptrdiff_t k = 0; k < m; ++k) {
            const Segment<CommonPoint> edge(common(get(k)), common(get(k + 1)));
            if (edge.contains(v)) {
                const PosLine across = perpendicularAt(v, k);
                const Event event{k, across, component,
                                  common(other.get(i - 1)), v,
                                  common(other.get(i + 1)), true};
                events.push_back(event);
                events.push_back(event);
                return;
            }
        }
    };
 
    bool previous = false;  // in[start] is false
    for (std::ptrdiff_t i = start + 1; i <= start + n; ++i) {
        const bool current = in[static_cast<std::size_t>(i % n)];
        const auto u = other.get(i - 1);
        const auto v = other.get(i);
        if (current && !previous) {
            ++component;
            if (const auto contact = extremeContact(u, v, true, true)) {
                events.push_back({contact->first, contact->second, component,
                                  common(u), common(v), common(v), false});
            }
        } else if (!current && previous) {
            if (const auto contact = extremeContact(u, v, false, true)) {
                events.push_back({contact->first, contact->second, component,
                                  common(u), common(v), common(v), false});
            }
        } else if (!current && intersects(Segment<typename OtherPolygon::PointType>(u, v))) {
            // Both endpoints outside: the edge meets the convex body in a
            // single sub-segment or touch point, a component of its own.
            ++component;
            if (const auto contact = extremeContact(u, v, true, false)) {
                events.push_back({contact->first, contact->second, component,
                                  common(u), common(v), common(v), false});
            }
            if (const auto contact = extremeContact(u, v, false, false)) {
                events.push_back({contact->first, contact->second, component,
                                  common(u), common(v), common(v), false});
            }
        }
        if (current) {
            addVertexContact(i);
        }
        previous = current;
    }
 
    if (events.empty()) {
        return false;  // the boundaries never meet: nested or disjoint
    }
 
    std::sort(events.begin(), events.end(),
              [&](const Event& x, const Event& y) {
        if (x.cedge != y.cedge) {
            return x.cedge < y.cedge;
        }
        const OrientedLine<CommonPoint> edge(common(get(x.cedge)),
                                             common(get(x.cedge + 1)));
        return edge.crossingOrder(x.where, y.where) ==
               std::partial_ordering::less;
    });
 
    // Whether ∂C heads strictly into the interior of P as it leaves the
    // event's contact point, decided against the recorded ∂P feature there.
    const auto entersInterior = [&](const Event& event) {
        const std::ptrdiff_t k = event.cedge;
        const OrientedLine<CommonPoint> edgeLine(common(get(k)),
                                                 common(get(k + 1)));
        const bool atEdgeEnd =
            edgeLine.crossingOrder(event.where,
                                   perpendicularAt(common(get(k + 1)), k)) ==
            std::partial_ordering::equivalent;
        const CommonPoint from = common(get(k + (atEdgeEnd ? 1 : 0)));
        const CommonPoint to = common(get(k + (atEdgeEnd ? 2 : 1)));
        const CommonNumber dx = to.x() - from.x();
        const CommonNumber dy = to.y() - from.y();
        if (!event.atVertex) {
            // Contact interior to polygon edge a->b: P's interior is strictly
            // to its left, so the arc dives inside iff it heads left.
            return orientationSign(event.a, event.b,
                                   CommonPoint(event.a.x() + dx,
                                               event.a.y() + dy)) > 0;
        }
        // Contact at polygon vertex b: the direction must point strictly into
        // the interior wedge between the edges a->b and b->c.
        const CommonPoint probe(event.b.x() + dx, event.b.y() + dy);
        const auto sideIn = orientationSign(event.a, event.b, probe);
        const auto sideOut = orientationSign(event.b, event.c, probe);
        return orientationSign(event.a, event.b, event.c) >= 0
                   ? (sideIn > 0 && sideOut > 0)
                   : (sideIn > 0 || sideOut > 0);
    };
 
    const std::size_t count = events.size();
    for (std::size_t t = 0; t < count; ++t) {
        const Event& current = events[t];
        const Event& next = events[(t + 1) % count];
        if (current.component != next.component && entersInterior(current)) {
            return true;
        }
    }
    return false;
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool Disk<PointType, LabelType>::separates(const OtherPolygon& other) const {
    if (isDegenerate() || other.isDegenerate()) {
        return false;
    }
 
    // Convex fast path -- same arc count as Disk::separates(Convex): removing
    // the disk disconnects a convex polygon iff its boundary leaves the disk in
    // >= 2 separate arcs. Each vertex is inside or outside the disk (an exact
    // squared-distance test), and an edge with both endpoints outside still
    // contributes an arc if it touches the disk.
    if (other.isConvex()) {
        const std::ptrdiff_t m = static_cast<std::ptrdiff_t>(other.size());
        int arcs = 0;
        bool prev_in = contains(other[m - 1]);
        for (std::ptrdiff_t i = 0; i < m; ++i) {
            const bool cur_in = contains(other[i]);
            if (!prev_in && !cur_in) {
                const Segment<typename OtherPolygon::PointType> edge(other.get(i - 1), other[i]);
                if (intersects(edge)) {
                    ++arcs;
                }
            } else if (prev_in && !cur_in) {
                ++arcs;
            }
            if (arcs >= 2) {
                return true;
            }
            prev_in = cur_in;
        }
        return false;
    }
 
    // General simple polygon. The arc count above over-reports here: a reflex
    // polygon can wrap around the disk and reconnect two of its boundary arcs,
    // so connectivity of P \ D is decided directly -- still without ever
    // constructing the (irrational) circle crossings. Triangulate P and glue
    // the convex per-triangle answers:
    //   - Within one triangle T, the components of T \ D correspond one-to-one
    //     to the maximal runs of pieces of dT outside the disk that are not
    //     interrupted by a contact of the (closed) disk with dT: the disk's
    //     bite off the convex T disconnects it exactly at those contacts.
    //   - A triangulation edge keeps at most two pieces outside the disk (the
    //     whole edge, or the parts adjacent to each endpoint), and which of
    //     them exist depends only on the edge, never on the incident triangle.
    //     Both triangles sharing the edge therefore resolve a surviving piece
    //     to the same slot, and any crossing between adjacent triangles goes
    //     through such a piece (the complement of the closed disk is open
    //     along the edge, so a crossing point sits in a piece of positive
    //     length).
    // Union-find over the pieces then yields the components of P \ D: the disk
    // separates the polygon iff at least two classes remain. No pieces at all
    // means the polygon is swallowed by the disk, which does not separate it.
    //
    // Layering: Triangulation lives in algorithm/triangulation.hpp, which
    // pgl.hpp includes after this header. Reach it only through the dependent
    // call other.triangulation(), never by naming the class here.
    const auto mesh = other.triangulation();
    using Mesh = std::decay_t<decltype(mesh)>;
    using TriId = typename Mesh::TriId;
    using VertexId = typename Mesh::VertexId;
    using Seg = Segment<typename Mesh::PointType>;
 
    std::vector<std::size_t> parent;
    const auto findRoot = [&parent](std::size_t x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    };
    const auto unite = [&](std::size_t a, std::size_t b) {
        parent[findRoot(a)] = findRoot(b);
    };
 
    // Whether the disk swallows a vertex depends on the vertex alone, so the
    // test is paid once per vertex instead of once per incident triangle.
    std::vector<char> swallowed(mesh.vertexIndexBound(), 0);
    for (const VertexId v : mesh.vertexIds()) {
        swallowed[v.index()] = contains(mesh.getShape(v)) ? 1 : 0;
    }
 
    // A piece is named by the edge that carries it and which end it hangs from.
    // Handle order stands in for the lexicographic order of the endpoints: all
    // the shared slot needs is that both incident triangles name the ends the
    // same way, and either order does that.
    enum : int { WHOLE = 0, NEAR_LOW = 1, NEAR_HIGH = 2 };
    using EdgeKey = std::pair<VertexId, VertexId>;
    std::map<EdgeKey, std::array<std::ptrdiff_t, 3>> pieceSlots;
    const auto slotOf = [&](VertexId a, VertexId b, int piece) {
        const EdgeKey key = a < b ? EdgeKey{a, b} : EdgeKey{b, a};
        auto [it, inserted] =
            pieceSlots.try_emplace(key, std::array<std::ptrdiff_t, 3>{-1, -1, -1});
        std::ptrdiff_t& slot = it->second[static_cast<std::size_t>(piece)];
        if (slot < 0) {
            slot = static_cast<std::ptrdiff_t>(parent.size());
            parent.push_back(parent.size());
        }
        return static_cast<std::size_t>(slot);
    };
 
    constexpr std::size_t NO_PIECE = std::numeric_limits<std::size_t>::max();
    mesh.visitTriangles([&](TriId t) {
        // Walk dT once, emitting its outside pieces in boundary order and
        // recording every disk contact between them; gap-free consecutive
        // pieces (wrap-around included) bound the same component of T \ D.
        std::size_t firstPiece = NO_PIECE;
        std::size_t lastPiece = NO_PIECE;
        bool gapBeforeFirst = false;
        bool gap = false;  // contact seen since the last emitted piece
        const auto emit = [&](std::size_t slot) {
            if (lastPiece == NO_PIECE) {
                firstPiece = slot;
                gapBeforeFirst = gap;
            } else if (!gap) {
                unite(lastPiece, slot);
            }
            lastPiece = slot;
            gap = false;
        };
 
        // vertices(t)[i] is the vertex at getShape(t)[i], so side i still runs
        // from vertex i to vertex i + 1 as the value walk had it.
        const auto v = mesh.vertices(t);
        const bool in[3] = {swallowed[v[0].index()] != 0, swallowed[v[1].index()] != 0,
                            swallowed[v[2].index()] != 0};
        for (int i = 0; i < 3; ++i) {
            const VertexId a = v[i];
            const VertexId b = v[(i + 1) % 3];
            const int nearA = a < b ? NEAR_LOW : NEAR_HIGH;
            const int nearB = nearA == NEAR_LOW ? NEAR_HIGH : NEAR_LOW;
            if (in[i] && in[(i + 1) % 3]) {  // edge swallowed by the disk
                gap = true;
            } else if (in[i]) {              // walk leaves the disk mid-edge
                gap = true;
                emit(slotOf(a, b, nearB));
            } else if (in[(i + 1) % 3]) {    // walk enters the disk mid-edge
                emit(slotOf(a, b, nearA));
                gap = true;
            } else if (intersects(Seg(mesh.getShape(a), mesh.getShape(b)))) {
                emit(slotOf(a, b, nearA));   // contact strictly inside the edge
                gap = true;
                emit(slotOf(a, b, nearB));
            } else {                         // edge clear of the disk
                emit(slotOf(a, b, WHOLE));
            }
        }
        if (firstPiece != NO_PIECE && !gap && !gapBeforeFirst) {
            unite(lastPiece, firstPiece);
        }
    });
 
    std::size_t classes = 0;
    for (std::size_t i = 0; i < parent.size(); ++i) {
        if (parent[i] == i && ++classes >= 2) {
            return true;
        }
    }
    return false;
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return separates(other.asPolygon());
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherTriangle& other) const {
    return separates(other.asPolygon());
}
 
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherDisk& other) const {
    if (other.isDegenerate()) {
        return false;
    }
    if (isDegenerate()) {
        // What this degenerate shape spans does the cutting: a segment answers
        // as one, and a point disconnects nothing.
        if (const auto spanned = getIfSegment()) {
            return spanned->separates(other);
        }
        return false;
    }
 
    // Convex fast path -- same crossing count as Convex::separates(Disk):
    // removing a convex polygon disconnects the disk when the polygon boundary
    // crosses the circle >= 4 times (>= 2 in/out arcs). Per edge:
    //   - one endpoint inside, one outside -> 1 crossing,
    //   - both endpoints outside and the edge carries a full chord -> 2,
    //   - otherwise 0.
    if (isConvex()) {
        int crossings = 0;
        for (const auto& edge : edgesView()) {
            const bool min_inside = other.interiorContains(edge.min());
            const bool max_inside = other.interiorContains(edge.max());
            if (min_inside != max_inside) {
                ++crossings;
            } else if (!min_inside && edge.separates(other)) {
                crossings += 2;
            }
            if (crossings >= 4) {
                return true;
            }
        }
        return false;
    }
 
    // General simple polygon P. The crossing count above over-counts here: a
    // non-convex P can take several separate bites out of the disk D, and the
    // bites' crossings do not disconnect D \ P. Connectivity is decided
    // directly instead -- still without constructing the (irrational) circle
    // crossings:
    //   - Enclose both shapes in a box polygon and build its constrained
    //     Delaunay triangulation with P's edges as constraints. Every triangle
    //     then lies entirely inside or entirely outside P; the sides are told
    //     apart by a flood fill from a box corner that flips whenever it steps
    //     over a constrained edge.
    //   - For an outside triangle T, the piece D ∩ T \ P is convex minus a
    //     part of its own boundary, hence connected. Two pieces adjacent
    //     across a triangulation edge e belong to the same component of D \ P
    //     iff (D ∩ e) \ P is nonempty. Only e's endpoints can lie in P (a P
    //     vertex; box corners are outside D by construction and e's relative
    //     interior never meets P), so that reduces to: the open disk meets e,
    //     or the closed disk touches e in a single interior point.
    // Union-find over the outside triangles then yields the components of
    // D \ P; the polygon separates the disk iff pieces of D \ P fall into at
    // least two classes. Pieces pinched at a P vertex on the circle stay
    // apart, because the pinch point itself belongs to P.
    using Common = std::common_type_t<typename PointType::NumberType,
                                      typename OtherDisk::NumberType>;
    using CPoint = Point<Common>;
    using CSeg = Segment<CPoint>;
 
    // Bounding box of both shapes, inflated so the box strictly contains them.
    const auto dbox = other.bbox();
    Common xlo = static_cast<Common>(dbox.min().x());
    Common ylo = static_cast<Common>(dbox.min().y());
    Common xhi = static_cast<Common>(dbox.max().x());
    Common yhi = static_cast<Common>(dbox.max().y());
    for (std::size_t i = 0; i < size(); ++i) {
        const auto v = (*this)[i];
        const Common vx = static_cast<Common>(v.x());
        const Common vy = static_cast<Common>(v.y());
        xlo = vx < xlo ? vx : xlo;
        ylo = vy < ylo ? vy : ylo;
        xhi = vx > xhi ? vx : xhi;
        yhi = vy > yhi ? vy : yhi;
    }
    xlo = xlo - Common{1};
    ylo = ylo - Common{1};
    xhi = xhi + Common{1};
    yhi = yhi + Common{1};
    const Polygon<CPoint> box(std::vector<CPoint>{
        CPoint(xlo, ylo), CPoint(xhi, ylo), CPoint(xhi, yhi), CPoint(xlo, yhi)});
 
    // P's edges become the constraints. The mesh keeps them flagged as such, and
    // that flag is what tells the flood fill below where the inside/outside
    // classification flips -- no separate set of P's edges to consult.
    std::vector<CSeg> constraints;
    constraints.reserve(size());
    for (std::size_t i = 0; i < size(); ++i) {
        const auto u = (*this)[i];
        const auto w = get(static_cast<std::ptrdiff_t>(i) + 1);
        constraints.emplace_back(
            CPoint(detail::asNumber<Common>(u.x()), detail::asNumber<Common>(u.y())),
            CPoint(detail::asNumber<Common>(w.x()), detail::asNumber<Common>(w.y())));
    }
 
    // Layering: Triangulation lives in algorithm/triangulation.hpp, which
    // pgl.hpp includes after this header. Reach it only through the dependent
    // call box.triangulation(constraints), never by naming the class here.
    const auto mesh = box.triangulation(constraints);
    using TriId = typename std::decay_t<decltype(mesh)>::TriId;
    // A handle's index is a slot of a table this size, so the per-triangle
    // bookkeeping below is plain vectors and crossing an edge is an array read.
    // The bound counts the storage rather than the triangles: the slots it holds
    // beyond them stay untouched, their state left at "unvisited".
    const std::size_t n = mesh.triangleIndexBound();
 
    // Flood fill from a box corner: -1 unvisited, 0 inside P, 1 outside P.
    std::vector<signed char> state(n, -1);
    std::vector<TriId> pending;
    for (const TriId t : mesh.incidentTriangles(mesh.getId(CPoint(xlo, ylo)))) {
        if (state[t.index()] < 0) {
            state[t.index()] = 1;
            pending.push_back(t);
        }
    }
    while (!pending.empty()) {
        const TriId t = pending.back();
        pending.pop_back();
        const std::size_t i = t.index();
        for (int k = 0; k < 3; ++k) {
            const auto nb = mesh.otherTriangle(t, k);
            if (!nb) {
                continue;
            }
            const std::size_t j = nb->index();
            if (state[j] < 0) {
                // Side k is constrained exactly when it is an edge of P, the
                // only place the classification flips. The box boundary is
                // constrained too, but it has no neighbor to cross to.
                state[j] = mesh.isConstrained(t, k)
                               ? static_cast<signed char>(1 - state[i])
                               : state[i];
                pending.push_back(*nb);
            }
        }
    }
 
    std::vector<std::size_t> parent(n);
    for (std::size_t i = 0; i < n; ++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;
    };
 
    mesh.visitTriangles([&](TriId t) {
        const std::size_t i = t.index();
        if (state[i] != 1) {
            return;
        }
        const auto tri = mesh.getShape(t);
        for (int k = 0; k < 3; ++k) {
            const auto nb = mesh.otherTriangle(t, k);
            if (!nb) {
                continue;
            }
            const std::size_t j = nb->index();
            if (j < i || state[j] != 1) {
                continue;  // shared edge handled from the smaller index only
            }
            const CSeg e(tri.get(k), tri.get(k + 1));
            if (other.intersects(e) &&
                (other.interiorsIntersect(e) ||
                 (!other.contains(e.min()) && !other.contains(e.max())))) {
                parent[findRoot(i)] = findRoot(j);
            }
        }
    });
 
    // D \ P has no isolated points and no pinches, so every component holds a
    // triangle whose open interior meets the open disk; count those classes.
    std::size_t classes = 0;
    std::vector<char> counted(n, 0);
    return mesh.visitTriangles([&](TriId t) {
        const std::size_t i = t.index();
        if (state[i] != 1 || !other.interiorsIntersect(mesh.getShape(t))) {
            return false;
        }
        const std::size_t root = findRoot(i);
        if (counted[root]) {
            return false;
        }
        counted[root] = 1;
        return ++classes >= 2;
    });
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherConvex& other) const {
    return separates(other.asPolygon());
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherPolygon& other) const {
    // Let A = *this and B = other. A.separates(B) asks whether B \ A is
    // disconnected. For a (topological) disk B this happens exactly when some
    // connected component of A ∩ B meets ∂B in two or more pieces -- a
    // "crosscut". A single bite (one contact arc) or an interior island (a hole,
    // zero contact arcs) leaves B connected, and two disjoint bites belong to
    // different components, so neither cuts. This is the same criterion the
    // Convex overload uses, generalized to a non-convex A: A may dip into B
    // through several separate pockets, so the count is per A ∩ B component, not
    // a raw boundary-arc tally.
    if (other.isDegenerate()) {
        // A polygon that collapses to a segment is cut as that segment is; a
        // point is never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
    if (isDegenerate()) {
        // What this degenerate shape spans does the cutting: a segment answers
        // as one, and a point disconnects nothing.
        if (const auto spanned = getIfSegment()) {
            return spanned->separates(other);
        }
        return false;
    }
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
 
    // A convex remover admits the event walk of Convex::separates(Polygon),
    // which stays in the shapes' native arithmetic; only a genuinely reflex A
    // needs the exact intersection machinery below.
    if (isConvex()) {
        return Convex<PointType>(vertices()).separates(other);
    }
 
    // Exact rationals: the components of A ∩ B are polygons whose vertices are
    // the intersections of the two boundaries, which are rational in general.
    using Number = Rational<BigInt>;  // exact, overflow-free (ERational's backing)
    using RPoint = Point<Number, typename OtherPolygon::PointType::LabelType>;
    using RPolygon = Polygon<RPoint>;
    using RSegment = Segment<RPoint>;
 
    const auto toR = [](const auto& p) {
        return RPoint(Number(p.x()), Number(p.y()));
    };
 
    // The 2-D components of A ∩ B. Lower-dimensional pieces (shared boundary
    // segments or touch points) cannot disconnect the 2-D body of B, so only the
    // polygon pieces matter.
    std::vector<RPolygon> pieces;
    for (const auto& piece : other.template intersection<Number>(*this)) {
        if (std::holds_alternative<RPolygon>(piece)) {
            pieces.push_back(std::get<RPolygon>(piece));
        }
    }
    if (pieces.empty()) {
        return false;
    }
 
    // The intersection returns pinched faces (two faces meeting at a single
    // vertex) as separate polygons, but such faces form one connected component
    // of the closed set A ∩ B. Union those that share a vertex so a bowtie
    // crosscut is counted as one component touching ∂B twice.
    const int k = static_cast<int>(pieces.size());
    std::vector<int> parent(k);
    for (int i = 0; i < k; ++i) {
        parent[i] = i;
    }
    const auto find = [&parent](int x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    };
    std::vector<std::pair<RPoint, int>> vertexOwner;
    for (int i = 0; i < k; ++i) {
        for (const auto& v : pieces[i].vertices()) {
            vertexOwner.emplace_back(v, i);
        }
    }
    std::sort(vertexOwner.begin(), vertexOwner.end(),
              [](const auto& a, const auto& b) { return a.first < b.first; });
    for (std::size_t i = 1; i < vertexOwner.size(); ++i) {
        if (vertexOwner[i].first == vertexOwner[i - 1].first) {
            parent[find(vertexOwner[i].second)] = find(vertexOwner[i - 1].second);
        }
    }
 
    // Component roots present at a point q of ∂B: every A ∩ B component whose
    // (closed) region contains q.
    const auto rootsAt = [&](const RPoint& q) {
        std::vector<int> roots;
        for (int i = 0; i < k; ++i) {
            if (pieces[i].contains(q)) {
                const int r = find(i);
                if (std::find(roots.begin(), roots.end(), r) == roots.end()) {
                    roots.push_back(r);
                }
            }
        }
        return roots;
    };
 
    // All piece vertices: these are exactly the points on ∂B where component
    // membership can change (crossings of ∂A with ∂B and shared-boundary ends),
    // so subdividing each edge of B at them makes every sub-edge uniformly
    // inside or outside a given component.
    std::vector<RPoint> breakpoints;
    for (const auto& piece : pieces) {
        for (const auto& v : piece.vertices()) {
            breakpoints.push_back(v);
        }
    }
 
    // Walk ∂B once, cyclically, sampling: each vertex, each interior crossing,
    // and the midpoint of each resulting sub-edge. Membership is constant between
    // consecutive breakpoints, so these samples label every arc (including
    // isolated point contacts, caught at the vertex/crossing samples).
    std::vector<std::vector<int>> elements;
    const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(other.size());
    for (std::ptrdiff_t i = 0; i < n; ++i) {
        const RPoint a = toR(other.get(i));
        const RPoint b = toR(other.get(i + 1));
 
        std::vector<RPoint> cuts;
        const RSegment edge(a, b);
        for (const auto& bp : breakpoints) {
            if (edge.interiorContains(bp)) {
                cuts.push_back(bp);
            }
        }
        std::sort(cuts.begin(), cuts.end(), [&](const RPoint& u, const RPoint& v) {
            return dotSign(u, v, a, b) > 0;
        });
        cuts.erase(std::unique(cuts.begin(), cuts.end()), cuts.end());
 
        std::vector<RPoint> chain;
        chain.push_back(a);
        chain.insert(chain.end(), cuts.begin(), cuts.end());
        chain.push_back(b);
 
        elements.push_back(rootsAt(a));  // the vertex of B
        for (std::size_t j = 0; j + 1 < chain.size(); ++j) {
            const RPoint mid((chain[j].x() + chain[j + 1].x()) / Number(2),
                             (chain[j].y() + chain[j + 1].y()) / Number(2));
            elements.push_back(rootsAt(mid));                  // the sub-edge
            if (j + 1 + 1 < chain.size()) {
                elements.push_back(rootsAt(chain[j + 1]));     // an interior crossing
            }
        }
    }
 
    // A component cuts B iff its contacts with ∂B fall into >= 2 maximal cyclic
    // runs. Count, per root, the runs as rising edges around the cycle.
    const std::size_t total = elements.size();
    const auto has = [](const std::vector<int>& s, int r) {
        return std::find(s.begin(), s.end(), r) != s.end();
    };
    std::vector<int> roots;
    for (const auto& s : elements) {
        for (int r : s) {
            if (std::find(roots.begin(), roots.end(), r) == roots.end()) {
                roots.push_back(r);
            }
        }
    }
    for (int r : roots) {
        int runs = 0;
        for (std::size_t i = 0; i < total; ++i) {
            const std::size_t prev = (i + total - 1) % total;
            if (has(elements[i], r) && !has(elements[prev], r)) {
                if (++runs >= 2) {
                    return true;
                }
            }
        }
    }
    return false;
}

 
template <class PointType, class LabelType, class Storage>
template <class OtherShape, class TouchesBoundary>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separatesOneDimensional(
    const OtherShape& other, TouchesBoundary touchesBoundary) const {
    if (points_.empty()) {
        return false;
    }
    if (points_.size() == 1) {
        // A single vertex disconnects `other` iff it is an interior point of it.
        return other.interiorContains((*this)[0]);
    }
    // Walk the connected components of (this ∩ other): consecutive edge pieces
    // belong to one component iff the shared chain vertex lies on `other`
    // (`other` is convex, so each edge meets it in one connected piece). A
    // component that avoids every boundary point of `other` disconnects it.
    bool active = false;
    bool touched = false;
    for (std::size_t i = 0; i + 1 < points_.size(); ++i) {
        const auto edge = this->template boundaryAt<false>(i);
        if (edge.intersects(other)) {
            const bool connected = active && other.contains((*this)[i]);
            if (!connected) {
                if (active && !touched) {
                    return true;
                }
                touched = false;
            }
            active = true;
            if (touchesBoundary(edge)) {
                touched = true;
            }
        } else {
            if (active && !touched) {
                return true;
            }
            active = false;
            touched = false;
        }
    }
    return active && !touched;
}
 
template <class PointType, class LabelType, class Storage>
template <bool OtherIsConvex, class OtherShape>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separatesTwoDimensional(const OtherShape& other) const {
    if (points_.size() <= 1) {
        return false;
    }
 
    // Removing the chain disconnects the region iff it removes a crosscut:
    // either a single edge already disconnects the region, or the chain has
    // points p < x < q in chain order with x in the region's interior and
    // p, q both non-interior (boundary or beyond) -- the piece of the chain
    // around x then runs boundary-to-boundary through the interior.
    //
    // The scan looks for the p/x/q pattern at the vertices. That is complete
    // for a convex region: an edge between two interior vertices stays
    // interior, so a mid-edge p, x, or q always has a vertex witness on the
    // correct side (an edge whose non-interior endpoints sandwich an interior
    // point is caught as a separating edge). A non-convex region breaks the
    // first rule -- an edge can leave the interior between two interior
    // vertices -- so OtherIsConvex == false additionally tests such edges,
    // each excursion acting as one more non-interior stop between its
    // endpoints. The last vertex pre-seeds the final non-interior stop, so
    // the scan can stop as soon as the interior stop is confirmed.
    bool a_check = false;  // saw a non-interior stop ...
    bool b_check = false;  // ... then an interior stop ...
    bool c_check = !other.interiorContains((*this)[size() - 1]);  // ... then a non-interior stop
    const auto pattern_found = [&](bool stop_is_interior) {
        if (!a_check) {
            a_check = !stop_is_interior;
        }
        else if (!b_check) {
            b_check = stop_is_interior;
        }
        else if (!c_check) {
            c_check = !stop_is_interior;
        }
        return a_check && b_check && c_check;
    };
 
    bool prev_interior = other.interiorContains((*this)[0]);
    pattern_found(prev_interior);
    for (std::size_t i = 1; i < size(); ++i) {
        const Segment<PointType> edge((*this)[i - 1], (*this)[i]);
        if (edge.separates(other)) {
            return true;
        }
        const bool cur_interior = other.interiorContains((*this)[i]);
        if constexpr (!OtherIsConvex) {
            if (prev_interior && cur_interior && !other.interiorContains(edge) &&
                pattern_found(false)) {
                return true;
            }
        }
        if (pattern_found(cur_interior)) {
            return true;
        }
        prev_interior = cur_interior;
    }
    return false;
}
 
template <class PointType, class LabelType, class Storage>
template<SegmentConcept OtherSegment>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherSegment& other) const {
    return separatesOneDimensional(other, [&other](const auto& edge) {
        return edge.contains(other.min()) || edge.contains(other.max());
    });
}
 
template <class PointType, class LabelType, class Storage>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherOrientedSegment& other) const {
    return separates(other.asSegment());
}
 
template <class PointType, class LabelType, class Storage>
template<LineConcept OtherLine>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherLine& other) const {
    // A line has no boundary, so any nonempty intersection with the bounded
    // chain disconnects it.
    return intersects(other);
}
 
template <class PointType, class LabelType, class Storage>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherOrientedLine& other) const {
    return intersects(other);
}
 
template <class PointType, class LabelType, class Storage>
template<RayConcept OtherRay>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherRay& other) const {
    return separatesOneDimensional(other, [&other](const auto& edge) {
        return edge.contains(other.source());
    });
}
 
template <class PointType, class LabelType, class Storage>
template<HalfplaneConcept OtherHalfplane>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherHalfplane& other) const {
    // No straight edge can disconnect a halfplane (Segment::separates is
    // always false there), but the chain can: bending through the interior
    // between two non-interior stops seals off a pocket against the boundary
    // line, which is exactly the crosscut pattern of the vertex scan.
    return separatesTwoDimensional(other);
}
 
template <class PointType, class LabelType, class Storage>
template<RectangleConcept OtherRectangle>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    if (!bbox().intersects(other)) {
        return false;
    }
    if (bbox().separates(other)) {
        return true;
    }
    return separatesTwoDimensional(other);
}
 
template <class PointType, class LabelType, class Storage>
template<TriangleConcept OtherTriangle>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherTriangle& other) const {
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    if (bbox().separates(other.bbox())) {
        return true;
    }
    return separatesTwoDimensional(other);
}
 
template <class PointType, class LabelType, class Storage>
template<DiskConcept OtherDisk>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherDisk& other) const {
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    return separatesTwoDimensional(other);
}
 
template <class PointType, class LabelType, class Storage>
template<ConvexConcept OtherConvex>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherConvex& other) const {
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    if (bbox().separates(other.bbox())) {
        return true;
    }
    return separatesTwoDimensional(other);
}
 
template <class PointType, class LabelType, class Storage>
template<PolygonConcept OtherPolygon>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherPolygon& other) const {
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    if (bbox().separates(other.bbox())) {
        return true;
    }
    // The polygon may be non-convex, so the scan must also spot excursion
    // edges (see separatesTwoDimensional).
    return separatesTwoDimensional<false>(other);
}
 
template <class PointType, class LabelType, class Storage>
template<MonotoneChainConcept OtherChain>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherChain& other) const {
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    if (bbox().separates(other.bbox())) {
        return true;
    }
 
    // Polyline this separates polyline other if there exist
    // three ordered points a,b,c in other, such that
    // !this->contains(a), this->contains(b), !this->contains(c)
    // These points a,b,c may not be vertices, but if they are on
    // the same edge, than the edge is separated. If they are not
    // on the same edge, then there are two edges that are not contained
    // in this with an intermediate point b that is contained
 
    bool a_check = false;
    bool b_check = false;
    bool c_check = !contains(other.get(-1));
 
    for (std::size_t i = 1; i < other.size(); ++i) {
        Segment<typename OtherChain::PointType> edge(other[i-1], other[i]);
        if (separates(edge)) {
            return true; // If edge is separated, then a,b,c are in edge
        }
        if (!a_check) { // We did not find point a st !this->contains(a) yet
            if (!contains(edge)) {
                a_check = true; // point of edge not contained in this is a
                b_check = contains(other[i]); // maybe b st this->contains(b) is in edge
            }
        }
        else if (!b_check) { // We did not find point b > a st this->contains(b) yet
            b_check = intersects(edge); // If the edge intersects this, we found b
        }
        else if (!c_check) { // We did not find point c > b st !this->contains(c) yet
            c_check = !contains(edge); // If edge is not contained, we found c
        }
 
        if (a_check && b_check && c_check) { // All points found
            return true;
        }
    }
 
    // There are no such a,b,c
    return false;
}
 
template <class PointType, class LabelType, class Storage>
template<PointConcept OtherPoint>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const Shape<OtherPoint>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}
 
template <class Number, class Label>
template<MonotoneChainConcept OtherChain>
constexpr bool Point<Number, Label>::separates(const OtherChain& other) const {
    // Removing a point disconnects the arc iff it is a relative-interior point.
    return other.interiorContains(*this);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool Segment<PointType, LabelType>::separates(const OtherChain& other) const {
    return detail::separatesChain(*this, other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherChain& other) const {
    return detail::separatesChain(asSegment(), other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool Line<PointType, LabelType>::separates(const OtherChain& other) const {
    return detail::separatesChain(*this, other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherChain& other) const {
    return detail::separatesChain(*this, other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool Ray<PointType, LabelType>::separates(const OtherChain& other) const {
    return detail::separatesChain(*this, other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherChain& other) const {
    return detail::separatesChain(*this, other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherChain& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return detail::separatesChain(*this, other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherChain& other) const {
    return detail::separatesChain(*this, other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool Disk<PointType, LabelType>::separates(const OtherChain& other) const {
    return detail::separatesChain(*this, other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool Convex<PointType, LabelType>::separates(const OtherChain& other) const {
    return detail::separatesChain(*this, other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherChain& other) const {
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    if (bbox().separates(other.bbox())) {
        return true;
    }
 
    // Polygon this separates polyline other if there exist
    // three ordered points a,b,c in other, such that
    // !this->contains(a), this->contains(b), !this->contains(c)
    // These points a,b,c may not be vertices, but if they are on
    // the same edge, than the edge is separated. If they are not
    // on the same edge, then there are two edges that are not contained
    // in this with an intermediate point b that is contained
 
    bool a_check = false;
    bool b_check = false;
    bool c_check = !contains(other.get(-1));
 
    for (std::size_t i = 1; i < other.size(); ++i) {
        Segment<typename OtherChain::PointType> edge(other[i-1], other[i]);
        if (separates(edge)) {
            return true; // If edge is separated, then a,b,c are in edge
        }
        if (!a_check) { // We did not find point a st !this->contains(a) yet
            if (!contains(edge)) {
                a_check = true; // point of edge not contained in this is a
                b_check = contains(other[i]); // maybe b st this->contains(b) is in edge
            }
        }
        else if (!b_check) { // We did not find point b > a st this->contains(b) yet
            b_check = intersects(edge); // If the edge intersects this, we found b
        }
        else if (!c_check) { // We did not find point c > b st !this->contains(c) yet
            c_check = !contains(edge); // If edge is not contained, we found c
        }
 
        if (a_check && b_check && c_check) { // All points found
            return true;
        }
    }
 
    // There are no such a,b,c
    return false;
}

 
namespace detail {

template <class Remover, class TargetEdgeRange>
bool separates1DSet(const Remover& remover, const TargetEdgeRange& targetEdges) {
    using TargetEdge = std::remove_cvref_t<std::ranges::range_value_t<TargetEdgeRange>>;
    using TargetNumber = typename TargetEdge::PointType::NumberType;
    using RemoverNumber = typename Remover::NumberType;
    constexpr bool approximate = std::is_floating_point_v<TargetNumber> ||
                                 std::is_floating_point_v<RemoverNumber>;
    using ExactNumber = std::conditional_t<approximate, double, Rational<BigInt>>;
    using ExactPoint = Point<ExactNumber>;
    using ExactSegment = Segment<ExactPoint>;
 
    // A free piece: a maximal subsegment of one target edge that survives the
    // removal. An endpoint marked open lies on the remover and is excluded
    // from the piece.
    struct Piece {
        ExactPoint lo;
        ExactPoint hi;
        bool loOpen;
        bool hiOpen;
    };
    std::vector<Piece> pieces;
 
    for (const auto& edge : targetEdges) {
        if (edge.min() == edge.max()) {
            continue;  // zero-length edge: adds no connectivity of its own
        }
        // Removed intervals along the edge, as closed intervals in the
        // lexicographic point order (the linear order along the edge).
        std::vector<std::pair<ExactPoint, ExactPoint>> removed;
        const auto addRemoved = [&removed](const auto& optionalPiece) {
            if (!optionalPiece) {
                return;
            }
            if (const auto* point = std::get_if<0>(&*optionalPiece)) {
                removed.emplace_back(ExactPoint(*point), ExactPoint(*point));
            } else {
                const auto& overlap = std::get<1>(*optionalPiece);
                removed.emplace_back(ExactPoint(overlap.min()), ExactPoint(overlap.max()));
            }
        };
        if constexpr (is_point_v<Remover>) {
            if (edge.contains(remover)) {
                removed.emplace_back(ExactPoint(remover), ExactPoint(remover));
            }
        } else if constexpr (is_segment_v<Remover>) {
            addRemoved(edge.template intersection<ExactNumber>(remover));
        } else if constexpr (is_halfplane_intersection_v<Remover>) {
            // A convex region meets each edge in one connected piece; its
            // segment intersection reports the point or subsegment directly.
            addRemoved(remover.template intersection<ExactNumber>(edge));
        } else if constexpr (is_polygon_v<Remover>) {
            // A (possibly reflex) polygon can bite several intervals out of
            // one edge; its segment intersection lists them all as pieces.
            for (const auto& piece : remover.template intersection<ExactNumber>(edge)) {
                if (const auto* point = std::get_if<0>(&piece)) {
                    removed.emplace_back(ExactPoint(*point), ExactPoint(*point));
                } else {
                    const auto& overlap = std::get<1>(piece);
                    removed.emplace_back(ExactPoint(overlap.min()), ExactPoint(overlap.max()));
                }
            }
        } else {
            static_assert(is_polyline_v<Remover> || is_monotone_chain_v<Remover>,
                          "unsupported remover kind");
            if (remover.size() == 1) {
                if (edge.contains(remover[0])) {
                    removed.emplace_back(ExactPoint(remover[0]), ExactPoint(remover[0]));
                }
            } else {
                for (const auto& removerEdge : remover.edgesView()) {
                    addRemoved(edge.template intersection<ExactNumber>(removerEdge));
                }
            }
        }
        std::sort(removed.begin(), removed.end());
        // Sweep the merged removed blocks and emit the complementary free
        // pieces. A free endpoint abutting a removed block is open; the edge
        // extremes are closed unless removed.
        ExactPoint cursor(edge.min());
        bool cursorOpen = false;
        std::size_t i = 0;
        while (i < removed.size()) {
            const ExactPoint blockLo = removed[i].first;
            ExactPoint blockHi = removed[i].second;
            for (++i; i < removed.size() && !(blockHi < removed[i].first); ++i) {
                if (blockHi < removed[i].second) {
                    blockHi = removed[i].second;
                }
            }
            if (cursor < blockLo) {
                pieces.push_back({cursor, blockLo, cursorOpen, true});
            }
            cursor = blockHi;
            cursorOpen = true;
        }
        const ExactPoint edgeHi(edge.max());
        if (cursor < edgeHi) {
            pieces.push_back({cursor, edgeHi, cursorOpen, false});
        }
    }
 
    if (pieces.size() < 2) {
        return false;  // zero or one piece is never disconnected
    }
 
    // Union-find over the pieces; two pieces join iff they share a point that
    // survives the removal.
    std::vector<std::size_t> parent(pieces.size());
    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;
    };
    const auto connected = [](const Piece& a, const Piece& b) {
        const ExactSegment sa(a.lo, a.hi);
        const ExactSegment sb(b.lo, b.hi);
        if (!sa.intersects(sb)) {
            return false;
        }
        if (sa.collinear(sb) && sa.min() < sb.max() && sb.min() < sa.max()) {
            // A positive-length overlap survives losing at most four
            // excluded endpoints.
            return true;
        }
        // Single common point: it joins the pieces unless it coincides with
        // an excluded (open) endpoint of either piece (the unique common
        // point equals x iff x lies on the other piece's segment).
        if (a.loOpen && sb.contains(a.lo)) {
            return false;
        }
        if (a.hiOpen && sb.contains(a.hi)) {
            return false;
        }
        if (b.loOpen && sa.contains(b.lo)) {
            return false;
        }
        if (b.hiOpen && sa.contains(b.hi)) {
            return false;
        }
        return true;
    };
    for (std::size_t i = 0; i < pieces.size(); ++i) {
        for (std::size_t j = i + 1; j < pieces.size(); ++j) {
            const std::size_t ri = findRoot(i);
            const std::size_t rj = findRoot(j);
            if (ri != rj && connected(pieces[i], pieces[j])) {
                parent[ri] = rj;
            }
        }
    }
    std::size_t components = 0;
    for (std::size_t i = 0; i < pieces.size(); ++i) {
        if (findRoot(i) == i) {
            ++components;
        }
    }
    return components >= 2;
}
 
}  // namespace detail
 
template <class PointType, class LabelType>
template<SegmentConcept OtherSegment>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherSegment& other) const {
    if (empty()) {
        return false;
    }
    if (size() == 1) {
        // A single vertex disconnects the segment iff it is an interior point of it.
        return other.interiorContains((*this)[0]);
    }
    return detail::separates1DSet(*this, std::array<OtherSegment, 1>{other});
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherPolyline& other) const {
    if (empty() || other.size() < 2) {
        return false;
    }
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    if (size() == 1) {
        return (*this)[0].separates(other);
    }
    return detail::separates1DSet(*this, other.edgesView());
}
 
template <class Number, class Label>
template<PolylineConcept OtherPolyline>
constexpr bool Point<Number, Label>::separates(const OtherPolyline& other) const {
    // Unlike a monotone chain, the polyline may reconnect around the removed
    // point through a self-crossing or a revisited vertex, so an interior
    // point is not necessarily a cut point.
    if (other.size() < 2) {
        return false;
    }
    return detail::separates1DSet(*this, other.edgesView());
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool Segment<PointType, LabelType>::separates(const OtherPolyline& other) const {
    if (other.size() < 2) {
        return false;
    }
    return detail::separates1DSet(*this, other.edgesView());
}
 
namespace detail {

template <class ExactPoint, class PolylineType, class CutterRange>
std::vector<Segment<ExactPoint>> arrangedPolylineEdges(const PolylineType& polyline,
                                                       const CutterRange& cutters) {
    using ExactSegment = Segment<ExactPoint>;
    using ExactNumber = typename ExactPoint::NumberType;
 
    std::vector<ExactSegment> edges;
    for (const auto& edge : polyline.edgesView()) {
        if (edge.min() != edge.max()) {
            edges.emplace_back(ExactPoint(edge.min()), ExactPoint(edge.max()));
        }
    }
 
    std::vector<ExactSegment> result;
    std::vector<ExactPoint> cuts;
    for (std::size_t i = 0; i < edges.size(); ++i) {
        cuts.clear();
        cuts.push_back(edges[i].min());
        cuts.push_back(edges[i].max());
        const auto addCuts = [&](const auto& s) {
            const auto piece = edges[i].template intersection<ExactNumber>(s);
            if (!piece) {
                return;
            }
            if (const auto* point = std::get_if<0>(&*piece)) {
                cuts.push_back(ExactPoint(*point));
            } else {
                const auto& overlap = std::get<1>(*piece);
                cuts.push_back(ExactPoint(overlap.min()));
                cuts.push_back(ExactPoint(overlap.max()));
            }
        };
        for (std::size_t j = 0; j < edges.size(); ++j) {
            if (j != i) {
                addCuts(edges[j]);
            }
        }
        for (const auto& cutter : cutters) {
            addCuts(cutter);
        }
        // All cut points lie on edge i, so the lexicographic point order is
        // the linear order along the edge.
        std::sort(cuts.begin(), cuts.end());
        cuts.erase(std::unique(cuts.begin(), cuts.end()), cuts.end());
        for (std::size_t k = 0; k + 1 < cuts.size(); ++k) {
            result.emplace_back(cuts[k], cuts[k + 1]);
        }
    }
    std::sort(result.begin(), result.end());
    result.erase(std::unique(result.begin(), result.end()), result.end());
    return result;
}
 
// Exact1DNumber, the exact coordinate type for a mixed pair, lives in
// predicates_helpers.hpp: the containment predicates need it too, and that
// header is parsed well before this one.

template <class Remover, class PolylineType>
bool separatesPolylineSet(const Remover& remover, const PolylineType& target) {
    if (target.size() < 2) {
        // Removing anything from at most one point cannot disconnect it.
        return false;
    }
    // A polyline is connected, so a remover that misses it removes nothing that
    // matters — the same shortcut @ref cellSeparates takes, and worth more here,
    // because the arrangement below is of the *target* alone and so is rebuilt
    // in full for every remover it is asked about. Only bounded removers can be
    // ruled out this way: a line, a ray and a half-plane have no bounding box,
    // which is exactly why they have no `bbox()` to call.
    if constexpr (requires { remover.bbox(); }) {
        if (!remover.bbox().intersects(target.bbox())) {
            return false;
        }
    }
    using ExactNumber = Exact1DNumber<typename PolylineType::NumberType,
                                      typename Remover::NumberType>;
    using ExactPoint = Point<ExactNumber>;
    const auto subEdges = arrangedPolylineEdges<ExactPoint>(
        target, std::array<Segment<ExactPoint>, 0>{});
 
    std::map<ExactPoint, std::size_t> nodes;
    std::vector<std::size_t> parent;
    const auto findRoot = [&parent](std::size_t x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    };
    const auto nodeOf = [&](const ExactPoint& p) {
        const auto [it, inserted] = nodes.try_emplace(p, parent.size());
        if (inserted) {
            parent.push_back(parent.size());
        }
        return it->second;
    };
 
    for (const auto& s : subEdges) {
        // contains(endpoint) implies intersects(s), so a sub-edge missed by
        // the remover always has both endpoints surviving.
        const bool minSurvives = !remover.contains(s.min());
        const bool maxSurvives = !remover.contains(s.max());
        std::optional<std::size_t> a;
        std::optional<std::size_t> b;
        if (minSurvives) {
            a = nodeOf(s.min());
        }
        if (maxSurvives) {
            b = nodeOf(s.max());
        }
        if (a && b && !remover.intersects(s)) {
            parent[findRoot(*a)] = findRoot(*b);
        }
    }
 
    std::size_t components = 0;
    for (std::size_t i = 0; i < parent.size(); ++i) {
        if (findRoot(i) == i && ++components >= 2) {
            return true;
        }
    }
    return false;
}

template <class PolylineType, class Region>
bool polylineSeparatesConvexRegion(const PolylineType& polyline, const Region& other) {
    if (polyline.size() < 2) {
        return false;
    }
    // Both operands are connected, so shapes that miss each other neither cut
    // nor are cut — as in @ref cellSeparates. A half-plane operand has no
    // bounding box to test, and no `bbox()` to call.
    if constexpr (requires { other.bbox(); }) {
        if (!polyline.bbox().intersects(other.bbox())) {
            return false;
        }
    }
    using ExactNumber = Exact1DNumber<typename PolylineType::NumberType,
                                      typename Region::NumberType>;
    using ExactPoint = Point<ExactNumber>;
    const auto subEdges = arrangedPolylineEdges<ExactPoint>(
        polyline, std::array<Segment<ExactPoint>, 0>{});
 
    std::map<ExactPoint, std::size_t> nodes;
    std::vector<std::size_t> parent(1, 0);  // slot 0: the contracted boundary
    const auto findRoot = [&parent](std::size_t x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    };
    const auto nodeOf = [&](const ExactPoint& p) {
        const auto [it, inserted] = nodes.try_emplace(p, parent.size());
        if (inserted) {
            parent.push_back(parent.size());
        }
        return it->second;
    };
 
    for (const auto& s : subEdges) {
        if (!other.interiorsIntersect(s)) {
            continue;  // no piece through the interior (skips boundary runs)
        }
        const std::size_t a = other.interiorContains(s.min()) ? nodeOf(s.min()) : 0;
        const std::size_t b = other.interiorContains(s.max()) ? nodeOf(s.max()) : 0;
        const std::size_t ra = findRoot(a);
        const std::size_t rb = findRoot(b);
        if (ra == rb) {
            return true;  // the sub-edge closes a cycle through the interior
        }
        parent[ra] = rb;
    }
    return false;
}

template <class PolylineType, class PolygonType>
bool polylineSeparatesPolygon(const PolylineType& polyline, const PolygonType& other) {
    if (polyline.size() < 2) {
        return false;
    }
    using ExactNumber = Exact1DNumber<typename PolylineType::NumberType,
                                      typename PolygonType::NumberType>;
    using ExactPoint = Point<ExactNumber>;
    using ExactSegment = Segment<ExactPoint>;
 
    std::vector<ExactSegment> cutters;
    for (const auto& edge : other.edgesView()) {
        cutters.emplace_back(ExactPoint(edge.min()), ExactPoint(edge.max()));
    }
    const auto subEdges = arrangedPolylineEdges<ExactPoint>(polyline, cutters);
 
    std::map<ExactPoint, std::size_t> nodes;
    std::vector<std::size_t> parent(1, 0);  // slot 0: the contracted boundary
    const auto findRoot = [&parent](std::size_t x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    };
    const auto nodeOf = [&](const ExactPoint& p) {
        const auto [it, inserted] = nodes.try_emplace(p, parent.size());
        if (inserted) {
            parent.push_back(parent.size());
        }
        return it->second;
    };
 
    for (const auto& s : subEdges) {
        const ExactPoint mid((s.min().x() + s.max().x()) / ExactNumber(2),
                             (s.min().y() + s.max().y()) / ExactNumber(2));
        if (!other.interiorContains(mid)) {
            continue;  // outside or along the boundary
        }
        // The sub-edge never crosses the boundary, so its endpoints are in the
        // closed polygon: on the boundary (contracted) or strictly inside.
        const std::size_t a = other.boundaryContains(s.min()) ? 0 : nodeOf(s.min());
        const std::size_t b = other.boundaryContains(s.max()) ? 0 : nodeOf(s.max());
        const std::size_t ra = findRoot(a);
        const std::size_t rb = findRoot(b);
        if (ra == rb) {
            return true;  // the sub-edge closes a cycle through the interior
        }
        parent[ra] = rb;
    }
    return false;
}
 
}  // namespace detail
 
template <class PointType, class LabelType>
template<OrientedSegmentConcept OtherOrientedSegment>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return separates(other.asSegment());
}
 
template <class PointType, class LabelType>
template<LineConcept OtherLine>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherLine& other) const {
    // A line has no boundary, so any nonempty intersection with the bounded
    // polyline disconnects it.
    return intersects(other);
}
 
template <class PointType, class LabelType>
template<OrientedLineConcept OtherOrientedLine>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return intersects(other);
}
 
template <class PointType, class LabelType>
template<RayConcept OtherRay>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherRay& other) const {
    if (empty()) {
        return false;
    }
    if (size() == 1) {
        // A single vertex disconnects the ray iff it is an interior point of it.
        return other.interiorContains((*this)[0]);
    }
    using ExactNumber = detail::Exact1DNumber<NumberType, typename OtherRay::PointType::NumberType>;
    using ExactPoint = Point<ExactNumber>;
    // Clip the ray to a segment leaving the polyline's bounding box: every
    // removed point lies inside the box, so the clipped segment has the same
    // free pieces as the ray, its far end standing in for the infinite tail.
    const ExactPoint source(other.source());
    const ExactPoint target(other.target());
    const ExactNumber dx = target.x() - source.x();
    const ExactNumber dy = target.y() - source.y();
    if (dx == ExactNumber{} && dy == ExactNumber{}) {
        return false;  // degenerate ray
    }
    const auto& box = bbox();
    ExactNumber scale(1);
    const auto extend = [&scale](const ExactNumber& s, const ExactNumber& d,
                                 const ExactNumber& lo, const ExactNumber& hi) {
        if (d > ExactNumber{}) {
            const ExactNumber needed = (hi - s) / d + ExactNumber(1);
            if (scale < needed) {
                scale = needed;
            }
        } else if (d < ExactNumber{}) {
            const ExactNumber needed = (s - lo) / -d + ExactNumber(1);
            if (scale < needed) {
                scale = needed;
            }
        }
    };
    extend(source.x(), dx, ExactNumber(box.min().x()), ExactNumber(box.max().x()));
    extend(source.y(), dy, ExactNumber(box.min().y()), ExactNumber(box.max().y()));
    const ExactPoint beyond(source.x() + scale * dx, source.y() + scale * dy);
    return detail::separates1DSet(
        *this, std::array<Segment<ExactPoint>, 1>{Segment<ExactPoint>(source, beyond)});
}
 
template <class PointType, class LabelType>
template<HalfplaneConcept OtherHalfplane>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherHalfplane& other) const {
    return detail::polylineSeparatesConvexRegion(*this, other);
}
 
template <class PointType, class LabelType>
template<RectangleConcept OtherRectangle>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    if (size() < 2) {
        return false;
    }
    if (!bbox().intersects(other)) {
        return false;
    }
    return detail::polylineSeparatesConvexRegion(*this, other);
}
 
template <class PointType, class LabelType>
template<TriangleConcept OtherTriangle>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherTriangle& other) const {
    if (size() < 2) {
        return false;
    }
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    return detail::polylineSeparatesConvexRegion(*this, other);
}
 
template <class PointType, class LabelType>
template<DiskConcept OtherDisk>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherDisk& other) const {
    if (size() < 2) {
        return false;
    }
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    return detail::polylineSeparatesConvexRegion(*this, other);
}
 
template <class PointType, class LabelType>
template<ConvexConcept OtherConvex>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherConvex& other) const {
    if (size() < 2) {
        return false;
    }
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    return detail::polylineSeparatesConvexRegion(*this, other);
}
 
template <class PointType, class LabelType>
template<PolygonConcept OtherPolygon>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherPolygon& other) const {
    if (size() < 2) {
        return false;
    }
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    return detail::polylineSeparatesPolygon(*this, other);
}
 
template <class PointType, class LabelType>
template<MonotoneChainConcept OtherChain>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherChain& other) const {
    if (empty() || other.size() < 2) {
        return false;
    }
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    if (size() == 1) {
        return (*this)[0].separates(other);
    }
    return detail::separates1DSet(*this, other.edgesView());
}
 
template <class PointType, class LabelType>
template<PointConcept OtherPoint>
constexpr bool Polyline<PointType, LabelType>::separates(const Shape<OtherPoint>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherPolyline& other) const {
    return asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool Line<PointType, LabelType>::separates(const OtherPolyline& other) const {
    return detail::separatesPolylineSet(*this, other);
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherPolyline& other) const {
    return asLine().separates(other);
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool Ray<PointType, LabelType>::separates(const OtherPolyline& other) const {
    return detail::separatesPolylineSet(*this, other);
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherPolyline& other) const {
    return detail::separatesPolylineSet(*this, other);
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherPolyline& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return detail::separatesPolylineSet(*this, other);
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherPolyline& other) const {
    return detail::separatesPolylineSet(*this, other);
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool Disk<PointType, LabelType>::separates(const OtherPolyline& other) const {
    return detail::separatesPolylineSet(*this, other);
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool Convex<PointType, LabelType>::separates(const OtherPolyline& other) const {
    return detail::separatesPolylineSet(*this, other);
}
 
template <class PointType, class LabelType, class Storage>
template<PolylineConcept OtherPolyline>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherPolyline& other) const {
    if (empty() || other.size() < 2) {
        return false;
    }
    if (size() == 1) {
        return (*this)[0].separates(other);
    }
    // The chain is neither convex nor a single interval on the polyline's
    // edges, so its removed pieces are computed exactly.
    return detail::separates1DSet(*this, other.edgesView());
}
 
template <class PointType, class LabelType>
template<PolylineConcept OtherPolyline>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherPolyline& other) const {
    if (other.size() < 2) {
        return false;
    }
    if (!bbox().intersects(other.bbox())) {
        return false;
    }
    // The polygon may be reflex, so it can bite several intervals out of one
    // polyline edge; the free pieces are joined geometrically.
    return detail::separates1DSet(*this, other.edgesView());
}
 
 
// ---------------------------------------------------------------------------
// HalfplaneIntersection
 
template <class PointType, class LabelType>
template <PointConcept OtherPoint>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherPoint&) const {
    // Removing anything from a single point leaves the point or nothing,
    // never two components.
    return false;
}
 
template <class PointType, class LabelType>
template <SegmentConcept OtherSegment>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherSegment& other) const {
    // The region meets the segment in one closed sub-interval; the remainder
    // is disconnected exactly when that interval is nonempty and touches
    // neither endpoint.
    return intersects(other) && !contains(other[0]) && !contains(other[1]);
}
 
template <class PointType, class LabelType>
template <OrientedSegmentConcept OtherOrientedSegment>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return separates(Segment<typename OtherOrientedSegment::PointType>(other[0], other[1]));
}
 
template <class PointType, class LabelType>
template <LineConcept OtherLine>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherLine& other) const {
    // The region meets the line in one closed interval; the remainder is
    // disconnected exactly when the interval is nonempty and bounded on both
    // sides (an unbounded interval swallows an end of the line, leaving a
    // single ray, and a full-line interval leaves nothing).
    if (empty() || halfplanes_.empty()) {
        return false;
    }
    const Halfplane<typename OtherLine::PointType> along(other[0], other[1]);
    const auto clip = clipLine(along);
    return !clip.empty && clip.entry >= 0 && clip.exit >= 0;
}
 
template <class PointType, class LabelType>
template <OrientedLineConcept OtherOrientedLine>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return separates(other.asLine());
}
 
template <class PointType, class LabelType>
template <RayConcept OtherRay>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherRay& other) const {
    // As for a line, but the interval additionally must not contain the
    // source (which would leave a single piece ahead of it), and only the
    // forward end has to be bounded.
    if (empty() || halfplanes_.empty() || contains(other.source())) {
        return false;
    }
    const Halfplane<typename OtherRay::PointType> along(other.source(), other.target());
    const auto clip = clipLine(along);
    if (clip.empty || clip.exit < 0) {
        return false;
    }
    // The interval must reach into the ray: the line may not leave the region
    // before the source.
    return constraintSide(static_cast<std::size_t>(clip.exit), other.source()) >= 0;
}
 
template <class PointType, class LabelType>
template <HalfplaneConcept OtherHalfplane>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherHalfplane& other) const {
    // Removing the region disconnects the half-plane exactly when the
    // region's closed contact set with the half-plane's extended boundary —
    // the boundary line plus the half-plane's arc of ideal directions — has at
    // least two connected components while the region also meets the
    // half-plane's interior. A bounded region only ever contributes the
    // boundary-line contact (one component), so only unbounded regions can
    // separate a half-plane. Complexity: O(n).
    if (empty() || halfplanes_.empty() || halfplanes_.size() == 1) {
        // The whole plane leaves nothing behind, and removing a single
        // half-plane leaves an open half-plane, which is connected.
        return false;
    }
    // The region must reach the half-plane's interior: with contact only on
    // the boundary line, the remainder stays connected through the interior.
    const SupStatus infimum = supStatus(other.opposite());
    if (infimum == SupStatus::below || infimum == SupStatus::on) {
        return false;
    }
    // Contact on the boundary line itself (one component when present,
    // possibly extended to the ideal endpoints of the line).
    const bool lineTouch = intersects(other.asLine());
    const auto reversed = other.opposite();
    int pieces = lineTouch ? 1 : 0;
    // Ideal contact: recession directions of the region that stay weakly
    // inside the half-plane. Each recession arc meets the half-plane's closed
    // ideal half-circle in at most one sub-arc.
    for (const auto& arc : recessionArcs()) {
        // Does the arc reach the line's ideal endpoints or the open ideal arc
        // strictly inside the half-plane?
        const bool atLeft = detail::arcContainsDirection(arc.first, arc.second, reversed);
        const bool atRight = detail::arcContainsDirection(arc.first, arc.second, other);
        const bool strictlyInside = detail::directionCross(other, arc.first) > 0 ||
                                    detail::directionCross(other, arc.second) > 0;
        if (!(atLeft || atRight || strictlyInside)) {
            continue;  // this recession arc never touches the half-plane's ideal boundary
        }
        if (lineTouch && (atLeft || atRight)) {
            continue;  // joined to the boundary-line contact at an ideal line endpoint
        }
        ++pieces;
    }
    return pieces >= 2;
}
 
// Forward area targets: removing the region from a bounded 2D shape B is
// decided on the part of the region near B (B is bounded), so clip the region
// to a box around B and defer to Convex::separates. A degenerate region
// reduces to its carrier (a line through B still cuts it).
 
template <class PointType, class LabelType>
template <RectangleConcept OtherRectangle>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    if (empty()) {
        return false;
    }
    using E = detail::region_exact_number_t<NumberType>;
    if (isDegenerate()) {
        return std::visit([&other](const auto& carrier) { return carrier.separates(other); },
                          detail::degenerateRegionCarrier(*this));
    }
    const auto clipped = detail::regionClippedToBox(*this, other.bbox());
    if (!clipped.isBounded()) {
        return false;
    }
    return clipped.template asConvex<E>().separates(other);
}
 
template <class PointType, class LabelType>
template <TriangleConcept OtherTriangle>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherTriangle& other) const {
    if (empty()) {
        return false;
    }
    using E = detail::region_exact_number_t<NumberType>;
    if (isDegenerate()) {
        return std::visit([&other](const auto& carrier) { return carrier.separates(other); },
                          detail::degenerateRegionCarrier(*this));
    }
    const auto clipped = detail::regionClippedToBox(*this, other.bbox());
    if (!clipped.isBounded()) {
        return false;
    }
    return clipped.template asConvex<E>().separates(other);
}
 
template <class PointType, class LabelType>
template <ConvexConcept OtherConvex>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherConvex& other) const {
    if (empty()) {
        return false;
    }
    if (other.isDegenerate()) {
        // A convex with fewer than three vertices is the segment or point it
        // spans: the segment overload answers for the former, and a point is
        // never disconnected.
        if (const auto spanned = other.getIfSegment()) {
            return separates(*spanned);
        }
        return false;
    }
    using E = detail::region_exact_number_t<NumberType>;
    if (isDegenerate()) {
        return std::visit([&other](const auto& carrier) { return carrier.separates(other); },
                          detail::degenerateRegionCarrier(*this));
    }
    const auto clipped = detail::regionClippedToBox(*this, other.bbox());
    if (!clipped.isBounded()) {
        return false;
    }
    return clipped.template asConvex<E>().separates(other);
}
 
template <class PointType, class LabelType>
template <DiskConcept OtherDisk>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherDisk& other) const {
    if (empty() || other.isDegenerate()) {
        return false;
    }
    using E = detail::region_exact_number_t<NumberType>;
    if (isDegenerate()) {
        return std::visit([&other](const auto& carrier) { return carrier.separates(other); },
                          detail::degenerateRegionCarrier(*this));
    }
    const auto clipped = detail::regionClippedToBox(*this, other.bbox());
    if (!clipped.isBounded()) {
        return false;
    }
    return clipped.template asConvex<E>().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonConcept OtherPolygon>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherPolygon& other) const {
    if (empty() || other.size() < 3) {
        return false;
    }
    using E = detail::region_exact_number_t<NumberType>;
    if (isDegenerate()) {
        return std::visit([&other](const auto& carrier) { return carrier.separates(other); },
                          detail::degenerateRegionCarrier(*this));
    }
    const auto clipped = detail::regionClippedToBox(*this, other.bbox());
    if (!clipped.isBounded()) {
        return false;
    }
    return clipped.template asConvex<E>().separates(other);
}
 
template <class PointType, class LabelType>
template <MonotoneChainConcept OtherChain>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherChain& other) const {
    // The chain is a one-dimensional arc; the convex region removes one
    // connected piece from each edge, so the arc-component test applies with
    // the region as the (convex) remover.
    if (empty()) {
        return false;
    }
    return detail::separatesChain(*this, other);
}
 
template <class PointType, class LabelType>
template <PolylineConcept OtherPolyline>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherPolyline& other) const {
    if (empty() || other.size() < 2) {
        return false;
    }
    return detail::separates1DSet(*this, other.edgesView());
}
 
 
// ---------------------------------------------------------------------------
// Reverse direction: removing a lower-ranked shape from the region.
 
namespace detail {

template <class Region, class Remover>
bool boundedRemoverSeparatesRegion(const Region& region, const Remover& remover) {
    using N = typename Region::NumberType;
    using E = region_exact_number_t<N>;
    using EPoint = Point<E, typename Region::PointType::LabelType>;
    const auto bounds = remover.bbox();
    const N margin(2);
    const E boxLoX(N(bounds.min().x()) - margin);
    const E boxLoY(N(bounds.min().y()) - margin);
    const E boxHiX(N(bounds.max().x()) + margin);
    const E boxHiY(N(bounds.max().y()) + margin);
    const auto clipped = regionClippedToBox(region, bounds);
    if (!clipped.isBounded()) {
        return false;  // the remover's neighborhood misses the region
    }
    const auto poly = clipped.template asConvex<E>();
    const std::size_t n = poly.size();
    if (n < 3) {
        return false;
    }
    // The boundary is walked as a cyclic sequence of pieces, each free
    // (surviving the removal) or removed. Box (ideal) edges are always free.
    std::vector<bool> free;
    for (std::size_t i = 0; i < n; ++i) {
        const EPoint u(poly[i]);
        const EPoint v(poly[(i + 1) % n]);
        const bool ideal =
            (u.x() == boxLoX && v.x() == boxLoX) || (u.x() == boxHiX && v.x() == boxHiX) ||
            (u.y() == boxLoY && v.y() == boxLoY) || (u.y() == boxHiY && v.y() == boxHiY);
        if (ideal) {
            free.push_back(true);
            continue;
        }
        const Segment<EPoint> edge(u, v);
        if constexpr (is_polygon_v<Remover>) {
            // A reflex remover may bite several blocks out of one edge; list
            // them by their parameter along u -> v and alternate free/removed.
            const EPoint dir = v - u;
            const E edgeLen = dir.x() * dir.x() + dir.y() * dir.y();
            std::vector<std::pair<E, E>> blocks;
            for (const auto& piece : remover.template intersection<E>(edge)) {
                if (const auto* overlap = std::get_if<1>(&piece)) {
                    const EPoint lo(overlap->min());
                    const EPoint hi(overlap->max());
                    E tLo = (lo - u).x() * dir.x() + (lo - u).y() * dir.y();
                    E tHi = (hi - u).x() * dir.x() + (hi - u).y() * dir.y();
                    if (tHi < tLo) {
                        std::swap(tLo, tHi);
                    }
                    if (tLo < tHi) {
                        blocks.emplace_back(tLo, tHi);
                    }
                }
            }
            std::sort(blocks.begin(), blocks.end());
            E cursor{};
            bool wroteAny = false;
            for (const auto& block : blocks) {
                if (cursor < block.first) {
                    free.push_back(true);
                }
                free.push_back(false);
                cursor = block.second;
                wroteAny = true;
            }
            if (!wroteAny || cursor < edgeLen) {
                free.push_back(true);
            }
        } else {
            // A convex remover meets the edge in one connected block.
            const bool uIn = remover.contains(u);
            const bool vIn = remover.contains(v);
            if (uIn && vIn) {
                free.push_back(false);
            } else if (uIn) {
                free.push_back(false);
                free.push_back(true);
            } else if (vIn) {
                free.push_back(true);
                free.push_back(false);
            } else if (remover.interiorsIntersect(edge)) {
                free.push_back(true);
                free.push_back(false);
                free.push_back(true);
            } else {
                free.push_back(true);
            }
        }
    }
    const std::size_t m = free.size();
    bool anyRemoved = false;
    for (const bool f : free) {
        anyRemoved = anyRemoved || !f;
    }
    if (!anyRemoved) {
        return false;  // the whole boundary survives: one component
    }
    // Count maximal cyclic runs of free pieces: starting the scan at a removed
    // piece keeps the runs from wrapping. Two runs means two components.
    std::size_t start = 0;
    while (free[start]) {
        ++start;
    }
    int components = 0;
    std::size_t i = 0;
    while (i < m) {
        if (!free[(start + i) % m]) {
            ++i;
            continue;
        }
        while (i < m && free[(start + i) % m]) {
            ++i;
        }
        ++components;
    }
    return components >= 2;
}
 
// Removing a bounded shape from the region: empty region -> false; degenerate
// region -> the remover cuts its carrier; full-dimensional region -> the
// boundary-component test (bounded regions delegate to Convex::separates for
// speed and reuse, which is the same predicate).
template <class Shape2, class Region>
bool boundedShapeSeparatesRegion(const Shape2& remover, const Region& region) {
    if (region.empty()) {
        return false;
    }
    if (region.isDegenerate()) {
        return std::visit([&remover](const auto& carrier) { return remover.separates(carrier); },
                          degenerateRegionCarrier(region));
    }
    using E = region_exact_number_t<typename Region::NumberType>;
    if (region.isBounded()) {
        return remover.separates(region.template asConvex<E>());
    }
    return boundedRemoverSeparatesRegion(region, remover);
}

template <class RemoverRegion, class Region>
bool regionSeparatesRegion(const RemoverRegion& remover, const Region& region) {
    using E = region_exact_number_t<typename RemoverRegion::NumberType>;
    using EPoint = Point<E, typename RemoverRegion::PointType::LabelType>;
    const HalfplaneIntersection<EPoint> removerExact(remover);
    const HalfplaneIntersection<EPoint> regionExact(region);
 
    bool started = false;
    E loX{}, loY{}, hiX{}, hiY{};
    const auto include = [&](const EPoint& point) {
        if (!started) {
            loX = hiX = point.x();
            loY = hiY = point.y();
            started = true;
            return;
        }
        loX = std::min(loX, point.x());
        loY = std::min(loY, point.y());
        hiX = std::max(hiX, point.x());
        hiY = std::max(hiY, point.y());
    };
 
    std::vector<Line<EPoint>> lines;
    for (const auto& regionSide : {&removerExact, &regionExact}) {
        for (const auto& halfplane : *regionSide) {
            include(halfplane.source());
            include(halfplane.target());
            lines.push_back(halfplane.asLine());
        }
    }
    for (std::size_t i = 0; i < lines.size(); ++i) {
        for (std::size_t j = i + 1; j < lines.size(); ++j) {
            const auto crossing = lines[i].template intersection<E>(lines[j]);
            if (crossing && crossing->index() == 0) {
                include(std::get<0>(*crossing));
            }
        }
    }
 
    // Both regions are full-dimensional and nonempty, so each stores at least
    // one half-plane and the coordinate pool is nonempty.
    const Rectangle<EPoint> bounds(EPoint(loX, loY), EPoint(hiX, hiY));
    const auto removerPoly = regionClippedToBox(removerExact, bounds).template asConvex<E>();
    const auto regionPoly = regionClippedToBox(regionExact, bounds).template asConvex<E>();
    return removerPoly.separates(regionPoly);
}
 
}  // namespace detail
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherRegion& other) const {
    // A degenerate side reduces to its carrier shape; the full-dimensional
    // pair takes the shared-box boundary-component count. Both are computed
    // over the exact coordinate type, since exact carriers and boxes are
    // generally not representable in the inputs' own coordinate type.
    if (empty() || other.empty()) {
        return false;
    }
    using E = detail::region_exact_number_t<NumberType>;
    using ERegion = HalfplaneIntersection<Point<E, typename PointType::LabelType>>;
    if (isDegenerate()) {
        const ERegion otherExact(other);
        return std::visit(
            [&otherExact](const auto& carrier) { return carrier.separates(otherExact); },
            detail::degenerateRegionCarrier(*this));
    }
    if (other.isDegenerate()) {
        const ERegion selfExact(*this);
        return std::visit(
            [&selfExact](const auto& carrier) { return selfExact.separates(carrier); },
            detail::degenerateRegionCarrier(other));
    }
    return detail::regionSeparatesRegion(*this, other);
}
 
template <class PointType, class LabelType>
template <PointConcept OtherPoint>
constexpr bool HalfplaneIntersection<PointType, LabelType>::separates(const Shape<OtherPoint>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}
 
template <class Number, class Label>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Point<Number, Label>::separates(const OtherRegion& other) const {
    // A point has empty interior; it can only sever a one-dimensional carrier.
    if (other.empty() || !other.isDegenerate()) {
        return false;
    }
    return std::visit([this](const auto& carrier) { return this->separates(carrier); },
                      detail::degenerateRegionCarrier(other));
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Segment<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (other.empty()) {
        return false;
    }
    if (other.isDegenerate()) {
        // A one-dimensional region is cut along its carrier.
        return std::visit([this](const auto& carrier) { return this->separates(carrier); },
                          detail::degenerateRegionCarrier(other));
    }
    // The same chord criterion as Segment::separates(Convex), which holds for
    // every convex region, bounded or not. Meeting the interior forces the
    // region to have interior strictly on both sides of the supporting line,
    // so the region's chord on that line is the line's interior trace; both
    // endpoints staying out of the interior then forces the segment to cover
    // that whole chord (a segment ending midway leaves a slit, not a split),
    // and a covered chord is necessarily bounded.
    return !isDegenerate() && !other.interiorContains(min()) &&
           !other.interiorContains(max()) && other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool OrientedSegment<PointType, LabelType>::separates(const OtherRegion& other) const {
    return asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Line<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (other.empty()) {
        return false;
    }
    if (other.isDegenerate()) {
        return std::visit([this](const auto& carrier) { return this->separates(carrier); },
                          detail::degenerateRegionCarrier(other));
    }
    // A full line severs a full-dimensional convex region exactly when the
    // region has interior on both sides of it.
    return detail::regionStrictlyOnBothSides(other, (*this)[0], (*this)[1]);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool OrientedLine<PointType, LabelType>::separates(const OtherRegion& other) const {
    return asLine().separates(other);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Ray<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (other.empty()) {
        return false;
    }
    if (other.isDegenerate()) {
        return std::visit([this](const auto& carrier) { return this->separates(carrier); },
                          detail::degenerateRegionCarrier(other));
    }
    // The same chord criterion as Ray::separates(Convex) with the far end at
    // infinity: only the source can leave the chord uncovered, since a source
    // outside the region sits before the whole chord and a source on the
    // boundary sits at its near end. An interior source leaves a slit.
    return !isDegenerate() && !other.interiorContains(source()) &&
           other.interiorsIntersect(*this);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Halfplane<PointType, LabelType>::separates(const OtherRegion& other) const {
    // Removing a closed half-plane from any convex set leaves a convex set,
    // which is connected; so a half-plane never disconnects the region.
    (void)other;
    return false;
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Rectangle<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    if (isDegenerate()) {
        return false;  // a degenerate rectangle is one-dimensional
    }
    return detail::boundedShapeSeparatesRegion(*this, other);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Triangle<PointType, LabelType>::separates(const OtherRegion& other) const {
    return detail::boundedShapeSeparatesRegion(*this, other);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Disk<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (isDegenerate()) {
        return false;
    }
    return detail::boundedShapeSeparatesRegion(*this, other);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Convex<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (isDegenerate()) {
        return false;
    }
    return detail::boundedShapeSeparatesRegion(*this, other);
}
 
template <class PointType, class LabelType, class Storage>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherRegion& other) const {
    if (other.empty()) {
        return false;
    }
    if (other.isDegenerate()) {
        return std::visit([this](const auto& carrier) { return this->separates(carrier); },
                          detail::degenerateRegionCarrier(other));
    }
    // The region is convex, so the crosscut scan applies as it does for a
    // Convex target; the bbox guards used there are skipped because an
    // unbounded region has no bounding box. An x-monotone chain never closes
    // a loop, so no sealed pocket escapes the scan.
    return separatesTwoDimensional(other);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Polyline<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (other.empty() || size() < 2) {
        return false;
    }
    if (other.isDegenerate()) {
        return std::visit([this](const auto& carrier) { return this->separates(carrier); },
                          detail::degenerateRegionCarrier(other));
    }
    // The region is convex, so the same cycle criterion as for a Convex
    // target applies; it is predicate-based and never asks for a bounding box.
    return detail::polylineSeparatesConvexRegion(*this, other);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherRegion>
constexpr bool Polygon<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (isDegenerate()) {
        return false;
    }
    return detail::boundedShapeSeparatesRegion(*this, other);
}
 
// ---------------------------------------------------------------------------
// PolygonWithHoles
//
// Every cut predicate a region takes part in — in either direction — is the
// same question: does `target ∖ remover` have at least two connected
// components? For a region neither operand's shape can be leaned on the way
// the polygon overloads above lean on theirs. A region's interior is not
// simply connected, so a crosscut of its boundary is not a cut; that interior
// need not even be connected — a band hole leaves two slabs — so a
// boundary-arc tally does not see all of its pieces; and the region carries
// pinch points and slits (decision (b) of the design), so it is not a surface
// and no Euler-characteristic count applies to it either.
//
// What survives all of that is the plainest reading of the question. Both
// operands are polygonal, so the arrangement of their boundaries cuts the
// plane into cells — open faces, open edges and vertices — on each of which
// membership in *both* operands is constant, and `target ∖ remover` is exactly
// the union of the kept cells. Two kept cells are adjacent in the remainder
// exactly when one is a face of the other (a cell's closure meets another kept
// cell only along shared faces, and `c ∪ f` is connected for a face `f` of
// `c`), so a union-find over the incidences counts the components with nothing
// assumed about either shape. Vertices and edges are cells in their own right,
// which is what keeps slits — region material with no area beside it — and
// pinch points from being lost.
//
// The construction is the @ref pgl::Arrangement of every boundary edge of both
// operands, which is that decomposition and carries the incidences directly.
// Coordinates are exact rationals, since the crossings are. A region without
// holes *is* its outer polygon, and forwards to it: that keeps the tested
// polygon implementations in charge of everything they already settle, and
// leaves the engine for the cases only a region reaches.
 
namespace detail {

template <class ExactPoint, class Shape>
void appendCutSegments(const Shape& shape, std::vector<Segment<ExactPoint>>& out) {
    const auto add = [&out](const auto& edge) {
        const ExactPoint lo(edge.min());
        const ExactPoint hi(edge.max());
        if (lo != hi) {
            out.emplace_back(lo, hi);
        }
    };
    if constexpr (is_point_v<Shape>) {
        (void)shape;
        (void)add;
    } else if constexpr (is_segment_v<Shape>) {
        add(shape);
    } else if constexpr (requires { shape.edgesView(); }) {
        for (const auto& edge : shape.edgesView()) {
            add(edge);
        }
    } else {
        // Everything else materializes its edges: a region and a set assemble
        // theirs from their rings, and a triangle and a rectangle return a
        // fixed-size array, which is why neither carries a lazy view.
        for (const auto& edge : shape.edges()) {
            add(edge);
        }
    }
}

template <class ExactPoint, class Shape>
void appendCutPoints(const Shape& shape, std::vector<ExactPoint>& out) {
    if constexpr (is_point_v<Shape>) {
        out.emplace_back(shape);
    } else if constexpr (is_segment_v<Shape>) {
        if (shape.min() == shape.max()) {
            out.emplace_back(shape.min());
        }
    } else if constexpr (is_monotone_chain_v<Shape> || is_polyline_v<Shape>) {
        if (shape.size() == 1) {
            out.emplace_back(shape[0]);
        }
    } else {
        (void)shape;
        (void)out;
    }
}

template <class Target>
bool disconnectedOnItsOwn(const Target& target) {
    if constexpr (is_polygon_set_v<Target>) {
        return !target.isConnected();
    } else {
        (void)target;
        return false;
    }
}

template <class Target, class Remover>
bool cellSeparates(const Target& target, const Remover& remover) {
    using ExactNumber = Exact1DNumber<typename Target::NumberType, typename Remover::NumberType>;
    using ExactPoint = Point<ExactNumber>;
    using ExactSegment = Segment<ExactPoint>;
 
    // Every operand the engine is called with is connected — a region included,
    // see @ref regionsAreConnected — so a remover that misses the target
    // removes nothing that matters. A set of regions is the exception: it may
    // come apart on its own, and then no remover has to touch it at all.
    if (!target.bbox().intersects(remover.bbox())) {
        return disconnectedOnItsOwn(target);
    }
 
    std::vector<ExactSegment> cuts;
    std::vector<ExactPoint> cutPoints;
    appendCutSegments<ExactPoint>(target, cuts);
    appendCutSegments<ExactPoint>(remover, cuts);
    appendCutPoints<ExactPoint>(target, cutPoints);
    appendCutPoints<ExactPoint>(remover, cutPoints);
    if (cuts.empty()) {
        return false;  // a target with no extent holds at most one component
    }
 
    // Layering: Arrangement lives in algorithm/arrangement.hpp, which pgl.hpp
    // includes after this header. Only the name is available here — see the
    // declaration in core/forward.hpp — and that is enough, since the type
    // depends on the operands and so is instantiated at the call site.
    const Arrangement<ExactPoint> arrangement(cuts, cutPoints);
    using FaceHandle = typename Arrangement<ExactPoint>::FaceId;
    using HalfedgeHandle = typename Arrangement<ExactPoint>::HalfedgeId;
 
    // One union-find node per kept cell. A cell is kept when its relative
    // interior lies in target ∖ remover, which one witness point decides.
    std::vector<std::size_t> parent;
    const auto findRoot = [&parent](std::size_t x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    };
    constexpr std::size_t dropped = std::numeric_limits<std::size_t>::max();
    const auto unite = [&](std::size_t a, std::size_t b) {
        if (a != dropped && b != dropped) {
            parent[findRoot(a)] = findRoot(b);
        }
    };
    const auto cellOf = [&](const ExactPoint& witness) {
        if (!target.contains(witness) || remover.contains(witness)) {
            return dropped;
        }
        parent.push_back(parent.size());
        return parent.size() - 1;
    };
 
    std::vector<std::size_t> faceCell(arrangement.faceCount(), dropped);
    for (std::uint32_t i = 0; i < arrangement.faceCount(); ++i) {
        const FaceHandle f(i);
        if (!arrangement.isUnbounded(f)) {
            faceCell[f.index()] = cellOf(arrangement.template witness<ExactNumber>(f));
        }
    }
    // A vertex the input left isolated — a cut point lying on no cut segment —
    // gets a cell of its own that the loop below never joins to anything, which
    // is right because it only arises when an operand is a single point: as the
    // target it is then the one cell that can survive at all, and as the remover
    // it is dropped along with every other point of the remover.
    std::vector<std::size_t> vertexCell(arrangement.vertexCount(), dropped);
    for (std::size_t i = 0; i < arrangement.vertexCount(); ++i) {
        vertexCell[i] = cellOf(arrangement.vertices()[i]);  // a vertex is its own witness
    }
 
    std::vector<std::size_t> edgeCell(arrangement.edgeCount(), dropped);
    for (std::uint32_t i = 0; i < arrangement.halfedgeCount(); ++i) {
        const HalfedgeHandle h(i);
        const std::size_t edge = h.index() / 2;
        if (h.index() % 2 == 0) {
            edgeCell[edge] = cellOf(arrangement.template witness<ExactNumber>(h));
        }
        const std::size_t side = edgeCell[edge];
        const std::size_t corner = vertexCell[arrangement.source(h).index()];
        const std::size_t face = faceCell[arrangement.face(h).index()];
        unite(side, corner);
        unite(side, face);
        unite(face, corner);
    }
 
    std::size_t components = 0;
    for (std::size_t i = 0; i < parent.size(); ++i) {
        if (findRoot(i) == i && ++components >= 2) {
            return true;
        }
    }
    return false;
}



template <class ExactPoint, class Linear, class Region>
std::optional<Segment<ExactPoint>> clipLinearToRegion(const Linear& linear, const Region& region) {
    using ExactNumber = typename ExactPoint::NumberType;
    const auto bounds = region.bbox();
    const ExactNumber margin(1);
    const Rectangle<ExactPoint> box(
        ExactPoint(ExactNumber(bounds.min().x()) - margin, ExactNumber(bounds.min().y()) - margin),
        ExactPoint(ExactNumber(bounds.max().x()) + margin, ExactNumber(bounds.max().y()) + margin));
    const auto piece = box.template intersection<ExactNumber>(linear);
    if (!piece) {
        return std::nullopt;
    }
    const auto* overlap = std::get_if<1>(&*piece);
    if (!overlap) {
        return std::nullopt;  // a single corner touch stays clear of the region
    }
    return *overlap;
}


template <bool LinearIsTarget, class Linear, class Region>
bool linearAndRegionSeparate(const Linear& linear, const Region& region) {
    using ExactNumber = Exact1DNumber<typename Linear::NumberType, typename Region::NumberType>;
    using ExactPoint = Point<ExactNumber>;
    const auto clip = clipLinearToRegion<ExactPoint>(linear, region);
    if (!clip) {
        // The operand stays clear of the region; only a region that was already
        // in several pieces answers anything but false.
        return LinearIsTarget ? false : disconnectedOnItsOwn(region);
    }
    if constexpr (LinearIsTarget) {
        return cellSeparates(*clip, region);
    } else {
        return cellSeparates(region, *clip);
    }
}

template <bool ConvexIsTarget, class ConvexOperand, class Region>
bool convexAndRegionSeparate(const ConvexOperand& convex, const Region& region) {
    using ExactNumber = Exact1DNumber<typename ConvexOperand::NumberType, typename Region::NumberType>;
    using ExactPoint = Point<ExactNumber>;
    HalfplaneIntersection<ExactPoint> clipped(convex);
    if (!clipped.isBounded()) {
        clipped = regionClippedToBox(clipped, region.bbox());
    }
    const auto run = [&region](const auto& operand) {
        if constexpr (ConvexIsTarget) {
            return cellSeparates(operand, region);
        } else {
            return cellSeparates(region, operand);
        }
    };
    if (clipped.empty()) {
        // As in @ref linearAndRegionSeparate: the operand cannot reach the
        // region, so only a region already in several pieces answers true.
        return ConvexIsTarget ? false : disconnectedOnItsOwn(region);
    }
    if (clipped.isDegenerate()) {
        return std::visit(
            [&run](const auto& carrier) {
                using Carrier = std::remove_cvref_t<decltype(carrier)>;
                if constexpr (is_point_v<Carrier> || is_segment_v<Carrier>) {
                    return run(carrier);
                } else {
                    return false;  // a bounded region carries no ray or line
                }
            },
            degenerateRegionCarrier(clipped));
    }
    return run(clipped.template asConvex<ExactNumber>().asPolygon());
}
 
// ---------------------------------------------------------------------------
// The disk pair.
//
// A circle is not polygonal, so the cell engine does not take it: it has no
// boundary to add to the arrangement, being neither a segment nor a source of
// rational crossings. What replaces it in each direction is the one property a
// disk has that a polygon does not — convexity — and each direction uses it
// differently. Neither ever constructs a circle crossing: every test below is a
// predicate on the disk, exact in the operands' own arithmetic.
//
// Convexity is also why the direction below keeps a triangulation where
// @ref cellSeparates dropped one. Its cells must each meet the disk in a
// connected set, and a general arrangement face need not: an L-shaped face and
// a disk covering both of its ends but not its corner meet in two pieces. A
// triangle cannot do that, so the arrangement's faces are cut into triangles
// before the union-find walks them, and only the split of the constraints —
// which the arrangement does by sweep — is shared with the polygonal engine.

template <class ExactPoint>
std::vector<Segment<ExactPoint>> splitCutSegments(const std::vector<Segment<ExactPoint>>& cuts) {
    // Layering: see the note in @ref cellSeparates.
    const Arrangement<ExactPoint> arrangement(cuts);
    return arrangement.boundedEdges();
}

template <class Region, class OtherDisk>
bool regionSeparatesDisk(const Region& region, const OtherDisk& disk) {
    using ExactNumber = Exact1DNumber<typename Region::NumberType, typename OtherDisk::NumberType>;
    using ExactPoint = Point<ExactNumber>;
    using ExactSegment = Segment<ExactPoint>;
 
    // A disk is connected, so a region that misses it removes nothing that
    // matters — the same shortcut @ref cellSeparates takes, and worth as much
    // here, since everything below is one arrangement and one triangulation.
    if (!region.bbox().intersects(disk.bbox())) {
        return false;
    }
 
    std::vector<ExactSegment> cuts;
    appendCutSegments<ExactPoint>(region, cuts);
    if (cuts.empty()) {
        return false;  // an empty region removes nothing
    }
 
    // A box strictly containing both shapes. Holding the whole disk is what
    // keeps every piece of `D ∖ A` inside the arrangement, and keeping the disk
    // off the box boundary is what leaves every vertex it reaches with a
    // complete star of triangles for the union-find to walk.
    const auto regionBox = region.bbox();
    const auto diskBox = disk.bbox();
    const ExactNumber margin(1);
    const ExactNumber loX =
        std::min(ExactNumber(regionBox.min().x()), ExactNumber(diskBox.min().x())) - margin;
    const ExactNumber loY =
        std::min(ExactNumber(regionBox.min().y()), ExactNumber(diskBox.min().y())) - margin;
    const ExactNumber hiX =
        std::max(ExactNumber(regionBox.max().x()), ExactNumber(diskBox.max().x())) + margin;
    const ExactNumber hiY =
        std::max(ExactNumber(regionBox.max().y()), ExactNumber(diskBox.max().y())) + margin;
    const Polygon<ExactPoint> box(std::vector<ExactPoint>{
        ExactPoint(loX, loY), ExactPoint(hiX, loY), ExactPoint(hiX, hiY), ExactPoint(loX, hiY)});
    const auto mesh = box.triangulation(splitCutSegments(cuts));
 
    // The disk in the arrangement's own arithmetic, so every test below is a
    // like-typed exact predicate.
    const Disk<ExactPoint> exactDisk(disk);
 
    std::vector<std::size_t> parent;
    const auto findRoot = [&parent](std::size_t x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    };
    const auto unite = [&](std::size_t a, std::size_t b) {
        parent[findRoot(a)] = findRoot(b);
    };
    // `dropped` marks a cell the disk does not reach and `unset` one not yet
    // classified. Real cell numbers are positions in `parent`, so neither
    // collides with one.
    constexpr std::size_t dropped = std::numeric_limits<std::size_t>::max();
    constexpr std::size_t unset = dropped - 1;
    const auto cellOf = [&](bool kept) {
        if (!kept) {
            return dropped;
        }
        parent.push_back(parent.size());
        return parent.size() - 1;
    };
 
    // Cells are keyed by handle rather than by position: a vertex indexes a
    // plain table, and an edge is named by its two vertex handles. Keyed by
    // coordinates instead, as the arrangement's own arithmetic makes them, every
    // step of every lookup would compare a pair of exact rationals.
    using TriId = typename std::decay_t<decltype(mesh)>::TriId;
    using VertexId = typename std::decay_t<decltype(mesh)>::VertexId;
    using EdgeKey = std::pair<VertexId, VertexId>;
    std::map<EdgeKey, std::size_t> edgeCells;
    std::vector<std::size_t> vertexCells(mesh.vertexIndexBound(), unset);
    const auto vertexCell = [&](VertexId v) {
        std::size_t& cell = vertexCells[v.index()];
        if (cell == unset) {
            const auto& vertex = mesh.getShape(v);
            cell = cellOf(exactDisk.contains(vertex) && !region.contains(vertex));
        }
        return cell;
    };
    const auto edgeCell = [&](VertexId a, VertexId b) {
        const EdgeKey key = a < b ? EdgeKey{a, b} : EdgeKey{b, a};
        const auto [it, inserted] = edgeCells.try_emplace(key, dropped);
        if (inserted) {
            // The disk reaches the open edge either through its own interior or
            // by touching it at a single point, and that point is interior to
            // the edge exactly when neither endpoint carries it. Membership in
            // the region is constant along the open edge, so its midpoint
            // decides it — an edge along a slit is region material, not
            // complement.
            const ExactSegment edge(mesh.getShape(a), mesh.getShape(b));
            const ExactPoint mid((edge.min().x() + edge.max().x()) / ExactNumber(2),
                                 (edge.min().y() + edge.max().y()) / ExactNumber(2));
            const bool met = exactDisk.interiorsIntersect(edge) ||
                             (exactDisk.intersects(edge) && !exactDisk.contains(edge.min()) &&
                              !exactDisk.contains(edge.max()));
            it->second = cellOf(met && !region.contains(mid));
        }
        return it->second;
    };
 
    mesh.visitTriangles([&](TriId t) {
        const auto triangle = mesh.getShape(t);
        const auto v = mesh.vertices(t);
        const std::size_t face =
            cellOf(exactDisk.interiorsIntersect(triangle) &&
                   !region.contains(triangle.template pointInside<ExactNumber>()));
        for (int k = 0; k < 3; ++k) {
            // vertices(t)[k] is the vertex at getShape(t)[k], so side k joins
            // the same two corners the value walk paired.
            const std::size_t side = edgeCell(v[k], v[(k + 1) % 3]);
            const std::size_t corner = vertexCell(v[k]);
            if (face != dropped && side != dropped) {
                unite(face, side);
            }
            if (face != dropped && corner != dropped) {
                unite(face, corner);
            }
            if (side != dropped && corner != dropped) {
                unite(side, corner);
            }
            const std::size_t next = vertexCell(v[(k + 1) % 3]);
            if (side != dropped && next != dropped) {
                unite(side, next);
            }
        }
    });
 
    std::size_t components = 0;
    for (std::size_t i = 0; i < parent.size(); ++i) {
        if (findRoot(i) == i && ++components >= 2) {
            return true;
        }
    }
    return false;
}

template <class Region>
std::vector<Segment<typename Region::PointType>> regionSlits(const Region& region) {
    using RegionPoint = typename Region::PointType;
    using RegionSegment = Segment<RegionPoint>;
 
    const auto edges = region.edges();
    std::vector<RegionPoint> vertices = region.vertices();
    std::sort(vertices.begin(), vertices.end());
    vertices.erase(std::unique(vertices.begin(), vertices.end()), vertices.end());
 
    std::vector<RegionSegment> slits;
    std::vector<RegionPoint> cuts;
    for (const auto& edge : edges) {
        cuts.clear();
        cuts.push_back(edge.min());
        cuts.push_back(edge.max());
        for (const auto& vertex : vertices) {
            if (edge.contains(vertex)) {
                cuts.push_back(vertex);
            }
        }
        // Every cut lies on the edge, so the lexicographic point order is the
        // linear order along it.
        std::sort(cuts.begin(), cuts.end());
        cuts.erase(std::unique(cuts.begin(), cuts.end()), cuts.end());
        for (std::size_t k = 0; k + 1 < cuts.size(); ++k) {
            const RegionSegment piece(cuts[k], cuts[k + 1]);
            int covers = 0;
            for (const auto& cover : edges) {
                if (cover.contains(piece.min()) && cover.contains(piece.max()) && ++covers == 2) {
                    slits.push_back(piece);
                    break;
                }
            }
        }
    }
    std::sort(slits.begin(), slits.end());
    slits.erase(std::unique(slits.begin(), slits.end()), slits.end());
    return slits;
}

template <class OtherDisk, class Region>
bool diskSeparatesRegion(const OtherDisk& disk, const Region& region) {
    using RegionPoint = typename Region::PointType;
    using RegionSegment = Segment<RegionPoint>;
 
    // A region is connected (@ref regionsAreConnected), so a disk that misses
    // it removes nothing that matters, and the triangulation below is not worth
    // building. The same shortcut @ref cellSeparates takes.
    if (!disk.bbox().intersects(region.bbox())) {
        return disconnectedOnItsOwn(region);
    }
 
    std::vector<std::size_t> parent;
    const auto findRoot = [&parent](std::size_t x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    };
    const auto unite = [&](std::size_t a, std::size_t b) {
        parent[findRoot(a)] = findRoot(b);
    };
 
    // Layering: Triangulation lives in algorithm/triangulation.hpp, which
    // pgl.hpp includes after this header. Reach it only through the dependent
    // call region.triangulation(), never by naming the class here.
    const auto mesh = region.triangulation();
    using TriId = typename std::decay_t<decltype(mesh)>::TriId;
    using VertexId = typename std::decay_t<decltype(mesh)>::VertexId;
 
    // One node per vertex the disk leaves behind; a vertex it swallows carries
    // nothing to connect, and is never handed to @c nodeOf. Nodes are kept by
    // vertex handle, so the walk indexes a table instead of searching a map
    // keyed by coordinates.
    constexpr std::size_t none = std::numeric_limits<std::size_t>::max();
    std::vector<std::size_t> nodes(mesh.vertexIndexBound(), none);
    const auto nodeOf = [&](VertexId v) {
        std::size_t& node = nodes[v.index()];
        if (node == none) {
            node = parent.size();
            parent.push_back(parent.size());
        }
        return node;
    };
 
    // Whether the disk swallows a vertex is a property of the vertex, so it is
    // decided once rather than once per incident triangle.
    std::vector<char> eatenAt(mesh.vertexIndexBound(), 0);
    for (const VertexId v : mesh.vertexIds()) {
        eatenAt[v.index()] = disk.contains(mesh.getShape(v)) ? 1 : 0;
    }
 
    mesh.visitTriangles([&](TriId t) {
        // vertices(t)[i] is the vertex at getShape(t)[i], so side i runs from
        // vertex i to vertex i + 1 as the value walk had it.
        const auto v = mesh.vertices(t);
        const bool eaten[3] = {eatenAt[v[0].index()] != 0, eatenAt[v[1].index()] != 0,
                               eatenAt[v[2].index()] != 0};
        std::size_t first = none;
        std::size_t last = none;
        bool gapBeforeFirst = false;
        bool gap = false;  // a contact with the disk since the last vertex kept
        const auto keep = [&](VertexId vertex) {
            const std::size_t node = nodeOf(vertex);
            if (last == none) {
                first = node;
                gapBeforeFirst = gap;
            } else if (!gap) {
                unite(last, node);
            }
            last = node;
            gap = false;
        };
        for (int i = 0; i < 3; ++i) {
            const VertexId from = v[i];
            const VertexId to = v[(i + 1) % 3];
            if (eaten[i] && eaten[(i + 1) % 3]) {  // the disk is convex: the edge is gone
                gap = true;
            } else if (eaten[i]) {                 // the walk leaves the disk mid-edge
                gap = true;
                keep(to);
            } else if (eaten[(i + 1) % 3]) {       // the walk enters the disk mid-edge
                keep(from);
                gap = true;
            } else if (disk.intersects(
                           RegionSegment(mesh.getShape(from), mesh.getShape(to)))) {
                keep(from);                        // a contact inside the edge
                gap = true;
                keep(to);
            } else {                               // the edge is clear of the disk
                keep(from);
                keep(to);
            }
        }
        if (first != none && !gap && !gapBeforeFirst) {
            unite(last, first);
        }
    });
 
    const auto slits = regionSlits(region);
    if (!slits.empty()) {
        // A slit endpoint is a ring vertex, hence a vertex of the mesh, but the
        // tip of a spur carries no triangle at all -- locating it would not find
        // it. So the handles come from a lookup over the stored vertices, built
        // once and only for a region that has slits in the first place.
        std::map<RegionPoint, VertexId> vertexAt;
        for (const VertexId v : mesh.vertexIds()) {
            vertexAt.emplace(mesh.getShape(v), v);
        }
        const auto nodeAtPoint = [&](const RegionPoint& p) {
            const auto it = vertexAt.find(p);
            return it == vertexAt.end() ? none : nodeOf(it->second);
        };
        for (const auto& slit : slits) {
            // Whatever the disk takes out of a slit, what is left hangs from an
            // endpoint — it meets the segment in a single subsegment — so each end
            // it spares is a piece of its own, and the two stay joined only when it
            // misses the slit altogether.
            const std::size_t low = disk.contains(slit.min()) ? none : nodeAtPoint(slit.min());
            const std::size_t high = disk.contains(slit.max()) ? none : nodeAtPoint(slit.max());
            if (low != none && high != none && !disk.intersects(slit)) {
                unite(low, high);
            }
        }
    }
 
    std::size_t components = 0;
    for (std::size_t i = 0; i < parent.size(); ++i) {
        if (findRoot(i) == i && ++components >= 2) {
            return true;
        }
    }
    return false;
}
 
}  // namespace detail

 
template <class PointType, class LabelType>
template <PointConcept OtherPoint>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherPoint&) const {
    return false;  // removing anything from a point leaves at most one piece
}
 
template <class PointType, class LabelType>
template <SegmentConcept OtherSegment>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherSegment& other) const {
    if (!hasHoles()) {
        return outer_.separates(other);
    }
    if (other.isDegenerate()) {
        return false;
    }
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <OrientedSegmentConcept OtherOrientedSegment>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherOrientedSegment& other) const {
    return separates(other.asSegment());
}
 
template <class PointType, class LabelType>
template <LineConcept OtherLine>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherLine& other) const {
    if (!hasHoles()) {
        return outer_.separates(other);
    }
    if (other.isDegenerate()) {
        return false;
    }
    return detail::linearAndRegionSeparate<true>(other, *this);
}
 
template <class PointType, class LabelType>
template <OrientedLineConcept OtherOrientedLine>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherOrientedLine& other) const {
    return separates(other.asLine());
}
 
template <class PointType, class LabelType>
template <RayConcept OtherRay>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherRay& other) const {
    if (!hasHoles()) {
        return outer_.separates(other);
    }
    if (other.isDegenerate()) {
        return false;
    }
    return detail::linearAndRegionSeparate<true>(other, *this);
}
 
template <class PointType, class LabelType>
template <HalfplaneConcept OtherHalfplane>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherHalfplane& other) const {
    if (!hasHoles()) {
        return outer_.separates(other);
    }
    if (other.isUndefined()) {
        return false;
    }
    return detail::convexAndRegionSeparate<true>(other, *this);
}
 
template <class PointType, class LabelType>
template <RectangleConcept OtherRectangle>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherRectangle& other) const {
    if (other.empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return separates(other.asPolygon());
}
 
template <class PointType, class LabelType>
template <TriangleConcept OtherTriangle>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherTriangle& other) const {
    return separates(other.asPolygon());
}
 
template <class PointType, class LabelType>
template <ConvexConcept OtherConvex>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherConvex& other) const {
    return separates(other.asPolygon());
}
 
template <class PointType, class LabelType>
template <PolygonConcept OtherPolygon>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherPolygon& other) const {
    if (!hasHoles()) {
        return outer_.separates(other);
    }
    if (other.isDegenerate()) {
        return false;
    }
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (!hasHoles()) {
        return outer_.separates(other);
    }
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    if (other.isDegenerate()) {
        return false;
    }
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <MonotoneChainConcept OtherChain>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherChain& other) const {
    if (!hasHoles()) {
        return outer_.separates(other);
    }
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PolylineConcept OtherPolyline>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherPolyline& other) const {
    if (!hasHoles()) {
        return outer_.separates(other);
    }
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <DiskConcept OtherDisk>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherDisk& other) const {
    if (!hasHoles()) {
        return outer_.separates(other);
    }
    if (other.isDegenerate()) {
        return false;  // a disk of radius zero is a point, and an undefined one has no circle
    }
    return detail::regionSeparatesDisk(*this, other);
}
 
template <class PointType, class LabelType>
template <HalfplaneIntersectionConcept OtherIntersection>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherIntersection& other) const {
    if (!hasHoles()) {
        return outer_.separates(other);
    }
    return detail::convexAndRegionSeparate<true>(other, *this);
}
 
// ---------------------------------------------------------------------------
// Reverse direction: removing a lower-ranked shape from the region.
//
// Every one of these is a real computation: a region is cut by shapes that
// cannot cut a simply connected target — a single point at a pinch, or a
// segment run from one hole to another — and no operand can be dismissed for
// its shape alone the way a half-plane can be dismissed against a convex
// region.
 
template <class Number, class Label>
template <PolygonWithHolesConcept OtherRegion>
bool Point<Number, Label>::separates(const OtherRegion& other) const {
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    // No single point cuts a region, however pinched. By the same duality as
    // @ref regionsAreConnected, the components of `A ∖ {p}` are counted by the
    // first cohomology of `(S² ∖ A) ∪ {p}`; each complement component is a
    // Jordan domain, so adding one of its boundary points keeps it simply
    // connected, and the union is then a wedge of simply connected spaces at
    // p. It is trivial, so what is left of the region is one piece — a pinch
    // point always has some other way around it, along the ring boundaries
    // that meet there.
    return false;
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool Segment<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    if (isDegenerate()) {
        return min().separates(other);
    }
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool OrientedSegment<PointType, LabelType>::separates(const OtherRegion& other) const {
    return asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool Line<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    if (isDegenerate()) {
        return false;
    }
    return detail::linearAndRegionSeparate<false>(*this, other);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool OrientedLine<PointType, LabelType>::separates(const OtherRegion& other) const {
    return asLine().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool Ray<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    if (isDegenerate()) {
        return false;
    }
    return detail::linearAndRegionSeparate<false>(*this, other);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool Halfplane<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    if (isUndefined()) {
        return false;
    }
    return detail::convexAndRegionSeparate<false>(*this, other);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool Rectangle<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (empty()) {
        // The empty set meets nothing and disconnects nothing.
        return false;
    }
    return asPolygon().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool Triangle<PointType, LabelType>::separates(const OtherRegion& other) const {
    return asPolygon().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool Convex<PointType, LabelType>::separates(const OtherRegion& other) const {
    return asPolygon().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool Polygon<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    if (isDegenerate()) {
        return false;
    }
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType, class Storage>
template <PolygonWithHolesConcept OtherRegion>
bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherRegion& other) const {
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool Polyline<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherRegion>
bool Disk<PointType, LabelType>::separates(const OtherRegion& other) const {
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    if (isDegenerate()) {
        // A disk of radius zero is a point, which never cuts a region (see
        // Point::separates above); an undefined disk determines no circle.
        return false;
    }
    return detail::diskSeparatesRegion(*this, other);
}
 
template <class PointType, class LabelType>
template <PolygonWithHolesConcept OtherHoledRegion>
bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherHoledRegion& other) const {
    if (!other.hasHoles()) {
        return separates(other.outer());
    }
    if (empty()) {
        return false;  // nothing is removed from a connected region
    }
    return detail::convexAndRegionSeparate<false>(*this, other);
}
 
// ---------------------------------------------------------------------------
// Runtime Shape argument: unwrap the stored alternative and re-dispatch. Every
// alternative has a per-shape overload above, so no fallback is needed.
 
template <class PointType, class LabelType>
template <PointConcept OtherPoint>
constexpr bool PolygonWithHoles<PointType, LabelType>::separates(const Shape<OtherPoint>& other) const {
    return std::visit(
        [this](const auto& value) {
            return this->separates(value);
        },
        other.variant());
}
 
 
template <class PointType, class LabelType>
bool PolygonSet<PointType, LabelType>::isConnected() const {
    // Each component is connected on its own, so the set is connected exactly
    // when the graph joining components that meet is — a union-find over the
    // pairs whose boxes overlap.
    if (components_.size() < 2) {
        return true;
    }
    std::vector<std::size_t> parent(components_.size());
    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::size_t pieces = components_.size();
    for (std::size_t i = 0; i < components_.size(); ++i) {
        for (std::size_t j = i + 1; j < components_.size(); ++j) {
            if (findRoot(i) == findRoot(j) ||
                !components_[i].bbox().intersects(components_[j].bbox())) {
                continue;
            }
            if (components_[i].intersects(components_[j])) {
                parent[findRoot(i)] = findRoot(j);
                --pieces;
            }
        }
    }
    return pieces == 1;
}
 
// ---------------------------------------------------------------------------
// The set of regions, in both directions.
//
// A set is the first shape in the library that may be disconnected, and that is
// what shapes these. As a **target** it can already be in several pieces, so
// none of the shortcuts that lean on a connected target apply and the engine has
// to count for real — including for a remover that misses it entirely, which
// still leaves it disconnected. As a **remover** it behaves like any other
// polygonal shape: its rings are cut segments like a region's.
//
// Neither direction folds over the components. Removing one component may leave
// the target whole while removing all of them cuts it, and a target's pieces are
// counted across the whole set at once.
 
template <class PointType, class LabelType>
template <detail::SetOperandConcept OtherShape>
bool PolygonSet<PointType, LabelType>::separates(const OtherShape& other) const {
    if (empty()) {
        return false;  // nothing removed, and every operand here is connected
    }
 
    if constexpr (PointConcept<OtherShape>) {
        return false;  // removing anything from a point leaves at most one piece
    } else if constexpr (OrientedSegmentConcept<OtherShape>) {
        return separates(other.asSegment());
    } else if constexpr (OrientedLineConcept<OtherShape>) {
        return separates(other.asLine());
    } else if constexpr (LineConcept<OtherShape> || RayConcept<OtherShape>) {
        if (other.isDegenerate()) {
            return false;
        }
        return detail::linearAndRegionSeparate<true>(other, *this);
    } else if constexpr (HalfplaneConcept<OtherShape>) {
        if (other.isUndefined()) {
            return false;
        }
        return detail::convexAndRegionSeparate<true>(other, *this);
    } else if constexpr (HalfplaneIntersectionConcept<OtherShape>) {
        return detail::convexAndRegionSeparate<true>(other, *this);
    } else if constexpr (DiskConcept<OtherShape>) {
        if (other.isDegenerate()) {
            return false;  // a disk of radius zero is a point, an undefined one has no circle
        }
        return detail::regionSeparatesDisk(*this, other);
    } else if constexpr (RectangleConcept<OtherShape> || TriangleConcept<OtherShape> ||
                         ConvexConcept<OtherShape>) {
        return separates(other.asPolygon());
    } else if constexpr (MonotoneChainConcept<OtherShape> || PolylineConcept<OtherShape>) {
        return detail::cellSeparates(other, *this);
    } else {
        if (other.isDegenerate()) {
            return false;
        }
        return detail::cellSeparates(other, *this);
    }
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool PolygonSet<PointType, LabelType>::separates(const OtherSet& other) const {
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PointConcept OtherPoint>
bool PolygonSet<PointType, LabelType>::separates(const Shape<OtherPoint>& other) const {
    return std::visit([this](const auto& value) { return this->separates(value); },
                      other.variant());
}
 
// Reverse direction: removing a lower-ranked shape from a set. Every one of
// these is a real computation — even a point's, since a set that is already in
// two pieces stays in two however little is taken out of it.
 
template <class Number, class Label>
template <PolygonSetConcept OtherSet>
bool Point<Number, Label>::separates(const OtherSet& other) const {
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool Segment<PointType, LabelType>::separates(const OtherSet& other) const {
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool OrientedSegment<PointType, LabelType>::separates(const OtherSet& other) const {
    return asSegment().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool Line<PointType, LabelType>::separates(const OtherSet& other) const {
    if (isDegenerate()) {
        return min().separates(other);
    }
    return detail::linearAndRegionSeparate<false>(*this, other);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool OrientedLine<PointType, LabelType>::separates(const OtherSet& other) const {
    return asLine().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool Ray<PointType, LabelType>::separates(const OtherSet& other) const {
    if (isDegenerate()) {
        return source().separates(other);
    }
    return detail::linearAndRegionSeparate<false>(*this, other);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool Halfplane<PointType, LabelType>::separates(const OtherSet& other) const {
    if (isUndefined()) {
        // A collapsed half-plane is the whole plane or none of it, so it removes
        // everything or nothing; either way only a set already in several pieces
        // comes apart.
        return detail::disconnectedOnItsOwn(other);
    }
    return detail::convexAndRegionSeparate<false>(*this, other);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool Rectangle<PointType, LabelType>::separates(const OtherSet& other) const {
    // No empty short-circuit here: a set of regions may already be in several
    // pieces, and removing the empty set from one of those still leaves it
    // disconnected. The cell engine decides, reached through the empty polygon.
    return asPolygon().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool Triangle<PointType, LabelType>::separates(const OtherSet& other) const {
    return asPolygon().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool Convex<PointType, LabelType>::separates(const OtherSet& other) const {
    return asPolygon().separates(other);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool Polygon<PointType, LabelType>::separates(const OtherSet& other) const {
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool PolygonWithHoles<PointType, LabelType>::separates(const OtherSet& other) const {
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType, class Storage>
template <PolygonSetConcept OtherSet>
bool MonotoneChain<PointType, LabelType, Storage>::separates(const OtherSet& other) const {
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool Polyline<PointType, LabelType>::separates(const OtherSet& other) const {
    return detail::cellSeparates(other, *this);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool Disk<PointType, LabelType>::separates(const OtherSet& other) const {
    if (isDegenerate()) {
        return detail::cellSeparates(other, a());
    }
    return detail::diskSeparatesRegion(*this, other);
}
 
template <class PointType, class LabelType>
template <PolygonSetConcept OtherSet>
bool HalfplaneIntersection<PointType, LabelType>::separates(const OtherSet& other) const {
    if (empty()) {
        return detail::disconnectedOnItsOwn(other);  // nothing is removed
    }
    return detail::convexAndRegionSeparate<false>(*this, other);
}
 
}  // namespace pgl