Another way to find the nearest neighbor points in a 2D plane I was inspired by another question to post another method of finding the two points in a plane that are closest to each other1.

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


Sort the points based on their X coordinates.
Take two left-most points. They give us our first guess at the shortest distance (call it D).
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.
Consider the next point to the right of those currently in the "band" as the current point.
Trim the band to remove points more than D units away in the X dimension from the current point.
Find the points in the band that are vertically within D units of the Y coordinate of the current point.
Look through the points in that rectangle (maximum of 6) to see if any is closer than D units from the current point.
If so, record the points and distance.
Repeat from step 4 for remaining points.


Here's the code:

#include &lt;iostream&gt;
#include &lt;cmath&gt;
#include &lt;algorithm&gt;
#include &lt;vector&gt;
#include &lt;set&gt;
#include &lt;cassert&gt;

struct point {
    double x, y;

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

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

double dist(point const &amp;a, point const &amp;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&lt;point, point&gt; min_dist(std::vector&lt;point&gt; points) {

    std::sort(points.begin(), points.end(),
        [](point const &amp;a, point const &amp;b) {
        return a.x &lt; 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&lt;point&gt; band{ *first, *last };

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

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

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

        assert(std::distance(begin, end) &lt;= 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&lt;point&gt; points{
        {1, 1},
        {17, 9},
        {23, 23},
        {3, 3},
        {100, 100},
        {200, 200},
        {24, 24},
        {300, 300}
    };

    auto r = min_dist(points);

    std::cout &lt;&lt; "Closest points: " &lt;&lt; r.first &lt;&lt; ", " &lt;&lt; r.second 
            &lt;&lt; ". Distance = "&lt;&lt; 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.
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 &lt; 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 &gt; 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}