分治法-最接近點對問題
背景:
計算機應用中經常採用點、圓等簡單的幾何物件表示物理實體,常需要了解其鄰域中其他幾何物件的信息
- 例如:在空中交通控制中,若將飛機作為空間中的一個點來處理,則具有最大碰撞危險的兩架飛機所處的點對,就是這個空間中距離最接近的一對點。
這類問題是計算幾何學中研究的基本問題之一。
我們著重考慮平面上的最接近點對問題。
最接近點對問題:
給定平面上的n個點,找出其中的一對點,使得在n個點組成的所有點對中,該點對的距離最小。
直觀解法:
- 將每一個點與其他n-1個點的距離算出,找出最小距離。
- 時間複雜度:T(n)=n(n-1)/2+n=O(n2),ps:求出n(n-1)/2個點對距離的時間+求出最短距離的時間
分治法:
- 分解:將n個點的集合分成大小近似相等的兩個子集。
- 求解:遞迴地求解兩個子集內部的最接近點對。
- 合併(關鍵問題):從子空間內部最接近點對,和兩個子空間之間的最接近點對中,選擇最接近點對。
1、一維最接近點對問題
演算法思路:
考慮將所給的n個點的集合S分成2個子集S1和S2,每個子集中約有n/2個點,然後在每個子集中遞迴地求其最接近的點對。
關鍵的問題是分治法中的合併步驟:
- 由S1和S2的最接近點對,還要求得原集合S中的最接近點對,因為S1和S2的最接近點對未必就是S的最接近點對。
- 如果這2個點分別在S1和S2中,則對於S1中任一點p,S2中最多隻有n/2個點
- 仍需做n^2/4(n/2*n/2)次計算和比較才能確定S的最接近點對。
- 合併步驟耗時為O(n^2),整個演算法所需計算時間T(n)應滿足: T(n)=2T(n/2)+O(n^2)。
- 因此,問題出在合併步驟耗時太多。
用x軸上某個點m將S劃分為2個子集S1和S2,使得S1={x∈S|x≤m};S2={x∈S|x>m}。
因此,所有p∈S1和q∈S2有p<q,遞迴地在S1和S2上找出其最接近點對{p1,p2}和{q1,q2},並設d=min{|p1-p2|,|q1-q2|}。
S中的最接近點對或者是{p1,p2},或者是{q1,q2},或者是某個{p3,q3},其中p3∈S1且q3∈S2。
如果S的最接近點對是{p3,q3},即|p3-q3|<d,則p3和q3兩者與m的距離不超過d,即|p3-m|<d,|q3-m|<d。
可得,p3∈(m-d,m],q3∈(m,m+d]。
由於在S1中,每個長度為d的半閉區間至多包含一個點(否則必有兩點距離小於d)。
並且m是S1和S2的分割點,因此(m-d,m]中至多包含S中的一個點,且為S1中最大點。
同理,(m,m+d]中也至多包含S中的一個點,且為S2中最小點。
因此,我們用線性時間就能找到區間(m-d,m]和(m,m+d]中所有點,即p3和q3。
按這種分治策略,合併步可在O(n)時間內完成。
還需考慮的一個問題:分割點m的選取,及S1和S2的劃分。
- 選取分割點m的一個基本要求是由此匯出集合S的一個線性分割。
- 即S=S1∪S2 ,S1∩S2=Φ,且S1={x|x≤m};S2={x|x>m}。
- 容易看出,如果選取m=[max(S)+min(S)]/2,可以滿足線性分割的要求。
- 選取分割點後,再用O(n)時間即可將S劃分成S1={x∈S|x≤m}和S2={x∈S|x>m}。
劃分可能出現的問題:造成劃分出的子集S1和S2的不平衡
- 例如在最壞情況下,|S1|=1,|S2|=n-1。
- 由此,產生的最壞情況下所需的計算時間T(n)=T(n-1)+O(n)。
- 該遞迴方程的解為:T(n)=O(n^2)
- 可用分治法中“平衡子問題”的方法來解決,如用S的n個點的座標的中位數來作分割點。
- 用選取中位數的線性時間演算法,可在O(n)時間內確定一個平衡的分割點m。
確定平衡點採用m=[max(S)+min(S)]/2的方法,程式如下:
#include <ctime>
#include <iostream>
using namespace std;
//點對結構體
struct Pair {
float d;//點對距離
float d1, d2;//點對座標
};
float Max(float s[], int p, int q);
float Min(float s[], int p, int q);
template<class Type>
void Swap(Type &x, Type &y);
template<class Type>
int Partition(Type s[], Type x, int l, int r);
Pair Cpair(float s[], int l, int r);
int main() {
srand((unsigned) time(NULL));
float s[] = {20.14, 3.04, 8.85, 31.72, 40.97, 81.27, 41.15, 45.13, 25.5, 81.84, 3.96, 5.18, 30.82, 72.23, 21.13,
57.59, 76.39, 60.28, 87.88, 13.67, 1.22, 7.82, 9.23, 29.09, 30.15};
Pair d;
d = Cpair(s, 0, 24);
cout << endl << "最近點對座標為:(d1:" << d.d1 << ",d2:" << d.d2 << ")";
cout << endl << "這兩點距離為:" << d.d << endl;
return 0;
}
float Max(float s[], int l, int r)//返回s[]中的最大值
{
float s_max = s[l];
for (int i = l + 1; i <= r; i++)
if (s_max < s[i])
s_max = s[i];
return s_max;
}
float Min(float s[], int l, int r)//返回s[]中的最小值
{
float s_min = s[l];
for (int i = l + 1; i <= r; i++)
if (s_min > s[i])
s_min = s[i];
return s_min;
}
template<class Type>
void Swap(Type &x, Type &y) {
Type temp = x;
x = y;
y = temp;
}
template<class Type>
int Partition(Type s[], Type x, int l, int r) {
int i = l - 1, j = r + 1;
while (true) {
while (s[++i] < x && i < r);
while (s[--j] > x);
if (i >= j) {
break;
}
Swap(s[i], s[j]);
}
return j;
}
//返回s[]中的具有最近距離的點對及其距離
Pair Cpair(float s[], int l, int r) {
Pair min_d = {99999, 0, 0};//最短距離
if (r - l < 1) return min_d;
float m1 = Max(s, l, r), m2 = Min(s, l, r);
float m = (m1 + m2) / 2;//找出點集中的中位數
//將點集中的各元素按與m的大小關係分組
int j = Partition(s, m, l, r);
Pair d1 = Cpair(s, l, j), d2 = Cpair(s, j + 1, r);//遞迴
float p = Max(s, l, j), q = Min(s, j + 1, r);
//返回s[]中的具有最近距離的點對及其距離
if (d1.d < d2.d) {
if ((q - p) < d1.d) {
min_d.d = (q - p);
min_d.d1 = q;
min_d.d2 = p;
return min_d;
} else return d1;
} else {
if ((q - p) < d2.d) {
min_d.d = (q - p);
min_d.d1 = q;
min_d.d2 = p;
return min_d;
} else return d2;
}
}
該演算法的分割步驟和合併步驟總共耗時O(n)。因此,演算法耗費的計算時間T(n)滿足遞迴方程:
解得:T(n)=O(nlogn)
2、二維最接近點對問題
①選取二維平面的一條垂直平分線L:x=m作為分割線。
━其中m為S中各點x座標的中位數,由此將S分割為S1={p∈S|px≤m}和S2={p∈S|px>m}。
②遞迴地在S1和S2上找出其最小距離d1和d2。
━現設d=min(d1,d2)。若S的最接近點對(p,q)之間的距離d(p,q)<d,則p必∈S1和q必∈S2。
③如果用符號P1和P2分別表示直線 L 的左右兩邊寬為d的區域,則p∈P1,q∈P2。
問題分析:
- 在一維的情形,距分割點距離為d的2個區間(m-d,m](m,m+d]中最多各有S中一個點。
- 因而這2個點成為唯一的末檢查過的最接近點對候選者。
- 在二維的情形,在最壞情況下有n2/4對這樣的候選者。
- 但是,P1和P2中的點具有稀疏性質,它使我們不必檢查所有這n^2/4對候選者。
問題的優化方案:
①考慮P1中的任意一點,它若與P2中的點q構成最接近點對的候選者,則必有:distance(p,q)<d。
②P2中滿足條件的點一定落在矩形R中,矩形R的大小為:d×2d。
③由d的定義可知:P2中任何2個點(qi∈S)的距離都不小於d,由此可以推出矩形R中最多隻有6個S中的點。
④因此,在分治法的合併步驟中最多隻需要檢查6×n/2=3n個候選者。
數學證明:
(一)合並時最多隻需檢查6×n/2=3n個候選者(矩形R中最多隻有6個S中的點)
補充知識:
抽屜原理:把m個元素任意放入n(n<m)個集合,則一定有一個集合呈至少要有k個元素。
其中 k= [m/n]([]表示向上取整)。
- 將矩形R長為2d的邊3等分,將長為d的邊2等分。
- 由此匯出6個(d/2)×(2d/3)的小矩形,如圖。
- 若矩形R中有多於6個S中的點,由抽屜原理易知,至少有一個小矩形中有2個以上S中的點
- 設u,v是位於同一小矩形中的2個點,則,
- 即:distance(u,v)<d,這與d的意義相矛盾(P2內不可能再出現距離<d的兩個點),命題得證。
⑤如何確定需要檢查的6個點?
- 可以將p和P2中所有S2的點投影到垂直線L上。
- 由於能與p點一起構成最接近點對候選者的S2中的點一定在矩形R中。
- 所以它們在直線 L 上的投影點與 p 在 L 上投影點的距離小於d。
- 根據上述分析,這種投影點最多隻有6個。因此,若將區域P1和P2中所有S中的點按其y座標排好序。
- 則:對P1中的所有點,只需一次掃描就可以找出所有候選者。
- 對排好序的點作一次掃描,可以找出所有最接近點對的候選者。
- 對P1中每個點,最多只需檢查P2中排好序的相繼6個點。
#include<time.h>
#include<iostream>
#include<cmath>
using namespace std;
const int M = 50;
//用類PointX和PointY表示依x座標和y座標排好序的點
class PointX {
public:
int operator<=(PointX a) const { return (x <= a.x); }
int ID; //點編號
float x, y; //點座標
};
class PointY {
public:
int operator<=(PointY a) const { return (y <= a.y); }
int p; //同一點在陣列x中的座標
float x, y; //點座標
};
float Random();
template<class Type>
float dis(const Type &u, const Type &v);
bool Cpair2(PointX X[], int n, PointX &a, PointX &b, float &d);
void closest(PointX X[], PointY Y[], PointY Z[], int l, int r, PointX &a, PointX &b, float &d);
template<typename Type>
void Copy(Type a[], Type b[], int left, int right);
template<class Type>
void Merge(Type c[], Type d[], int l, int m, int r);
template<class Type>
void MergeSort(Type a[], Type b[], int left, int right);
int main() {
srand((unsigned) time(NULL));
int length;
cout << "請輸入點對數:";
cin >> length;
PointX X[M];
cout << "隨機生成的二維點對為:" << endl;
for (int i = 0; i < length; i++) {
X[i].ID = i;
X[i].x = Random();
X[i].y = Random();
cout << "(" << X[i].x << "," << X[i].y << ") ";
}
PointX a;
PointX b;
float d;
Cpair2(X, length, a, b, d);
cout << endl;
cout << "最鄰近點對為:(" << a.x << "," << a.y << ")和(" << b.x << "," << b.y << ") " << endl;
cout << "最鄰近距離為: " << d << endl;
return 0;
}
float Random() {
float result = rand() % 10000;
return result * 0.01;
}
//平面上任意兩點u和v之間的距離可計算如下
template<class Type>
inline float dis(const Type &u, const Type &v) {
float dx = u.x - v.x;
float dy = u.y - v.y;
return sqrt(dx * dx + dy * dy);
}
bool Cpair2(PointX X[], int n, PointX &a, PointX &b, float &d) {
if (n < 2) return false;
PointX *tmpX = new PointX[n];
MergeSort(X, tmpX, 0, n - 1);
PointY *Y = new PointY[n];
for (int i = 0; i < n; i++) //將陣列X中的點複製到陣列Y中
{
Y[i].p = i;
Y[i].x = X[i].x;
Y[i].y = X[i].y;
}
PointY *tmpY = new PointY[n];
MergeSort(Y, tmpY, 0, n - 1);
PointY *Z = new PointY[n];
closest(X, Y, Z, 0, n - 1, a, b, d);
delete[]Y;
delete[]Z;
delete[]tmpX;
delete[]tmpY;
return true;
}
void closest(PointX X[], PointY Y[], PointY Z[], int l, int r, PointX &a, PointX &b, float &d) {
if (r - l == 1) //兩點的情形
{
a = X[l];
b = X[r];
d = dis(X[l], X[r]);
return;
}
if (r - l == 2) //3點的情形
{
float d1 = dis(X[l], X[l + 1]);
float d2 = dis(X[l + 1], X[r]);
float d3 = dis(X[l], X[r]);
if (d1 <= d2 && d1 <= d3) {
a = X[l];
b = X[l + 1];
d = d1;
return;
}
if (d2 <= d3) {
a = X[l + 1];
b = X[r];
d = d2;
} else {
a = X[l];
b = X[r];
d = d3;
}
return;
}
//多於3點的情形,用分治法
int m = (l + r) / 2;
int f = l, g = m + 1;
//在演算法預處理階段,將陣列X中的點依x座標排序,將陣列Y中的點依y座標排序
//演算法分割階段,將子陣列X[l:r]均勻劃分成兩個不想交的子集,取m=(l+r)/2
//X[l:m]和X[m+1:r]就是滿足要求的分割。
for (int i = l; i <= r; i++) {
if (Y[i].p > m) Z[g++] = Y[i];
else Z[f++] = Y[i];
}
closest(X, Z, Y, l, m, a, b, d);
float dr;
PointX ar, br;
closest(X, Z, Y, m + 1, r, ar, br, dr);
if (dr < d) {
a = ar;
b = br;
d = dr;
}
Merge(Z, Y, l, m, r);//重構陣列Y
//d矩形條內的點置於Z中
int k = l;
for (int i = l; i <= r; i++) {
if (fabs(X[m].x - Y[i].x) < d) {
Z[k++] = Y[i];
}
}
//搜尋Z[l:k-1]
for (int i = l; i < k; i++) {
for (int j = i + 1; j < k && Z[j].y - Z[i].y < d; j++) {
float dp = dis(Z[i], Z[j]);
if (dp < d) {
d = dp;
a = X[Z[i].p];
b = X[Z[j].p];
}
}
}
}
template<class Type>
void Merge(Type c[], Type d[], int l, int m, int r) {
int i = l, j = m + 1, k = l;
while ((i <= m) && (j <= r)) {
if (c[i] <= c[j]) {
d[k++] = c[i++];
} else {
d[k++] = c[j++];
}
}
if (i > m) {
for (int q = j; q <= r; q++) {
d[k++] = c[q];
}
} else {
for (int q = i; q <= m; q++) {
d[k++] = c[q];
}
}
}
template<class Type>
void MergeSort(Type a[], Type b[], int left, int right) {
if (left < right) {
int i = (left + right) / 2;
MergeSort(a, b, left, i);
MergeSort(a, b, i + 1, right);
Merge(a, b, left, i, right);//合併到陣列b
Copy(a, b, left, right);//複製回陣列a
}
}
template<typename Type>
void Copy(Type a[], Type b[], int left, int right) {
for (int i = left; i <= right; i++)
a[i] = b[i];
}