distance.hpp

Go to the documentation of this file.

#pragma once
 
#include "implementation/intersection.hpp"
#include "shape/convex.hpp"

 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <ranges>
 
 
namespace pgl {
 
namespace detail {

template <class ResultNumber, class Self, class OtherShape>
constexpr ResultNumber maxVertexSquaredDistance(const Self& self, const OtherShape& other) {
    const auto self_vertices = self.vertices();
    ResultNumber worst = other.template squaredDistance<ResultNumber>(self_vertices[0]);
    for (std::size_t index = 1; index < self_vertices.size(); ++index) {
        const ResultNumber current = other.template squaredDistance<ResultNumber>(self_vertices[index]);
        if (worst < current) {
            worst = current;
        }
    }
    return worst;
}
 
}  // namespace detail
 
// -----------------------------------------------------------------------------
// Point
 
template <class Number, class Label>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Point<Number, Label>::squaredDistance(const OtherPoint& other) const {
    const ResultNumber dx = detail::asNumber<ResultNumber>(x()) - detail::asNumber<ResultNumber>(other.x());
    const ResultNumber dy = detail::asNumber<ResultNumber>(y()) - detail::asNumber<ResultNumber>(other.y());
    return dx * dx + dy * dy;
}
 
template <class Number, class Label>
template <class ApproximateNumber, PointConcept OtherPoint>
ApproximateNumber Point<Number, Label>::distance(const OtherPoint& other) const {
    return std::sqrt(this->template squaredDistance<ApproximateNumber>(other));
}
 
template <class Number, class Label>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Point<Number, Label>::distanceL1(const OtherPoint& other) const {
    const ResultNumber dx = detail::asNumber<ResultNumber>(x()) - detail::asNumber<ResultNumber>(other.x());
    const ResultNumber dy = detail::asNumber<ResultNumber>(y()) - detail::asNumber<ResultNumber>(other.y());
    return pgl::detail::abs(dx) + pgl::detail::abs(dy);
}
 
template <class Number, class Label>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Point<Number, Label>::distanceLInf(const OtherPoint& other) const {
    const ResultNumber dx = detail::asNumber<ResultNumber>(x()) - detail::asNumber<ResultNumber>(other.x());
    const ResultNumber dy = detail::asNumber<ResultNumber>(y()) - detail::asNumber<ResultNumber>(other.y());
    return std::max(pgl::detail::abs(dx), pgl::detail::abs(dy));
}
 
template <class Number, class Label>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Point<Number, Label>::squaredHausdorffDistance(const OtherPoint& other) const {
    return squaredDistance<ResultNumber>(other);
}
 
// -----------------------------------------------------------------------------
// Segment
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Segment<PointType, LabelType>::squaredDistance(const OtherPoint& point) const {
    const Point<ResultNumber> a(min());
    const Point<ResultNumber> b(max());
    const Point<ResultNumber> p(point);
 
    const auto ab = b - a;
    const auto ap = p - a;
    const auto bp = p - b;
    const auto squared_length = ab * ab;
 
    if (squared_length == 0) {
        return a.template squaredDistance<ResultNumber>(p);
    }
 
    // Which of the three pieces of the segment's Voronoi diagram p falls in.
    // Both clamps are sign tests on a projection, so dotSign settles them in
    // the promoted coordinate type rather than in the caller's ResultNumber —
    // the same reason Triangle::isRectangle prefers it to the plain dot
    // product operator. `ap . ab >= ab . ab` is `bp . ab >= 0`, which saves
    // comparing against the squared length as well.
    if (dotSign(ap, ab) <= 0) {
        return static_cast<ResultNumber>(ap * ap);
    }
 
    if (dotSign(bp, ab) >= 0) {
        return static_cast<ResultNumber>(bp * bp);
    }
 
    const ResultNumber twice_triangle_area = ab.x() * ap.y() - ab.y() * ap.x();
    return static_cast<ResultNumber>((twice_triangle_area * twice_triangle_area) / squared_length);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Segment<PointType, LabelType>::squaredDistance(const OtherSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    const auto this_min_to_other = other.template squaredDistance<ResultNumber>(min());
    const auto this_max_to_other = other.template squaredDistance<ResultNumber>(max());
    const auto other_min_to_this = this->template squaredDistance<ResultNumber>(other.min());
    const auto other_max_to_this = this->template squaredDistance<ResultNumber>(other.max());
 
    const auto best_from_this = this_min_to_other < this_max_to_other ? this_min_to_other : this_max_to_other;
    const auto best_from_other = other_min_to_this < other_max_to_this ? other_min_to_this : other_max_to_this;
 
    return best_from_this < best_from_other ? best_from_this : best_from_other;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Segment<PointType, LabelType>::squaredHausdorffDistance(const OtherSegment& other) const {
    const auto this_min_to_other = other.template squaredDistance<ResultNumber>(min());
    const auto this_max_to_other = other.template squaredDistance<ResultNumber>(max());
    const auto other_min_to_this = this->template squaredDistance<ResultNumber>(other.min());
    const auto other_max_to_this = this->template squaredDistance<ResultNumber>(other.max());
 
    const auto worst_from_this = this_min_to_other > this_max_to_other ? this_min_to_other : this_max_to_other;
    const auto worst_from_other = other_min_to_this > other_max_to_this ? other_min_to_this : other_max_to_this;
 
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Segment<PointType, LabelType>::squaredHausdorffDistance(const OtherPoint& point) const {
    // The point-side supremum (point-to-segment distance) is never larger than
    // the segment-side supremum (a min is never above a max), so the segment's
    // farthest endpoint from the point is the whole answer.
    return detail::maxVertexSquaredDistance<ResultNumber>(*this, point);
}
 
// -----------------------------------------------------------------------------
// OrientedSegment
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto OrientedSegment<PointType, LabelType>::squaredDistance(const OtherPoint& point) const {
    return static_cast<Segment<PointType>>(*this).template squaredDistance<ResultNumber>(point);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto OrientedSegment<PointType, LabelType>::squaredDistance(const OtherSegment& other) const {
    return static_cast<Segment<PointType>>(*this).template squaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto OrientedSegment<PointType, LabelType>::squaredDistance(const OtherOrientedSegment& other) const {
    return static_cast<Segment<PointType>>(*this).template squaredDistance<ResultNumber>(
        static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto OrientedSegment<PointType, LabelType>::squaredHausdorffDistance(const OtherSegment& other) const {
    return static_cast<Segment<PointType>>(*this).template squaredHausdorffDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto OrientedSegment<PointType, LabelType>::squaredHausdorffDistance(const OtherOrientedSegment& other) const {
    return static_cast<Segment<PointType>>(*this).template squaredHausdorffDistance<ResultNumber>(
        static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto OrientedSegment<PointType, LabelType>::squaredHausdorffDistance(const OtherPoint& point) const {
    return static_cast<Segment<PointType>>(*this).template squaredHausdorffDistance<ResultNumber>(point);
}
 
// -----------------------------------------------------------------------------
// Line
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Line<PointType, LabelType>::squaredDistance(const OtherPoint& point) const {
    if (isDegenerate()) {
        return min().template squaredDistance<ResultNumber>(point);
    }
 
    const Point<ResultNumber> a(min());
    const Point<ResultNumber> b(max());
    const Point<ResultNumber> p(point);
    const auto ab = b - a;
    const auto ap = p - a;
    const ResultNumber twice_triangle_area = ab.x() * ap.y() - ab.y() * ap.x();
    const ResultNumber squared_length = ab * ab;
    return static_cast<ResultNumber>((twice_triangle_area * twice_triangle_area) / squared_length);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto Line<PointType, LabelType>::squaredDistance(const OtherLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    return other.template squaredDistance<ResultNumber>(min());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Line<PointType, LabelType>::squaredDistance(const OtherSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    const auto source_distance = this->template squaredDistance<ResultNumber>(other.min());
    const auto target_distance = this->template squaredDistance<ResultNumber>(other.max());
    return source_distance < target_distance ? source_distance : target_distance;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Line<PointType, LabelType>::squaredDistance(const OtherOrientedSegment& other) const {
    return this->template squaredDistance<ResultNumber>(static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
// -----------------------------------------------------------------------------
// OrientedLine
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto OrientedLine<PointType, LabelType>::squaredDistance(const OtherPoint& point) const {
    return this->asLine().template squaredDistance<ResultNumber>(point);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto OrientedLine<PointType, LabelType>::squaredDistance(const OtherLine& other) const {
    return this->asLine().template squaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr auto OrientedLine<PointType, LabelType>::squaredDistance(const OtherOrientedLine& other) const {
    return this->asLine().template squaredDistance<ResultNumber>(other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto OrientedLine<PointType, LabelType>::squaredDistance(const OtherSegment& other) const {
    return this->asLine().template squaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto OrientedLine<PointType, LabelType>::squaredDistance(const OtherOrientedSegment& other) const {
    return this->asLine().template squaredDistance<ResultNumber>(other);
}
 
// -----------------------------------------------------------------------------
// Ray
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Ray<PointType, LabelType>::squaredDistance(const OtherPoint& point) const {
    if (isDegenerate()) {
        return source().template squaredDistance<ResultNumber>(point);
    }
 
    const Point<ResultNumber> a(source());
    const Point<ResultNumber> b(target());
    const Point<ResultNumber> p(point);
    const auto ab = b - a;
    const auto ap = p - a;
    const auto squared_length = ab * ab;
 
    if (squared_length == 0) {
        return static_cast<ResultNumber>(ap * ap);
    }
 
    // Behind the source, so the source is the nearest point; the projection's
    // sign is promoted for the same reason as the segment's clamps above.
    if (dotSign(ap, ab) <= 0) {
        return static_cast<ResultNumber>(ap * ap);
    }
 
    const ResultNumber twice_triangle_area = ab.x() * ap.y() - ab.y() * ap.x();
    return static_cast<ResultNumber>((twice_triangle_area * twice_triangle_area) / squared_length);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto Ray<PointType, LabelType>::squaredDistance(const OtherLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    return other.template squaredDistance<ResultNumber>(source());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr auto Ray<PointType, LabelType>::squaredDistance(const OtherOrientedLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    return other.template squaredDistance<ResultNumber>(source());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Ray<PointType, LabelType>::squaredDistance(const OtherSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    const auto source_to_other = other.template squaredDistance<ResultNumber>(source());
    const auto other_min_to_this = this->template squaredDistance<ResultNumber>(other.min());
    const auto other_max_to_this = this->template squaredDistance<ResultNumber>(other.max());
    const auto best_from_segment = other_min_to_this < other_max_to_this ? other_min_to_this : other_max_to_this;
    return source_to_other < best_from_segment ? source_to_other : best_from_segment;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Ray<PointType, LabelType>::squaredDistance(const OtherOrientedSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    const auto source_to_other = other.template squaredDistance<ResultNumber>(source());
    const auto other_source_to_this = this->template squaredDistance<ResultNumber>(other.source());
    const auto other_target_to_this = this->template squaredDistance<ResultNumber>(other.target());
    const auto best_from_segment = other_source_to_this < other_target_to_this ? other_source_to_this : other_target_to_this;
    return source_to_other < best_from_segment ? source_to_other : best_from_segment;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr auto Ray<PointType, LabelType>::squaredDistance(const OtherRay& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    const auto this_source_to_other = other.template squaredDistance<ResultNumber>(source());
    const auto other_source_to_this = this->template squaredDistance<ResultNumber>(other.source());
    return this_source_to_other < other_source_to_this ? this_source_to_other : other_source_to_this;
}
 
// -----------------------------------------------------------------------------
// Halfplane
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Halfplane<PointType, LabelType>::squaredDistance(const OtherPoint& point) const {
    if (intersects(point)) {
        return ResultNumber{};
    }
    return asLine().template squaredDistance<ResultNumber>(point);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Halfplane<PointType, LabelType>::squaredDistance(const OtherSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return asLine().template squaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Halfplane<PointType, LabelType>::squaredDistance(const OtherOrientedSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return asLine().template squaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto Halfplane<PointType, LabelType>::squaredDistance(const OtherLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return asLine().template squaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr auto Halfplane<PointType, LabelType>::squaredDistance(const OtherOrientedLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return asLine().template squaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr auto Halfplane<PointType, LabelType>::squaredDistance(const OtherRay& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return asLine().template squaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr auto Halfplane<PointType, LabelType>::squaredDistance(const OtherHalfplane& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return asLine().template squaredDistance<ResultNumber>(other.asLine());
}
 
// -----------------------------------------------------------------------------
// Rectangle
 
template <class PointType, class LabelType>
template <class Left, class Right>
constexpr auto Rectangle<PointType, LabelType>::axisDistance(const Left& first_min, const Left& first_max, const Right& second_min, const Right& second_max)
    -> std::common_type_t<Left, Right> {
    using Distance = std::common_type_t<Left, Right>;
 
    if (first_max < second_min) {
        return static_cast<Distance>(second_min) - static_cast<Distance>(first_max);
    }
 
    if (second_max < first_min) {
        return static_cast<Distance>(first_min) - static_cast<Distance>(second_max);
    }
 
    return Distance{};
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Rectangle<PointType, LabelType>::squaredDistance(const OtherPoint& point) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    const ResultNumber dx = static_cast<ResultNumber>(axisDistance(min().x(), max().x(), point.x(), point.x()));
    const ResultNumber dy = static_cast<ResultNumber>(axisDistance(min().y(), max().y(), point.y(), point.y()));
    return dx * dx + dy * dy;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto Rectangle<PointType, LabelType>::squaredDistance(const OtherLine& other) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    const auto rectangle_vertices = vertices();
    auto best_distance = other.template squaredDistance<ResultNumber>(rectangle_vertices[0]);
    for (std::size_t index = 1; index < rectangle_vertices.size(); ++index) {
        const auto current_distance = other.template squaredDistance<ResultNumber>(rectangle_vertices[index]);
        if (current_distance < best_distance) {
            best_distance = current_distance;
        }
    }
 
    return best_distance;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr auto Rectangle<PointType, LabelType>::squaredDistance(const OtherOrientedLine& other) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    return this->template squaredDistance<ResultNumber>(other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Rectangle<PointType, LabelType>::squaredDistance(const OtherSegment& other) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    auto best_distance = this->template squaredDistance<ResultNumber>(other.min());
    const auto other_max_distance = this->template squaredDistance<ResultNumber>(other.max());
    if (other_max_distance < best_distance) {
        best_distance = other_max_distance;
    }
 
    const auto rectangle_edges = edges();
    for (const auto& edge : rectangle_edges) {
        const auto current_distance = edge.template squaredDistance<ResultNumber>(other);
        if (current_distance < best_distance) {
            best_distance = current_distance;
        }
    }
 
    return best_distance;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Rectangle<PointType, LabelType>::squaredDistance(const OtherOrientedSegment& other) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    return this->template squaredDistance<ResultNumber>(static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr auto Rectangle<PointType, LabelType>::squaredDistance(const OtherRay& other) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    auto best_distance = this->template squaredDistance<ResultNumber>(other.source());
    const auto rectangle_edges = edges();
    for (const auto& edge : rectangle_edges) {
        const auto current_distance = other.template squaredDistance<ResultNumber>(edge);
        if (current_distance < best_distance) {
            best_distance = current_distance;
        }
    }
 
    return best_distance;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr auto Rectangle<PointType, LabelType>::squaredDistance(const OtherHalfplane& other) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template squaredDistance<ResultNumber>(other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto Rectangle<PointType, LabelType>::squaredDistance(const OtherRectangle& other) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    const ResultNumber dx = static_cast<ResultNumber>(axisDistance(min().x(), max().x(), other.min().x(), other.max().x()));
    const ResultNumber dy = static_cast<ResultNumber>(axisDistance(min().y(), max().y(), other.min().y(), other.max().y()));
    return dx * dx + dy * dy;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto Rectangle<PointType, LabelType>::squaredHausdorffDistance(const OtherRectangle& other) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Rectangle<PointType, LabelType>::squaredHausdorffDistance(const OtherPoint& point) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    return detail::maxVertexSquaredDistance<ResultNumber>(*this, point);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Rectangle<PointType, LabelType>::squaredHausdorffDistance(const OtherSegment& other) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Rectangle<PointType, LabelType>::squaredHausdorffDistance(const OtherOrientedSegment& other) const {
    // There is no nearest point of the empty set, so it has no distance to
    // anything; asking for one is a precondition violation.
    assert(!empty());
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
// -----------------------------------------------------------------------------
// Triangle
 
template <class PointType, class LabelType>
template <class ResultNumber, class OtherShape>
constexpr ResultNumber Triangle<PointType, LabelType>::edgeMinSquaredDistance(const OtherShape& other) const {
    const auto triangle_edges = edges();
    auto best = triangle_edges[0].template squaredDistance<ResultNumber>(other);
    for (std::size_t index = 1; index < triangle_edges.size(); ++index) {
        const auto current = triangle_edges[index].template squaredDistance<ResultNumber>(other);
        if (current < best) {
            best = current;
        }
    }
    return best;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, class OtherShape>
constexpr ResultNumber Triangle<PointType, LabelType>::vertexMinSquaredDistance(const OtherShape& other) const {
    const auto triangle_vertices = vertices();
    auto best = other.template squaredDistance<ResultNumber>(triangle_vertices[0]);
    for (std::size_t index = 1; index < triangle_vertices.size(); ++index) {
        const auto current = other.template squaredDistance<ResultNumber>(triangle_vertices[index]);
        if (current < best) {
            best = current;
        }
    }
    return best;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Triangle<PointType, LabelType>::squaredDistance(const OtherPoint& point) const {
    if (intersects(point)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(point);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Triangle<PointType, LabelType>::squaredDistance(const OtherSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Triangle<PointType, LabelType>::squaredDistance(const OtherOrientedSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto Triangle<PointType, LabelType>::squaredDistance(const OtherLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template vertexMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr auto Triangle<PointType, LabelType>::squaredDistance(const OtherOrientedLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template vertexMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr auto Triangle<PointType, LabelType>::squaredDistance(const OtherRay& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr auto Triangle<PointType, LabelType>::squaredDistance(const OtherHalfplane& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template squaredDistance<ResultNumber>(other.asLine());
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto Triangle<PointType, LabelType>::squaredDistance(const OtherRectangle& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr auto Triangle<PointType, LabelType>::squaredDistance(const OtherTriangle& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Triangle<PointType, LabelType>::squaredHausdorffDistance(const OtherPoint& point) const {
    return detail::maxVertexSquaredDistance<ResultNumber>(*this, point);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Triangle<PointType, LabelType>::squaredHausdorffDistance(const OtherSegment& other) const {
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Triangle<PointType, LabelType>::squaredHausdorffDistance(const OtherOrientedSegment& other) const {
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto Triangle<PointType, LabelType>::squaredHausdorffDistance(const OtherRectangle& other) const {
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr auto Triangle<PointType, LabelType>::squaredHausdorffDistance(const OtherTriangle& other) const {
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
// Convex
 
template <class PointType_, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Convex<PointType_, LabelType>::squaredDistance(const OtherPoint& point) const {
    // Promoted so a query point of a different coordinate type (e.g. a Disk's
    // inherently-double center) combines safely with this polygon's own vertex
    // type, matching the promotion orientationDeterminant already applies below.
    using Num = detail::dot_coordinate_t<NumberType, typename OtherPoint::NumberType>;
    if (contains(point)) {
        return ResultNumber{};
    }
    if (size() <= 32) {
        auto edgeVector = edges();
        ResultNumber best = edgeVector[0].template squaredDistance<ResultNumber>(point);
        for (auto & e : edgeVector) {
            const ResultNumber current = e.template squaredDistance<ResultNumber>(point);
            if (current < best) {
                best = current;
            }
        }
        return best;
    }
 
    const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(size());
    const auto V = [this](std::ptrdiff_t i) { return get(i); };
    const auto facing = [&](std::ptrdiff_t i) {
        return orientationDeterminant(V(i), V(i + 1), point) < 0;
    };
 
    // Reference point c = (v0 + v1 + v2) / 3, strictly interior. To stay exact we
    // work with vectors scaled by 3 (a positive factor, so signs are preserved).
    const Num sx = detail::asNumber<Num>(V(0).x()) + detail::asNumber<Num>(V(1).x()) + detail::asNumber<Num>(V(2).x());
    const Num sy = detail::asNumber<Num>(V(0).y()) + detail::asNumber<Num>(V(1).y()) + detail::asNumber<Num>(V(2).y());
    const Num bx = Num{3} * detail::asNumber<Num>(point.x()) - sx;   // 3 * (q - c)
    const Num by = Num{3} * detail::asNumber<Num>(point.y()) - sy;
 
    // Sign of (V(k) - c) x (q - c). Around the boundary this changes sign exactly
    // twice: at the edge the ray c->q exits through (which faces q) and at the
    // opposite edge (which does not).
    const auto dSign = [&](std::ptrdiff_t k) -> int {
        const Num ax = Num{3} * detail::asNumber<Num>(V(k).x()) - sx;   // 3 * (V(k) - c)
        const Num ay = Num{3} * detail::asNumber<Num>(V(k).y()) - sy;
        const Num cross = ax * by - ay * bx;
        return (cross > Num{0}) - (cross < Num{0});
    };
    // D(k) above is linear in V(k) with gradient (by, -bx); its argmax/argmin are
    // the genuinely unimodal support function, so they bracket the two sign
    // changes and are found in O(log n) with the existing cyclic search.
    const auto gDot = [&](std::ptrdiff_t k) {
        return detail::asNumber<Num>(V(k).x()) * by - detail::asNumber<Num>(V(k).y()) * bx;
    };
    const auto indices = std::views::iota(std::ptrdiff_t{0}, n);
    const std::ptrdiff_t mPos = *detail::cyclicMax(
        indices.begin(), indices.end(), [&](std::ptrdiff_t k) { return gDot(k); });
    const std::ptrdiff_t mNeg = *detail::cyclicMax(
        indices.begin(), indices.end(), [&](std::ptrdiff_t k) { return -gDot(k); });
 
    const auto wrap = [n](std::ptrdiff_t a) { return ((a % n) + n) % n; };
    // Largest offset in [0, len] for which pred(from + step * offset) holds; the
    // predicate is true...false over the bracket, so a binary search locates the
    // transition.
    const auto lastWhile =
        [&](std::ptrdiff_t from, std::ptrdiff_t step, std::ptrdiff_t len, auto pred) {
            std::ptrdiff_t lo = 0, hi = len;
            while (lo < hi) {
                const std::ptrdiff_t mid = (lo + hi + 1) / 2;
                if (pred(from + step * mid)) {
                    lo = mid;
                } else {
                    hi = mid - 1;
                }
            }
            return lo;
        };
 
    // Exit edge (faces q): the + -> - sign change between the two D extrema.
    const std::ptrdiff_t anchor = mPos + lastWhile(mPos, 1, wrap(mNeg - mPos),
        [&](std::ptrdiff_t k) { return dSign(k) > 0; });
    // Back edge (does not face q): the - -> + sign change on the other side.
    const std::ptrdiff_t back = mNeg + lastWhile(mNeg, 1, wrap(mPos - mNeg),
        [&](std::ptrdiff_t k) { return dSign(k) < 0; });
 
    // Grow the contiguous chain of q-facing edges around the anchor.
    const std::ptrdiff_t hiEdge = anchor + lastWhile(anchor, 1, wrap(back - anchor), facing);
    const std::ptrdiff_t loEdge = anchor - lastWhile(anchor, -1, wrap(anchor - back), facing);
 
    // Over the facing chain the per-edge distance is unimodal: pick the closest.
    const auto edgeDist = [&](std::ptrdiff_t i) {
        return Segment<PointType>(V(i), V(i + 1)).template squaredDistance<ResultNumber>(point);
    };
    std::ptrdiff_t lo = 0, hi = hiEdge - loEdge;
    while (lo < hi) {
        const std::ptrdiff_t mid = (lo + hi) / 2;
        if (edgeDist(loEdge + mid) < edgeDist(loEdge + mid + 1)) {
            hi = mid;
        } else {
            lo = mid + 1;
        }
    }
    return edgeDist(loEdge + lo);
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Convex<PointType_, LabelType>::squaredDistance(const OtherSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    if (size() <= 32) {
        auto edgeVector = edges();
        ResultNumber best = edgeVector[0].template squaredDistance<ResultNumber>(other);
        for (auto& e : edgeVector) {
            const ResultNumber current = e.template squaredDistance<ResultNumber>(other);
            if (current < best) {
                best = current;
            }
        }
        return best;
    }
 
    // Disjoint, large polygon. Two witnesses come from the segment endpoints; the
    // third (polygon vertex closest to the segment interior) lies on a vertex
    // extremal along the segment normal, found in O(log n).
    const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(size());
 
    ResultNumber best = this->template squaredDistance<ResultNumber>(other.min());
    const ResultNumber from_max = this->template squaredDistance<ResultNumber>(other.max());
    if (from_max < best) {
        best = from_max;
    }
 
    // Support functional along the segment normal. The orientation determinant of
    // (segment.min, segment.max, vertex) is the signed-distance functional — it is
    // linear, so it is cyclic-unimodal over the CCW vertices and cyclicMax locates
    // its extrema; those two vertices are the polygon's nearest to the segment's
    // supporting line from either side, the only interior-witness candidates. The
    // orientation primitive keeps the coordinate promotion overflow-safe.
    const PointType a(other.min());
    const PointType b(other.max());
    const auto support = [&](std::ptrdiff_t i) {
        return orientationDeterminant(a, b, get(i));
    };
    const auto indices = std::views::iota(std::ptrdiff_t{0}, n);
    const std::ptrdiff_t hi = *detail::cyclicMax(
        indices.begin(), indices.end(), [&](std::ptrdiff_t k) { return support(k); });
    const std::ptrdiff_t lo = *detail::cyclicMax(
        indices.begin(), indices.end(), [&](std::ptrdiff_t k) { return -support(k); });
 
    const ResultNumber from_hi = other.template squaredDistance<ResultNumber>(get(hi));
    if (from_hi < best) {
        best = from_hi;
    }
    const ResultNumber from_lo = other.template squaredDistance<ResultNumber>(get(lo));
    if (from_lo < best) {
        best = from_lo;
    }
 
    return best;
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Convex<PointType_, LabelType>::squaredDistance(const OtherOrientedSegment& other) const {
    return this->template squaredDistance<ResultNumber>(
        static_cast<Segment<typename OtherOrientedSegment::PointType>>(other));
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr auto Convex<PointType_, LabelType>::squaredDistance(const OtherConvex& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    // Disjoint convex polygons: the closest point of one lies on its boundary, so
    // the minimum over the edges of either polygon of the edge-to-polygon distance
    // is the polygon-to-polygon distance. Query the edges of the polygon with
    // fewer vertices to minimize the number of O(log n) segment searches.
    const auto minOverEdges = [](const auto& source, const auto& target) {
        const auto edgeVector = source.edges();
        ResultNumber best = target.template squaredDistance<ResultNumber>(edgeVector[0]);
        for (const auto& edge : edgeVector) {
            const ResultNumber current = target.template squaredDistance<ResultNumber>(edge);
            if (current < best) {
                best = current;
            }
        }
        return best;
    };
 
    if (size() <= other.size()) {
        return minOverEdges(*this, other);
    }
    return minOverEdges(other, *this);
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr auto Convex<PointType_, LabelType>::squaredDistance(const OtherTriangle& other) const {
    return this->template squaredDistance<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto Convex<PointType_, LabelType>::squaredDistance(const OtherRectangle& other) const {
    return this->template squaredDistance<ResultNumber>(other.asConvex());
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto Convex<PointType_, LabelType>::squaredDistance(const OtherLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    if (size() <= 8) {
        ResultNumber best = other.template squaredDistance<ResultNumber>(get(0));
        for (std::ptrdiff_t i = 1; i < static_cast<std::ptrdiff_t>(size()); ++i) {
            const ResultNumber current = other.template squaredDistance<ResultNumber>(get(i));
            if (current < best) {
                best = current;
            }
        }
        return best;
    }
 
    // Disjoint line: the whole polygon lies on one side, so its closest point to
    // the line is the vertex extremal along the line normal toward the line. The
    // orientation determinant of (line.min, line.max, vertex) is the signed
    // distance functional — linear, hence cyclic-unimodal over the CCW vertices.
    // cyclicMax locates its two extrema (the polygon's nearest and farthest
    // vertices from the line), the only closest-vertex candidates. Using the
    // orientation primitive keeps the coordinate promotion overflow-safe.
    const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(size());
    const PointType a(other.min());
    const PointType b(other.max());
    const auto support = [&](std::ptrdiff_t i) {
        return orientationDeterminant(a, b, get(i));
    };
    const auto indices = std::views::iota(std::ptrdiff_t{0}, n);
    const std::ptrdiff_t hi = *detail::cyclicMax(
        indices.begin(), indices.end(), [&](std::ptrdiff_t k) { return support(k); });
    const std::ptrdiff_t lo = *detail::cyclicMax(
        indices.begin(), indices.end(), [&](std::ptrdiff_t k) { return -support(k); });
 
    const ResultNumber from_hi = other.template squaredDistance<ResultNumber>(get(hi));
    const ResultNumber from_lo = other.template squaredDistance<ResultNumber>(get(lo));
    return from_hi < from_lo ? from_hi : from_lo;
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr auto Convex<PointType_, LabelType>::squaredDistance(const OtherOrientedLine& other) const {
    return this->template squaredDistance<ResultNumber>(other.asLine());
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr auto Convex<PointType_, LabelType>::squaredDistance(const OtherRay& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    if (size() <= 16) {
        auto edgeVector = edges();
        ResultNumber best = other.template squaredDistance<ResultNumber>(edgeVector[0]);
        for (const auto& e : edgeVector) {
            const ResultNumber current = other.template squaredDistance<ResultNumber>(e);
            if (current < best) {
                best = current;
            }
        }
        return best;
    }
 
    // Disjoint, large polygon. One witness comes from the ray's source; the other
    // (polygon vertex closest to the ray interior) lies on a vertex extremal along
    // the ray normal, found in O(log n). The orientation determinant of
    // (ray.source, ray.target, vertex) is the signed-distance functional — linear,
    // hence cyclic-unimodal over the CCW vertices — and stays overflow-safe.
    const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(size());
 
    ResultNumber best = this->template squaredDistance<ResultNumber>(other.source());
 
    const PointType a(other.source());
    const PointType b(other.target());
    const auto support = [&](std::ptrdiff_t i) {
        return orientationDeterminant(a, b, get(i));
    };
    const auto indices = std::views::iota(std::ptrdiff_t{0}, n);
    const std::ptrdiff_t hi = *detail::cyclicMax(
        indices.begin(), indices.end(), [&](std::ptrdiff_t k) { return support(k); });
    const std::ptrdiff_t lo = *detail::cyclicMax(
        indices.begin(), indices.end(), [&](std::ptrdiff_t k) { return -support(k); });
 
    const ResultNumber from_hi = other.template squaredDistance<ResultNumber>(get(hi));
    if (from_hi < best) {
        best = from_hi;
    }
    const ResultNumber from_lo = other.template squaredDistance<ResultNumber>(get(lo));
    if (from_lo < best) {
        best = from_lo;
    }
 
    return best;
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr auto Convex<PointType_, LabelType>::squaredDistance(const OtherHalfplane& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
 
    // Disjoint: the polygon lies entirely on the far side of the boundary line,
    // so the closest point of the half-plane is on its boundary.
    return this->template squaredDistance<ResultNumber>(other.asLine());
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Convex<PointType_, LabelType>::squaredHausdorffDistance(const OtherPoint& point) const {
    return detail::maxVertexSquaredDistance<ResultNumber>(*this, point);
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Convex<PointType_, LabelType>::squaredHausdorffDistance(const OtherSegment& other) const {
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Convex<PointType_, LabelType>::squaredHausdorffDistance(const OtherOrientedSegment& other) const {
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto Convex<PointType_, LabelType>::squaredHausdorffDistance(const OtherRectangle& other) const {
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr auto Convex<PointType_, LabelType>::squaredHausdorffDistance(const OtherTriangle& other) const {
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
template <class PointType_, class LabelType>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr auto Convex<PointType_, LabelType>::squaredHausdorffDistance(const OtherConvex& other) const {
    const auto worst_from_this = detail::maxVertexSquaredDistance<ResultNumber>(*this, other);
    const auto worst_from_other = detail::maxVertexSquaredDistance<ResultNumber>(other, *this);
    return worst_from_this > worst_from_other ? worst_from_this : worst_from_other;
}
 
// -----------------------------------------------------------------------------
// Disk
 
namespace detail {

template <class Float = double, class DiskType, class OtherShape>
Float diskExteriorSquaredDistance(const DiskType& disk, const OtherShape& other) {
    const Float gap =
        std::sqrt(other.template squaredDistance<Float>(disk.template center<Float>()))
        - disk.template radius<Float>();
    return gap * gap;
}
 
}  // namespace detail
 
template <class PointType_, class LabelType>
template <class ResultNumber, DiskConcept OtherDisk>
detail::floating_result_t<ResultNumber> Convex<PointType_, LabelType>::squaredDistance(
    const OtherDisk& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(other, *this);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PointConcept OtherPoint>
detail::floating_result_t<ResultNumber> Disk<PointType_, TLabel>::squaredDistance(
    const OtherPoint& point) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (contains(point)) {
        return Float{0};
    }
 
    // Exterior point: the nearest point of the disk is on the circle, so the
    // distance is |point - center| - radius. The squared distance is generally
    // irrational, hence the floating result. Evaluate the gap directly as
    // sqrt(dc2) - sqrt(r2); this well-conditioned form avoids the catastrophic
    // cancellation of the algebraically equal dc2 + r2 - 2*sqrt(dc2 * r2).
    const Float gap = center<Float>().template distance<Float>(point) - radius<Float>();
    return gap * gap;
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, SegmentConcept OtherSegment>
detail::floating_result_t<ResultNumber> Disk<PointType_, TLabel>::squaredDistance(
    const OtherSegment& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
detail::floating_result_t<ResultNumber> Disk<PointType_, TLabel>::squaredDistance(
    const OtherOrientedSegment& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, LineConcept OtherLine>
detail::floating_result_t<ResultNumber> Disk<PointType_, TLabel>::squaredDistance(
    const OtherLine& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
detail::floating_result_t<ResultNumber> Disk<PointType_, TLabel>::squaredDistance(
    const OtherOrientedLine& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RayConcept OtherRay>
detail::floating_result_t<ResultNumber> Disk<PointType_, TLabel>::squaredDistance(
    const OtherRay& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
detail::floating_result_t<ResultNumber> Disk<PointType_, TLabel>::squaredDistance(
    const OtherHalfplane& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
detail::floating_result_t<ResultNumber> Disk<PointType_, TLabel>::squaredDistance(
    const OtherRectangle& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
detail::floating_result_t<ResultNumber> Disk<PointType_, TLabel>::squaredDistance(
    const OtherTriangle& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(*this, other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, DiskConcept OtherDisk>
detail::floating_result_t<ResultNumber> Disk<PointType_, TLabel>::squaredDistance(
    const OtherDisk& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
 
    // Disjoint disks: the nearest points lie on the two circles along the line
    // through the centers, so the distance is the center separation minus both
    // radii. Subtracting the radii from the directly computed center distance is
    // well-conditioned (no catastrophic cancellation in the squared form).
    const Float gap = center<Float>().template distance<Float>(other.template center<Float>())
                      - radius<Float>() - other.template radius<Float>();
    return gap * gap;
}
 
// -----------------------------------------------------------------------------
// Polygon
 
template <class PointType_, class TLabel>
template <class ResultNumber, class OtherShape>
constexpr ResultNumber Polygon<PointType_, TLabel>::edgeMinSquaredDistance(const OtherShape& other) const {
    const auto boundaryEdges = edges();
    ResultNumber best = boundaryEdges[0].template squaredDistance<ResultNumber>(other);
    for (std::size_t index = 1; index < boundaryEdges.size(); ++index) {
        const ResultNumber current = boundaryEdges[index].template squaredDistance<ResultNumber>(other);
        if (current < best) {
            best = current;
        }
    }
    return best;
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherPoint& point) const {
    if (intersects(point)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(point);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherOrientedSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherOrientedLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RayConcept OtherRay>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherRay& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherHalfplane& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherRectangle& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherTriangle& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherConvex& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolygonConcept OtherPolygon>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherPolygon& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, class DiskPointType, class DiskLabel>
detail::floating_result_t<ResultNumber> Polygon<PointType_, TLabel>::squaredDistance(
    const Disk<DiskPointType, DiskLabel>& disk) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(disk)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(disk, *this);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, MonotoneChainConcept OtherChain>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherChain& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
// -----------------------------------------------------------------------------
// MonotoneChain
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, class OtherShape>
constexpr ResultNumber MonotoneChain<PointType, LabelType, Storage>::edgeMinSquaredDistance(const OtherShape& other) const {
    assert(size() >= 1);
    if (size() == 1) {
        // A chain collapsed to a single vertex has no edge to measure from;
        // its distance is the distance from that vertex.
        return (*this)[0].template squaredDistance<ResultNumber>(other);
    }
    ResultNumber best = this->template boundaryAt<false>(0).template squaredDistance<ResultNumber>(other);
    for (std::size_t index = 1; index + 1 < size(); ++index) {
        const ResultNumber current =
            this->template boundaryAt<false>(index).template squaredDistance<ResultNumber>(other);
        if (current < best) {
            best = current;
        }
    }
    return best;
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherPoint& point) const {
    if (intersects(point)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(point);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherOrientedSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherOrientedLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, RayConcept OtherRay>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherRay& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherHalfplane& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherRectangle& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherTriangle& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherConvex& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, MonotoneChainConcept OtherChain>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::squaredDistance(const OtherChain& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType, class Storage>
template <class ResultNumber, class DiskPointType, class DiskLabel>
detail::floating_result_t<ResultNumber> MonotoneChain<PointType, LabelType, Storage>::squaredDistance(
    const Disk<DiskPointType, DiskLabel>& disk) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(disk)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(disk, *this);
}
 
// -----------------------------------------------------------------------------
// Polyline
 
template <class PointType, class LabelType>
template <class ResultNumber, class OtherShape>
constexpr ResultNumber Polyline<PointType, LabelType>::edgeMinSquaredDistance(const OtherShape& other) const {
    assert(size() >= 1);
    if (size() == 1) {
        // A polyline collapsed to a single vertex has no edge to measure from;
        // its distance is the distance from that vertex.
        return (*this)[0].template squaredDistance<ResultNumber>(other);
    }
    ResultNumber best = this->template boundaryAt<false>(0).template squaredDistance<ResultNumber>(other);
    for (std::size_t index = 1; index + 1 < size(); ++index) {
        const ResultNumber current =
            this->template boundaryAt<false>(index).template squaredDistance<ResultNumber>(other);
        if (current < best) {
            best = current;
        }
    }
    return best;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherPoint& point) const {
    if (intersects(point)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(point);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PolylineConcept OtherPolyline>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherPolyline& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherOrientedSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherOrientedLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherRay& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherHalfplane& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherRectangle& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherTriangle& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherConvex& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, MonotoneChainConcept OtherChain>
constexpr auto Polyline<PointType, LabelType>::squaredDistance(const OtherChain& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, class DiskPointType, class DiskLabel>
detail::floating_result_t<ResultNumber> Polyline<PointType, LabelType>::squaredDistance(
    const Disk<DiskPointType, DiskLabel>& disk) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(disk)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(disk, *this);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, PolylineConcept OtherPolyline>
constexpr auto Polygon<PointType_, TLabel>::squaredDistance(const OtherPolyline& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
 
// -----------------------------------------------------------------------------
// HalfplaneIntersection
//
// When the region and the other shape are disjoint, the closest point of the
// convex region lies on its boundary, so the distance is the minimum over the
// region's boundary edges (segments, rays, or lines) of the edge-to-shape
// distance. Edges are built with an exact coordinate type; the caller's
// ResultNumber governs only the final cast. Querying the empty region is
// undefined behavior.
 
namespace detail {
 
template <class ResultNumber, class Region, class Other>
constexpr ResultNumber regionEdgesSquaredDistance(const Region& region, const Other& other) {
    ResultNumber best{};
    bool has = false;
    for (std::size_t i = 0; i < region.size(); ++i) {
        const ResultNumber current = std::visit(
            [&other](const auto& piece) {
                return static_cast<ResultNumber>(piece.template squaredDistance<ResultNumber>(other));
            },
            region.template edge<ResultNumber>(i));
        if (!has || current < best) {
            best = current;
            has = true;
        }
    }
    return best;
}
 
}  // namespace detail
 
#define PGL_HPI_SQUARED_DISTANCE(ConceptName, ArgType)                                        \
    template <class PointType, class LabelType>                                                \
    template <class ResultNumber, ConceptName ArgType>                                         \
    constexpr auto HalfplaneIntersection<PointType, LabelType>::squaredDistance(               \
        const ArgType& other) const {                                                          \
        if (intersects(other)) {                                                                \
            return ResultNumber{};                                                              \
        }                                                                                       \
        return detail::regionEdgesSquaredDistance<ResultNumber>(*this, other);                  \
    }
 
PGL_HPI_SQUARED_DISTANCE(PointConcept, OtherPoint)
PGL_HPI_SQUARED_DISTANCE(SegmentConcept, OtherSegment)
PGL_HPI_SQUARED_DISTANCE(OrientedSegmentConcept, OtherOrientedSegment)
PGL_HPI_SQUARED_DISTANCE(LineConcept, OtherLine)
PGL_HPI_SQUARED_DISTANCE(OrientedLineConcept, OtherOrientedLine)
PGL_HPI_SQUARED_DISTANCE(RayConcept, OtherRay)
PGL_HPI_SQUARED_DISTANCE(HalfplaneConcept, OtherHalfplane)
PGL_HPI_SQUARED_DISTANCE(RectangleConcept, OtherRectangle)
PGL_HPI_SQUARED_DISTANCE(TriangleConcept, OtherTriangle)
PGL_HPI_SQUARED_DISTANCE(ConvexConcept, OtherConvex)
PGL_HPI_SQUARED_DISTANCE(MonotoneChainConcept, OtherChain)
PGL_HPI_SQUARED_DISTANCE(PolylineConcept, OtherPolyline)
PGL_HPI_SQUARED_DISTANCE(PolygonConcept, OtherPolygon)
 
#undef PGL_HPI_SQUARED_DISTANCE
 
template <class PointType, class LabelType>
template <class ResultNumber, DiskConcept OtherDisk>
detail::floating_result_t<ResultNumber>
HalfplaneIntersection<PointType, LabelType>::squaredDistance(const OtherDisk& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (intersects(other)) {
        return Float{0};
    }
    // A distance to a disk is irrational, so the whole scan runs in Float; the
    // edges themselves stay exact, hence the separate E.
    using E = detail::region_exact_number_t<NumberType>;
    Float best{0};
    bool has = false;
    for (std::size_t i = 0; i < size(); ++i) {
        const Float current = std::visit(
            [&other](const auto& piece) { return other.template squaredDistance<Float>(piece); },
            this->template edge<E>(i));
        if (!has || current < best) {
            best = current;
            has = true;
        }
    }
    return best;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneIntersectionConcept OtherRegion>
constexpr auto HalfplaneIntersection<PointType, LabelType>::squaredDistance(const OtherRegion& other) const {
    // Disjoint convex regions realize their distance on this region's
    // boundary; each edge re-dispatches into the other region's own edge scan.
    if (intersects(other)) {
        return ResultNumber{};
    }
    return detail::regionEdgesSquaredDistance<ResultNumber>(*this, other);
}
 
 
// -----------------------------------------------------------------------------
// PolygonWithHoles
//
// The region is closed, so a shape it misses is nearest to a point of ∂A — the
// edges of the outer ring and of every hole. One scan serves every operand.
 
template <class PointType, class LabelType>
template <class ResultNumber, class OtherShape>
constexpr ResultNumber PolygonWithHoles<PointType, LabelType>::edgeMinSquaredDistance(const OtherShape& other) const {
    ResultNumber best{};
    bool seeded = false;
    anyBoundaryEdge([&](const auto& edge) {
        const ResultNumber current = edge.template squaredDistance<ResultNumber>(other);
        if (!seeded || current < best) {
            best = current;
            seeded = true;
        }
        return false;
    });
    return best;
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PointConcept OtherPoint>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherPoint& point) const {
    if (intersects(point)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(point);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, SegmentConcept OtherSegment>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedSegmentConcept OtherOrientedSegment>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherOrientedSegment& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, LineConcept OtherLine>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, OrientedLineConcept OtherOrientedLine>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherOrientedLine& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RayConcept OtherRay>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherRay& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherHalfplane& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, RectangleConcept OtherRectangle>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherRectangle& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, TriangleConcept OtherTriangle>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherTriangle& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, ConvexConcept OtherConvex>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherConvex& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PolygonConcept OtherPolygon>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherPolygon& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
// Against another region the scan over this region's ring edges is still
// enough: the other region's own ring edges are what its edge-to-shape distance
// minimizes over, so the pair is measured boundary against boundary either way.
template <class PointType, class LabelType>
template <class ResultNumber, PolygonWithHolesConcept OtherRegion>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherRegion& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, MonotoneChainConcept OtherChain>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherChain& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PolylineConcept OtherPolyline>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherPolyline& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
// The operand may be unbounded, which changes nothing: the region is bounded and
// closed, so a disjoint pair still realizes its distance on ∂A, and each ring
// edge hands the query on to the other region's own edge scan.
template <class PointType, class LabelType>
template <class ResultNumber, HalfplaneIntersectionConcept OtherIntersection>
constexpr auto PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherIntersection& other) const {
    if (intersects(other)) {
        return ResultNumber{};
    }
    return this->template edgeMinSquaredDistance<ResultNumber>(other);
}
 
// The nearest point of a disjoint disk lies on its circle, so the gap is the
// distance to its center less its radius — generally irrational, hence the
// floating result every other distance to a Disk returns.
template <class PointType, class LabelType>
template <class ResultNumber, DiskConcept OtherDisk>
detail::floating_result_t<ResultNumber>
PolygonWithHoles<PointType, LabelType>::squaredDistance(const OtherDisk& other) const {
    using Float = detail::floating_result_t<ResultNumber>;
    if (other.isDegenerate()) {
        // A degenerate disk has no centre or radius to measure from: it is the
        // point a() when its radius is zero, and undefined otherwise.
        return this->template squaredDistance<Float>(other.a());
    }
    if (intersects(other)) {
        return Float{0};
    }
    return detail::diskExteriorSquaredDistance<Float>(other, *this);
}
 
 
// ---------------------------------------------------------------------------
// PolygonSet
//
// The distance to a union is the minimum of the distances, with no caveat and
// no operand it fails for. An empty set has no component to measure from and
// reports zero.
 
template <class PointType, class LabelType>
template <class ResultNumber, detail::SetOperandConcept OtherShape>
auto PolygonSet<PointType, LabelType>::squaredDistance(const OtherShape& other) const {
    return minOverComponents([&other](const ComponentType& component) {
        return component.template squaredDistance<ResultNumber>(other);
    });
}
 
template <class PointType, class LabelType>
template <class ResultNumber, PolygonSetConcept OtherSet>
auto PolygonSet<PointType, LabelType>::squaredDistance(const OtherSet& other) const {
    ResultNumber best{};
    bool seeded = false;
    for (const auto& component : other) {
        const ResultNumber current = this->template squaredDistance<ResultNumber>(component);
        if (!seeded || current < best) {
            best = current;
            seeded = true;
        }
    }
    return best;
}
 
}  // namespace pgl