polygonset.hpp

Go to the documentation of this file.

#pragma once
 
#include "shape/polygonwithholes.hpp"
 
#include <algorithm>
#include <cassert>
#include <compare>
#include <concepts>
#include <cstddef>
#include <iterator>
#include <optional>
#include <ostream>
#include <ranges>
#include <type_traits>
#include <utility>
#include <variant>
#include <vector>
 
 
namespace pgl {
 
namespace detail {

template <class T>
concept SetOperandConcept =
    shapeRank<std::remove_cvref_t<T>> >= 0 &&
    shapeRank<std::remove_cvref_t<T>> < shapeRank<PolygonSet<Point<>, NoLabel>> &&
    !is_empty_shape_v<T>;

template <class T>
concept SetBooleanOperandConcept = PolygonalRegionConcept<T>;

template <class T>
concept SetMinkowskiOperandConcept =
    SetBooleanOperandConcept<T> || is_segment_v<T> || is_oriented_segment_v<T> ||
    is_polyline_v<T> || is_monotone_chain_v<T>;

template <class ResultNumber, class Component, class Operand>
concept ComponentDistanceL1Concept =
    requires(const Component& component, const Operand& operand) {
        component.template distanceL1<ResultNumber>(operand);
    };

template <class ResultNumber, class Component, class Operand>
concept ComponentDistanceLInfConcept =
    requires(const Component& component, const Operand& operand) {
        component.template distanceLInf<ResultNumber>(operand);
    };
 
}  // namespace detail
 
template <class PointType = Point<>, class Label>
struct PolygonSet;
 
// Explicit deduction guides, rather than the implicit ones the constrained
// constructors below would otherwise synthesize: clang 18 mishandles the
// constraints on such a guide, and CI compiles with it.
PolygonSet() -> PolygonSet<Point<>, NoLabel>;
 
template <PolygonWithHolesConcept Component>
PolygonSet(Component&&) -> PolygonSet<typename std::remove_cvref_t<Component>::PointType, NoLabel>;
 
template <std::ranges::input_range ComponentRange>
    requires detail::is_polygon_with_holes_v<std::ranges::range_value_t<ComponentRange>>
PolygonSet(ComponentRange&&)
    -> PolygonSet<typename std::ranges::range_value_t<ComponentRange>::PointType, NoLabel>;
 
template <std::ranges::input_range ComponentRange>
    requires detail::is_polygon_with_holes_v<std::ranges::range_value_t<ComponentRange>>
PolygonSet(ComponentRange&&, bool)
    -> PolygonSet<typename std::ranges::range_value_t<ComponentRange>::PointType, NoLabel>;
 

template <class PointType_, class TLabel>
struct PolygonSet {
    using PointType = PointType_;
    using NumberType = typename PointType::NumberType;
    using LabelType = TLabel;
    using ComponentType = PolygonWithHoles<PointType>;
 
    class VertexIterator;
    using PolygonType = Polygon<PointType>;
    using EdgeType = Segment<PointType>;
    static_assert(detail::is_point_v<PointType>, "PolygonSet requires pgl::Point vertices");

    constexpr PolygonSet() = default;

    constexpr explicit PolygonSet(ComponentType component) {
        if (!component.isDegenerate()) {
            components_.push_back(std::move(component));
        }
    }

    template <std::ranges::input_range ComponentRange>
        requires detail::is_polygon_with_holes_v<std::ranges::range_value_t<ComponentRange>>
    constexpr PolygonSet(ComponentRange&& components, bool trusted = false) {
        for (const auto& component : components) {
            components_.emplace_back(component);
        }
        if (!trusted) {
            normalize();
        }
    }

    template <PointConcept OtherPointType, class OtherLabelType>
        requires(std::constructible_from<PointType, const OtherPointType&>)
    constexpr PolygonSet(const PolygonSet<OtherPointType, OtherLabelType>& other) {
        components_.reserve(other.componentCount());
        for (const auto& component : other.components()) {
            components_.emplace_back(component);
        }
    }

    template <class A = LabelType>
        requires(detail::has_label_v<A>)
    constexpr A& label() const {
        return label_;
    }
 
    // -------------------------------------------------------------------------
    // Component access
    //
    // Deliberately not `size()` / `operator[]`: `size()` counts defining points
    // on Polygon, Convex, Polyline and MonotoneChain, and a name whose meaning
    // differs per shape is a trap in generic code. @ref PolygonWithHoles made
    // the same call for its holes.

    [[nodiscard]] constexpr std::size_t componentCount() const {
        return components_.size();
    }

    [[nodiscard]] constexpr const ComponentType& component(std::size_t index) const {
        assert(index < components_.size());
        return components_[index];
    }

    [[nodiscard]] constexpr const std::vector<ComponentType>& components() const {
        return components_;
    }

    [[nodiscard]] constexpr auto begin() const { return components_.begin(); }

    [[nodiscard]] constexpr auto cbegin() const { return components_.cbegin(); }

    [[nodiscard]] constexpr auto end() const { return components_.end(); }

    [[nodiscard]] constexpr auto cend() const { return components_.cend(); }

    constexpr void addComponent(ComponentType component) {
        if (component.isDegenerate()) {
            return;
        }
        const auto position = std::lower_bound(components_.begin(), components_.end(), component);
        if (position != components_.end() && *position == component) {
            return;
        }
        components_.insert(position, std::move(component));
        resetCache();
    }

    constexpr void eraseComponent(std::size_t index) {
        assert(index < components_.size());
        components_.erase(components_.begin() + static_cast<std::ptrdiff_t>(index));
        resetCache();
    }

    constexpr bool eraseComponent(const ComponentType& component) {
        const auto position = std::lower_bound(components_.begin(), components_.end(), component);
        if (position == components_.end() || !(*position == component)) {
            return false;
        }
        components_.erase(position);
        resetCache();
        return true;
    }

    [[nodiscard]] constexpr std::size_t holeCount() const {
        std::size_t total = 0;
        for (const auto& component : components_) {
            total += component.holeCount();
        }
        return total;
    }

    [[nodiscard]] constexpr bool hasHoles() const {
        for (const auto& component : components_) {
            if (component.hasHoles()) {
                return true;
            }
        }
        return false;
    }

    [[nodiscard]] constexpr std::size_t vertexCount() const {
        std::size_t total = 0;
        for (const auto& component : components_) {
            total += component.vertexCount();
        }
        return total;
    }

    [[nodiscard]] constexpr std::vector<PointType> vertices() const {
        std::vector<PointType> result;
        result.reserve(vertexCount());
        for (const auto& component : components_) {
            for (const auto& vertex : component.vertices()) {
                result.push_back(vertex);
            }
        }
        return result;
    }

    [[nodiscard]] constexpr auto verticesView() const {
        return std::ranges::subrange(verticesBegin(), verticesEnd());
    }

    [[nodiscard]] constexpr VertexIterator verticesBegin() const {
        return VertexIterator(this, 0);
    }

    [[nodiscard]] constexpr VertexIterator verticesEnd() const {
        return VertexIterator(this, components_.size());
    }

    [[nodiscard]] constexpr std::vector<EdgeType> edges() const {
        std::vector<EdgeType> result;
        result.reserve(vertexCount());
        for (const auto& component : components_) {
            for (const auto& edge : component.edges()) {
                result.push_back(edge);
            }
        }
        return result;
    }

    [[nodiscard]] constexpr std::vector<OrientedSegment<PointType>> orientedEdges() const {
        std::vector<OrientedSegment<PointType>> result;
        result.reserve(vertexCount());
        for (const auto& component : components_) {
            for (const auto& edge : component.orientedEdges()) {
                result.push_back(edge);
            }
        }
        return result;
    }
 
    // -------------------------------------------------------------------------
    // Value semantics

    [[nodiscard]] constexpr auto operator<=>(const PolygonSet& other) const {
        if (auto cmp = components_.size() <=> other.components_.size(); cmp != 0) {
            return cmp;
        }
        for (std::size_t i = 0; i < components_.size(); ++i) {
            if (auto cmp = components_[i] <=> other.components_[i]; cmp != 0) {
                return cmp;
            }
        }
        return std::strong_ordering::equal;
    }

    [[nodiscard]] constexpr bool operator==(const PolygonSet& other) const {
        if (components_.size() != other.components_.size()) {
            return false;
        }
        for (std::size_t i = 0; i < components_.size(); ++i) {
            if (!(components_[i] == other.components_[i])) {
                return false;
            }
        }
        return true;
    }

    template<AnyShapeConcept OtherShape>
    [[nodiscard]] constexpr bool samePointSet(const OtherShape& other) const;
 
    // -------------------------------------------------------------------------
    // State queries

    [[nodiscard]] constexpr bool empty() const {
        return components_.empty();
    }

    [[nodiscard]] constexpr bool isDegenerate() const {
        using Exact = detail::promoted_number_t<NumberType>;
        return twiceArea<Exact>() == Exact(0);
    }

    [[nodiscard]] constexpr bool isPoint() const {
        return components_.size() == 1 && components_[0].isPoint();
    }

    [[nodiscard]] constexpr bool isSegment() const {
        return components_.size() == 1 && components_[0].isSegment();
    }

    [[nodiscard]] constexpr bool isUndefined() const {
        return !isPoint() && !isSegment() && isDegenerate();
    }

    template <class Rational = pgl::Rational<pgl::BigInt>>
    [[nodiscard]] bool isSimple() const {
        for (const auto& component : components_) {
            if (!component.template isSimple<Rational>()) {
                return false;
            }
        }
        return true;
    }

    template <class Rational = pgl::Rational<pgl::BigInt>>
    [[nodiscard]] bool isValid() const;

    [[nodiscard]] bool isRegular() const {
        for (const auto& component : components_) {
            if (!component.isRegular()) {
                return false;
            }
        }
        return true;
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType::LabelType>> regularized() const;
 
    // -------------------------------------------------------------------------
    // Measures

    template <class ResultNumber = NumberType>
    [[nodiscard]] constexpr ResultNumber twiceArea() const {
        ResultNumber total{};
        for (const auto& component : components_) {
            total += component.template twiceArea<ResultNumber>();
        }
        return total;
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr auto area() const {
        ResultNumber result = static_cast<ResultNumber>(twiceArea());
        return result / ResultNumber(2);
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr Point<ResultNumber> centroid() const;

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr Point<ResultNumber> verticesCentroid() const;

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] Point<ResultNumber> pointInside() const;

    [[nodiscard]] constexpr Segment<PointType> diameter() const {
        std::vector<PointType> hullPoints;
        hullPoints.reserve(vertexCount());
        for (const auto& component : components_) {
            for (const auto& vertex : component.outer()) {
                hullPoints.push_back(vertex);
            }
        }
        return Convex<PointType>(std::move(hullPoints)).diameter();
    }

    [[nodiscard]] constexpr Convex<PointType> convexHull() const {
        std::vector<PointType> hullPoints;
        hullPoints.reserve(vertexCount());
        for (const auto& component : components_) {
            for (const auto& vertex : component.outer()) {
                hullPoints.push_back(vertex);
            }
        }
        return Convex<PointType>(std::move(hullPoints));
    }

    [[nodiscard]] constexpr const Rectangle<PointType>& bbox() const;

    template <class ResultNumber = grid_number_t<typename PointType_::NumberType>>
        requires(detail::extended_integral<ResultNumber> || std::same_as<ResultNumber, BigInt>)
    [[nodiscard]] std::vector<Point<ResultNumber, typename PointType::LabelType>>
    latticePoints() const;

    template <std::floating_point ResultNumber = double>
    [[nodiscard]] constexpr Rectangle<Point<ResultNumber>> fbox() const;
 
    // -------------------------------------------------------------------------
    // Decompositions

    auto triangulation() const;

    template <class SegmentRange>
    auto triangulation(const SegmentRange& segments) const;

    [[nodiscard]] std::vector<Convex<PointType>> convexPartition() const;

    [[nodiscard]] std::vector<Convex<PointType>> convexCovering() const;

    template <class ResultNumber = grid_number_t<typename PointType_::NumberType>>
        requires(std::signed_integral<ResultNumber>)
    [[nodiscard]] auto asBitMatrix() const;
 
    // -------------------------------------------------------------------------
    // Boolean operations
    //
    // This is what the shape is for. `Polygon` and `PolygonWithHoles` already
    // produce a set of regions from a difference, a union or a symmetric
    // difference, and now say so in their return type; a set does the same and
    // takes one back, so the four operations are **closed** and a result can be
    // fed straight into the next one.
    //
    // The engine does not care that a receiver is a set: it is the arrangement
    // of both operands' boundaries with one witness test per cell, and a set
    // contributes its components' rings the way a region contributes its own.
    // In particular the operands go in together rather than being folded over
    // one component at a time, so one arrangement settles the whole answer.

    template <class ResultNumber = division_result_t<NumberType>, detail::SetBooleanOperandConcept OtherShape>
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType::LabelType>>
    difference(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneIntersectionConcept OtherIntersection>
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType::LabelType>>
    difference(const OtherIntersection& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneConcept OtherHalfplane>
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType::LabelType>>
    difference(const OtherHalfplane& other) const;

    template <class ResultNumber = division_result_t<NumberType>, detail::SetBooleanOperandConcept OtherShape>
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType::LabelType>>
    regularizedUnion(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, detail::SetBooleanOperandConcept OtherShape>
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType::LabelType>>
    regularizedIntersection(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, detail::SetBooleanOperandConcept OtherShape>
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType::LabelType>>
    symmetricDifference(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneIntersectionConcept OtherIntersection>
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType::LabelType>>
    regularizedIntersection(const OtherIntersection& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneConcept OtherHalfplane>
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType::LabelType>>
    regularizedIntersection(const OtherHalfplane& other) const;

    template <class ResultNumber = division_result_t<NumberType>, detail::SetBooleanOperandConcept OtherShape>
    [[nodiscard]] std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                           Polyline<Point<ResultNumber, typename PointType::LabelType>>,
                                           PolygonWithHoles<Point<ResultNumber, typename PointType::LabelType>>>>
    intersection(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneIntersectionConcept OtherIntersection>
    [[nodiscard]] std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                           Polyline<Point<ResultNumber, typename PointType::LabelType>>,
                                           PolygonWithHoles<Point<ResultNumber, typename PointType::LabelType>>>>
    intersection(const OtherIntersection& other) const;

    template <class ResultNumber = division_result_t<NumberType>, HalfplaneConcept OtherHalfplane>
    [[nodiscard]] std::vector<std::variant<Point<ResultNumber, typename PointType::LabelType>,
                                           Polyline<Point<ResultNumber, typename PointType::LabelType>>,
                                           PolygonWithHoles<Point<ResultNumber, typename PointType::LabelType>>>>
    intersection(const OtherHalfplane& other) const;

    template <class ResultNumber = NumberType, class EmptyPoint>
    [[nodiscard]] constexpr EmptyShape<EmptyPoint> intersection(const EmptyShape<EmptyPoint>&) const {
        return {};
    }
 
    // -------------------------------------------------------------------------
    // Predicates
    //
    // With `A = ⋃ Aᵢ`, interiors pairwise disjoint and contacts finite, four of
    // the five relations fold over the components outright:
    //
    //   intersects(x)          ∃i Aᵢ.intersects(x)          — A is the union
    //   interiorsIntersect(x)  ∃i Aᵢ.interiorsIntersect(x)   — A° = ⋃ Aᵢ°
    //   interiorContains(x)    ∃i Aᵢ.interiorContains(x)     — see below
    //   contains(point)        ∃i Aᵢ.contains(point)
    //
    // `A° = ⋃ Aᵢ°` is what the no-shared-edge clause of @ref isValid buys, and
    // it settles `interiorContains` for every operand: the `Aᵢ°` are open and
    // pairwise disjoint, so a connected operand inside their union is inside one
    // of them, and every shape but another set is connected.
    //
    // `contains` and `boundaryContains` are the two that do not fold, and they
    // fail for exactly one kind of operand. `∃i` is still necessary and
    // sufficient whenever the operand stays connected after finitely many points
    // are removed — every operand with area, and every point — because the
    // components meet at finitely many points at most. A **one-dimensional**
    // operand does not: two unit squares meeting corner to corner at the origin
    // contain the segment from (-1,-1) to (1,1) between them and neither
    // contains it alone. Those operands are settled by splitting them at every
    // component-boundary contact and classifying each piece, which is exact and
    // costs nothing when the components do not touch at all (@ref isPinched).

    template <detail::SetOperandConcept OtherShape>
    [[nodiscard]] bool intersects(const OtherShape& other) const;

    template <detail::SetOperandConcept OtherShape>
    [[nodiscard]] bool interiorsIntersect(const OtherShape& other) const;

    template <detail::SetOperandConcept OtherShape>
    [[nodiscard]] bool contains(const OtherShape& other) const;

    template <detail::SetOperandConcept OtherShape>
    [[nodiscard]] bool interiorContains(const OtherShape& other) const;

    template <SegmentConcept OtherSegment>
    [[nodiscard]] bool interiorContainsInterior(const OtherSegment& other) const;

    template <detail::SetOperandConcept OtherShape>
    [[nodiscard]] bool boundaryContains(const OtherShape& other) const;

    template <detail::SetOperandConcept OtherShape>
    [[nodiscard]] bool separates(const OtherShape& other) const;

    template <detail::SetOperandConcept OtherShape>
    [[nodiscard]] bool crosses(const OtherShape& other) const;
 
    // -------------------------------------------------------------------------
    // The self pair. A set operand is the one operand that need not be
    // connected, so the relations that lean on the operand's connectedness fold
    // over *its* components instead of assuming it lies in one piece.

    template <PolygonSetConcept OtherSet>
    [[nodiscard]] bool intersects(const OtherSet& other) const;

    template <PolygonSetConcept OtherSet>
    [[nodiscard]] bool interiorsIntersect(const OtherSet& other) const;

    template <PolygonSetConcept OtherSet>
    [[nodiscard]] bool contains(const OtherSet& other) const;

    template <PolygonSetConcept OtherSet>
    [[nodiscard]] bool interiorContains(const OtherSet& other) const;

    template <PolygonSetConcept OtherSet>
    [[nodiscard]] bool boundaryContains(const OtherSet& other) const;

    template <PolygonSetConcept OtherSet>
    [[nodiscard]] bool separates(const OtherSet& other) const;

    template <PolygonSetConcept OtherSet>
    [[nodiscard]] bool crosses(const OtherSet& other) const;
 
    // -------------------------------------------------------------------------
    // The empty set is a subset of every shape, so its containment relations are
    // true and its intersection relations are false.

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool contains(const EmptyShape<EmptyPoint>&) const {
        return true;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool interiorContains(const EmptyShape<EmptyPoint>&) const {
        return true;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool boundaryContains(const EmptyShape<EmptyPoint>&) const {
        return true;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool intersects(const EmptyShape<EmptyPoint>&) const {
        return false;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool interiorsIntersect(const EmptyShape<EmptyPoint>&) const {
        return false;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool separates(const EmptyShape<EmptyPoint>&) const {
        return false;
    }

    template <class EmptyPoint>
    [[nodiscard]] constexpr bool crosses(const EmptyShape<EmptyPoint>&) const {
        return false;
    }
 
    // -------------------------------------------------------------------------
    // Runtime Shape argument: visit the wrapped alternative and re-dispatch to
    // the matching per-shape overload (defined in the implementation layer).

    template <PointConcept OtherPoint>
    [[nodiscard]] bool contains(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] bool interiorContains(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] bool boundaryContains(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] bool intersects(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] bool interiorsIntersect(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] bool separates(const Shape<OtherPoint>& other) const;

    template <PointConcept OtherPoint>
    [[nodiscard]] bool crosses(const Shape<OtherPoint>& other) const;
 
    // -------------------------------------------------------------------------
    // Distances
    //
    // The distance to a union is the minimum of the distances, with no caveat at
    // all: every one of these is exactly the minimum over the components.

    template <class ResultNumber = division_result_t<NumberType>, detail::SetOperandConcept OtherShape>
    [[nodiscard]] auto squaredDistance(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PolygonSetConcept OtherSet>
    [[nodiscard]] auto squaredDistance(const OtherSet& other) const;

    template <class ResultNumber = NumberType, BoundedPolygonalConcept OtherShape>
        requires detail::ClosestPairConcept<PolygonSet<PointType_, TLabel>, OtherShape>
    [[nodiscard]] auto closestSegments(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, class OtherShape>
        requires detail::ClosestPointsPairConcept<PolygonSet<PointType_, TLabel>, OtherShape>
    [[nodiscard]] auto closestPoints(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, detail::SetOperandConcept OtherShape>
        requires detail::ComponentDistanceL1Concept<ResultNumber, PolygonWithHoles<PointType_>, OtherShape>
    [[nodiscard]] auto distanceL1(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PolygonSetConcept OtherSet>
    [[nodiscard]] auto distanceL1(const OtherSet& other) const;

    template <class ResultNumber = division_result_t<NumberType>, detail::SetOperandConcept OtherShape>
        requires detail::ComponentDistanceLInfConcept<ResultNumber, PolygonWithHoles<PointType_>, OtherShape>
    [[nodiscard]] auto distanceLInf(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PolygonSetConcept OtherSet>
    [[nodiscard]] auto distanceLInf(const OtherSet& other) const;

    template <class ResultNumber = division_result_t<NumberType>, PointConcept OtherPoint>
    [[nodiscard]] constexpr auto intersection(const Shape<OtherPoint>& other) const {
        return other.template intersection<ResultNumber>(*this);
    }

    template <class ResultNumber = division_result_t<NumberType>, PointConcept OtherPoint>
    [[nodiscard]] auto regularizedIntersection(const Shape<OtherPoint>& other) const {
        return other.template regularizedIntersection<ResultNumber>(*this);
    }

    template <class ResultNumber = division_result_t<NumberType>, PointConcept OtherPoint>
    [[nodiscard]] auto regularizedUnion(const Shape<OtherPoint>& other) const {
        return other.template regularizedUnion<ResultNumber>(*this);
    }

    template <class ResultNumber = division_result_t<NumberType>, PointConcept OtherPoint>
    [[nodiscard]] auto difference(const Shape<OtherPoint>& other) const {
        return Shape<OtherPoint>(*this).template difference<ResultNumber>(other);
    }

    template <class ResultNumber = division_result_t<NumberType>, PointConcept OtherPoint>
    [[nodiscard]] auto symmetricDifference(const Shape<OtherPoint>& other) const {
        return other.template symmetricDifference<ResultNumber>(*this);
    }

    template <class ResultNumber = double, PointConcept OtherPoint>
    [[nodiscard]] constexpr auto distanceL1(const Shape<OtherPoint>& other) const {
        return other.template distanceL1<ResultNumber>(*this);
    }

    template <class ResultNumber = double, PointConcept OtherPoint>
    [[nodiscard]] constexpr auto distanceLInf(const Shape<OtherPoint>& other) const {
        return other.template distanceLInf<ResultNumber>(*this);
    }

    [[nodiscard]] bool isPinched() const;

    [[nodiscard]] bool isConnected() const;

    template <class OtherShape>
        requires MinkowskiSummableConcept<PolygonSet<PointType_, TLabel>, OtherShape>
    [[nodiscard]] constexpr auto minkowskiSum(const OtherShape& other) const;

    template <class OtherShape>
        requires MinkowskiSummableConcept<PolygonSet<PointType_, TLabel>, OtherShape>
    [[nodiscard]] constexpr auto minkowskiErosion(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>, class OtherShape>
        requires (!MinkowskiSummableConcept<PolygonSet<PointType_, TLabel>, OtherShape> &&
                  BoundedPolygonalConcept<OtherShape>)
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType_::LabelType>>
    minkowskiErosion(const OtherShape& other) const;

    template <class ResultNumber = division_result_t<NumberType>,
              detail::SetMinkowskiOperandConcept OtherShape>
    [[nodiscard]] PolygonSet<Point<ResultNumber, typename PointType::LabelType>>
    minkowskiSum(const OtherShape& other) const;
 
    // -------------------------------------------------------------------------
    // Transformations

    template <class TranslationNumber, class TranslationLabel>
    constexpr PolygonSet& operator+=(const Point<TranslationNumber, TranslationLabel>& translation) {
        for (auto& component : components_) {
            component += translation;
        }
        resetCache();
        return *this;
    }

    template <class TranslationNumber, class TranslationLabel>
    constexpr PolygonSet& operator-=(const Point<TranslationNumber, TranslationLabel>& translation) {
        return *this += (-translation);
    }

    template <class Scalar>
        requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
    constexpr PolygonSet& operator*=(const Scalar& scalar) {
        for (auto& component : components_) {
            component *= scalar;
        }
        normalize();
        return *this;
    }

    template <class Scalar>
        requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
    constexpr PolygonSet& operator/=(const Scalar& scalar) {
        for (auto& component : components_) {
            component /= scalar;
        }
        normalize();
        return *this;
    }

    [[nodiscard]] constexpr PolygonSet rotated90(int k) const {
        return mappedComponents([k](const ComponentType& component) { return component.rotated90(k); });
    }

    constexpr void rotate90(int k) {
        auto saved = label_;
        *this = rotated90(k);
        label_ = std::move(saved);
    }

    template <class OtherNumber>
    [[nodiscard]] constexpr PolygonSet scaledUpX(const OtherNumber scalar) const {
        return mappedComponents(
            [scalar](const ComponentType& component) { return component.scaledUpX(scalar); });
    }

    template <class OtherNumber>
    constexpr void scaleUpX(const OtherNumber scalar) {
        auto saved = label_;
        *this = scaledUpX(scalar);
        label_ = std::move(saved);
    }

    template <class OtherNumber>
    [[nodiscard]] constexpr PolygonSet scaledUpY(const OtherNumber scalar) const {
        return mappedComponents(
            [scalar](const ComponentType& component) { return component.scaledUpY(scalar); });
    }

    template <class OtherNumber>
    constexpr void scaleUpY(const OtherNumber scalar) {
        auto saved = label_;
        *this = scaledUpY(scalar);
        label_ = std::move(saved);
    }

    template <class OtherNumber>
    [[nodiscard]] constexpr PolygonSet scaledDownX(const OtherNumber scalar) const {
        return mappedComponents(
            [scalar](const ComponentType& component) { return component.scaledDownX(scalar); });
    }

    template <class OtherNumber>
    constexpr void scaleDownX(const OtherNumber scalar) {
        auto saved = label_;
        *this = scaledDownX(scalar);
        label_ = std::move(saved);
    }

    template <class OtherNumber>
    [[nodiscard]] constexpr PolygonSet scaledDownY(const OtherNumber scalar) const {
        return mappedComponents(
            [scalar](const ComponentType& component) { return component.scaledDownY(scalar); });
    }

    template <class OtherNumber>
    constexpr void scaleDownY(const OtherNumber scalar) {
        auto saved = label_;
        *this = scaledDownY(scalar);
        label_ = std::move(saved);
    }

    class VertexIterator {
      public:
        using iterator_category = std::forward_iterator_tag;
        using iterator_concept = std::forward_iterator_tag;
        using value_type = PointType;
        using difference_type = std::ptrdiff_t;
        using reference = value_type;
 
        constexpr VertexIterator() = default;
 
        constexpr value_type operator*() const {
            assert(set != nullptr);
            return *inner;
        }
 
        constexpr VertexIterator& operator++() {
            ++inner;
            skipExhausted();
            return *this;
        }
 
        constexpr VertexIterator operator++(int) {
            VertexIterator copy(*this);
            ++(*this);
            return copy;
        }
 
        constexpr bool operator==(const VertexIterator& other) const = default;
 
      private:
        friend struct PolygonSet;
 
        constexpr VertexIterator(const PolygonSet* set_arg, std::size_t component_arg)
            : set(set_arg), component(component_arg) {
            enterComponent();
            skipExhausted();
        }
 
        // Seeds `inner` from the current component, or clears it at the end so
        // that a walked-to-the-end iterator compares equal to verticesEnd().
        constexpr void enterComponent() {
            const std::size_t count = set == nullptr ? 0 : set->components_.size();
            inner = component < count ? set->components_[component].verticesBegin()
                                      : typename ComponentType::VertexIterator{};
        }
 
        // Steps over components holding no vertices. A canonical set has none —
        // addComponent() drops the ones without area — but this keeps the end
        // state reachable by increment alone for a set built any other way.
        constexpr void skipExhausted() {
            const std::size_t count = set == nullptr ? 0 : set->components_.size();
            while (component < count && inner == set->components_[component].verticesEnd()) {
                ++component;
                enterComponent();
            }
        }
 
        const PolygonSet* set = nullptr;
        std::size_t component = 0;
        typename ComponentType::VertexIterator inner{};
    };
 
  private:
    template <class ComponentRelation>
    constexpr bool anyComponent(ComponentRelation&& relation) const {
        for (const auto& component : components_) {
            if (relation(component)) {
                return true;
            }
        }
        return false;
    }

    template <class OtherSegment>
    bool segmentIn(const OtherSegment& segment, bool boundaryOnly) const;

    template <class OtherChain>
    bool chainIn(const OtherChain& chain, bool boundaryOnly) const;

    template <class OtherRegion>
    bool regionIn(const OtherRegion& region) const;

    template <class ComponentDistance>
    auto minOverComponents(ComponentDistance&& distance) const {
        using ResultNumber = std::decay_t<decltype(distance(std::declval<const ComponentType&>()))>;
        ResultNumber best{};
        bool seeded = false;
        for (const auto& component : components_) {
            const ResultNumber current = distance(component);
            if (!seeded || current < best) {
                best = current;
                seeded = true;
            }
        }
        return best;
    }

    template <class ComponentTransform>
    constexpr PolygonSet mappedComponents(ComponentTransform&& transform) const {
        PolygonSet result;
        result.components_.reserve(components_.size());
        for (const auto& component : components_) {
            result.components_.push_back(transform(component));
        }
        result.normalize();
        return result;
    }
 
    std::vector<ComponentType> components_{};
    [[no_unique_address]] mutable LabelType label_{};
 
    // Cached bounding box. A region reads its own off its outer ring, which
    // caches one already; a set has no single ring to ask, so it caches the
    // union of the components' boxes itself.
    mutable Rectangle<PointType> bbox_{};
 
    // Tri-state cache for @ref isPinched: -1 not yet computed, 0 no two
    // components touch, 1 some two do.
    mutable signed char pinched_ = -1;
 
    // Memoized hash, computed lazily by std::hash<PolygonSet>, with the same
    // sentinel scheme as Polygon and PolygonWithHoles: hashUnset_ means "not yet
    // computed", and the one true hash colliding with it is remapped so the
    // sentinel is never stored as a real value.
    static constexpr std::size_t hashUnset_ = pgl::detail::numeric_limits<std::size_t>::max();
    mutable std::size_t hash_ = hashUnset_;
    friend struct std::hash<PolygonSet>;
 
    template <class OtherPointType, class OtherLabelType>
    friend struct PolygonSet;
 
    constexpr void resetCache() const {
        hash_ = hashUnset_;
        bbox_ = {};
        pinched_ = -1;
    }

    constexpr void normalize() {
        std::erase_if(components_,
                      [](const ComponentType& component) { return component.isDegenerate(); });
        std::sort(components_.begin(), components_.end());
        components_.erase(std::unique(components_.begin(), components_.end()), components_.end());
        resetCache();
    }
};
 
// `set + point` is the translating Minkowski sum, spelled by the generic
// operator+ in implementation/minkowski.hpp like every other shape's: writing a
// second one here would shadow it and, holding the set's own point type, would
// silently truncate a translation that does not fit it.

template <class PointType, class LabelType, class TranslationNumber, class TranslationLabel>
constexpr auto operator-(const PolygonSet<PointType, LabelType>& set,
                         const Point<TranslationNumber, TranslationLabel>& translation) {
    return set + (-translation);
}
 
template <class PointType, class LabelType, class Scalar>
    requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
constexpr auto operator*(const PolygonSet<PointType, LabelType>& set, const Scalar& scalar) {
    using ResultPointType = Point<decltype(std::declval<PointType>().x() * scalar), typename PointType::LabelType>;
    PolygonSet<ResultPointType, LabelType> result(set);
    result *= scalar;
    if constexpr (detail::has_label_v<LabelType>) {
        result.label() = LabelType{};
    }
    return result;
}
 
template <class Scalar, class PointType, class LabelType>
    requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
constexpr auto operator*(const Scalar& scalar, const PolygonSet<PointType, LabelType>& set) {
    return set * scalar;
}
 
template <class PointType, class LabelType, class Scalar>
    requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
constexpr auto operator/(const PolygonSet<PointType, LabelType>& set, const Scalar& scalar) {
    using ResultPointType = Point<decltype(std::declval<PointType>().x() / scalar), typename PointType::LabelType>;
    PolygonSet<ResultPointType, LabelType> result(set);
    result /= scalar;
    if constexpr (detail::has_label_v<LabelType>) {
        result.label() = LabelType{};
    }
    return result;
}
 
template <class PointType, class LabelType>
std::ostream& operator<<(std::ostream& stream, const PolygonSet<PointType, LabelType>& set);
 
}  // namespace pgl