bounding.hpp

Go to the documentation of this file.

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

 
#include "shape/rectangle.hpp"
 
namespace pgl {
 
// -----------------------------------------------------------------------------
// Point
 
template <class Number, class Label>
constexpr Rectangle<Point<Number, Label>> Point<Number, Label>::bbox() const {
    return Rectangle<Point<Number, Label>>(*this, *this, true);
}
 
template <class Number, class Label>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> Point<Number, Label>::fbox() const {
    return Rectangle<Point<ResultNumber>>(
        detail::lowerFloatingBound<ResultNumber>(x()),
        detail::lowerFloatingBound<ResultNumber>(y()),
        detail::upperFloatingBound<ResultNumber>(x()),
        detail::upperFloatingBound<ResultNumber>(y()),
        true);
}
 
template <class Number, class Label>
constexpr std::array<Point<Number, Label>, 1> Point<Number, Label>::vertices() const {
    return {*this};
}
 
template <class Number, class Label>
constexpr std::array<Segment<Point<Number, Label>>, 0> Point<Number, Label>::edges() const {
    return {};
}
 
template <class Number, class Label>
constexpr std::array<OrientedSegment<Point<Number, Label>>, 0> Point<Number, Label>::orientedEdges() const {
    return {};
}
 
// -----------------------------------------------------------------------------
// Segment
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber, class Value>
constexpr ResultNumber Segment<PointType, LabelType>::lowerCoordinateBound(const Value& value) {
    if constexpr (requires { value.template lowerBound<ResultNumber>(); }) {
        return value.template lowerBound<ResultNumber>();
    } else {
        return static_cast<ResultNumber>(value);
    }
}
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber, class Value>
constexpr ResultNumber Segment<PointType, LabelType>::upperCoordinateBound(const Value& value) {
    if constexpr (requires { value.template upperBound<ResultNumber>(); }) {
        return value.template upperBound<ResultNumber>();
    } else {
        return static_cast<ResultNumber>(value);
    }
}
 
template <class PointType, class LabelType>
constexpr Rectangle<PointType> Segment<PointType, LabelType>::bbox() const {
    if (min().y() < max().y()) {
        return Rectangle<PointType>(min().x(), min().y(), max().x(), max().y(), true);
    }
    return Rectangle<PointType>(min().x(), max().y(), max().x(), min().y(), true);
}
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> Segment<PointType, LabelType>::fbox() const {
    ResultNumber xmin = lowerCoordinateBound<ResultNumber>(min().x());
    ResultNumber ymin = lowerCoordinateBound<ResultNumber>(min().y());
    ResultNumber xmax = upperCoordinateBound<ResultNumber>(max().x());
    ResultNumber ymax = upperCoordinateBound<ResultNumber>(max().y());
    if (ymin > ymax) {
        std::swap(ymin, ymax);
    }
    return Rectangle<Point<ResultNumber>>(xmin, ymin, xmax, ymax, true);
}
 
template <class PointType, class LabelType>
constexpr std::array<PointType, 2> Segment<PointType, LabelType>::vertices() const {
    return {min(), max()};
}
 
template <class PointType, class LabelType>
constexpr std::array<Segment<PointType, LabelType>, 1> Segment<PointType, LabelType>::edges() const {
    return {*this};
}
 
template <class PointType, class LabelType>
constexpr std::array<OrientedSegment<PointType>, 1> Segment<PointType, LabelType>::orientedEdges() const {
    return {OrientedSegment<PointType>(min(), max())};
}
 
// -----------------------------------------------------------------------------
// OrientedSegment
 
template <class PointType, class LabelType>
constexpr Rectangle<PointType> OrientedSegment<PointType, LabelType>::bbox() const {
    return static_cast<Segment<PointType>>(*this).bbox();
}
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> OrientedSegment<PointType, LabelType>::fbox() const {
    return static_cast<Segment<PointType>>(*this).template fbox<ResultNumber>();
}
 
template <class PointType, class LabelType>
constexpr std::array<PointType, 2> OrientedSegment<PointType, LabelType>::vertices() const {
    return {source(), target()};
}
 
template <class PointType, class LabelType>
constexpr std::array<Segment<PointType>, 1> OrientedSegment<PointType, LabelType>::edges() const {
    return {static_cast<Segment<PointType>>(*this)};
}
 
template <class PointType, class LabelType>
constexpr std::array<OrientedSegment<PointType, LabelType>, 1> OrientedSegment<PointType, LabelType>::orientedEdges() const {
    return {*this};
}
 
// -----------------------------------------------------------------------------
// Rectangle
 
template <class PointType, class LabelType>
constexpr Rectangle<PointType, LabelType> Rectangle<PointType, LabelType>::bbox() const {
    return *this;
}
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> Rectangle<PointType, LabelType>::fbox() const {
    return Rectangle<Point<ResultNumber>>(
        detail::lowerFloatingBound<ResultNumber>(min().x()),
        detail::lowerFloatingBound<ResultNumber>(min().y()),
        detail::upperFloatingBound<ResultNumber>(max().x()),
        detail::upperFloatingBound<ResultNumber>(max().y()),
        true);
}
 
template <class PointType, class LabelType>
constexpr typename Rectangle<PointType, LabelType>::PointType Rectangle<PointType, LabelType>::bottomRight() const {
    return makeCorner(max().x(), min().y());
}
 
template <class PointType, class LabelType>
constexpr typename Rectangle<PointType, LabelType>::PointType Rectangle<PointType, LabelType>::topLeft() const {
    return makeCorner(min().x(), max().y());
}
 
template <class PointType, class LabelType>
template <bool Oriented>
constexpr typename Rectangle<PointType, LabelType>::template BoundaryType<Oriented> Rectangle<PointType, LabelType>::boundaryAt(std::size_t index) const {
    assert(index < size());
 
    const auto bottom_left = min();
    const auto bottom_right = bottomRight();
    const auto top_right = max();
    const auto top_left = topLeft();
 
    switch (index) {
    case 0:
        return BoundaryType<Oriented>(bottom_left, bottom_right);
    case 1:
        return BoundaryType<Oriented>(bottom_right, top_right);
    case 2:
        return BoundaryType<Oriented>(top_right, top_left);
    default:
        return BoundaryType<Oriented>(top_left, bottom_left);
    }
}
 
template <class PointType, class LabelType>
constexpr std::array<typename Rectangle<PointType, LabelType>::PointType, 4> Rectangle<PointType, LabelType>::vertices() const {
    assert(!empty());
    return {
        min(),
        bottomRight(),
        max(),
        topLeft(),
    };
}
 
template <class PointType, class LabelType>
constexpr std::array<Segment<PointType>, 4> Rectangle<PointType, LabelType>::edges() const {
    assert(!empty());
    return {
        boundaryAt<false>(0),
        boundaryAt<false>(1),
        boundaryAt<false>(2),
        boundaryAt<false>(3),
    };
}
 
template <class PointType, class LabelType>
constexpr std::array<OrientedSegment<PointType>, 4> Rectangle<PointType, LabelType>::orientedEdges() const {
    assert(!empty());
    return {
        boundaryAt<true>(0),
        boundaryAt<true>(1),
        boundaryAt<true>(2),
        boundaryAt<true>(3),
    };
}
 
// -----------------------------------------------------------------------------
// Triangle
 
template <class PointType, class LabelType>
constexpr Rectangle<PointType> Triangle<PointType, LabelType>::bbox() const {
    return Rectangle<PointType>(points_);
}
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> Triangle<PointType, LabelType>::fbox() const {
    return bbox().template fbox<ResultNumber>();
}
 
template <class PointType, class LabelType>
constexpr std::array<PointType, 3> Triangle<PointType, LabelType>::vertices() const {
    return points_;
}
 
template <class PointType, class LabelType>
constexpr std::array<Segment<PointType>, 3> Triangle<PointType, LabelType>::edges() const {
    return {
        Segment<PointType>(a(), b()),
        Segment<PointType>(b(), c()),
        Segment<PointType>(c(), a()),
    };
}
 
template <class PointType, class LabelType>
constexpr std::array<OrientedSegment<PointType>, 3> Triangle<PointType, LabelType>::orientedEdges() const {
    return {
        OrientedSegment<PointType>(a(), b()),
        OrientedSegment<PointType>(b(), c()),
        OrientedSegment<PointType>(c(), a()),
    };
}
 
// template <class PointType>
// template <std::ranges::input_range Range>
//     requires std::constructible_from<PointType, std::ranges::range_reference_t<Range>>
// constexpr void Rectangle<PointType, LabelType>::assignBoundingBox(Range&& points) {
//     auto iterator = std::ranges::begin(points);
//     const auto sentinel = std::ranges::end(points);
//
//     if (iterator == sentinel) {
//         throw std::invalid_argument("Rectangle bounding box requires at least one point");
//     }
//
//     PointType min_corner(*iterator);
//     PointType max_corner(*iterator);
//     std::optional<PointType> exact_min_corner(min_corner);
//     std::optional<PointType> exact_max_corner(max_corner);
//     ++iterator;
//
//     for (; iterator != sentinel; ++iterator) {
//         const PointType point(*iterator);
//         const PointType previous_min = min_corner;
//         const PointType previous_max = max_corner;
//
//         if (point.x() < min_corner.x()) {
//             min_corner = makeCorner(point.x(), min_corner.y());
//         }
//         if (point.y() < min_corner.y()) {
//             min_corner = makeCorner(min_corner.x(), point.y());
//         }
//         if (max_corner.x() < point.x()) {
//             max_corner = makeCorner(point.x(), max_corner.y());
//         }
//         if (max_corner.y() < point.y()) {
//             max_corner = makeCorner(max_corner.x(), point.y());
//         }
//
//         if (point.x() == min_corner.x() && point.y() == min_corner.y()) {
//             exact_min_corner = point;
//         } else if (previous_min != min_corner) {
//             exact_min_corner.reset();
//         }
//
//         if (point.x() == max_corner.x() && point.y() == max_corner.y()) {
//             exact_max_corner = point;
//         } else if (previous_max != max_corner) {
//             exact_max_corner.reset();
//         }
//     }
//
//     points_[0] = exact_min_corner.has_value() ? std::move(*exact_min_corner) : std::move(min_corner);
//     points_[1] = exact_max_corner.has_value() ? std::move(*exact_max_corner) : std::move(max_corner);
// }
 
template <class PointType, class LabelType>
template <PointConcept OtherPoint>
constexpr void Rectangle<PointType, LabelType>::insert(const OtherPoint& point) {
    const NumberType x = static_cast<NumberType>(point.x());
    const NumberType y = static_cast<NumberType>(point.y());
 
    if (empty()) {
        // The empty rectangle bounds nothing, so it cannot be grown: it becomes
        // the inserted point outright.
        points_[0] = points_[1] = makeCorner(x, y);
        return;
    }
 
    const PointType old_min = min();
    const PointType old_max = max();
 
    NumberType min_x = old_min.x();
    NumberType min_y = old_min.y();
    NumberType max_x = old_max.x();
    NumberType max_y = old_max.y();
 
    if (x < min_x) {
        min_x = x;
    }
    if (y < min_y) {
        min_y = y;
    }
    if (max_x < x) {
        max_x = x;
    }
    if (max_y < y) {
        max_y = y;
    }
 
    points_[0] = (min_x == old_min.x() && min_y == old_min.y())
        ? old_min
        : makeCorner(min_x, min_y);
    points_[1] = (max_x == old_max.x() && max_y == old_max.y())
        ? old_max
        : makeCorner(max_x, max_y);
}
 
template <class PointType, class LabelType>
template <RectangleConcept OtherRectangle>
constexpr void Rectangle<PointType, LabelType>::insert(const OtherRectangle& other) {
    if (other.empty()) {
        // Its corners are inverted placeholders, not points to enclose.
        return;
    }
    insert(other.min());
    insert(other.max());
}
 
template <class PointType, class LabelType>
template <class TShape>
    requires(!detail::is_point_v<TShape> && !RectangleConcept<TShape> && requires(const TShape& shape) { shape.bbox(); })
constexpr void Rectangle<PointType, LabelType>::insert(const TShape& shape) {
    insert(shape.bbox());
}
 
 
// ---------------------------------------------------------------------------
// Convex
 
template <class PointType, class LabelType>
constexpr const Rectangle<PointType>& Convex<PointType, LabelType>::bbox() const {
    if (!bbox_.empty()) {
        return bbox_;
    }
    if (points_.empty()) {
        return bbox_ = {};
    }
    if (points_.size() <= 6) {
        return bbox_ = Rectangle<PointType>(points_) + translation_;
    }
 
    // Find the peak of a plain unimodal range [first, last): the keys ascend
    // weakly to a single maximum, then descend weakly. O(log n) binary search.
    const auto unimodalMax = [](auto first, auto last, auto key) {
        auto lo = first;
        auto hi = last - 1;
        while (lo < hi) {
            const auto mid = lo + (hi - lo) / 2;
            if (key(*mid) < key(*(mid + 1)))
                lo = mid + 1;
            else
                hi = mid;
        }
        return lo;
    };
 
    // The boundary splits at its leftmost (points_[0]) and rightmost (maxIndex())
    // vertices into two x-monotone chains. On each chain y is plain unimodal (the
    // split removes the rotation cyclicMax would otherwise have to handle). The
    // lower chain [0, mi] is contiguous; the upper chain [mi, n) wraps back to
    // points_[0], so that vertex is folded into the max_y candidate explicitly.
    const NumberType min_x = points_[0].x();
    const auto mi = maxIndex();
    const NumberType max_x = points_[mi].x();
    const NumberType min_y = unimodalMax(points_.begin(), points_.begin() + mi + 1,
        [](const PointType& p) { return -p.y(); })->y();
    const NumberType max_y = std::max(
        unimodalMax(points_.begin() + mi, points_.end(),
            [](const PointType& p) { return p.y(); })->y(),
        points_[0].y());
 
    return bbox_ = Rectangle<PointType>(min_x, min_y, max_x, max_y, true) + translation_;
}
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> Convex<PointType, LabelType>::fbox() const {
    return bbox().template fbox<ResultNumber>();
}
 
template <class PointType, class LabelType>
template <PointConcept OtherPoint>
constexpr void Convex<PointType, LabelType>::insert(const OtherPoint& point) {
    const PointType vertex = static_cast<PointType>(point);
    if (contains(vertex)) {
        return;
    }
    std::vector<PointType> points = vertices();
    points.push_back(vertex);
    rebuildHull(points);
}
 
template <class PointType, class LabelType>
template <class TShape>
    requires(!detail::is_point_v<TShape> && requires(const TShape& shape) { shape.vertices(); })
constexpr void Convex<PointType, LabelType>::insert(const TShape& shape) {
    // The hull of a shape is the hull of its vertices, so the point-range
    // overload does the work.
    insert(shape.vertices());
}
 
// ---------------------------------------------------------------------------
// Polygon
 
template <class PointType, class LabelType>
constexpr const Rectangle<PointType>& Polygon<PointType, LabelType>::bbox() const {
    if (!bbox_.empty()) {
        return bbox_;
    }
    if (points_.empty()) {
        return bbox_ = {};
    }
    return bbox_ = Rectangle<PointType>(points_) + translation_;
}
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> Polygon<PointType, LabelType>::fbox() const {
    return bbox().template fbox<ResultNumber>();
}
 
// ---------------------------------------------------------------------------
// PolygonSet
 
template <class PointType, class LabelType>
constexpr const Rectangle<PointType>& PolygonSet<PointType, LabelType>::bbox() const {
    if (!bbox_.empty()) {
        return bbox_;
    }
    if (components_.empty()) {
        return bbox_ = {};
    }
    // Every hole lies inside its own component's outer ring, so the components'
    // own boxes — each of them its outer ring's — already cover every vertex.
    Rectangle<PointType> box = components_.front().bbox();
    for (std::size_t i = 1; i < components_.size(); ++i) {
        box.insert(components_[i].bbox());
    }
    return bbox_ = box;
}
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> PolygonSet<PointType, LabelType>::fbox() const {
    return bbox().template fbox<ResultNumber>();
}
 
// ---------------------------------------------------------------------------
// MonotoneChain
 
template <class PointType, class LabelType, class Storage>
constexpr const Rectangle<PointType>& MonotoneChain<PointType, LabelType, Storage>::bbox() const {
    if (!bbox_.empty()) {
        return bbox_;
    }
    if (points_.empty()) {
        return bbox_ = {};
    }
    return bbox_ = Rectangle<PointType>(points_) + translation_;
}
 
template <class PointType, class LabelType, class Storage>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> MonotoneChain<PointType, LabelType, Storage>::fbox() const {
    return bbox().template fbox<ResultNumber>();
}
 
// ---------------------------------------------------------------------------
// Polyline
 
template <class PointType, class LabelType>
constexpr const Rectangle<PointType>& Polyline<PointType, LabelType>::bbox() const {
    if (!bbox_.empty()) {
        return bbox_;
    }
    if (points_.empty()) {
        return bbox_ = {};
    }
    return bbox_ = Rectangle<PointType>(points_) + translation_;
}
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> Polyline<PointType, LabelType>::fbox() const {
    return bbox().template fbox<ResultNumber>();
}
 
 
// ---------------------------------------------------------------------------
// HalfplaneIntersection
 
template <class PointType, class LabelType>
template <class ResultNumber>
constexpr Rectangle<Point<ResultNumber, typename PointType::LabelType>>
HalfplaneIntersection<PointType, LabelType>::bbox() const {
    using ResultPoint = Point<ResultNumber, typename PointType::LabelType>;
    if (empty() || !isBounded()) {
        throw std::logic_error("HalfplaneIntersection::bbox is only defined for a nonempty bounded region");
    }
    using C = detail::promoted_number_t<detail::promoted_number_t<NumberType>>;
    struct Fraction {
        C num;
        C den;  // positive
    };
    // The extreme point in an outward normal direction a = (dy, -dx) sits
    // where the stored boundary directions bracket the query direction
    // (dx, dy) — the same bracketing as supStatus, but producing the extreme
    // coordinate as an exact fraction. O(log n) per direction.
    const auto extreme = [&](NumberType dx, NumberType dy, bool wantX) -> Fraction {
        const HalfplaneType query(PointType(NumberType(0), NumberType(0)), PointType(dx, dy));
        const std::size_t n = halfplanes_.size();
        const std::size_t pos = linearUpperBound(query);
        const std::size_t predIdx = (pos == 0 ? n : pos) - 1;
        const auto& pred = halfplanes_[predIdx];
        if (detail::directionEqual(pred, query)) {
            // Attained along an axis-parallel edge, where the wanted
            // coordinate is constant.
            return {static_cast<C>(wantX ? pred.source().x() : pred.source().y()), C(1)};
        }
        const auto& succ = halfplanes_[pos == n ? 0 : pos];
        // Vertex of the two boundary lines a_i · v = b_i.
        const auto row = [](const HalfplaneType& h) {
            const C sx = static_cast<C>(h.source().x());
            const C sy = static_cast<C>(h.source().y());
            const C tx = static_cast<C>(h.target().x());
            const C ty = static_cast<C>(h.target().y());
            struct Row { C ax, ay, b; };
            return Row{ty - sy, sx - tx, sx * ty - sy * tx};
        };
        const auto r1 = row(pred);
        const auto r2 = row(succ);
        const C den = r1.ax * r2.ay - r1.ay * r2.ax;  // positive: bounded gaps are below pi
        assert(den > C(0));
        if (wantX) {
            return {r1.b * r2.ay - r2.b * r1.ay, den};
        }
        return {r1.ax * r2.b - r2.ax * r1.b, den};
    };
    // Outward normals +x, -x, +y, -y correspond to boundary directions
    // (0, 1), (0, -1), (-1, 0), (1, 0).
    const Fraction xMax = extreme(NumberType(0), NumberType(1), true);
    const Fraction xMin = extreme(NumberType(0), NumberType(-1), true);
    const Fraction yMax = extreme(NumberType(-1), NumberType(0), false);
    const Fraction yMin = extreme(NumberType(1), NumberType(0), false);
    const auto convert = [](const Fraction& f, bool roundUp) -> ResultNumber {
        if constexpr (std::is_floating_point_v<ResultNumber>) {
            (void)roundUp;
            return detail::asNumber<ResultNumber>(f.num) / detail::asNumber<ResultNumber>(f.den);
        } else if constexpr (is_Rational_v<ResultNumber>) {
            // Exact when the fraction fits the rational's storage type.
            (void)roundUp;
            return detail::asNumber<ResultNumber>(f.num) / detail::asNumber<ResultNumber>(f.den);
        } else if constexpr (is_Rational_v<C> || std::is_floating_point_v<C>) {
            // Exact (Rational) or approximate (floating) coordinates: divide
            // in the working type and convert; no outward rounding applies.
            (void)roundUp;
            return static_cast<ResultNumber>(f.num / f.den);
        } else {
            // Integer coordinates and result: round outward so the box
            // encloses the region.
            C quotient = f.num / f.den;
            const C remainder = f.num % f.den;
            if (remainder != C(0)) {
                if (roundUp && f.num > C(0)) {
                    ++quotient;
                }
                if (!roundUp && f.num < C(0)) {
                    --quotient;
                }
            }
            return static_cast<ResultNumber>(quotient);
        }
    };
    return Rectangle<ResultPoint>(ResultPoint(convert(xMin, false), convert(yMin, false)),
                                  ResultPoint(convert(xMax, true), convert(yMax, true)));
}
 
template <class PointType, class LabelType>
template <std::floating_point ResultNumber>
constexpr Rectangle<Point<ResultNumber>> HalfplaneIntersection<PointType, LabelType>::fbox() const {
    const auto box = bbox<ResultNumber>();
    return Rectangle<Point<ResultNumber>>(Point<ResultNumber>(box.min().x(), box.min().y()),
                                          Point<ResultNumber>(box.max().x(), box.max().y()));
}
 
}  // namespace pgl