sortpoints.hpp

Go to the documentation of this file.

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

 
#include <algorithm>
#include <type_traits>
#include <vector>
 
 
namespace pgl {

template <class Number, class Label, class CenterNumber, class CenterLabel>
void sortAround(std::vector<Point<Number, Label>>& points,
                const Point<CenterNumber, CenterLabel>& p) {
    if (points.size() < 2)
        return;
 
    // A point equal to p has no direction around it, and would compare tied
    // with every other point. Park those at the end and sort the rest.
    const auto first = points.begin();
    const auto last = std::partition(first, points.end(),
                                     [&p](const auto& q) { return q != p; });
    if (last - first < 2)
        return;
 
    // The smallest point defines the reference direction (angle zero).
    const auto reference = *std::min_element(first, last);
 
    // Splits off the directions in the half turn [0, pi) measured
    // counterclockwise from the +x direction, that is the points above the
    // horizontal line through p, ties on it broken by the side of p they fall
    // on. The coordinates are promoted the way the orientation predicate
    // promotes them, so the split and the sort agree on every direction.
    using Compare = std::common_type_t<Number, CenterNumber>;
    const auto firstHalf = [&p](const auto& q) {
        const auto vertical = detail::strongOrder(detail::asNumber<Compare>(q.y()),
                                                  detail::asNumber<Compare>(p.y()));
        if (vertical != 0)
            return vertical > 0;
        return detail::strongOrder(detail::asNumber<Compare>(q.x()),
                                   detail::asNumber<Compare>(p.x())) > 0;
    };
    const auto middle = std::partition(first, last, firstHalf);
 
    // The center is an operand of every orientation sign the sort takes —
    // O(n log n) of them — so its coordinates are converted for the filter once
    // here rather than once per comparison. The two points being compared are
    // different ones each time and are converted as they come.
    using SignCoordinate = detail::orientation_coordinate_t<CenterNumber, Number, Number>;
    const auto center = detail::filtered<SignCoordinate>(p);
 
    // Neither part spans more than half a turn, so within one of them the
    // orientation sign is a consistent order on the directions.
    const auto less = [&p, &center](const auto& a, const auto& b) {
        const auto turn = detail::orientationSignOf(
            center, detail::filtered<SignCoordinate>(a), detail::filtered<SignCoordinate>(b))
            .value();
        if (turn > 0)
            return true;   // b is counterclockwise of a, so a has the smaller angle.
        if (turn < 0)
            return false;
 
        // Same direction from p: the farther point comes first. The distance
        // is taken in the promoted type, so points carrying more precision
        // than the center still compare exactly.
        return p.template squaredDistance<Compare>(a) >
               p.template squaredDistance<Compare>(b);
    };
    std::sort(first, middle, less);
    std::sort(middle, last, less);
 
    // The points now run counterclockwise from the +x direction: rotate the
    // block sharing the reference direction to the front. The search compares
    // directions alone, so it lands on the first point of that block rather
    // than on the reference itself, which the distance tie-break may have put
    // behind others pointing the same way.
    const auto above = firstHalf(reference);
    const auto start = std::lower_bound(
        above ? first : middle, above ? middle : last, reference,
        [&center](const auto& a, const auto& b) {
            return detail::orientationSignOf(center,
                                             detail::filtered<SignCoordinate>(a),
                                             detail::filtered<SignCoordinate>(b))
                       .value() > 0;
        });
    std::rotate(first, start, last);
}
 
namespace detail {
 
// Recursive median Hilbert sort (CGAL's median policy): split the range into the
// four quadrants the Hilbert curve visits, using nested medians, then recurse
// into each quadrant with the rotated/reflected curve state. `xAxis` selects the
// primary split axis; `upX`/`upY` give the curve's direction along each axis.
// Only the two coordinate comparators are used, so it is exact for any numeric
// type and never constructs intermediate coordinates.
template <class RandomIt, class LessX, class LessY>
void hilbertSortMedian(RandomIt begin, RandomIt end, bool xAxis, bool upX, bool upY,
                       const LessX& lessX, const LessY& lessY) {
    if (end - begin <= 1) {
        return;
    }
 
    const RandomIt m0 = begin, m4 = end;
    const RandomIt m2 = m0 + (m4 - m0) / 2;  // primary-axis median of the whole range
    const RandomIt m1 = m0 + (m2 - m0) / 2;  // secondary-axis median of the lower half
    const RandomIt m3 = m2 + (m4 - m2) / 2;  // secondary-axis median of the upper half
 
    // Rearranges [first,last) so *mid is its median under `less`, ascending when
    // `up`, descending otherwise.
    const auto split = [](RandomIt first, RandomIt mid, RandomIt last, const auto& less,
                          bool up) {
        if (up) {
            std::nth_element(first, mid, last, less);
        } else {
            std::nth_element(first, mid, last,
                             [&less](const auto& a, const auto& b) { return less(b, a); });
        }
    };
 
    if (xAxis) {
        split(m0, m2, m4, lessX, upX);
        split(m0, m1, m2, lessY, upY);
        split(m2, m3, m4, lessY, !upY);
    } else {
        split(m0, m2, m4, lessY, upY);
        split(m0, m1, m2, lessX, upX);
        split(m2, m3, m4, lessX, !upX);
    }
 
    hilbertSortMedian(m0, m1, !xAxis, upY, upX, lessX, lessY);
    hilbertSortMedian(m1, m2, xAxis, upX, upY, lessX, lessY);
    hilbertSortMedian(m2, m3, xAxis, upX, upY, lessX, lessY);
    hilbertSortMedian(m3, m4, !xAxis, !upY, !upX, lessX, lessY);
}
 
}  // namespace detail

template <class Number, class Label>
void hilbertSort(std::vector<Point<Number, Label>>& points) {
    const auto lessX = [](const Point<Number, Label>& a, const Point<Number, Label>& b) {
        return a.x() < b.x();
    };
    const auto lessY = [](const Point<Number, Label>& a, const Point<Number, Label>& b) {
        return a.y() < b.y();
    };
    detail::hilbertSortMedian(points.begin(), points.end(), true, false, false, lessX, lessY);
}
 
}  // namespace pgl