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

Question

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;
}

Answer 1

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.

Answer 2

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]]