Another way to find the nearest neighbor points in a 2D plane

Question

I was inspired by another question to post another method of finding the two points in a plane that are closest to each other¹.

This one is a line-sweep algorithm. It works roughly like this:

1. Sort the points based on their X coordinates.
2. Take two left-most points. They give us our first guess at the shortest distance (call it D).
3. Insert those two points into a set that's sorted based on Y coordinates. This collection forms a vertical "band" of points whose X coordinates are within D units of the X coordinate of the current point.
4. Consider the next point to the right of those currently in the "band" as the current point.
5. Trim the band to remove points more than D units away in the X dimension from the current point.
6. Find the points in the band that are vertically within D units of the Y coordinate of the current point.
7. Look through the points in that rectangle (maximum of 6) to see if any is closer than D units from the current point.
8. If so, record the points and distance.
9. Repeat from step 4 for remaining points.

Here's the code:

#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <set>
#include <cassert>

struct point {
    double x, y;

    // Used by the `set<point>` to keep the points in the band
    // sorted by Y coordinates.
    bool operator<(point const &other) const {
        return y < other.y;
    }

    friend std::ostream &operator<<(std::ostream &os, point const &p) {
        return os << "(" << p.x << ", " << p.y << ")";
    }
};

double dist(point const &a, point const &b) {
    return std::hypot(a.x - b.x, a.y - b.y);
}

// We're going to modify the input we receive, so we receive it by value.
// If we knew that the source was going to be modifiable, we could receive
// by non-const reference instead;
std::pair<point, point> min_dist(std::vector<point> points) {

    std::sort(points.begin(), points.end(),
        [](point const &a, point const &b) {
        return a.x < b.x;
    }
    );

    // First and last points from `point` that are currently in the "band".
    auto first = points.cbegin();
    auto last = first + 1;

    // The two closest points we've found so far:
    auto first_point = *first;
    auto second_point = *last;

    std::set<point> band{ *first, *last };

    double d = dist(*first, *last);

    while (++last != points.end()) {
        while (last->x - first->x > d) {
            band.erase(*first);
            ++first;
        }

        auto begin = band.lower_bound({ 0, last->y - d });
        auto end = band.upper_bound({ 0, last->y + d });

        assert(std::distance(begin, end) <= 6);

        for (auto p = begin; p != end; ++p) {
            d = std::min(d, dist(*p, *last));
            first_point = *p;
            second_point = *last;
        }

        band.insert(*last);
    }
    return std::make_pair(first_point, second_point);
}

int main() {
    std::vector<point> points{
        {1, 1},
        {17, 9},
        {23, 23},
        {3, 3},
        {100, 100},
        {200, 200},
        {24, 24},
        {300, 300}
    };

    auto r = min_dist(points);

    std::cout << "Closest points: " << r.first << ", " << r.second 
            << ". Distance = "<< dist(r.first, r.second);
}

---

^{1. In this case, I've used Euclidian distance, but another metric such as Manhattan distance could be used as well.}

Answer 1

Your algorithm looks (and works) great, but I found two problems with part of code responsible for finding closest pair within band.

Wrong pair of points is returned as result

The closest distance in d variable is correctly updated, but first_point and second_points are updated even if d < dist(*p, *last), which leads to returning wrong pair of points. You can reproduce the problem on following set of points: {-1, 10}, {0, 10}, {0, 9}, {0, 8}, {0, 7}, {0, 6}, {0, 5}, {0, 4}, {0, 3}, {0, 2}, {0, 1}, {0, 0}, {0.99, 0.99}

You can fix that by replacing your for loop with the following one:

for (auto p = begin; p != end; ++p) {
    if (d > dist(*p, *last)) {
        first_point = *p;
        second_point = *last;
        d = dist(first_point, second_point);
    }
}

for loop is not executed when there is only one point in band

If there is only one element in band, then both lower_bound and upper_bound return iterator pointing to the same element in band, which means that for loop will not be executed, because the condition p != end is not met. You can reproduce the problem on following set of points: {-1, 0}, {0 ,0}, {0.1, 0.1}