跳转至

平面最近點對

引入

給定 \(n\) 個二維平面上的點，求一組歐幾里得距離最近的點對。

下面我們介紹一種時間複雜度為 \(O(n\log n)\) 的分治算法來解決這個問題。該算法在 1975 年由 Franco P. Preparata 提出，Preparata 和 Michael Ian Shamos 證明了該算法在決策樹模型下是最優的。

過程

與常規的分治算法一樣，我們將這個有 \(n\) 個點的集合拆分成兩個大小相同的集合 \(S_1, S_2\)，並不斷遞歸下去。但是我們遇到了一個難題：如何合併？即如何求出一個點在 \(S_1\) 中，另一個點在 \(S_2\) 中的最近點對？這裏我們先假設合併操作的時間複雜度為 \(O(n)\)，可知算法總複雜度為 \(T(n) = 2T(\frac{n}{2}) + O(n) = O(n\log n)\)。

我們先將所有點按照 \(x_i\) 為第一關鍵字、\(y_i\) 為第二關鍵字排序，並以點 \(p_m (m = \lfloor \frac{n}{2} \rfloor)\) 為分界點，拆分點集為 \(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} \]

並遞歸下去，求出兩點集各自內部的最近點對，設距離為 \(h_1,h_2\)，取較小值設為 \(h\)。

現在該合併了！我們試圖找到這樣的一組點對，其中一個屬於 \(A_1\)，另一個屬於 \(A_2\)，且二者距離小於 \(h\)。因此我們將所有橫座標與 \(x_m\) 的差小於 \(h\) 的點放入集合 \(B\)：

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

結合圖像，直線 \(m\) 將點分成了兩部分。\(m\) 左側為 \(A_1\) 點集，右側為為 \(A_2\) 點集。

再根據 \(B = \{ p_i \ \big | \ \lvert x_i - x_m \rvert < h \}\) 規則，得到綠色點組成的 \(B\) 點集。nearest-points1

對於 \(B\) 中的每個點 \(p_i\)，我們當前目標是找到一個同樣在 \(B\) 中、且到其距離小於 \(h\) 的點。為了避免兩個點之間互相考慮，我們只考慮那些縱座標小於 \(y_i\) 的點。顯然對於一個合法的點 \(p_j\)，\(y_i - y_j\) 必須小於 \(h\)。於是我們獲得了一個集合 \(C(p_i)\)：

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

在點集 \(B\) 中選一點 \(p_i\)，根據 \(C(p_i) = \{ p_j\ \big |\ p_j \in B,\ y_i - h < y_j \le y_i \}\) 的規則，得到了由紅色方框內的黃色點組成的 \(C\) 點集。

nearest-points2

如果我們將 \(B\) 中的點按照 \(y_i\) 排序，\(C(p_i)\) 將很容易得到，即緊鄰 \(p_i\) 的連續幾個點。

由此我們得到了合併的步驟：

  1. 構建集合 \(B\)。
  2. 將 \(B\) 中的點按照 \(y_i\) 排序。通常做法是 \(O(n\log n)\)，但是我們可以改變策略優化到 \(O(n)\)（下文講解）。
  3. 對於每個 \(p_i \in B\) 考慮 \(p_j \in C(p_i)\)，對於每對 \((p_i,p_j)\) 計算距離並更新答案（當前所處集合的最近點對）。

注意到我們上文提到了兩次排序，因為點座標全程不變，第一次排序可以只在分治開始前進行一次。我們令每次遞歸返回當前點集按 \(y_i\) 排序的結果，對於第二次排序，上層直接使用下層的兩個分別排序過的點集歸併即可。

似乎這個算法仍然不優，\(|C(p_i)|\) 將處於 \(O(n)\) 數量級，導致總複雜度不對。其實不然，其最大大小為 \(7\)，我們給出它的證明：

複雜度證明

我們已經瞭解到，\(C(p_i)\) 中的所有點的縱座標都在 \((y_i-h,y_i]\) 範圍內；且 \(C(p_i)\) 中的所有點，和 \(p_i\) 本身，橫座標都在 \((x_m-h,x_m+h)\) 範圍內。這構成了一個 \(2h \times h\) 的矩形。

我們再將這個矩形拆分為兩個 \(h \times h\) 的正方形，不考慮 \(p_i\)，其中一個正方形中的點為 \(C(p_i) \cap A_1\)，另一個為 \(C(p_i) \cap A_2\)，且兩個正方形內的任意兩點間距離大於 \(h\)。（因為它們來自同一下層遞歸）

我們將一個 \(h \times h\) 的正方形拆分為四個 \(\frac{h}{2} \times \frac{h}{2}\) 的小正方形。可以發現，每個小正方形中最多有 \(1\) 個點：因為該小正方形中任意兩點最大距離是對角線的長度，即 \(\frac{h}{\sqrt 2}\)，該數小於 \(h\)。

nearest-points3

由此，每個正方形中最多有 \(4\) 個點，矩形中最多有 \(8\) 個點，去掉 \(p_i\) 本身，\(\max(C(p_i))=7\)。

實現

我們使用一個結構體來存儲點，並定義用於排序的函數對象：

結構體定義 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
struct pt {
  int x, y, id;
};

struct cmp_x {
  bool operator()(const pt& a, const pt& b) const {
    return a.x < b.x || (a.x == b.x && a.y < b.y);
  }
};

struct cmp_y {
  bool operator()(const pt& a, const pt& b) const { return a.y < b.y; }
};

int n;
vector<pt> a;

為了方便實現遞歸，我們引入 upd_ans() 輔助函數來計算兩點間距離並嘗試更新答案：

答案更新函數 
1
2
3
4
5
6
7
8
double mindist;
int ansa, ansb;

void upd_ans(const pt& a, const pt& b) {
  double dist =
      sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + .0);
  if (dist < mindist) mindist = dist, ansa = a.id, ansb = b.id;
}

下面是遞歸本身：假設在調用前 a[] 已按 \(x_i\) 排序。如果 \(r-l\) 過小，使用暴力算法計算 \(h\)，終止遞歸。

我們使用 std::inplace_merge() 來執行歸併排序，並創建輔助緩衝區 t[]，\(B\) 存儲在其中。

主體函數 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
void rec(int l, int r) {
  if (r - l <= 3) {
    for (int i = l; i <= r; ++i)
      for (int j = i + 1; j <= r; ++j) upd_ans(a[i], a[j]);
    sort(a + l, a + r + 1, &cmp_y);
    return;
  }

  int m = (l + r) >> 1;
  int midx = a[m].x;
  rec(l, m), rec(m + 1, r);
  inplace_merge(a + l, a + m + 1, a + r + 1, &cmp_y);

  static pt t[MAXN];
  int tsz = 0;
  for (int i = l; i <= r; ++i)
    if (abs(a[i].x - midx) < mindist) {
      for (int j = tsz - 1; j >= 0 && a[i].y - t[j].y < mindist; --j)
        upd_ans(a[i], t[j]);
      t[tsz++] = a[i];
    }
}

在主函數中，這樣開始遞歸即可：

調用接口 
1
2
3
sort(a, a + n, &cmp_x);
mindist = 1E20;
rec(0, n - 1);

推廣：平面最小周長三角形

上述算法有趣地推廣到這個問題：在給定的一組點中，選擇三個點，使得它們兩兩的距離之和最小。

算法大體保持不變，每次嘗試找到一個比當前答案周長 \(d\) 更小的三角形，將所有橫座標與 \(x_m\) 的差小於 \(\frac{d}{2}\) 的點放入集合 \(B\)，嘗試更新答案。（周長為 \(d\) 的三角形的最長邊小於 \(\frac{d}{2}\)）

非分治算法

過程

其實，除了上面提到的分治算法，還有另一種時間複雜度同樣是 \(O(n \log n)\) 的非分治算法。

我們可以考慮一種常見的統計序列的思想：對於每一個元素，將它和它的左邊所有元素的貢獻加入到答案中。平面最近點對問題同樣可以使用這種思想。

具體地，我們把所有點按照 \(x_i\) 為第一關鍵字、\(y_i\) 為第二關鍵字排序，並建立一個以 \(y_i\) 為第一關鍵字、\(x_i\) 為第二關鍵字排序的 multiset。對於每一個位置 \(i\)，我們執行以下操作：

  1. 將所有滿足 \(x_i - x_j >= d\) 的點從集合中刪除。它們不會再對答案有貢獻。
  2. 對於集合內滿足 \(\lvert y_i - y_j \rvert < d\) 的所有點，統計它們和 \(p_i\) 的距離。
  3. 將 \(p_i\) 插入到集合中。

由於每個點最多會被插入和刪除一次，所以插入和刪除點的時間複雜度為 \(O(n \log n)\)，而統計答案部分的時間複雜度證明與分治算法的時間複雜度證明方法類似，讀者不妨一試。

實現

參考代碼 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <set>
const int N = 200005;
int n;
double ans = 1e20;

struct point {
  double x, y;

  point(double x = 0, double y = 0) : x(x), y(y) {}
};

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

struct cmp_y {
  bool operator()(const point &a, const point &b) const { return a.y < b.y; }
};

void upd_ans(const point &a, const point &b) {
  double dist = sqrt(pow((a.x - b.x), 2) + pow((a.y - b.y), 2));
  if (ans > dist) ans = dist;
}

point a[N];
std::multiset<point, cmp_y> s;

int main() {
  scanf("%d", &n);
  for (int i = 0; i < n; i++) scanf("%lf%lf", &a[i].x, &a[i].y);
  std::sort(a, a + n, cmp_x());
  for (int i = 0, l = 0; i < n; i++) {
    while (l < i && a[i].x - a[l].x >= ans) s.erase(s.find(a[l++]));
    for (auto it = s.lower_bound(point(a[i].x, a[i].y - ans));
         it != s.end() && it->y - a[i].y < ans; it++)
      upd_ans(*it, a[i]);
    s.insert(a[i]);
  }
  printf("%.4lf", ans);
  return 0;
}

期望線性做法

其實，除了上面提到的時間複雜度為 \(O(n \log n)\) 的做法，還有一種 期望 複雜度為 \(O(n)\) 的算法。

首先將點對 隨機打亂，我們將維護前綴點集的答案。考慮從前 \(i - 1\) 個點求出第 \(i\) 個點的答案。

記前 \(i - 1\) 個點的最近點對距離為 \(s\)，我們將平面以 \(s\) 為邊長劃分為若干個網格，並存下每個網格內的點（使用 哈希表）， 然後檢查第 \(i\) 個點所在網格的周圍九個網格中的所有點，並更新答案。注意到需檢查的點的個數是 \(O(1)\) 的，因為前 \(i - 1\) 個點的最近點對距離為 \(s\)， 從而每個網格不超過 4 個點。

如果這一過程中，答案被更新，我們就重構網格圖，否則不重構。在前 \(i\) 個點中，最近點對包含 \(i\) 的概率為 \(O\left(\frac{1}{i}\right)\)， 而重構網格的代價為 \(O(i)\)，從而第 \(i\) 個點的期望代價為 \(O(1)\)。於是對於 \(n\) 個點，該算法期望為 \(O(n)\)。

習題

參考資料與拓展閲讀

本頁面中的分治算法部分主要譯自博文 Нахождение пары ближайших точек 與其英文翻譯版 Finding the nearest pair of points。其中俄文版版權協議為 Public Domain + Leave a Link；英文版版權協議為 CC-BY-SA 4.0。

知乎專欄：計算幾何 - 最近點對問題

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用