Given a collection of points on a 2D plane, find the pair that is closest to each other

Full disclosure: I'm working on this for an online course. However, my goal is really just to get a pointer to where the issue is.

The goal is to implement the closest points problem, that is, given a set of points on a 2D plane, find the shortest distance between two points. After lots of stress-testing and debugging, I am confident the algorithm is correct. However, it is not fast enough, which is the problem I want to solve.

My algorithm is the implementation of what is described in that page, and includes the following optimisations:

1. Sort the arrays before passing them to the function;
2. Keep track of minimum distance so far and stop if it equals 0;
3. When considering the strip in the middle, only calculate the distance between a point and the next one if the distance between their x-coordinates and the distance between their y-coordinates is smaller than the minimum distance found so far;

The main parts of the code are

double min_distance(const vector<Point> points_x, const vector<Point> points_y, double current_delta) {
// Base case - sets of 3 points or fewer
if (points_x.size() <= 3) return brute_force(points_x);
if (current_delta == 0) return current_delta;

int mid = points_x.size() / 2;
int mid_x = points_x[mid].x;

vector<Point> x_left;
vector<Point> x_right;
vector<Point> y_left;
vector<Point> y_right;

// Creates the x_left and x_right arrays
for (int i = 0; i < points_x.size(); i++) {
if (i < mid) {
x_left.push_back(points_x[i]);
} else if (i == mid) {
x_left.push_back(points_x[i]);
x_right.push_back(points_x[i]);
} else {
x_right.push_back(points_x[i]);
}
}
// Creates the y_left and y_right arrays
for (int i = 0; i < points_y.size(); i++) {
if (points_y[i].x < mid_x) {
y_left.push_back(points_y[i]);
} else if (points_y[i].x == mid_x) {
y_left.push_back(points_y[i]);
y_right.push_back(points_y[i]);
} else {
y_right.push_back(points_y[i]);
}
}

// Recursively solve left
double min_left = min_distance(x_left, y_left, current_delta);
if (min_left == 0) return min_left;
// Recursively solve left
double min_right = min_distance(x_right, y_right, current_delta);
if (min_right == 0) return min_right;

double delta = 0;
if (min_left == -1 && min_right != -1) {
delta = min_right;
} else if (min_left != -1 && min_right == -1) {
delta = min_left;
} else {
delta = min(min_left, min_right);
}

// Find values of the x's that determine the borders on the strip
double min_x = mid_x - delta;
double max_x = mid_x + delta;
// Create array y with all points in the strip
vector<Point> y_strip;
// Creates the y_strip sorted by its y coordinate
for (int i = 0; i < points_y.size(); i++) {
// Add all the points inside the strip
if (points_y[i].x >= min_x && points_y[i].x <= max_x) {
y_strip.push_back(points_y[i]);
}
}
// Find the minimum distance in the strip
double min_strip = mid_min_distance(y_strip, delta);
return min(delta, min_strip);
}


for recursively calculating the smallest distance on the left and right sides, and

double mid_min_distance(const vector<Point> y_strip, double delta) {
// If mid_region is empty or contains just 1 point
if (y_strip.size() < 2) return delta;

double mid_min = minimal_distance(y_strip[0], y_strip[1]);
double mid_min_distance = mid_min;
// Brute force to inner points
for (int i = 0; i < y_strip.size(); i++) {
for (int j = i + 1; j < y_strip.size(); j++) {
// If the vertical/horizontal distance between the points
// is greater than delta, break the loop
if (std::abs(y_strip[i].y - y_strip[j].y) < delta &&
std::abs(y_strip[i].x - y_strip[j].x) < delta) {
mid_min_distance = minimal_distance(y_strip[i], y_strip[j]);
mid_min = min(mid_min_distance, mid_min);
} else {
break;
}
}
}
return mid_min;
}


to calculate the distance between points that lie on the strip surrounding the middle x value. Despite seeming O(n^2), that should not be the case since there is a small limit on the number of times the inner loop actually runs. Tracing the execution of the program with the debugger, this seems to be working correctly.

This should run in less than 2 seconds for the input range I must consider but my code is taking 4 seconds. This means I'm missing something big and thus I just wanted to get a sense at where to look.

EDIT: Taking into account some answers here, I re-wrote parts of the code (thoroughly stress-tested the new version - appears to be correct)

double min_distance(const vector<Point> &points_x, int x_left, int x_right,
const vector<Point> &points_y, double &current_delta) {
// No distance is lower than zero so halt
if (current_delta == 0) return current_delta;
// Base case
int current_size = x_right - x_left;
if (current_size < 4) return brute_force(points_x, x_left, x_right);

int mid = x_left + ((x_right - x_left) / 2);
int mid_x = points_x[mid].x;
// Create the y vectors from previously sorted points_y
vector<Point> y_left;
vector<Point> y_right;
int i = 0;
int y_size = points_y.size();
for (i = 0; i < y_size; i++) {
if (points_y[i].x < mid_x) {
y_left.push_back(points_y[i]);
} else if (points_y[i].x == mid_x) {
y_left.push_back(points_y[i]);
y_right.push_back(points_y[i]);
} else {
y_right.push_back(points_y[i]);
}
}
// Recursively solve left and right
double min_left = min_distance(points_x, x_left, mid, y_left, current_delta);
if (min_left == 0) return min_left;
double min_right = min_distance(points_x, mid, x_right, y_right, current_delta);
if (min_right == 0) return min_right;

if (min_left == -1 && min_right != -1) {
current_delta = min_right;
} else if (min_left != -1 && min_right == -1) {
current_delta = min_left;
} else if (min_left != -1 && min_right != -1) {
current_delta = min(min_left, min_right);
}
// Create y-sorted array with all points in the strip
vector<Point> y_strip;
for (Point p : points_y) {
if (abs(p.x - mid_x) <= current_delta) {
y_strip.push_back(p);
}
}
// Find the minimum distance in the strip region
// and return the lowest of the two
return strip_min_distance(y_strip, current_delta);
}

double strip_min_distance(const vector<Point> &y_strip, double &current_delta) {
int size = y_strip.size();
if (size < 2) return current_delta;

double mid_min = euclidean_distance(y_strip[0], y_strip[1]);
double mid_min_distance = mid_min;
// Brute for to inner points
// Points whose vertical distance exceeds current_delta can
// be ignored
int i = 0, j = 0;
for (i = 0; i < size; i++) {
for (j = i + 1; j < size; j++) {
if (abs(y_strip[i].y - y_strip[j].y) < current_delta) {
mid_min_distance = euclidean_distance(y_strip[i], y_strip[j]);
mid_min = min(mid_min_distance, mid_min);
} else {
break;
}
}
}
current_delta = min(current_delta, mid_min);
return current_delta;
}

• Welcome to Code Review! Please avoid using unqualified word "problem", since our code reviews only correctly working code. I believe your concern is performance, so you could mention that instead of problem. Also, I see "in that page", but there is no link. May be you forgot to embed it? – Incomputable Mar 20 '17 at 21:45
• A brief description of what the code is supposed to be doing would be useful. That is, what is the code supposed to do? "The closest points problem" is not specific enough for us to be able to guess. – Edward Mar 20 '17 at 21:48
• I believe what you're trying to say is, "given a collection of points on a 2D plane, find the pair that is closest to each other"? – 200_success Mar 20 '17 at 21:50

Consider using range-based for loop

If you don't need information about index of element, then you can iterate over the vector of points with following code:

for (Point p : points)


By using range-based for loop and std::abs to calculate the borders of the strip you can get the following code for finding all points in the strip (which in my opinion is more readable):

vector<Point> y_strip;

for (Point p : points_y) {
if (std::abs(p.x - mid_x) <= delta) {
y_strip.push_back(p)
}
}


Additionally, by using std::abs you don't have to create temporary variables min_x and max_x.

Performance of push_back

At the beginning you create 4 vectors and fill them with push_back function which might result in memory reallocation. Consider using reserve function after initialization of each vector to avoid unnecessary reallocation.

Algorith correctness

I think that the following condition in mid_min_distance function is unnecessary:

std::abs(y_strip[i].x - y_strip[j].x) < delta


It prevents your algorith from checking all possible pairs of points in the strip as it breaks nested for loop too early. Test your algorithm for the following sets of points(in all examples correct answer is $\sqrt{13} \approx 3.60555$):

1. $(-8, 1), (-1, 1), (0, -6), (1, 4), (4, 1)$
2. $(-4, 1), (-1, 4), (0, -6), (1, 1), (8, 1)$
3. $(-8, -1), (-1, -1), (0, 6), (1, -4), (4, -1)$
4. $(-4, -1), (-1, -4), (0, 6), (1, -1), (8, -1)$

Additionally, the following code won't work properly when all points have the same x coordinate:

// Creates the y_left and y_right arrays
for (int i = 0; i < points_y.size(); i++) {
if (points_y[i].x < mid_x) {
y_left.push_back(points_y[i]);
} else if (points_y[i].x == mid_x) {
y_left.push_back(points_y[i]);
y_right.push_back(points_y[i]);
} else {
y_right.push_back(points_y[i]);
}
}


Both y_left and y_right will contain all the points from points_y vector.

• I re-wrote the code from scratch to take into account the feedback here (some things I haven't yet implemented, i.e the push_back / reserve). The 2 tests you suggested are returning the correct value. The point regarding the x coordinates all being the same makes sense, I have to think about how to go around that. Still, the algorithm is 3.98 seconds to run in a test where the max is 2 seconds so I feel like something big is still going unnoticed as the suggestions given here seem like marginal improvements. – mcansado Mar 23 '17 at 0:12
• The last point you made was precisely the issue. A case with 100000 points, where they all share the x-coordinate. I've changed the code now and it worked. Thank you! – mcansado Mar 23 '17 at 7:41
• You're welcome, I'm glad to see that my remarks were helpful! – adasikow Mar 23 '17 at 19:43
• One more remark about algorithm correctness - in my first two examples I wrongly assumed that points in y_strip are sorted by y coordinate ascending. I've added two more examples in my answer to show that original algorithm (with std::abs(y_strip[i].x - y_strip[j].x) < delta condition in mid_min_distance function) will return incorrect value for one of these four examples. – adasikow Mar 23 '17 at 20:06
• But they are sorted by y-ascending. In the main function (not shown as it didn't seem strictly necessary for the question), 2 different arrays are created: points_x and points_y. Each is sorted in ascending order of their respective coordinate. Since y_strip is being created from points_y, it'll be sorted in y_ascending order. – mcansado Mar 24 '17 at 10:03

Comments out of sync with code

  // Recursively solve left
double min_left = min_distance(x_left, y_left, current_delta);
if (min_left == 0) return min_left;

// Recursively solve left
double min_right = min_distance(x_right, y_right, current_delta);
if (min_right == 0) return min_right;


Pretty sure the second comment should be "recursively solve right"--or omitted entirely, since it's pretty easily deduced from the code.

Parameter passing

Right now, you're passing a couple of std::vector parameters by value. You almost certainly want to pass by reference instead. I'd also at least consider passing either iterators or (perhaps) something range-like, such as a pair of iterators signifying a range in a collection.

Using this, creating the left and right "arrays" might look something like this:

range x_left{0, x_mid};
range x_right{x_mid, x_end};

range y_left{0, y_mid};
range y_right{mid, y_end};


The big difference here is that we're just keeping track of positions in the original vectors, rather than copying the entire contents of each original vector every time we make a recursive call (which is compounded when you pass by value, so you end up making not just one, but two copies of each sub-array to do your recursive call).

[[probably more to follow]]

• Absolutely right, will change the code and question. – mcansado Mar 20 '17 at 23:04
• @mcansado , please don't do that. It is an established practice that once you posted the question and get an answer, your code in the question is set in stone. Otherwise it would invalidate answers. You can ask a follow up question if you want to get review on improved version. – Incomputable Mar 20 '17 at 23:26
• Undid it, that makes sense actually - first question here. I'll try and implement those changes and see how that works. – mcansado Mar 21 '17 at 7:16