xysweep.hpp

Go to the documentation of this file.

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

 
#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
 
namespace pgl::detail {

inline constexpr std::size_t xySweepMinSegments = 128;

template <SegmentConcept Segment, class Visitor>
bool visitXYSweepPairs(const std::vector<Segment>& segments, Visitor visit) {
    using PointType = typename Segment::PointType;
    using NumberType = typename PointType::NumberType;
    using IndexedSegment = pgl::Segment<PointType, std::uint32_t>;
 
    const std::size_t count = segments.size();
    if (count < 2) {
        return false;  // no pair to visit
    }
    if (count > static_cast<std::size_t>(std::numeric_limits<std::uint32_t>::max())) {
        throw std::length_error("xy sweep exceeds its 32-bit segment capacity");
    }
 
    // The index travels with the segment as its label, so a segment reported by
    // the tree names its position in the input.
    std::vector<IndexedSegment> indexed;
    indexed.reserve(count);
    for (std::size_t i = 0; i < count; ++i) {
        indexed.emplace_back(segments[i].min(), segments[i].max(),
                             static_cast<std::uint32_t>(i));
    }
 
    struct Event {
        NumberType x;
        bool closing;
        std::uint32_t index;
    };
    std::vector<Event> events;
    events.reserve(2 * count);
    for (const IndexedSegment& segment : indexed) {
        const auto box = segment.bbox();
        events.push_back({box.min().x(), false, segment.label()});
        events.push_back({box.max().x(), true, segment.label()});
    }
    // Openings come before closings at a shared abscissa, so two boxes meeting
    // in a single vertical line are still compared, and a vertical segment is
    // queried before the same abscissa closes it.
    std::sort(events.begin(), events.end(), [](const Event& a, const Event& b) {
        if (a.x < b.x) {
            return true;
        }
        if (b.x < a.x) {
            return false;
        }
        return !a.closing && b.closing;
    });
 
    pgl::IntervalTree<IndexedSegment, pgl::ProjectionAxis::y> active;
    for (const Event& event : events) {
        const IndexedSegment& segment = indexed[event.index];
        if (event.closing) {
            active.erase(segment);
            continue;
        }
        // Everything the tree holds already overlaps this segment in x, so the
        // shapes it reports for the y-extent are exactly the boxes that overlap.
        const bool stopped =
            active.visitProjectionsIntersecting(segment, [&](const IndexedSegment& other) {
                return visit(static_cast<std::size_t>(event.index),
                             static_cast<std::size_t>(other.label()));
            });
        if (stopped) {
            return true;
        }
        active.insert(segment);
    }
    return false;
}
 
}  // namespace pgl::detail
 
 
namespace pgl {

template<class Rational = pgl::Rational<pgl::BigInt>, class Container>
auto xyCrossings(const Container &segments) {
    using Point = Container::value_type::PointType;
    // Rebuilt carrying the input's own label type: defaulting it to NoLabel
    // here dropped a segment label that the caller had put on to identify the
    // pairs coming back, and dropped it silently, since the conversion is
    // well-formed either way.
    using Seg = pgl::Segment<Point, typename Container::value_type::LabelType>;
    std::vector<Seg> v;
    for (const auto &s : segments) {
        Seg converted = s;
        v.push_back(converted);
    }
 
    std::vector<std::array<Seg,2>> ret;
    pgl::detail::visitXYSweepPairs(v, [&v,&ret](std::size_t i, std::size_t j) {
        Seg s1 = v[i];
        Seg s2 = v[j];
        if (s1.crosses(s2)) {
            if (s2 < s1)
                std::swap(s1,s2);
            ret.push_back({s1,s2});
        }
    });
 
    return ret;
}

template<class Rational = pgl::Rational<pgl::BigInt>, class Container>
auto xyIntersections(const Container &segments) {
    using Point = Container::value_type::PointType;
    // Rebuilt carrying the input's own label type: defaulting it to NoLabel
    // here dropped a segment label that the caller had put on to identify the
    // pairs coming back, and dropped it silently, since the conversion is
    // well-formed either way.
    using Seg = pgl::Segment<Point, typename Container::value_type::LabelType>;
    std::vector<Seg> v;
    for (const auto &s : segments) {
        Seg converted = s;
        v.push_back(converted);
    }
 
    std::vector<std::array<Seg,2>> ret;
    pgl::detail::visitXYSweepPairs(v, [&v,&ret](std::size_t i, std::size_t j) {
        Seg s1 = v[i];
        Seg s2 = v[j];
        if (s1.intersects(s2)) {
            if (s2 < s1)
                std::swap(s1,s2);
            ret.push_back({s1,s2});
        }
    });
 
    return ret;
}

template <class PointType_, class LabelType>
template <class Rational>
bool Polygon<PointType_, LabelType>::isSimple() const {
    using Number = typename PointType::NumberType;
    const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(size());
    if (n < 3) {
        return false;
    }
 
    std::vector<pgl::Segment<PointType>> edges;
    edges.reserve(static_cast<std::size_t>(n));
    for (std::ptrdiff_t i = 0; i < n; ++i) {
        pgl::Segment<PointType> edge(get(i), get(i + 1));
        if (edge.isDegenerate()) {
            return false;  // zero-length edge / repeated vertex
        }
        edges.push_back(edge);
    }
 
    // Exact sweep for large integer/rational polygons; the bounding-box sweep
    // for large floating-point ones, which the exact sweep cannot handle; brute
    // force for the small ones either way.
    if constexpr (!std::is_floating_point_v<Number>) {
        if (n > 8) {
            pgl::detail::BentleyOttmann<Rational, pgl::Segment<PointType>> bo;
            return bo.testPolygon(edges);
        }
    } else if (edges.size() > detail::xySweepMinSegments) {
        const std::size_t last = edges.size() - 1;
        bool notSimple = false;
        detail::visitXYSweepPairs(edges, [&edges, last, &notSimple](std::size_t a, std::size_t b) {
            // The ring closes, so the first and last edges are consecutive too.
            const bool adjacent = (a + 1 == b) || (b + 1 == a) ||
                                  (a == 0 && b == last) || (b == 0 && a == last);
            notSimple = adjacent ? edges[a].interiorsIntersect(edges[b])
                                 : edges[a].intersects(edges[b]);
            return notSimple;  // the first violation ends the sweep
        });
        return !notSimple;
    }
 
    for (std::ptrdiff_t i = 0; i < n; ++i) {
        for (std::ptrdiff_t j = i + 1; j < n; ++j) {
            const bool adjacent = (j == i + 1) || (i == 0 && j == n - 1);
            if (adjacent) {
                if (edges[i].interiorsIntersect(edges[j])) {
                    return false;  // consecutive edges overlap beyond the shared vertex
                }
            } else if (edges[i].intersects(edges[j])) {
                return false;  // non-adjacent edges must be disjoint
            }
        }
    }
    return true;
}

template <class PointType_, class TLabel>
template <class Rational>
bool Polyline<PointType_, TLabel>::isSimple() const {
    using Number = typename PointType::NumberType;
    const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(size());
    if (n < 2) {
        return true;  // no edge: vacuously simple
    }
 
    std::vector<pgl::Segment<PointType>> edges;
    edges.reserve(static_cast<std::size_t>(n - 1));
    for (std::ptrdiff_t i = 0; i + 1 < n; ++i) {
        pgl::Segment<PointType> edge(get(i), get(i + 1));
        if (edge.isDegenerate()) {
            return false;  // zero-length edge / repeated consecutive vertex
        }
        edges.push_back(edge);
    }
 
    // Exact sweep for large integer/rational polylines; the bounding-box sweep
    // for large floating-point ones, which the exact sweep cannot handle; brute
    // force for the small ones either way.
    if constexpr (!std::is_floating_point_v<Number>) {
        if (edges.size() > 8) {
            pgl::detail::BentleyOttmann<Rational, pgl::Segment<PointType>> bo;
            return bo.testPolyLine(edges);
        }
    } else if (edges.size() > detail::xySweepMinSegments) {
        bool notSimple = false;
        detail::visitXYSweepPairs(edges, [&edges, &notSimple](std::size_t a, std::size_t b) {
            // In an open chain the first and last edges are NOT adjacent, so a
            // closed polyline (first vertex equal to the last) is not simple.
            const bool adjacent = (a + 1 == b) || (b + 1 == a);
            notSimple = adjacent ? edges[a].interiorsIntersect(edges[b])
                                 : edges[a].intersects(edges[b]);
            return notSimple;  // the first violation ends the sweep
        });
        return !notSimple;
    }
 
    const std::ptrdiff_t m = static_cast<std::ptrdiff_t>(edges.size());
    for (std::ptrdiff_t i = 0; i < m; ++i) {
        for (std::ptrdiff_t j = i + 1; j < m; ++j) {
            // In an open chain the first and last edges are NOT adjacent, so a
            // closed polyline (first vertex equal to the last) is not simple.
            const bool adjacent = (j == i + 1);
            if (adjacent) {
                if (edges[i].interiorsIntersect(edges[j])) {
                    return false;  // consecutive edges overlap beyond the shared vertex
                }
            } else if (edges[i].intersects(edges[j])) {
                return false;  // non-adjacent edges must be disjoint
            }
        }
    }
    return true;
}
 
}  // namespace pgl