minkowski.hpp

Go to the documentation of this file.

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

 
#include <algorithm>
#include <cstddef>
#include <optional>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
 
namespace pgl {
 
namespace detail {

template <class T> struct minkowskiPointOf { using type = typename T::PointType; };
template <class Number, class Label> struct minkowskiPointOf<Point<Number, Label>> { using type = Point<Number, Label>; };
template <class PointType> struct minkowskiPointOf<Shape<PointType>> { using type = PointType; };
template <class T> using minkowskiPointOf_t = typename minkowskiPointOf<std::remove_cvref_t<T>>::type;

template <class A, class B>
using minkowskiPoint_t = Point<
    std::common_type_t<typename minkowskiPointOf_t<A>::NumberType,
                       typename minkowskiPointOf_t<B>::NumberType>,
    typename minkowskiPointOf_t<A>::LabelType>;

template <class A, class B>
using minkowskiRegionPoint_t = Point<
    std::conditional_t<is_halfplane_intersection_v<A> || is_halfplane_intersection_v<B>,
                       division_result_t<typename minkowskiPoint_t<A, B>::NumberType>,
                       typename minkowskiPoint_t<A, B>::NumberType>,
    typename minkowskiPoint_t<A, B>::LabelType>;
 
template <class A, class B>
    requires MinkowskiSummableConcept<A, B>
constexpr auto minkowskiSumOf(const A& a, const B& b);

template <class ShapeT, class TranslationNumber, class TranslationLabel>
constexpr auto minkowskiTranslated(const ShapeT& shape,
                                   const Point<TranslationNumber, TranslationLabel>& translation) {
    using Translation = Point<TranslationNumber, TranslationLabel>;
    using ResultPoint = minkowskiPoint_t<ShapeT, Translation>;
    using ResultNumber = typename ResultPoint::NumberType;
 
    // Moves one vertex, promoting the coordinates exactly as Point + Point does.
    const auto moved = [&translation](const auto& vertex) {
        return ResultPoint(
            detail::asNumber<ResultNumber>(vertex.x()) + detail::asNumber<ResultNumber>(translation.x()),
            detail::asNumber<ResultNumber>(vertex.y()) + detail::asNumber<ResultNumber>(translation.y()));
    };
 
    if constexpr (is_point_v<ShapeT>) {
        return moved(shape);
    } else if constexpr (is_empty_shape_v<ShapeT>) {
        // The empty set has no point to move; only the point type promotes.
        return EmptyShape<ResultPoint>{};
    } else if constexpr (is_shape_v<ShapeT>) {
        return std::visit(
            [&translation](const auto& alternative) {
                return Shape<ResultPoint>(minkowskiTranslated(alternative, translation));
            },
            shape.variant());
    } else {
        using Label = typename ShapeT::LabelType;
 
        if constexpr (is_segment_v<ShapeT>) {
            return Segment<ResultPoint, Label>(moved(shape.min()), moved(shape.max()));
        } else if constexpr (is_oriented_segment_v<ShapeT>) {
            return OrientedSegment<ResultPoint, Label>(moved(shape.source()), moved(shape.target()));
        } else if constexpr (is_line_v<ShapeT>) {
            return Line<ResultPoint, Label>(moved(shape.min()), moved(shape.max()));
        } else if constexpr (is_oriented_line_v<ShapeT>) {
            return OrientedLine<ResultPoint, Label>(moved(shape.source()), moved(shape.target()));
        } else if constexpr (is_ray_v<ShapeT>) {
            return Ray<ResultPoint, Label>(moved(shape.source()), moved(shape.target()));
        } else if constexpr (is_halfplane_v<ShapeT>) {
            return Halfplane<ResultPoint, Label>(moved(shape.source()), moved(shape.target()));
        } else if constexpr (is_rectangle_v<ShapeT>) {
            // Translation preserves the corner order, so the result needs no
            // normalizing -- and an empty rectangle stays empty.
            return Rectangle<ResultPoint, Label>(moved(shape.min()), moved(shape.max()), true);
        } else if constexpr (is_triangle_v<ShapeT>) {
            return Triangle<ResultPoint, Label>(moved(shape.a()), moved(shape.b()), moved(shape.c()));
        } else if constexpr (is_disk_v<ShapeT>) {
            return Disk<ResultPoint, Label>(moved(shape.a()), moved(shape.b()), moved(shape.c()));
        } else {
            // Vector-backed shapes translate in place: their own operator+=
            // keeps whatever cached state and canonical order the shape
            // maintains, which rebuilding from vertices would throw away.
            const auto translate = [&translation](auto result) {
                result += translation;
                if constexpr (has_label_v<Label>) {
                    result.label() = Label{};
                }
                return result;
            };
 
            if constexpr (is_convex_v<ShapeT>) {
                return translate(Convex<ResultPoint, Label>(shape));
            } else if constexpr (is_polygon_v<ShapeT>) {
                return translate(Polygon<ResultPoint, Label>(shape));
            } else if constexpr (is_monotone_chain_v<ShapeT>) {
                return translate(MonotoneChain<ResultPoint, Label>(shape));
            } else if constexpr (is_polyline_v<ShapeT>) {
                return translate(Polyline<ResultPoint, Label>(shape));
            } else if constexpr (is_polygon_with_holes_v<ShapeT>) {
                return translate(PolygonWithHoles<ResultPoint, Label>(shape));
            } else if constexpr (is_polygon_set_v<ShapeT>) {
                return translate(PolygonSet<ResultPoint, Label>(shape));
            } else {
                static_assert(is_halfplane_intersection_v<ShapeT>,
                              "minkowskiTranslated has no branch for this shape kind: every "
                              "shape is closed under translation, so a new one needs one here");
                return translate(HalfplaneIntersection<ResultPoint, Label>(shape));
            }
        }
    }
}

template <class ResultPoint, class ShapeT>
constexpr std::vector<ResultPoint> minkowskiVertices(const ShapeT& shape) {
    using ResultNumber = typename ResultPoint::NumberType;
 
    std::vector<ResultPoint> vertices;
    const auto append = [&vertices](const auto& vertex) {
        vertices.emplace_back(detail::asNumber<ResultNumber>(vertex.x()),
                              detail::asNumber<ResultNumber>(vertex.y()));
    };
 
    if constexpr (is_convex_v<ShapeT>) {
        // Convex already stores its hull counterclockwise with the collinear
        // vertices pruned, so a scan here would be pure overhead.
        vertices.reserve(shape.size());
        for (const auto& vertex : shape) {
            append(vertex);
        }
        return vertices;
    } else {
        // At most four vertices, and any of them may be degenerate (a flat
        // triangle, a zero-width rectangle, a collapsed segment), so the scan
        // both orders them and drops what collapsed. An empty shape has no
        // vertices to read, and its caller absorbs the empty list.
        if (coversNoPoint(shape)) {
            return vertices;
        }
        for (const auto& vertex : shape.vertices()) {
            append(vertex);
        }
        return grahamScan(vertices);
    }
}

template <class HalfplaneT, class ShapeT>
constexpr auto minkowskiHalfplaneSum(const HalfplaneT& halfplane, const ShapeT& shape) {
    using ResultPoint = minkowskiPoint_t<HalfplaneT, ShapeT>;
    using ResultNumber = typename ResultPoint::NumberType;
    using ResultHalfplane = Halfplane<ResultPoint, typename HalfplaneT::LabelType>;
 
    const auto& source = halfplane.source();
    const auto& target = halfplane.target();
    const ResultNumber dx = detail::asNumber<ResultNumber>(target.x()) - detail::asNumber<ResultNumber>(source.x());
    const ResultNumber dy = detail::asNumber<ResultNumber>(target.y()) - detail::asNumber<ResultNumber>(source.y());
 
    bool found = false;
    ResultPoint support(ResultNumber{}, ResultNumber{});
    ResultNumber best{};
    for (const auto& vertex : shape.vertices()) {
        const ResultNumber x = static_cast<ResultNumber>(vertex.x());
        const ResultNumber y = static_cast<ResultNumber>(vertex.y());
        const ResultNumber side = dx * y - dy * x;  // cross(d, vertex)
        if (!found || side < best) {
            found = true;
            best = side;
            support = ResultPoint(x, y);
        }
    }
    if (!found) {
        return ResultHalfplane(ResultPoint(detail::asNumber<ResultNumber>(source.x()),
                                           detail::asNumber<ResultNumber>(source.y())),
                               ResultPoint(detail::asNumber<ResultNumber>(target.x()),
                                           detail::asNumber<ResultNumber>(target.y())));
    }
    const auto moved = [&support](const auto& point) {
        return ResultPoint(detail::asNumber<ResultNumber>(point.x()) + support.x(),
                           detail::asNumber<ResultNumber>(point.y()) + support.y());
    };
    return ResultHalfplane(moved(source), moved(target));
}

template <class ResultPoint>
struct MinkowskiPolyhedron {
    std::vector<ResultPoint> directions;
    std::vector<ResultPoint> anchors;
    std::vector<ResultPoint> recessions;
    bool empty = false;
};

template <class ResultPoint>
constexpr typename ResultPoint::NumberType minkowskiCross(const ResultPoint& u,
                                                          const ResultPoint& v) {
    return u.x() * v.y() - u.y() * v.x();
}

template <class ResultPoint, class ShapeT>
constexpr MinkowskiPolyhedron<ResultPoint> minkowskiPolyhedronOf(const ShapeT& shape) {
    using ResultNumber = typename ResultPoint::NumberType;
 
    MinkowskiPolyhedron<ResultPoint> polyhedron;
    const ResultNumber zero{};
    const auto cast = [](const auto& point) {
        return ResultPoint(detail::asNumber<ResultNumber>(point.x()),
                           detail::asNumber<ResultNumber>(point.y()));
    };
    const auto reversed = [](const ResultPoint& vector) {
        return ResultPoint(-vector.x(), -vector.y());
    };
 
    if constexpr (is_point_v<ShapeT>) {
        // A single point bounds every direction and recedes in none. The
        // dispatcher translates instead of coming here, so this is only for
        // completeness.
        polyhedron.anchors.push_back(cast(shape));
    } else if constexpr (is_line_v<ShapeT> || is_oriented_line_v<ShapeT> || is_ray_v<ShapeT> ||
                         is_halfplane_v<ShapeT>) {
        const ResultPoint source = cast(shape[0]);
        const ResultPoint target = cast(shape[1]);
        const ResultPoint forward(target.x() - source.x(), target.y() - source.y());
        polyhedron.anchors.push_back(source);
        if constexpr (is_halfplane_v<ShapeT>) {
            // One constraint, and a recession cone that is the half-plane
            // itself: it runs off along its boundary both ways and inward along
            // the normal, so `cross(d, ·) >= 0` on those three holds exactly for
            // the positive multiples of the boundary direction.
            polyhedron.directions.push_back(forward);
            polyhedron.recessions.push_back(forward);
            polyhedron.recessions.push_back(reversed(forward));
            polyhedron.recessions.emplace_back(-forward.y(), forward.x());
        } else if constexpr (is_ray_v<ShapeT>) {
            // Two constraints pin the ray to its supporting line and one caps it
            // at the source; the cap's boundary direction is the source-ward
            // normal, since `cross((dy, -dx), p)` is `dot(d, p)`.
            polyhedron.directions.push_back(forward);
            polyhedron.directions.push_back(reversed(forward));
            polyhedron.directions.emplace_back(forward.y(), -forward.x());
            polyhedron.recessions.push_back(forward);
        } else {
            // A line is capped nowhere and recedes both ways, so a direction
            // bounds it only when it is parallel to it.
            polyhedron.directions.push_back(forward);
            polyhedron.directions.push_back(reversed(forward));
            polyhedron.recessions.push_back(forward);
            polyhedron.recessions.push_back(reversed(forward));
        }
    } else if constexpr (is_halfplane_intersection_v<ShapeT>) {
        if (shape.empty()) {
            polyhedron.empty = true;
            return polyhedron;
        }
        const std::size_t n = shape.size();
        if (n == 0) {
            // The whole plane constrains nothing and recedes everywhere.
            const ResultNumber one = static_cast<ResultNumber>(1);
            polyhedron.anchors.emplace_back(zero, zero);
            polyhedron.recessions.emplace_back(one, zero);
            polyhedron.recessions.emplace_back(-one, zero);
            polyhedron.recessions.emplace_back(zero, one);
            polyhedron.recessions.emplace_back(zero, -one);
            return polyhedron;
        }
 
        std::vector<ResultPoint> sources;
        std::vector<ResultNumber> offsets;
        sources.reserve(n);
        offsets.reserve(n);
        polyhedron.directions.reserve(n);
        for (const auto& halfplane : shape) {
            const ResultPoint source = cast(halfplane.source());
            const ResultPoint target = cast(halfplane.target());
            const ResultPoint direction(target.x() - source.x(), target.y() - source.y());
            sources.push_back(source);
            polyhedron.directions.push_back(direction);
            offsets.push_back(minkowskiCross(direction, source));
        }
        if (n == 1) {
            // One stored constraint *is* a half-plane, and its recession cone is
            // two-dimensional, which no pair of edge directions would report.
            const ResultPoint& forward = polyhedron.directions.front();
            polyhedron.anchors.push_back(sources.front());
            polyhedron.recessions.push_back(forward);
            polyhedron.recessions.push_back(reversed(forward));
            polyhedron.recessions.emplace_back(-forward.y(), forward.x());
            return polyhedron;
        }
 
        // The stored half-planes are sorted by boundary direction, and a
        // cyclically consecutive pair meets in a vertex exactly when it turns
        // left. Where it does not, the region's boundary runs off to infinity
        // between them: forward along the first edge, backward along the second.
        for (std::size_t i = 0; i < n; ++i) {
            const std::size_t next = (i + 1) % n;
            const ResultPoint& here = polyhedron.directions[i];
            const ResultPoint& there = polyhedron.directions[next];
            const ResultNumber turn = minkowskiCross(here, there);
            if (turn > zero) {
                // Cramer's rule on the two boundary lines.
                polyhedron.anchors.emplace_back(
                    (offsets[i] * there.x() - offsets[next] * here.x()) / turn,
                    (offsets[i] * there.y() - offsets[next] * here.y()) / turn);
            } else {
                polyhedron.recessions.push_back(here);
                polyhedron.recessions.push_back(reversed(there));
            }
        }
        if (polyhedron.anchors.empty()) {
            // No vertex at all, with two or more constraints: the region is a
            // slab or a line, and each stored boundary line lies in it whole.
            polyhedron.anchors = std::move(sources);
        }
    } else {
        // Bounded and convex: the hull vertices, counterclockwise and with the
        // collinear ones dropped, are both the anchors and the edge directions.
        std::vector<ResultPoint> vertices = minkowskiVertices<ResultPoint>(shape);
        if (vertices.empty()) {
            polyhedron.empty = true;
            return polyhedron;
        }
        if (vertices.size() >= 2) {
            // A shape that collapsed to a segment contributes both directions,
            // as the two edges of the degenerate polygon it is.
            polyhedron.directions.reserve(vertices.size());
            for (std::size_t i = 0; i < vertices.size(); ++i) {
                const ResultPoint& from = vertices[i];
                const ResultPoint& to = vertices[(i + 1) % vertices.size()];
                polyhedron.directions.emplace_back(to.x() - from.x(), to.y() - from.y());
            }
        }
        polyhedron.anchors = std::move(vertices);
    }
    return polyhedron;
}

template <class ResultPoint>
constexpr std::optional<ResultPoint> minkowskiInfimumPoint(
    const MinkowskiPolyhedron<ResultPoint>& polyhedron, const ResultPoint& direction) {
    using ResultNumber = typename ResultPoint::NumberType;
 
    for (const auto& recession : polyhedron.recessions) {
        if (crossSign(direction, recession) < 0) {
            return std::nullopt;
        }
    }
    const ResultPoint* best = nullptr;
    ResultNumber least{};
    for (const auto& anchor : polyhedron.anchors) {
        const ResultNumber value = minkowskiCross(direction, anchor);
        if (best == nullptr || value < least) {
            best = &anchor;
            least = value;
        }
    }
    if (best == nullptr) {
        return std::nullopt;
    }
    return *best;
}

template <class A, class B>
constexpr auto minkowskiPolyhedralSum(const A& a, const B& b) {
    using ResultPoint = minkowskiRegionPoint_t<A, B>;
    using Region = HalfplaneIntersection<ResultPoint>;
 
    const MinkowskiPolyhedron<ResultPoint> left = minkowskiPolyhedronOf<ResultPoint>(a);
    const MinkowskiPolyhedron<ResultPoint> right = minkowskiPolyhedronOf<ResultPoint>(b);
    if (left.empty || right.empty) {
        // Nothing to add to anything: the empty convex polygon is the region
        // constructor that spells the empty region.
        return Region(Convex<ResultPoint>());
    }
 
    Region region;
    const auto clamp = [&region, &left, &right](const ResultPoint& direction) {
        const std::optional<ResultPoint> here = minkowskiInfimumPoint(left, direction);
        if (!here) {
            return;
        }
        const std::optional<ResultPoint> there = minkowskiInfimumPoint(right, direction);
        if (!there) {
            return;
        }
        const ResultPoint base(here->x() + there->x(), here->y() + there->y());
        region.insert(Halfplane<ResultPoint>(
            base, ResultPoint(base.x() + direction.x(), base.y() + direction.y())));
    };
    for (const auto& direction : left.directions) {
        clamp(direction);
    }
    for (const auto& direction : right.directions) {
        clamp(direction);
    }
    return region;
}

template <class PointType>
constexpr int minkowskiDirectionOrder(const PointType& u, const PointType& v) {
    using Number = typename PointType::NumberType;
 
    const auto half = [](const PointType& direction) {
        const Number origin{};
        return (direction.y() < origin || (direction.y() == origin && direction.x() < origin)) ? 1
                                                                                               : 0;
    };
 
    const int halfU = half(u);
    const int halfV = half(v);
    if (halfU != halfV) {
        return halfU < halfV ? -1 : 1;
    }
 
    const auto turn = crossSign(u, v);
    if (turn > 0) {
        return -1;  // v lies counterclockwise from u, so u has the smaller angle
    }
    if (turn < 0) {
        return 1;
    }
    return 0;
}

template <class A, class B>
constexpr auto minkowskiConvexSum(const A& a, const B& b) {
    using ResultPoint = minkowskiPoint_t<A, B>;
    using ResultNumber = typename ResultPoint::NumberType;
 
    std::vector<ResultPoint> left = minkowskiVertices<ResultPoint>(a);
    std::vector<ResultPoint> right = minkowskiVertices<ResultPoint>(b);
    if (left.empty() || right.empty()) {
        return Convex<ResultPoint>();  // an empty convex polygon absorbs
    }
 
    // Start each walk at its bottom-most (then leftmost) vertex: the outgoing
    // edge there has an angle in [0, π), so the edges that follow it are sorted
    // by angle over the whole [0, 2π) turn.
    const auto rotateToLowest = [](std::vector<ResultPoint>& vertices) {
        auto lowest = vertices.begin();
        for (auto it = vertices.begin() + 1; it != vertices.end(); ++it) {
            if (it->y() < lowest->y() || (it->y() == lowest->y() && it->x() < lowest->x())) {
                lowest = it;
            }
        }
        std::rotate(vertices.begin(), lowest, vertices.end());
    };
    rotateToLowest(left);
    rotateToLowest(right);
 
    // A one-vertex operand contributes no edge, which is what makes the merge
    // below degrade gracefully into a plain translation.
    const auto edgesOf = [](const std::vector<ResultPoint>& vertices) {
        std::vector<ResultPoint> edges;
        if (vertices.size() >= 2) {
            edges.reserve(vertices.size());
            for (std::size_t i = 0; i < vertices.size(); ++i) {
                const ResultPoint& from = vertices[i];
                const ResultPoint& to = vertices[(i + 1) % vertices.size()];
                edges.emplace_back(to.x() - from.x(), to.y() - from.y());
            }
        }
        return edges;
    };
    const std::vector<ResultPoint> leftEdges = edgesOf(left);
    const std::vector<ResultPoint> rightEdges = edgesOf(right);
 
    std::vector<ResultPoint> boundary;
    boundary.reserve(leftEdges.size() + rightEdges.size() + 1);
    ResultPoint current(left.front().x() + right.front().x(),
                        left.front().y() + right.front().y());
    boundary.push_back(current);
 
    std::size_t i = 0;
    std::size_t j = 0;
    while (i < leftEdges.size() || j < rightEdges.size()) {
        ResultPoint step(ResultNumber{}, ResultNumber{});
        if (j == rightEdges.size()) {
            step = leftEdges[i++];
        } else if (i == leftEdges.size()) {
            step = rightEdges[j++];
        } else {
            const int order = minkowskiDirectionOrder(leftEdges[i], rightEdges[j]);
            if (order < 0) {
                step = leftEdges[i++];
            } else if (order > 0) {
                step = rightEdges[j++];
            } else {
                // Equal directions fuse, which is what keeps the result free of
                // collinear vertices.
                step = ResultPoint(leftEdges[i].x() + rightEdges[j].x(),
                                   leftEdges[i].y() + rightEdges[j].y());
                ++i;
                ++j;
            }
        }
        current = ResultPoint(current.x() + step.x(), current.y() + step.y());
        boundary.push_back(current);
    }
    if (boundary.size() > 1) {
        boundary.pop_back();  // the walk closed back onto its first vertex
    }
 
    // Convex stores its hull from the lexicographically smallest vertex; the
    // walk above starts from the bottom-most one instead.
    std::rotate(boundary.begin(), std::min_element(boundary.begin(), boundary.end()),
                boundary.end());
    return Convex<ResultPoint>(std::move(boundary), true);
}

template <class A, class B>
    requires MinkowskiSummableConcept<A, B>
constexpr auto minkowskiSumOf(const A& a, const B& b) {
    using ResultPoint = minkowskiPoint_t<A, B>;
 
    if constexpr (is_shape_v<A> || is_shape_v<B>) {
        using ResultShape = Shape<ResultPoint>;
        // Only the pair of stored alternatives decides whether the sum exists,
        // and that is not known until run time.
        // The constructibility is part of the test, not a consequence of it: a
        // pair whose sum is a HalfplaneIntersection over rational coordinates
        // has an answer, but not one a wrapper over integral points can hold.
        const auto sum = [](const auto& left, const auto& right) -> ResultShape {
            if constexpr (requires { ResultShape(minkowskiSumOf(left, right)); }) {
                return ResultShape(minkowskiSumOf(left, right));
            } else {
                throw std::logic_error(
                    "Shape::minkowskiSum is not defined for this pair of alternatives, or its "
                    "result does not fit the wrapper's point type");
            }
        };
        if constexpr (is_shape_v<A> && is_shape_v<B>) {
            return std::visit(sum, a.variant(), b.variant());
        } else if constexpr (is_shape_v<A>) {
            return std::visit([&b, &sum](const auto& left) { return sum(left, b); }, a.variant());
        } else {
            return std::visit([&a, &sum](const auto& right) { return sum(a, right); }, b.variant());
        }
    } else if constexpr (is_empty_shape_v<A> || is_empty_shape_v<B>) {
        return EmptyShape<ResultPoint>{};
    } else if constexpr (is_point_v<B>) {
        return minkowskiTranslated(a, b);
    } else if constexpr (is_point_v<A>) {
        return minkowskiTranslated(b, a);
    } else if constexpr (is_rectangle_v<A> && is_rectangle_v<B>) {
        // Two axis-aligned rectangles are the one non-trivial pair closed under
        // the sum: opposite corners simply add, minima to minima and maxima to
        // maxima. Sweeping anything over the empty set covers nothing.
        if (a.empty() || b.empty()) {
            return Rectangle<ResultPoint>();
        }
        return Rectangle<ResultPoint>(minkowskiSumOf(a.min(), b.min()),
                                      minkowskiSumOf(a.max(), b.max()), true);
    } else if constexpr (is_halfplane_v<A> && !UnboundedConvexConcept<B>) {
        // A half-plane absorbs anything bounded, convex or not, and stays one.
        return minkowskiHalfplaneSum(a, b);
    } else if constexpr (is_halfplane_v<B> && !UnboundedConvexConcept<A>) {
        return minkowskiHalfplaneSum(b, a);
    } else if constexpr (UnboundedConvexConcept<A> || UnboundedConvexConcept<B>) {
        // One operand at least is stored as constraints rather than as
        // vertices, and so is their sum.
        return minkowskiPolyhedralSum(a, b);
    } else {
        return minkowskiConvexSum(a, b);
    }
}
 
}  // namespace detail
 
// -----------------------------------------------------------------------------
// Member entry points
//
// Each shape forwards to the same dispatcher; the pairs a shape actually
// accepts are fixed by MinkowskiSummableConcept on the declaration.
 
#define PGL_DEFINE_MINKOWSKI_SUM(SHAPE)                                                  \
    template <class PointType, class LabelType>                                          \
    template <class OtherShape>                                                          \
        requires MinkowskiSummableConcept<SHAPE<PointType, LabelType>, OtherShape>       \
    constexpr auto SHAPE<PointType, LabelType>::minkowskiSum(const OtherShape& other) const { \
        return detail::minkowskiSumOf(*this, other);                                     \
    }
 
PGL_DEFINE_MINKOWSKI_SUM(Segment)
PGL_DEFINE_MINKOWSKI_SUM(OrientedSegment)
PGL_DEFINE_MINKOWSKI_SUM(Line)
PGL_DEFINE_MINKOWSKI_SUM(OrientedLine)
PGL_DEFINE_MINKOWSKI_SUM(Ray)
PGL_DEFINE_MINKOWSKI_SUM(Halfplane)
PGL_DEFINE_MINKOWSKI_SUM(Rectangle)
PGL_DEFINE_MINKOWSKI_SUM(Triangle)
PGL_DEFINE_MINKOWSKI_SUM(Disk)
PGL_DEFINE_MINKOWSKI_SUM(Convex)
PGL_DEFINE_MINKOWSKI_SUM(Polygon)
PGL_DEFINE_MINKOWSKI_SUM(Polyline)
PGL_DEFINE_MINKOWSKI_SUM(PolygonWithHoles)
PGL_DEFINE_MINKOWSKI_SUM(PolygonSet)
PGL_DEFINE_MINKOWSKI_SUM(HalfplaneIntersection)
 
#undef PGL_DEFINE_MINKOWSKI_SUM
 
// The one curved pair with an answer: centres add, radii add. It carries a
// ResultNumber of its own because the radii are square roots of what a disk
// stores, so this sum leaves the lattice where every other single-shape sum
// stays on it.
template <class PointType_, class TLabel>
template <class ResultNumber, DiskConcept OtherDisk>
Disk<Point<ResultNumber, typename Disk<PointType_, TLabel>::PointLabelType>>
Disk<PointType_, TLabel>::minkowskiSum(const OtherDisk& other) const {
    using ResultPoint = Point<ResultNumber, PointLabelType>;
 
    const auto leftCenter = center<ResultNumber>();
    const auto rightCenter = other.template center<ResultNumber>();
    const ResultNumber radii =
        radius<ResultNumber>() + other.template radius<ResultNumber>();
    return Disk<ResultPoint>(ResultPoint(leftCenter.x() + rightCenter.x(),
                                         leftCenter.y() + rightCenter.y()),
                             radii);
}
 
// The other sum whose support point is not a vertex: a half-plane slides out by
// the disk's radius along its own unit normal, and both of those are square
// roots. Same ResultNumber convention as the disk pair.
template <class PointType_, class TLabel>
template <class ResultNumber, DiskConcept OtherDisk>
Halfplane<Point<ResultNumber, typename PointType_::LabelType>>
Halfplane<PointType_, TLabel>::minkowskiSum(const OtherDisk& other) const {
    using ResultPoint = Point<ResultNumber, typename PointType_::LabelType>;
 
    const ResultNumber dx = detail::asNumber<ResultNumber>(target().x()) - detail::asNumber<ResultNumber>(source().x());
    const ResultNumber dy = detail::asNumber<ResultNumber>(target().y()) - detail::asNumber<ResultNumber>(source().y());
    // Reported the way Disk::radius reports it: an exact result type has no
    // square root to offer, and says so rather than rounding silently.
    if constexpr (!requires(ResultNumber v) { std::sqrt(v); }) {
        throw std::runtime_error("std::sqrt is not available for the requested ResultNumber type");
    } else {
        const ResultNumber length = std::sqrt(dx * dx + dy * dy);
 
        // The outward normal of `{p : cross(d, p - s) >= 0}` is `(dy, -dx)/|d|`,
        // and the disk's support point in it is its centre moved r that way.
        const auto center = other.template center<ResultNumber>();
        const ResultNumber radius = other.template radius<ResultNumber>();
        const ResultNumber offsetX = center.x() + radius * dy / length;
        const ResultNumber offsetY = center.y() - radius * dx / length;
 
        const auto moved = [&offsetX, &offsetY](const auto& point) {
            return ResultPoint(detail::asNumber<ResultNumber>(point.x()) + offsetX,
                               detail::asNumber<ResultNumber>(point.y()) + offsetY);
        };
        return Halfplane<ResultPoint>(moved(source()), moved(target()));
    }
}
 
template <class PointType_, class TLabel>
template <class ResultNumber, HalfplaneConcept OtherHalfplane>
Halfplane<Point<ResultNumber, typename Disk<PointType_, TLabel>::PointLabelType>>
Disk<PointType_, TLabel>::minkowskiSum(const OtherHalfplane& other) const {
    return other.template minkowskiSum<ResultNumber>(*this);
}
 
template <class Number, class Label>
template <class OtherShape>
    requires MinkowskiSummableConcept<Point<Number, Label>, OtherShape>
constexpr auto Point<Number, Label>::minkowskiSum(const OtherShape& other) const {
    return detail::minkowskiSumOf(*this, other);
}
 
template <class PointType>
template <class OtherShape>
    requires MinkowskiSummableConcept<EmptyShape<PointType>, OtherShape>
constexpr auto EmptyShape<PointType>::minkowskiSum(const OtherShape& other) const {
    return detail::minkowskiSumOf(*this, other);
}
 
template <class PointType, class LabelType, class Storage>
template <class OtherShape>
    requires MinkowskiSummableConcept<MonotoneChain<PointType, LabelType, Storage>, OtherShape>
constexpr auto MonotoneChain<PointType, LabelType, Storage>::minkowskiSum(const OtherShape& other) const {
    return detail::minkowskiSumOf(*this, other);
}
 
template <class PointType>
template <class OtherShape>
    requires MinkowskiSummableConcept<Shape<PointType>, OtherShape>
constexpr auto Shape<PointType>::minkowskiSum(const OtherShape& other) const {
    return detail::minkowskiSumOf(*this, other);
}

template <class A, class B>
    requires MinkowskiSummableConcept<A, B>
[[nodiscard]] constexpr auto operator+(const A& a, const B& b) {
    return a.minkowskiSum(b);
}
 
}  // namespace pgl