redbluesweep.hpp

Go to the documentation of this file.

#pragma once
 
#include "algorithm/intersections.hpp"

 
#include <algorithm>
#include <cstddef>
#include <ranges>
#include <set>
#include <type_traits>
#include <utility>
#include <vector>
 
 
namespace pgl {

enum class BoundaryContact {
    Disjoint,  
    Touching,  
    Crossing   
};
 
namespace detail {

template <class Number>
class RedBlueSweeper {
  public:
    using SweepPoint = pgl::Point<Number>;
    using SweepSegment = pgl::Segment<SweepPoint>;
 
  private:
    struct Edge {
        SweepPoint lo;   
        SweepPoint hi;   
        bool blue = false;

        bool upright() const { return lo.x() == hi.x(); }
        SweepSegment segment() const { return SweepSegment(lo, hi); }
    };

    enum class Kind : unsigned char {
        Right = 0,       
        Upper = 1,       
        UprightLow = 2,  
        Left = 3         
    };
 
    struct Event {
        SweepPoint at;
        std::size_t edge = 0;
        Kind kind = Kind::Left;
 
        bool operator<(const Event& other) const {
            if (at < other.at) return true;
            if (other.at < at) return false;
            if (kind != other.kind) return kind < other.kind;
            return edge < other.edge;
        }
    };

    struct Below {
        const std::vector<Edge>* edges = nullptr;
        using is_transparent = void;

        static int sideOf(const Edge& one, const Edge& that) {
            int side = signOf(orientationSign(one.lo, one.hi, that.lo));
            if (side == 0) {
                side = signOf(orientationSign(one.lo, one.hi, that.hi));
            }
            // one.lo -> one.hi points right, so counterclockwise means above.
            return side > 0 ? -1 : (side < 0 ? 1 : 0);
        }
 
        int compare(std::size_t first, std::size_t second) const {
            if (first == second) {
                return 0;
            }
            const Edge& one = (*edges)[first];
            const Edge& other = (*edges)[second];
            const int side =
                one.lo.x() <= other.lo.x() ? sideOf(one, other) : -sideOf(other, one);
            if (side != 0) {
                return side;
            }
            return first < second ? -1 : 1;
        }
 
        bool operator()(std::size_t first, std::size_t second) const {
            return compare(first, second) < 0;
        }
 
        // A bare point probes the status at the sweep line: it sorts below every
        // edge it lies on, so lower_bound(p) is the first edge at or above p.
        bool operator()(std::size_t first, const SweepPoint& probe) const {
            const Edge& one = (*edges)[first];
            return signOf(orientationSign(one.lo, one.hi, probe)) > 0;
        }
        bool operator()(const SweepPoint& probe, std::size_t second) const {
            const Edge& other = (*edges)[second];
            return signOf(orientationSign(other.lo, other.hi, probe)) <= 0;
        }
    };
 
    std::vector<Edge> edges_;
    std::vector<Event> events_;
    std::set<std::size_t, Below> status_;
    bool touching_ = false;

    struct Extent {
        Number xlo{}, ylo{}, xhi{}, yhi{};
        bool any = false;
 
        void insert(const SweepPoint& point) {
            if (!any) {
                xlo = xhi = point.x();
                ylo = yhi = point.y();
                any = true;
                return;
            }
            if (point.x() < xlo) xlo = point.x();
            if (xhi < point.x()) xhi = point.x();
            if (point.y() < ylo) ylo = point.y();
            if (yhi < point.y()) yhi = point.y();
        }
 
        bool overlaps(const SweepPoint& lo, const SweepPoint& hi) const {
            return any && !(hi.x() < xlo) && !(xhi < lo.x()) &&
                   !(std::max(lo.y(), hi.y()) < ylo) && !(yhi < std::min(lo.y(), hi.y()));
        }
    };
 
    template <class Range>
    void gather(const Range& range, bool blue, Extent& extent) {
        for (const auto& edge : range) {
            Edge stored{SweepPoint(edge.min()), SweepPoint(edge.max()), blue};
            extent.insert(stored.lo);
            extent.insert(stored.hi);
            edges_.push_back(std::move(stored));
        }
    }

    bool adjacent(std::size_t first, std::size_t second) {
        const Edge& one = edges_[first];
        const Edge& other = edges_[second];
        if (one.blue == other.blue) {
            return false;
        }
        const SweepSegment mine = one.segment();
        const SweepSegment theirs = other.segment();
        if (mine.crosses(theirs)) {
            return true;
        }
        if (mine.intersects(theirs)) {
            touching_ = true;
        }
        return false;
    }

    bool scanUpright(std::size_t index) {
        const Edge& upright = edges_[index];
        const SweepSegment bar = upright.segment();
        for (auto it = status_.lower_bound(upright.lo); it != status_.end(); ++it) {
            const Edge& other = edges_[*it];
            if (signOf(orientationSign(other.lo, other.hi, upright.hi)) < 0) {
                break;  // past the top of the vertical
            }
            if (other.blue == upright.blue) {
                continue;
            }
            if (other.segment().crosses(bar)) {
                return true;
            }
            touching_ = true;
            break;
        }
        return false;
    }

    void uprightPairs(std::size_t begin, std::size_t end) {
        std::vector<std::size_t> red;
        std::vector<std::size_t> blue;
        for (std::size_t k = begin; k < end; ++k) {
            if (events_[k].kind != Kind::UprightLow) {
                continue;
            }
            (edges_[events_[k].edge].blue ? blue : red).push_back(events_[k].edge);
        }
        std::size_t i = 0;
        std::size_t j = 0;
        while (i < red.size() && j < blue.size()) {
            const Edge& one = edges_[red[i]];
            const Edge& other = edges_[blue[j]];
            if (one.hi.y() < other.lo.y()) {
                ++i;
            } else if (other.hi.y() < one.lo.y()) {
                ++j;
            } else {
                touching_ = true;
                return;
            }
        }
    }

    bool batch(std::size_t begin, std::size_t end) {
        // A point that is an endpoint of edges of both colours is a meeting the
        // status structure cannot see: an edge ending there is erased before one
        // starting there is inserted, so such a pair is never adjacent. It is
        // never a crossing either — a shared endpoint makes one of the four
        // orientations in Segment::crosses vanish — so recognizing it here, from
        // the events alone, closes the only gap in the adjacency argument.
        for (std::size_t k = begin; k < end;) {
            std::size_t run = k;
            bool red = false;
            bool blue = false;
            while (run < end && !(events_[k].at < events_[run].at)) {
                (edges_[events_[run].edge].blue ? blue : red) = true;
                ++run;
            }
            if (red && blue) {
                touching_ = true;
            }
            k = run;
        }
 
        // Verticals see the status twice, since neither snapshot alone holds
        // every edge spanning this x: before the deletions it still has the ones
        // ending here, after the insertions it has the ones starting here.
        for (std::size_t k = begin; k < end; ++k) {
            if (events_[k].kind == Kind::UprightLow && scanUpright(events_[k].edge)) {
                return true;
            }
        }
        uprightPairs(begin, end);
 
        for (std::size_t k = begin; k < end; ++k) {
            if (events_[k].kind != Kind::Right) {
                continue;
            }
            const auto it = status_.find(events_[k].edge);
            if (it == status_.end()) {
                continue;
            }
            const bool hasBelow = it != status_.begin();
            auto below = it;
            if (hasBelow) {
                --below;
            }
            const auto above = std::next(it);
            status_.erase(it);
            if (hasBelow && above != status_.end() && adjacent(*below, *above)) {
                return true;
            }
        }
 
        for (std::size_t k = begin; k < end; ++k) {
            if (events_[k].kind != Kind::Left) {
                continue;
            }
            const auto [it, fresh] = status_.insert(events_[k].edge);
            if (!fresh) {
                continue;
            }
            if (it != status_.begin() && adjacent(*std::prev(it), *it)) {
                return true;
            }
            const auto above = std::next(it);
            if (above != status_.end() && adjacent(*it, *above)) {
                return true;
            }
        }
 
        for (std::size_t k = begin; k < end; ++k) {
            if (events_[k].kind == Kind::UprightLow && scanUpright(events_[k].edge)) {
                return true;
            }
        }
        return false;
    }
 
  public:
    template <class RedRange, class BlueRange>
    RedBlueSweeper(const RedRange& red, const BlueRange& blue) : status_(Below{&edges_}) {
        Extent redExtent;
        Extent blueExtent;
        gather(red, false, redExtent);
        const std::size_t split = edges_.size();
        gather(blue, true, blueExtent);
 
        // An edge whose box misses the other colour's box cannot meet anything
        // of that colour, and dropping it can only remove intersections, never
        // create one. On two boundaries that overlap everywhere this costs one
        // extra pass; on two that barely touch it removes almost all the work.
        std::vector<Edge> kept;
        kept.reserve(edges_.size());
        for (std::size_t i = 0; i < edges_.size(); ++i) {
            const Extent& against = i < split ? blueExtent : redExtent;
            if (against.overlaps(edges_[i].lo, edges_[i].hi)) {
                kept.push_back(edges_[i]);
            }
        }
        edges_ = std::move(kept);
 
        events_.reserve(2 * edges_.size());
        for (std::size_t i = 0; i < edges_.size(); ++i) {
            const Edge& edge = edges_[i];
            if (edge.upright()) {
                events_.push_back(Event{edge.lo, i, Kind::UprightLow});
                events_.push_back(Event{edge.hi, i, Kind::Upper});
            } else {
                events_.push_back(Event{edge.lo, i, Kind::Left});
                events_.push_back(Event{edge.hi, i, Kind::Right});
            }
        }
    }
 
    BoundaryContact run() {
        std::sort(events_.begin(), events_.end());
        std::size_t i = 0;
        while (i < events_.size()) {
            std::size_t j = i;
            while (j < events_.size() && !(events_[i].at.x() < events_[j].at.x())) {
                ++j;
            }
            if (batch(i, j)) {
                return BoundaryContact::Crossing;
            }
            i = j;
        }
        return touching_ ? BoundaryContact::Touching : BoundaryContact::Disjoint;
    }
};
 
}  // namespace detail

template <class RedRange, class BlueRange>
BoundaryContact redBlueSweep(const RedRange& red, const BlueRange& blue) {
    using RedNumber = typename std::ranges::range_value_t<RedRange>::PointType::NumberType;
    using BlueNumber = typename std::ranges::range_value_t<BlueRange>::PointType::NumberType;
    using Number = std::common_type_t<RedNumber, BlueNumber>;
    return detail::RedBlueSweeper<Number>(red, blue).run();
}

template <class RedRange, class BlueRange>
bool boundariesCross(const RedRange& red, const BlueRange& blue) {
    return redBlueSweep(red, blue) == BoundaryContact::Crossing;
}

template <class RedRange, class BlueRange>
bool boundariesMeet(const RedRange& red, const BlueRange& blue) {
    return redBlueSweep(red, blue) != BoundaryContact::Disjoint;
}

struct SweepContact {
    bool crossed;
    bool met;
};
 
template <class RedRange, class BlueRange>
SweepContact boundaryContactBits(const RedRange& red, const BlueRange& blue) {
    const BoundaryContact contact = redBlueSweep(red, blue);
    return {contact == BoundaryContact::Crossing, contact != BoundaryContact::Disjoint};
}


#ifndef PGL_BOUNDARY_STRATEGY
#define PGL_BOUNDARY_STRATEGY 0
#endif
 
namespace detail {

template <class Shape>
constexpr std::size_t boundaryEdgeCount(const Shape& shape) {
    if constexpr (PolygonWithHolesConcept<Shape>) {
        return shape.vertexCount();
    } else {
        return shape.size();
    }
}

constexpr bool clearsEdgeFloor(std::size_t redEdges, std::size_t blueEdges) {
    constexpr std::size_t kMinEdges = 128;
    return redEdges >= kMinEdges && blueEdges >= kMinEdges;
}

constexpr bool sweepBeatsChains(std::size_t redEdges, std::size_t redChains,
                                std::size_t blueEdges, std::size_t blueChains) {
    const std::size_t edges = redEdges + blueEdges;
    return 3 * redChains * blueChains > edges;
}
 
}  // namespace detail

template <class RedShape, class BlueShape>
constexpr bool preferSweep([[maybe_unused]] const RedShape& red,
                           [[maybe_unused]] const BlueShape& blue) {
#if PGL_BOUNDARY_STRATEGY == 1
    return false;
#elif PGL_BOUNDARY_STRATEGY == 2
    return true;
#else
    const std::size_t redEdges = detail::boundaryEdgeCount(red);
    const std::size_t blueEdges = detail::boundaryEdgeCount(blue);
    // Settled on edge counts alone — O(1), read straight off size()/vertexCount()
    // — before either chainCount() is read, so an operand pair too small to be
    // worth a sweep never pays for that O(n) pass (see clearsEdgeFloor).
    if (!detail::clearsEdgeFloor(redEdges, blueEdges)) {
        return false;
    }
    return detail::sweepBeatsChains(redEdges, red.chainCount(), blueEdges, blue.chainCount());
#endif
}


template <PolygonConcept OuterPolygon, PolygonConcept InnerPolygon>
bool sweepContains(const OuterPolygon& outer, const InnerPolygon& inner) {
    using InnerPoint = typename InnerPolygon::PointType;
 
    if (inner.size() == 0) {
        return true;
    }
    if (!outer.bbox().contains(inner.bbox())) {
        return false;
    }
    if (inner.size() == 1) {
        return outer.contains(inner[0]);
    }
    if (const auto vertex = inner.getIfPoint()) {
        return outer.contains(*vertex);
    }
    if (outer.isPoint()) {
        // `inner` has two distinct vertices by the test above.
        return false;
    }
    if (!outer.contains(inner.get(0))) {
        return false;
    }
 
    const BoundaryContact contact = redBlueSweep(outer.edgesView(), inner.edgesView());
    if (contact == BoundaryContact::Crossing) {
        return false;
    }
    if (contact == BoundaryContact::Disjoint) {
        // One vertex inside and no boundary contact leaves the whole of `inner`
        // inside: its boundary cannot reach out without meeting `outer`'s.
        return true;
    }
 
    // The boundaries meet without a crossing having been found. Whether one is
    // hiding behind the touchings no longer matters: checking every edge settles
    // containment outright.
    for (std::size_t i = 0; i < inner.size(); ++i) {
        if (!outer.contains(Segment<InnerPoint>(inner[i], inner[(i + 1) % inner.size()]))) {
            return false;
        }
    }
    return true;
}
 
}  // namespace pgl