Another way to find the nearest neighbor points in a 2D plane
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}