disk.hpp

Go to the documentation of this file.

#pragma once
 
#include "shape/triangle.hpp"

 
#include <array>
#include <cassert>
#include <cmath>
#include <compare>
#include <concepts>
#include <cstddef>
#include <iterator>
#include <numbers>
#include <optional>
#include <vector>
#include <ostream>
#include <type_traits>
#include <utility>
#include <stdexcept>
 
 
namespace pgl {
 
template <class PointType = Point<>, class Label = NoLabel>
struct Disk;

Disk() -> Disk<Point<>, NoLabel>;

template <class PointType>
Disk(PointType, typename PointType::NumberType) -> Disk<PointType, NoLabel>;

template <PointConcept PointType>
Disk(PointType, PointType, PointType) -> Disk<PointType, NoLabel>;

template <class Number>
Disk(Number, Number, Number) -> Disk<Point<Number>, NoLabel>;

template <class PointType_, class TLabel>
struct Disk {
    using PointType = PointType_;
    using NumberType = PointType::NumberType;
    using PointLabelType = PointType::LabelType;
    using LabelType = TLabel;
 
    static_assert(detail::is_point_v<PointType>, "Disk requires pgl::Point boundary points");

    constexpr Disk() = default;

    constexpr Disk(PointType first, PointType second, PointType third)
        requires PointConcept<PointType>
        : points_(canonicalizePoints(std::move(first), std::move(second), std::move(third))) {}

    template <class A>
        requires(detail::has_label_v<LabelType> && std::constructible_from<LabelType, A&&>)
    constexpr Disk(PointType first, PointType second, PointType third, A&& label)
        requires PointConcept<PointType>
        : Disk(std::move(first), std::move(second), std::move(third)) {
        label_ = std::forward<A>(label);
    }

    constexpr Disk(PointType center, NumberType radius)
        requires PointConcept<PointType> {
        const NumberType zero{};
        const NumberType r = radius < zero ? -radius : radius;
        const NumberType cx = center.x();
        const NumberType cy = center.y();
        const auto point_label = detail::copyLabel<PointLabelType>(center);
 
        points_ = canonicalizePoints(
            Point<NumberType, PointLabelType>(cx - r, cy,     point_label),
            Point<NumberType, PointLabelType>(cx + r, cy,     point_label),
            Point<NumberType, PointLabelType>(cx,     cy + r, point_label)
        );
    }

    template <class A>
        requires(detail::has_label_v<LabelType> && std::constructible_from<LabelType, A&&>)
    constexpr Disk(PointType center, NumberType radius, A&& label)
        requires PointConcept<PointType> : Disk(std::move(center), radius) {
        label_ = std::forward<A>(label);
    }

    constexpr Disk(NumberType x, NumberType y, NumberType radius)
        : Disk(PointType(x, y), radius) {}

    template <class A>
        requires(detail::has_label_v<LabelType> && std::constructible_from<LabelType, A&&>)
    constexpr Disk(NumberType x, NumberType y, NumberType radius, A&& label)
        : Disk(PointType(x, y), radius, std::forward<A>(label)) {}

    template <PointConcept OtherPointType, class OtherLabelType>
        requires(std::constructible_from<PointType, const OtherPointType&>)
    constexpr Disk(const Disk<OtherPointType, OtherLabelType>& other)
        : points_(canonicalizePoints(PointType(other.a()), PointType(other.b()), PointType(other.c()))),
          label_(detail::copyLabel<LabelType>(other)) {}

    template <PointConcept OtherPointType, class OtherLabelType>
        requires(std::constructible_from<PointType, const OtherPointType&>)
    constexpr Disk& operator=(const Disk<OtherPointType, OtherLabelType>& other) {
        points_ = canonicalizePoints(PointType(other.a()), PointType(other.b()), PointType(other.c()));
        label_ = detail::copyLabel<LabelType>(other);
        return *this;
    }

    template <class A = LabelType>
        requires(detail::has_label_v<A>)
    constexpr A& label() const {
        return label_;
    }

    constexpr const PointType& operator[](std::size_t index) const {
        assert(index < size());
        return points_[index];
    }

    static constexpr std::size_t size() {
        return 3;
    }

    constexpr const PointType& get(std::ptrdiff_t index) const {
        const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(size());
        return (*this)[static_cast<std::size_t>(((index % n) + n) % n)];
    }

    constexpr std::ptrdiff_t index(const PointType& point) const {
        for (std::ptrdiff_t i = 0; i < static_cast<std::ptrdiff_t>(size()); ++i) {
            if ((*this)[static_cast<std::size_t>(i)] == point) {
                return i;
            }
        }
        return -1;
    }

    constexpr const PointType& a() const {
        return points_[0];
    }

    constexpr const PointType& b() const {
        return points_[1];
    }

    constexpr const PointType& c() const {
        return points_[2];
    }

    constexpr auto begin() const {
        return PointIterator(this, 0);
    }

    constexpr auto cbegin() const {
        return PointIterator(this, 0);
    }

    constexpr auto end() const {
        return PointIterator(this, size());
    }

    constexpr auto cend() const {
        return PointIterator(this, size());
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr Point<ResultNumber, PointLabelType> center() const {
        if (auto cr = centerAndRadius()) {
            return Point<ResultNumber, PointLabelType>(cr->first);
        }
 
        if (isDegenerate()) {
            return Point<ResultNumber, PointLabelType>(a().x(),a().y(), PointLabelType{});
        }
 
        using Coordinate = detail::promoted_number_t<std::common_type_t<NumberType, ResultNumber>>;
        const Coordinate ax = static_cast<Coordinate>(a().x()), ay = static_cast<Coordinate>(a().y());
        const Coordinate bx = static_cast<Coordinate>(b().x()), by = static_cast<Coordinate>(b().y());
        const Coordinate cx = static_cast<Coordinate>(c().x()), cy = static_cast<Coordinate>(c().y());
 
        const Coordinate aa = ax * ax + ay * ay;
        const Coordinate bb = bx * bx + by * by;
        const Coordinate cc = cx * cx + cy * cy;
        const Coordinate two = static_cast<Coordinate>(2);
 
        const Coordinate denominator = two * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by));
 
        return Point<ResultNumber, PointLabelType>(
            (aa * (by - cy) + bb * (cy - ay) + cc * (ay - by)) / denominator,
            (aa * (cx - bx) + bb * (ax - cx) + cc * (bx - ax)) / denominator,
            PointLabelType{}
        );
    }
 
    // /**
    //  * @brief Returns the center in this disk's exact coordinate type.
    //  *
    //  * For integral coordinates the circumcenter is generally rational, so the
    //  * result is `Rational<promoted NumberType>`; for floating-point coordinates
    //  * it is the promoted floating-point type.
    //  */
    // [[nodiscard]] constexpr auto center() const {
    //     using ResultNumber = std::conditional_t<
    //         detail::extended_integral<NumberType>,
    //         Rational<detail::promoted_number_t<NumberType>>,
    //         detail::promoted_number_t<NumberType>>;
    //     return center<ResultNumber>();
    // }

    template <class ResultNumber = double>
    [[nodiscard]] constexpr ResultNumber radius() const {
        if (auto cr = centerAndRadius()) {
            return static_cast<ResultNumber>(cr->second);
        }
 
        if constexpr (!requires(ResultNumber v) { std::sqrt(v); }) {
            throw std::runtime_error("std::sqrt is not available for the requested ResultNumber type");
        } else {
            return std::sqrt(squaredRadius<ResultNumber>());
        }
    }

    constexpr bool isDegenerate() const {
        return orientationSign(a(), b(), c()) == std::partial_ordering::equivalent;
    }

    [[nodiscard]] constexpr bool isPoint() const {
        return a() == b() && b() == c();
    }

    [[nodiscard]] constexpr std::optional<PointType> getIfPoint() const {
        if (!isPoint()) {
            return std::nullopt;
        }
        return a();
    }

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

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr ResultNumber squaredRadius() const {
        if (auto cr = centerAndRadius()) {
            return detail::asNumber<ResultNumber>(cr->second) * detail::asNumber<ResultNumber>(cr->second);
        }
 
        if (isDegenerate()) {
            return ResultNumber{};
        }
 
        // Circumradius formula:
        // R^2 = |AB|^2 |BC|^2 |CA|^2 / (4 * det(A,B,C)^2),
        // where det(A,B,C) is twice the signed area of triangle ABC
        const ResultNumber ab2 = a().template squaredDistance<ResultNumber>(b());
        const ResultNumber bc2 = b().template squaredDistance<ResultNumber>(c());
        const ResultNumber ca2 = c().template squaredDistance<ResultNumber>(a());
        const ResultNumber determinant = static_cast<ResultNumber>(orientationDeterminant(a(), b(), c()));
        const ResultNumber four = static_cast<ResultNumber>(4);
 
        return (ab2 * bc2 * ca2) / (four * determinant * determinant);
    }
 
    // /**
    //  * @brief Returns the squared radius in this disk's exact result type.
    //  *
    //  * Exact for integral coordinates, where the result is
    //  * `Rational<promoted NumberType>`; for floating-point coordinates it is the
    //  * promoted floating-point type.
    //  */
    // [[nodiscard]] constexpr auto squaredRadius() const {
    //     using ResultNumber = std::conditional_t<
    //         detail::extended_integral<NumberType>,
    //         Rational<detail::promoted_number_t<NumberType>>,
    //         detail::promoted_number_t<NumberType>>;
    //     return squaredRadius<ResultNumber>();
    // }

    template <std::floating_point ResultNumber = double>
    [[nodiscard]] constexpr ResultNumber area() const {
        return std::numbers::pi_v<ResultNumber> * squaredRadius<ResultNumber>();
    }

    [[nodiscard]] constexpr Rectangle<PointType> bbox() const {
        // Exact when built from a center and radius: center +/- (radius, radius).
        if (auto cr = centerAndRadius()) {
            const auto radius = cr->second;
            return Rectangle<PointType>(cr->first - PointType(radius, radius),
                                        cr->first + PointType(radius, radius));
        }
 
        // Degenerate disk (collinear boundary points): no finite circle, so the
        // bounding box of the three boundary points contains it.
        if (isDegenerate()) {
            const auto lo = [](NumberType u, NumberType v, NumberType w) { return u < v ? (u < w ? u : w) : (v < w ? v : w); };
            const auto hi = [](NumberType u, NumberType v, NumberType w) { return u > v ? (u > w ? u : w) : (v > w ? v : w); };
            return Rectangle<PointType>(
                PointType(lo(a().x(), b().x(), c().x()), lo(a().y(), b().y(), c().y())),
                PointType(hi(a().x(), b().x(), c().x()), hi(a().y(), b().y(), c().y())));
        }
 
        // General case: the circumcenter is the rational point (nx/d, ny/d) and
        // the radius is sqrt(s)/d, where d = 2*det(A,B,C) > 0 and
        // s = |AB|^2 |BC|^2 |CA|^2 = (radius*d)^2. The box spans
        // [cx +/- r] x [cy +/- r]; how it is rounded depends on the coordinate
        // type, but the result always contains the disk (and may be larger).
        if constexpr (std::floating_point<NumberType>) {
            // Real square root is available: center +/- radius, tight up to
            // floating-point rounding.
            const NumberType r = radius<NumberType>();
            const auto ctr = center<NumberType>();
            return Rectangle<PointType>(ctr - PointType(r, r), ctr + PointType(r, r));
        } else {
            // The circumcentre numerators are degree three in the coordinates,
            // and the radius bound below multiplies three side lengths, so a
            // single promotion overflows well inside the coordinate type's own
            // range — an `int` disk spanning a few million already lost its box.
            // Promote twice for an integral coordinate, as @ref inCircleDeterminant
            // does, which puts every `int` input safely inside the intermediate.
            // A rational coordinate promotes once: its storage grows as needed.
            using Wide = std::conditional_t<
                detail::extended_integral<NumberType>,
                detail::promoted_number_t<detail::promoted_number_t<NumberType>>,
                detail::promoted_number_t<NumberType>>;
            const Wide ax = a().x(), ay = a().y();
            const Wide bx = b().x(), by = b().y();
            const Wide cx = c().x(), cy = c().y();
            const Wide aa = ax * ax + ay * ay;
            const Wide bb = bx * bx + by * by;
            const Wide cc = cx * cx + cy * cy;
 
            Wide d  = Wide{2} * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by));
            Wide nx = aa * (by - cy) + bb * (cy - ay) + cc * (ay - by);
            Wide ny = aa * (cx - bx) + bb * (ax - cx) + cc * (bx - ax);
            if (d < Wide{0}) {           // normalize to d > 0
                d = -d; nx = -nx; ny = -ny;
            }
 
            const Wide ab2 = (ax - bx) * (ax - bx) + (ay - by) * (ay - by);
            const Wide bc2 = (bx - cx) * (bx - cx) + (by - cy) * (by - cy);
            const Wide ca2 = (cx - ax) * (cx - ax) + (cy - ay) * (cy - ay);
            // ceil(sqrt(n)) for a non-negative *integer* n: floor-sqrt (std::sqrt
            // when usable, else integer Newton) followed by an upward correction.
            // The argument is always an integer, so this terminates.
            const auto ceilIntSqrt = [](auto n) {
                using Int = decltype(n);
                Int r{};
                if constexpr (requires { std::sqrt(n); }) {
                    r = static_cast<Int>(std::sqrt(n));
                } else if (n >= Int{2}) {
                    Int x = n, y = (n + Int{1}) / Int{2};
                    while (y < x) { x = y; y = (x + n / x) / Int{2}; }
                    r = x;
                } else {
                    r = n;
                }
                while (r * r < n) { ++r; }
                return r;
            };
 
            if constexpr (pgl::is_Rational_v<Wide>) {
                // sqrt(s) = sqrt(sn*sd)/sd <= ceil(sqrt(sn*sd))/sd =: u, an exact
                // rational upper bound whose only square root is over the integer
                // sn*sd -- so it never iterates on the (irrational) sqrt itself.
                // Reduced once up front: both parts are wanted, and it also keeps
                // the integer square root below off a needlessly wide sn*sd.
                const Wide s = (ab2 * bc2 * ca2).simplified();   // (radius * d)^2
                using Int = std::remove_cvref_t<decltype(s.numerator())>;
                const Int sn = s.numerator();
                const Int sd = s.denominator();
                const Wide u(ceilIntSqrt(sn * sd), sd);   // u >= sqrt(s)
                return Rectangle<PointType>(
                    PointType(static_cast<NumberType>((nx - u) / d),
                              static_cast<NumberType>((ny - u) / d)),
                    PointType(static_cast<NumberType>((nx + u) / d),
                              static_cast<NumberType>((ny + u) / d)));
            } else {
                // Integer coordinates: round (nx +/- q)/d outward to integers,
                // with q >= radius*d = sqrt(ab2*bc2*ca2).
                //
                // Bound each factor's root separately rather than forming that
                // product: it is degree six in the coordinates, so it overflows
                // a fixed-width Wide at coordinates around a thousand — where
                // the box quietly stopped containing its disk. Each ceil-sqrt is
                // at least the real root, so their product is at least the root
                // of the product, which is all the radius needs; every
                // intermediate stays as wide as a single squared distance. The
                // box remains the outer bound promised above, wider by under one
                // unit per factor.
                const Wide q = ceilIntSqrt(ab2) * ceilIntSqrt(bc2) * ceilIntSqrt(ca2);
                const auto floorDiv = [](Wide p, Wide den) {   // den > 0
                    return p >= Wide{0} ? p / den : -((-p + den - Wide{1}) / den);
                };
                const auto ceilDiv = [](Wide p, Wide den) {    // den > 0
                    return p >= Wide{0} ? (p + den - Wide{1}) / den : -((-p) / den);
                };
                return Rectangle<PointType>(
                    PointType(static_cast<NumberType>(floorDiv(nx - q, d)),
                              static_cast<NumberType>(floorDiv(ny - q, d))),
                    PointType(static_cast<NumberType>(ceilDiv(nx + q, d)),
                              static_cast<NumberType>(ceilDiv(ny + q, d))));
            }
        }
    }

    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 {
        const ResultNumber r = radius<ResultNumber>();
        const auto center_point = center<ResultNumber>();
        const ResultNumber cx = center_point.x();
        const ResultNumber cy = center_point.y();
        return Rectangle<Point<ResultNumber>>(
            Point<ResultNumber>(cx - r, cy - r),
            Point<ResultNumber>(cx + r, cy + r));
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr Point<ResultNumber, PointLabelType> pointInside() const {
        const ResultNumber two = static_cast<ResultNumber>(2);
        return Point<ResultNumber, PointLabelType>(
            (detail::asNumber<ResultNumber>(a().x()) + detail::asNumber<ResultNumber>(b().x())) / two,
            (detail::asNumber<ResultNumber>(a().y()) + detail::asNumber<ResultNumber>(b().y())) / two,
            PointLabelType{});
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] constexpr Segment<Point<ResultNumber, PointLabelType>> diameter() const {
        using ResultPoint = Point<ResultNumber, PointLabelType>;
 
        if (auto cr = centerAndRadius()) {
            return Segment<ResultPoint>(a(),b());
        }
 
        if (isDegenerate()) {
            return Segment<ResultPoint>();
        }
 
        const ResultPoint center_point = center<ResultNumber>();
        ResultPoint anti_a = center_point + (center_point-a());
        return Segment<ResultPoint>(static_cast<ResultPoint>(a()), anti_a);
    }

    template <class ResultNumber = double, PointConcept OtherPoint>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherPoint& point) const;

    template <class ResultNumber = double, SegmentConcept OtherSegment>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherSegment& other) const;

    template <class ResultNumber = double, OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherOrientedSegment& other) const;

    template <class ResultNumber = double, LineConcept OtherLine>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherLine& other) const;

    template <class ResultNumber = double, OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherOrientedLine& other) const;

    template <class ResultNumber = double, RayConcept OtherRay>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherRay& other) const;

    template <class ResultNumber = double, HalfplaneConcept OtherHalfplane>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherHalfplane& other) const;

    template <class ResultNumber = double, RectangleConcept OtherRectangle>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherRectangle& other) const;

    template <class ResultNumber = double, TriangleConcept OtherTriangle>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherTriangle& other) const;

    template <class ResultNumber = double, DiskConcept OtherDisk>
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherDisk& other) const;

    template <class ResultNumber = double, typename OtherShape>
        requires (!PointConcept<OtherShape> && detail::shapeRank<OtherShape> > detail::shapeRank<Disk>
                  && requires(const OtherShape& o, const Disk& self) { o.squaredDistance(self); })
    [[nodiscard]] detail::floating_result_t<ResultNumber> squaredDistance(const OtherShape& other) const {
        return other.template squaredDistance<ResultNumber>(*this);
    }

    template <class ResultNumber = double, PointConcept OtherPoint>
    [[nodiscard]] detail::floating_result_t<ResultNumber> distanceL1(const OtherPoint& point) const;

    template <class ResultNumber = double, PointConcept OtherPoint>
    [[nodiscard]] detail::floating_result_t<ResultNumber> distanceLInf(const OtherPoint& point) 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 = double, PointConcept OtherPoint>
    [[nodiscard]] detail::floating_result_t<ResultNumber>
    distanceL1(const Shape<OtherPoint>& other) const {
        using FloatingResult = detail::floating_result_t<ResultNumber>;
        return other.template distanceL1<FloatingResult>(*this);
    }

    template <class ResultNumber = double, PointConcept OtherPoint>
    [[nodiscard]] detail::floating_result_t<ResultNumber>
    distanceLInf(const Shape<OtherPoint>& other) const {
        using FloatingResult = detail::floating_result_t<ResultNumber>;
        return other.template distanceLInf<FloatingResult>(*this);
    }

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

    template <class ResultNumber = NumberType, PointConcept OtherPoint>
    [[nodiscard]] constexpr std::optional<Point<ResultNumber, typename PointType::LabelType>>
    intersection(const OtherPoint& other) const;

    template <class ResultNumber = NumberType, class EmptyPoint>
    [[nodiscard]] constexpr EmptyShape<EmptyPoint> intersection(const EmptyShape<EmptyPoint>&) const {
        return {};
    }

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

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool contains(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool contains(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool contains(const OtherOrientedLine& other) const;

    template <RayConcept OtherRay>
    [[nodiscard]] constexpr bool contains(const OtherRay& other) const;

    template <HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool contains(const OtherHalfplane& other) const;

    template <TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool contains(const OtherTriangle& other) const;

    template <RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool contains(const OtherRectangle& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool contains(const OtherConvex& other) const;

    template<PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool contains(const OtherPolygon& other) const;

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool contains(const OtherDisk& other) const;

    [[nodiscard]] constexpr bool contains(const Shape<PointType>& other) const;
 
    // The empty set is a subset of every shape, so its containment relations are
    // true; the symmetric intersection/crossing predicates reach the empty set
    // through Disk's existing generic OtherShape fallbacks.
    template <class EmptyPoint>
    [[nodiscard]] constexpr bool contains(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 interiorContains(const EmptyShape<EmptyPoint>&) const {
        return true;
    }

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

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

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool interiorContains(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool interiorContains(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool interiorContains(const OtherOrientedLine& other) const;

    template <RayConcept OtherRay>
    [[nodiscard]] constexpr bool interiorContains(const OtherRay& other) const;

    template <HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool interiorContains(const OtherHalfplane& other) const;

    template <TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool interiorContains(const OtherTriangle& other) const;

    template <RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool interiorContains(const OtherRectangle& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool interiorContains(const OtherConvex& other) const;

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool interiorContains(const OtherDisk& other) const;

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

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool intersects(const OtherDisk& other) const;

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

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

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool boundaryContains(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool boundaryContains(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool boundaryContains(const OtherOrientedLine& other) const;

    template <RayConcept OtherRay>
    [[nodiscard]] constexpr bool boundaryContains(const OtherRay& other) const;

    template <HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool boundaryContains(const OtherHalfplane& other) const;

    template <TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool boundaryContains(const OtherTriangle& other) const;

    template <RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool boundaryContains(const OtherRectangle& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool boundaryContains(const OtherConvex& other) const;

    template<PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool boundaryContains(const OtherPolygon& other) const;

    template <DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool boundaryContains(const OtherDisk& other) const;

    [[nodiscard]] constexpr bool boundaryContains(const Shape<PointType>& other) const;
 
    // --- not-yet-implemented predicate pairs (throw); see implementation ---
    template<PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool interiorContains(const OtherPolygon& other) const;

    template<PolygonConcept OtherPolygon>
    [[nodiscard]] constexpr bool separates(const OtherPolygon& other) const;

    template<MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool contains(const OtherChain& other) const;

    template<MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool boundaryContains(const OtherChain& other) const;

    template<MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool interiorContains(const OtherChain& other) const;

    template<MonotoneChainConcept OtherChain>
    [[nodiscard]] constexpr bool separates(const OtherChain& other) const;

    template<PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool contains(const OtherPolyline& other) const;

    template<PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool boundaryContains(const OtherPolyline& other) const;

    template<PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool interiorContains(const OtherPolyline& other) const;

    template<PolylineConcept OtherPolyline>
    [[nodiscard]] constexpr bool separates(const OtherPolyline& other) const;

    template<HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool contains(const OtherRegion& other) const;

    template<HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool boundaryContains(const OtherRegion& other) const;

    template<HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool interiorContains(const OtherRegion& other) const;

    template<HalfplaneIntersectionConcept OtherRegion>
    [[nodiscard]] constexpr bool separates(const OtherRegion& other) const;

    template<PolygonWithHolesConcept OtherRegion>
    [[nodiscard]] constexpr bool contains(const OtherRegion& other) const;

    template<PolygonWithHolesConcept OtherRegion>
    [[nodiscard]] constexpr bool boundaryContains(const OtherRegion& other) const;

    template<PolygonWithHolesConcept OtherRegion>
    [[nodiscard]] constexpr bool interiorContains(const OtherRegion& other) const;

    template<PolygonWithHolesConcept OtherRegion>
    [[nodiscard]] bool separates(const OtherRegion& other) const;
 
    // -------------------------------------------------------------------------
    // A set of regions
    //
    // It outranks every other shape, so the symmetric relations reach it through
    // the rank-based forwarders and only the asymmetric ones are answered here.
    // A set is the union of its components, so it is contained exactly when
    // every component is — no matter what this shape is.

    template<PolygonSetConcept OtherSet>
    [[nodiscard]] constexpr bool contains(const OtherSet& other) const {
        for (const auto& component : other) {
            if (!contains(component)) {
                return false;
            }
        }
        return true;
    }

    template<PolygonSetConcept OtherSet>
    [[nodiscard]] constexpr bool boundaryContains(const OtherSet& other) const {
        for (const auto& component : other) {
            if (!boundaryContains(component)) {
                return false;
            }
        }
        return true;
    }

    template<PolygonSetConcept OtherSet>
    [[nodiscard]] constexpr bool interiorContains(const OtherSet& other) const {
        for (const auto& component : other) {
            if (!interiorContains(component)) {
                return false;
            }
        }
        return true;
    }

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

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

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

    template <OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool separates(const OtherOrientedSegment& other) const;

    template <LineConcept OtherLine>
    [[nodiscard]] constexpr bool separates(const OtherLine& other) const;

    template <OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool separates(const OtherOrientedLine& other) const;

    template <ConvexConcept OtherConvex>
    [[nodiscard]] constexpr bool separates(const OtherConvex& other) const;

    template<RayConcept OtherRay>
    [[nodiscard]] constexpr bool separates(const OtherRay& other) const;

    template<HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool separates(const OtherHalfplane& other) const;

    template<RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool separates(const OtherRectangle& other) const;

    template<TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool separates(const OtherTriangle& other) const;

    template<DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool separates(const OtherDisk& other) const;
 
    // --- Symmetric predicate dispatch: intersects / interiorsIntersect / crosses ---
    // For a symmetric predicate a.method(b) == b.method(a). The canonical
    // implementor of each pair is the shape that appears later in pgl.hpp; the
    // earlier shape forwards to it. Disk is therefore canonical against every
    // shape up to and including Disk, and forwards to Convex/Polygon (which come
    // later). Pairs whose disk-side geometry is not yet implemented throw, in the
    // style of the other "not implemented yet" stubs.

    template<PointConcept OtherPoint>
    [[nodiscard]] constexpr bool intersects(const OtherPoint& other) const;
    template<OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool intersects(const OtherOrientedSegment& other) const;
    template<LineConcept OtherLine>
    [[nodiscard]] constexpr bool intersects(const OtherLine& other) const;
    template<OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool intersects(const OtherOrientedLine& other) const;
    template<RayConcept OtherRay>
    [[nodiscard]] constexpr bool intersects(const OtherRay& other) const;
    template<HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool intersects(const OtherHalfplane& other) const;
    template<RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool intersects(const OtherRectangle& other) const;
    template<TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool intersects(const OtherTriangle& other) const;

    template<PointConcept OtherPoint>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherPoint& other) const;
    template<SegmentConcept OtherSegment>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherSegment& other) const;
    template<OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherOrientedSegment& other) const;
    template<LineConcept OtherLine>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherLine& other) const;
    template<OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherOrientedLine& other) const;
    template<RayConcept OtherRay>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherRay& other) const;
    template<HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherHalfplane& other) const;
    template<RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherRectangle& other) const;
    template<TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherTriangle& other) const;
    template<DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherDisk& other) const;

    template<PointConcept OtherPoint>
    [[nodiscard]] constexpr bool crosses(const OtherPoint& other) const;
    template<SegmentConcept OtherSegment>
    [[nodiscard]] constexpr bool crosses(const OtherSegment& other) const;
    template<OrientedSegmentConcept OtherOrientedSegment>
    [[nodiscard]] constexpr bool crosses(const OtherOrientedSegment& other) const;
    template<LineConcept OtherLine>
    [[nodiscard]] constexpr bool crosses(const OtherLine& other) const;
    template<OrientedLineConcept OtherOrientedLine>
    [[nodiscard]] constexpr bool crosses(const OtherOrientedLine& other) const;
    template<RayConcept OtherRay>
    [[nodiscard]] constexpr bool crosses(const OtherRay& other) const;
    template<HalfplaneConcept OtherHalfplane>
    [[nodiscard]] constexpr bool crosses(const OtherHalfplane& other) const;
    template<RectangleConcept OtherRectangle>
    [[nodiscard]] constexpr bool crosses(const OtherRectangle& other) const;
    template<TriangleConcept OtherTriangle>
    [[nodiscard]] constexpr bool crosses(const OtherTriangle& other) const;
    template<DiskConcept OtherDisk>
    [[nodiscard]] constexpr bool crosses(const OtherDisk& other) const;

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

    template<typename OtherShape>
        requires (!PointConcept<OtherShape> && detail::shapeRank<OtherShape> > detail::shapeRank<Disk>)
    [[nodiscard]] constexpr bool intersects(const OtherShape& other) const {
        return other.intersects(*this);
    }

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

    template<typename OtherShape>
        requires (!PointConcept<OtherShape> && detail::shapeRank<OtherShape> > detail::shapeRank<Disk>)
    [[nodiscard]] constexpr bool interiorsIntersect(const OtherShape& other) const {
        return other.interiorsIntersect(*this);
    }

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

    template<typename OtherShape>
        requires (!PointConcept<OtherShape> && detail::shapeRank<OtherShape> > detail::shapeRank<Disk>)
    [[nodiscard]] constexpr bool crosses(const OtherShape& other) const {
        return other.crosses(*this);
    }

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

    template <class Scalar>
        requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
    constexpr Disk& operator/=(const Scalar& scalar) {
        for (auto& point : points_) {
            point /= scalar;
        }
        points_ = canonicalizePoints(points_[0], points_[1], points_[2]);
        return *this;
    }

    constexpr bool operator==(const Disk& other) const {
        // Degenerate == non-degenerate: false
        // Degenerate == degenerate with same points: true
        // Degenerate == degenerate with different points: false
        if (isDegenerate() || other.isDegenerate()) {
            if (isDegenerate() != other.isDegenerate()) {
                return false;
            }
 
            return a() == other.a() &&
                   b() == other.b() &&
                   c() == other.c();
        }
 
        // For true circles, three non-collinear boundary points define the
        // circle uniquely. So the other disk is the same if its three stored points lies on the disk boundary.
        return inCircleSign(a(), b(), c(), other.a()) == std::partial_ordering::equivalent &&
               inCircleSign(a(), b(), c(), other.b()) == std::partial_ordering::equivalent &&
               inCircleSign(a(), b(), c(), other.c()) == std::partial_ordering::equivalent;
    }

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

    constexpr std::partial_ordering operator<=>(const Disk& other) const {
        // Keep degenerate disks in a separate ordering class from true circles.
        if (isDegenerate() != other.isDegenerate()) {
            return isDegenerate() ? std::partial_ordering::less : std::partial_ordering::greater;
        }
 
        // Order true circles from smaller to larger, then by center.
        if (!isDegenerate()) {
            using ExactNumber = division_result_t<NumberType>;
            const auto radius_order = pgl::detail::threeWay(squaredRadius<ExactNumber>(),
                                                             other.template squaredRadius<ExactNumber>());
            if (radius_order != 0) {
                return radius_order;
            }
            return center<ExactNumber>() <=> other.template center<ExactNumber>();
        }
 
        // Degenerate disks fall back to their canonical point representation.
        const auto first_point_order = a() <=> other.a();
        if (first_point_order != 0) {
            return first_point_order;
        }
 
        const auto second_point_order = b() <=> other.b();
        if (second_point_order != 0) {
            return second_point_order;
        }
 
        return c() <=> other.c();
    }

    [[nodiscard]] constexpr Disk rotated90(int k = 1) const;

    constexpr void rotate90(int k = 1);

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

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

    template <class ResultNumber = double, DiskConcept OtherDisk>
    [[nodiscard]] std::optional<Disk<Point<ResultNumber, PointLabelType>>>
    minkowskiErosion(const OtherDisk& other) const;

    template <class ResultNumber = double, HalfplaneConcept OtherHalfplane>
    [[nodiscard]] EmptyShape<Point<ResultNumber, PointLabelType>>
    minkowskiErosion(const OtherHalfplane& other) const;

    template <class ResultNumber = double, DiskConcept OtherDisk>
    [[nodiscard]] Disk<Point<ResultNumber, PointLabelType>>
    minkowskiSum(const OtherDisk& other) const;

    template <class ResultNumber = double, HalfplaneConcept OtherHalfplane>
    [[nodiscard]] Halfplane<Point<ResultNumber, PointLabelType>>
    minkowskiSum(const OtherHalfplane& other) const;

    template <PointConcept OtherPoint>
    constexpr Disk& operator+=(const OtherPoint& translation) {
        for (auto& point : points_) {
            point += translation;
        }
        return *this;
    }

    template <PointConcept OtherPoint>
    constexpr Disk& operator-=(const OtherPoint& translation) {
        for (auto& point : points_) {
            point -= translation;
        }
        return *this;
    }

    template <class Scalar>
        requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
    constexpr Disk& operator*=(const Scalar& scalar) {
        for (auto& point : points_) {
            point *= scalar;
        }
        points_ = canonicalizePoints(points_[0], points_[1], points_[2]);
        return *this;
    }

    class PointIterator {
      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 = const PointType&;

        constexpr PointIterator() = default;

        constexpr reference operator*() const {
            assert(disk != nullptr);
            return (*disk)[index];
        }

        constexpr PointIterator& operator++() {
            ++index;
            return *this;
        }

        constexpr PointIterator operator++(int) {
            PointIterator copy(*this);
            ++(*this);
            return copy;
        }

        constexpr bool operator==(const PointIterator& other) const = default;
 
      private:
        friend struct Disk;

        constexpr PointIterator(const Disk* disk_arg, std::size_t index_arg)
            : disk(disk_arg), index(index_arg) {}
 
        const Disk* disk = nullptr;  
        std::size_t index = 0;       
    };

    template <class OtherShape>
    [[nodiscard]] constexpr bool pointInsideInteriorContainedIn(const OtherShape& shape) const {
        const auto p = pointInside<NumberType>();
        if (interiorContains(p)) {
            return shape.interiorContains(p);
        }
        return (shape * 2).interiorContains((*this * 2).template pointInside<NumberType>());
    }
 
  private:
    static constexpr std::array<PointType, 3> canonicalizePoints(PointType first, PointType second, PointType third) {
        std::array<PointType, 3> points{
            std::move(first),
            std::move(second),
            std::move(third),
        };
 
        // Sort the three stored boundary points lexicographically
        if (points[1] < points[0]) {
            std::swap(points[0], points[1]);
        }
        if (points[2] < points[1]) {
            std::swap(points[1], points[2]);
        }
        if (points[1] < points[0]) {
            std::swap(points[0], points[1]);
        }
 
        // Keep non-degenerate disks in counterclockwise order
        if (orientationSign(points[0], points[1], points[2]) == std::partial_ordering::less) {
            std::swap(points[1], points[2]);
        }
 
        return points;
    }

    [[nodiscard]] constexpr std::optional<std::pair<PointType, NumberType>> centerAndRadius() const {
        if (a().y() == b().y()) {
            NumberType r = c().y() - a().y();
            NumberType center_x = a().x() + r;
            if (r >= 0 && b().x() == center_x + r && c().x() == center_x) {
                return std::make_pair(PointType(center_x, a().y()), r);
            }
        }
        return std::nullopt;
    }
 
    std::array<PointType, 3> points_{};        
    [[no_unique_address]] mutable LabelType label_{};  
};



template <class PointType, class LabelType, class TranslationNumber, class TranslationLabel>
constexpr auto operator-(const Disk<PointType, LabelType>& disk, const Point<TranslationNumber, TranslationLabel>& translation) {
    auto result = disk;
    result -= translation;
    if constexpr (detail::has_label_v<LabelType>) {
        result.label() = LabelType{};
    }
    return result;
}

template <class PointType, class LabelType, class Scalar>
    requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
constexpr auto operator*(const Disk<PointType, LabelType>& disk, const Scalar& scalar) {
    auto result = disk;
    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 Disk<PointType, LabelType>& disk) {
    return disk * scalar;
}

template <class PointType, class LabelType, class Scalar>
    requires(!detail::is_point_v<Scalar> && !TransformationConcept<Scalar>)
constexpr auto operator/(const Disk<PointType, LabelType>& disk, const Scalar& scalar) {
    auto result = disk;
    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 Disk<PointType, LabelType>& disk);
 
}  // namespace pgl