AIOJ

Introduction

Given $n$ points in the two-dimensional plane, find a pair of points with minimum Euclidean distance.

We introduce a divide-and-conquer algorithm with time complexity $O(n\log n)$ for this problem. The algorithm was proposed by Franco P. Preparata in 1975, and Preparata and Michael Ian Shamos proved that it is optimal in the decision-tree model.

Procedure

As in a conventional divide-and-conquer algorithm, we split the set of $n$ points into two equally sized sets $S_1, S_2$ and recurse. However, we face a difficulty: how should we merge? In other words, how do we find the closest pair with one point in $S_1$ and the other in $S_2$ ? For now, assume that merging takes $O(n)$ time; then the total complexity is $T(n) = 2T(\frac{n}{2}) + O(n) = O(n\log n)$ .

First sort all points by $x_i$ as the primary key and $y_i$ as the secondary key, use point $p_m (m = \lfloor \frac{n}{2} \rfloor)$ as the dividing point, and split the point set into $A_1,A_2$ :

\begin{aligned} A_1 &= \{p_i \ \big | \ i = 0 \ldots m \}\\ A_2 &= \{p_i \ \big | \ i = m + 1 \ldots n-1 \} \end{aligned}

Recurse to find the closest pair within each point set. Let their distances be $h_1,h_2$ , and let $h$ be the smaller one.

Now we merge. We seek a pair of points, one in $A_1$ and one in $A_2$ , whose distance is less than $h$ . Therefore, every point whose horizontal distance from $x_m$ is less than $h$ is put into set $B$ :

B = \{ p_i \ \big | \ \lvert x_i - x_m \rvert < h \}

As shown in the figure, line $m$ divides the points into two parts: the points to the left of $m$ form $A_1$ , and those to the right form $A_2$ .

Applying the rule

B = \{ p_i \ \big | \ \lvert x_i - x_m \rvert < h \}

gives the set

B

formed by the green points.

nearest-points1

For every point in $B$ , say $p_i$ , our goal is to find another point also in $B$ whose distance from it is less than $h$ . To avoid considering each pair twice, we consider only points whose vertical coordinate is less than $y_i$ . Clearly, for a valid point $p_j$ , $y_i - y_j$ must be less than $h$ . Thus we obtain the set $C(p_i)$ :

C(p_i) = \{ p_j\ \big |\ p_j \in B,\ y_i - h < y_j \le y_i \}

Choose a point in $B$ , say $p_i$ . Applying the rule $C(p_i) = \{ p_j\ \big |\ p_j \in B,\ y_i - h < y_j \le y_i \}$ gives the set $C$ consisting of the yellow points inside the red box.

nearest-points2

If we sort the points in $B$ by $y_i$ , then $C(p_i)$ is easy to obtain: it is just the consecutive points immediately adjacent to $p_i$ .

We therefore obtain the merge steps:

Construct set $B$ .
Sort the points in $B$ by $y_i$ . The usual approach takes $O(n\log n)$ , but we can change the strategy and optimize it to $O(n)$ (as explained below).
For each $p_i \in B$ , consider $p_j \in C(p_i)$ ; compute the distance for every pair $(p_i,p_j)$ and update the answer (the closest pair in the current set).

Notice that sorting was mentioned twice above. Since the point coordinates never change, the first sort needs to be performed only once, before the divide-and-conquer begins. Let each recursive call return its current point set sorted by $y_i$ ; for the second sort, the upper level can directly merge the two point sets already sorted by the lower level.

The algorithm still seems inefficient: $|C(p_i)|$ could be $O(n)$ , which would make the total complexity too large. In fact, this is not the case: its maximum size is $7$ . We prove this as follows:

Complexity proof

We know that every point in $C(p_i)$ has vertical coordinate in the range $(y_i-h,y_i]$ ; moreover, all points in $C(p_i)$ and $p_i$ itself have horizontal coordinates in $(x_m-h,x_m+h)$ . Together they form a rectangle of size $2h \times h$ .

Split this rectangle into two $h \times h$ squares. Excluding $p_i$ , the points in one square form $C(p_i) \cap A_1$ , and those in the other form $C(p_i) \cap A_2$ ; the distance between any two points within either square is greater than $h$ (because they come from the same lower-level recursive call).

Split an $h \times h$ square into four $\frac{h}{2} \times \frac{h}{2}$ small squares. Each small square contains at most $1$ point: the greatest possible distance between two points in it is the diagonal length, $\frac{h}{\sqrt 2}$ , which is less than $h$ .

nearest-points3

Thus each square contains at most $4$ points, the rectangle contains at most $8$ points, and after removing $p_i$ itself, $\max(C(p_i))=7$ .

Reference implementation

#include <algorithm>
#include <cstdio>
#include <vector>
using namespace std;

const int N = 500000 + 10;

struct point {
  int x, y, id;
};

int n, A, B;
point a[N];
// mindist is the squared closest distance
long long mindist;

// Update the answer
void upd_ans(const point& a, const point& b) {
  long long dist =
      1LL * (a.x - b.x) * (a.x - b.x) + 1LL * (a.y - b.y) * (a.y - b.y);
  if (dist < mindist) {
    mindist = dist;
    A = a.id;
    B = b.id;
  }
}

// Use [l, r) to represent the current divide-and-conquer interval
void DC(int l, int r) {
  // The current interval contains only one point; return directly
  if (l + 1 == r) return;

  int m = (l + r) >> 1;
  int midx = a[m].x;
  DC(l, m);
  DC(m, r);
  // Use std::inplace_merge() to merge-sort
  inplace_merge(a + l, a + m, a + r,
                [&](point a, point b) { return a.y < b.y; });

  vector<point> t;
  for (int i = l; i < r; i++)
    // Compare squared distances, and do not include equality
    if (1LL * (a[i].x - midx) * (a[i].x - midx) < mindist) t.push_back(a[i]);
  for (int i = 0; i < t.size(); i++)
    for (int j = i + 1; j < t.size(); j++) {
      if (1LL * (t[i].y - t[j].y) * (t[i].y - t[j].y) >= mindist) break;
      upd_ans(t[i], t[j]);
    }
}

void Solve() {
  scanf("%d", &n);
  for (int i = 0; i < n; i++) {
    scanf("%d %d", &a[i].x, &a[i].y);
    a[i].id = i;
  }
  // Sort by horizontal coordinate before calling
  sort(a, a + n, [&](point x, point y) { return x.x < y.x; });
  mindist = 9'000'000'000'000'000'000LL;
  DC(0, n);
  printf("%d %d\n", A, B);
}

int main() {
  int T;
  scanf("%d", &T);
  while (T--) Solve();
  return 0;
}

Extension: minimum-perimeter triangle in the plane

The algorithm above has an interesting extension: choose three points from a given set so that the sum of their three pairwise distances is minimized.

The algorithm is largely unchanged. At each step, try to find a triangle whose perimeter is smaller than the current answer $d$ . Put all points whose horizontal distance from $x_m$ is less than $\frac{d}{2}$ into $B$ and try to update the answer. (A triangle with perimeter $d$ has a longest side shorter than $\frac{d}{2}$ .)

Non-divide-and-conquer algorithm

Besides the divide-and-conquer algorithm described above, there is another algorithm with the same time complexity, $O(n \log n)$ , that does not use divide and conquer.

We can use a common idea from sequence processing: for each element, add its contribution with all elements to its left to the answer. The closest-pair problem in the plane can be handled in the same way.

Specifically, sort all points by $x_i$ as the primary key and $y_i$ as the secondary key, and build a multiset keyed by $y_i$ . At each position $i$ , perform the following operations:

Remove from the set every point satisfying $x_i - x_j \ge d$ . Such points can no longer contribute to the answer.
For every point in the set satisfying $\lvert y_i - y_j \rvert < d$ , compute its distance from $p_i$ .
Insert $p_i$ into the set.

Because each point is inserted and removed at most once, inserting and removing points takes $O(n \log n)$ time. The time-complexity proof for counting the answer is similar to the proof for the divide-and-conquer algorithm, and the reader may try it as an exercise.

Reference implementation

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <set>
using namespace std;

constexpr int N = 500000 + 10;

struct point {
  int x, y, id;

  point() {}

  bool operator<(const point& a) const {
    return x < a.x || (x == a.x && y < a.y);
  }
};

int n, A, B;
long long mindist;
point a[N];

void upd_ans(const point& a, const point& b) {
  long long dist =
      1LL * (a.x - b.x) * (a.x - b.x) + 1LL * (a.y - b.y) * (a.y - b.y);
  if (dist < mindist) {
    mindist = dist;
    A = a.id;
    B = b.id;
  }
}

void Solve() {
  scanf("%d", &n);
  for (int i = 0; i < n; i++) {
    scanf("%d %d", &a[i].x, &a[i].y);
    a[i].id = i;
  }
  sort(a, a + n);
  multiset<pair<double, point>> s;
  mindist = 9'000'000'000'000'000'000LL;
  for (int i = 0, l = 0; i < n; i++) {
    for (; l < i && 1LL * (a[i].x - a[l].x) * (a[i].x - a[l].x) >= mindist; l++)
      s.erase(s.find({a[l].y, a[l]}));
    // Be careful about floating-point errors
    for (auto it =
             s.lower_bound({(double)a[i].y - sqrt(mindist) + 1e-6, point()});
         it != s.end() &&
         1LL * (it->first - a[i].y) * (it->first - a[i].y) < mindist;
         it++)
      upd_ans(it->second, a[i]);
    s.insert({a[i].y, a[i]});
  }
  printf("%d %d\n", A, B);
}

int main() {
  int T;
  scanf("%d", &T);
  while (T--) Solve();
  return 0;
}

Expected-linear algorithm

Besides the $O(n \log n)$ method above, there is also an algorithm with expected complexity $O(n)$ .

First shuffle the points randomly, and maintain the answer for the prefix point set. From the first $i - 1$ points, consider finding the answer for point $i$ .

Let the closest-pair distance among the first $i - 1$ points be $s$ . Partition the plane into a grid of cells with side length $s$ , and store the points in each cell using a hash table. Then inspect all points in the nine cells surrounding the cell containing the $i$ -th point and update the answer. Only $O(1)$ points need to be inspected: because the closest-pair distance among the first $i - 1$ points is $s$ , each cell contains at most $4$ points.

If the answer is updated during this process, rebuild the grid; otherwise, do not rebuild it. Among the first $i$ points, the probability that the closest pair contains $i$ is $O\left(\frac{1}{i}\right)$ , while rebuilding the grid costs $O(i)$ , so the expected cost for the $i$ -th point is $O(1)$ . Thus for $n$ points, the expected complexity is $O(n)$ .

Exercises

References and further reading

The divide-and-conquer section of this page is primarily translated from the blog post Нахождение пары ближайших точек and its English translation Finding the nearest pair of points. The Russian version is licensed under Public Domain + Leave a Link; the English version is licensed under CC-BY-SA 4.0.

Zhihu column: Computational geometry — the closest-pair problem

Closest pair of points