closestpair.hpp

Go to the documentation of this file.

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

 
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <iterator>
#include <type_traits>
#include <utility>
#include <vector>
 
 
namespace pgl {
 
namespace detail {
 
template <class Container>
using closest_pair_input_point_t =
    std::remove_cvref_t<decltype(*std::begin(std::declval<const Container&>()))>;

template <class PointType>
using closest_pair_coordinate_t = promoted_number_t<typename PointType::NumberType>;
 
template <class PointType>
using closest_pair_distance_t =
    decltype(std::declval<const PointType&>()
                 .template squaredDistance<closest_pair_coordinate_t<PointType>>(
                     std::declval<const PointType&>()));

template <class PointType>
struct ClosestPairCandidate {
    PointType first;
    PointType second;
    closest_pair_distance_t<PointType> squaredDistance;
};

template <class Coordinate>
inline constexpr std::size_t closestPairBruteForceThreshold =
    arbitraryPrecision<Coordinate> ? 4 : (is_Rational_v<Coordinate> ? 12 : 6);

template <class PointType>
constexpr bool closestPairLessY(const PointType& a, const PointType& b) {
    if (a.y() < b.y())
        return true;
    if (b.y() < a.y())
        return false;
    return a.x() < b.x();
}

template <class PointType>
void closestPairBaseCase(const PointType* points, std::size_t count,
                         ClosestPairCandidate<PointType>& best) {
    using Coordinate = closest_pair_coordinate_t<PointType>;
    for (std::size_t i = 0; i < count; ++i) {
        for (std::size_t j = i + 1; j < count; ++j) {
            const auto squared = points[i].template squaredDistance<Coordinate>(points[j]);
            if (squared < best.squaredDistance) {
                best = {points[i], points[j], squared};
            }
        }
    }
}

template <class PointType>
std::size_t closestPairGatherStrip(const PointType* points, std::size_t count, std::size_t half,
                                   closest_pair_coordinate_t<PointType> splitX,
                                   const ClosestPairCandidate<PointType>& best, PointType* strip) {
    using Coordinate = closest_pair_coordinate_t<PointType>;
 
    // Squared throughout, so the half-width is never rooted and the test stays
    // exact: |x - splitX| < d becomes (x - splitX)^2 < d^2. A best of zero
    // admits nothing, which is right — no pair can beat a coincident one.
    const auto inside = [&](std::size_t index) {
        const Coordinate dx = detail::asNumber<Coordinate>(points[index].x()) - splitX;
        return dx * dx < best.squaredDistance;
    };
 
    std::size_t low = half;
    while (low > 0 && inside(low - 1)) {
        --low;
    }
    std::size_t high = half;
    while (high < count && inside(high)) {
        ++high;
    }
 
    std::size_t stripCount = 0;
    for (std::size_t i = low; i < high; ++i) {
        strip[stripCount++] = points[i];
    }
    return stripCount;
}

template <class PointType>
void closestPairScanStrip(const PointType* strip, std::size_t stripCount,
                          ClosestPairCandidate<PointType>& best) {
    using Coordinate = closest_pair_coordinate_t<PointType>;
    for (std::size_t i = 0; i < stripCount; ++i) {
        for (std::size_t j = i + 1; j < stripCount; ++j) {
            const Coordinate dy =
                detail::asNumber<Coordinate>(strip[j].y()) - detail::asNumber<Coordinate>(strip[i].y());
            if (dy * dy >= best.squaredDistance) {
                break;
            }
            const auto squared = strip[i].template squaredDistance<Coordinate>(strip[j]);
            if (squared < best.squaredDistance) {
                best = {strip[i], strip[j], squared};
            }
        }
    }
}

template <std::size_t Threshold, class PointType>
void closestPairRecursive(const PointType* points, std::size_t count, PointType* scratch,
                          ClosestPairCandidate<PointType>& best) {
    using Coordinate = closest_pair_coordinate_t<PointType>;
    static_assert(Threshold >= 3,
                  "a range above the threshold must split into halves of at least two points");
    assert(count >= 2);
 
    // Nothing above needs a y-order, so brute force leaves the range untouched.
    if (count <= Threshold) {
        closestPairBaseCase(points, count, best);
        return;
    }
 
    // The vertical line through points[half] separates the halves: everything
    // left of it has a smaller-or-equal abscissa, everything right of it a
    // greater-or-equal one.
    const std::size_t half = count / 2;
    const Coordinate splitX = static_cast<Coordinate>(points[half].x());
 
    closestPairRecursive<Threshold>(points, half, scratch, best);
    closestPairRecursive<Threshold>(points + half, count - half, scratch + half, best);
 
    // The range is in x-order, so the strip comes out in x-order too and has to
    // be sorted before it can be scanned.
    const std::size_t stripCount =
        closestPairGatherStrip(points, count, half, splitX, best, scratch);
    std::sort(scratch, scratch + stripCount, closestPairLessY<PointType>);
    closestPairScanStrip(scratch, stripCount, best);
}

template <std::size_t Threshold, class Container>
[[nodiscard]] Segment<closest_pair_input_point_t<Container>>
closestPairDriver(const Container& input) {
    using InputPoint = closest_pair_input_point_t<Container>;
    using Coordinate = closest_pair_coordinate_t<InputPoint>;
    static_assert(PointConcept<InputPoint>,
                  "closestPair requires a container of pgl::Point values");
 
    std::vector<InputPoint> points(std::begin(input), std::end(input));
    assert(points.size() >= 2);
 
    std::sort(points.begin(), points.end());
    std::vector<InputPoint> scratch(points.size());
 
    ClosestPairCandidate<InputPoint> best{
        points[0], points[1],
        points[0].template squaredDistance<Coordinate>(points[1])};
    closestPairRecursive<Threshold>(points.data(), points.size(), scratch.data(), best);
 
    return Segment<InputPoint>(best.first, best.second);
}
 
}  // namespace detail

template <class Container>
[[nodiscard]] Segment<detail::closest_pair_input_point_t<Container>>
closestPair(const Container& input) {
    using Coordinate =
        detail::closest_pair_coordinate_t<detail::closest_pair_input_point_t<Container>>;
    return detail::closestPairDriver<detail::closestPairBruteForceThreshold<Coordinate>>(input);
}
 
}  // namespace pgl