使用分治法和蠻力法求解最近點對
問題描述
對於平面上給定的N個點，給出所有點對的最短距離，即輸入是平面上的N個點，輸出是N點中具有最短距離的兩點。求解
建立點類，使之具有兩個屬性，x座標和y座標
class myPoint {
public:
int x; //x座標
int y; //y座標
};
比較函式
bool compare(myPoint a, myPoint b, int type) { //若type=1，表示比較x座標；若type=2，表示比較y座標。
if (type == 1) {
return a.x > b.x;
}
else 隨機生成點集，放在X集和Y集中。注意：X集和Y集中點的內容是一樣的，只是排序的方法不同而已，X集中的點按照x座標進行排序，Y集中的點按照y座標進行排序
void Initial(int N, myPoint* X, myPoint* Y) {
srand(N); //種子為資料的規模
int i;
for (i = 0; i < N; i++) { //產生隨機數，X集和Y集中點的內容是一樣的，只是排序的方法不同而已
X[i].x = rand();
X[i].y = rand 距離函式
double Distance(myPoint p1, myPoint p2) { //計算兩個點的距離
return sqrt((p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y));
}快排演算法(使用分治法前，必須先排序)
void quickSort(myPoint* P, int start, int end, int type) { //若type=1，表示比較x座標；若type=2 蠻力法
假設輸入規模為n，則應該求出每個點與其他n-1個點之間的距離，並將它與最小距離做比較，若小於最小距離，則更新最小距離。則總共應該執行操作n(n-1)次，故演算法的複雜度為O(n2)。
double force(int start, int end, myPoint* P, myPoint& P1, myPoint& P2) { //暴力法
int i, j;
if (end - start<1) { //只有一個點時，直接返回
return 0.0;
}
double minDis = Distance(P[start], P[start + 1]); //初始化
P1 = P[start];
P2 = P[start+1];
for (i = start; i <= end; i++) {
for (j = start; j <= end; j++) {
if (i != j && Distance(P[i], P[j])<minDis) {
minDis = Distance(P[i], P[j]);
P1 = P[i];
P2 = P[j];
}
}
}
return minDis;
}2.分治法
（1） 預排序：
將點集存放到陣列X中，陣列X中的點集按照x座標進行排序；將點集存放到陣列Y中，陣列Y中的點集按照y座標進行排序。陣列X和陣列Y中總體的內容是相同的，都包含整個點集的X座標和Y座標，只是排序不一樣而已。排序過程呼叫快速排序演算法，其演算法複雜度是O(n logn)
（2） 判斷：
在分治演算法中，若判斷到點集的數量小於等於3，則直接呼叫蠻力法進行排序；若點集格式大於3，則採用分治的思想來求解。
（3） 分解：
找出一條垂直線L,將點集按照x座標的順序的不同分成兩部分。由於之前做了預排序，所以取陣列的中位數，就可以直接點集近乎平分成兩半。
（4） 解決：
將上一步分成兩半的子點集，遞迴呼叫分治演算法對問題規模進行分解。當遞迴分解到符合點集的數目小於等於3後，通過蠻力法求解最小距離並且組成返回。此時可以得到左邊的最小距離minL，右邊的最小距離minR。選取min（minL，minR）作為當前的最小的距離s。
（5） 合併：
判斷左、右兩部分子點集中，是不是存在一個點a在左邊的子點集，另外一個點b在右邊的子點集，使得a,b之間的距離小於之前求得的最小距離s。那麼左邊的n/2 個點分別與右邊的n/2個點進行比較，就要比較n2/4次，使得演算法複雜度為O(n2)，這使得與蠻力法沒多大差別，我們尋找更好的解決方法。我們嘗試縮小搜尋的範圍，由於當前的最小距離是s，所以在距離中間線的寬度為s之外的點，不可能從對邊找到一個距離小於s的點。接著我們繼續縮小搜尋的範圍，對於一個點p，若要找到一個與它的距離不大於s的點q，則q的縱座標與p的縱座標的差值必須不大於s。所以我們得到了一個可能的滿足條件的矩形區域。
如上圖所示，若點p為圖中的粉紅色小點，且點a可能的存在範圍為藍線區域。則綠色的區域所圍成的邊長為2S的矩形為q點可能的存在範圍。由互相排斥的性質可得，極端情況下，這個區域裡面最多總共有12個點，圖中的每個紫色的位置最多有一個點，總共6個點；圖中的黃色點位置，每個位置可能有兩個點，一個屬於左邊，一個屬於右邊，總共6個點。所以加起來總共有12個點。現在，我們一開始對y座標進行排序的陣列Y可以派上用場了。遍歷陣列Y，若某個點落在到中線距離小於s的長條形區域上，也就是圖中的兩條紅線所夾區域，則加入到新陣列C中。所以，在合併的過程中，複雜度為O（n）。有因為迭代了logn層遞迴呼叫，所以分治過程演算法複雜度為O（nlogn）。又因為預排序的演算法複雜度為O（nlogn），所以整個分治演算法的複雜度依舊為O（nlogn）
我們做一條直線橫切綠色矩形，將綠色矩形分割為兩半。直線上面的點即屬於上半矩形又屬於下半矩形，所以每半個綠矩形現在最多有八個點。我們遍歷陣列C中的所有點，判斷每個點後面的7個點是否與自身的距離小於s，若小於s則更新s。只判斷後面7個點而不判斷前面的點是因為當使有的點都判斷後面7個點就足夠複雜所有情況了。
double divide_conquer(int start, int end, myPoint* X, myPoint* Y, myPoint& P1, myPoint& P2) {
if (end - start < 3) //當點的數量小於等於3時，用暴力法。由於是按照x座標來劃分區域，固X陣列的點的數量和位置君未發生改變，故用X陣列不用Y陣列
return force(start, end, X,P1,P2);
int mid = (start + end) / 2; //按照x座標的排序，從中間切開
int leftLen = 0, rightLen = 0, i, j; //leftLen，和rightLen分別指示左右兩邊點集的個數
myPoint* YL = new myPoint[(end - start + 1) ]; //儲存中線左邊的點集，按照Y座標排序
myPoint* YR = new myPoint[(end - start + 1) ]; //儲存中線右邊的點集，按照Y座標排序
for (i = 0; i <= end - start; i++) { //遍歷集合中每一個點，他們都是按照y的順序進行排序的，若點的x座標小於等於中線，則分割到左子集，否則分割到右子集
//通過此操作YL和YR中的點的座標依舊是按照y座標排序好的
if (Y[i].x <= X[mid].x) {
YL[leftLen++] = Y[i];
}
else {
YR[rightLen++] = Y[i];
}
}
double left = divide_conquer(start, mid, X, YL,P1,P2); //遞迴求左邊部分的最短距離
myPoint leftP1 = P1; //儲存左支距離最小的兩個點的座標
myPoint leftP2 = P2;
double right = divide_conquer(mid + 1, end, X, YR,P1,P2); //遞迴求右邊部分的最短距離
double minDis; //取小的一個
if (left < right) {
minDis = left;
P1 = leftP1;
P2 = leftP2;
}
else { //此時並不需要更新P1，P2的值，因為此時P1和P2的值就是右支的最小的一對點的座標
minDis = right;
}
myPoint* newY = new myPoint[(end - start + 1)]; //新建一個數組，用於儲存以中界線為中心，寬度為 2*minDis 的垂直帶形區域內的點
int newYLen = 0;
double leftBorder = X[mid].x - minDis; //區域的左邊界
double rightBorder = X[mid].x + minDis; //區域的右邊界
for (i = 0; i <= end - start; i++) { //遍歷Y集中所有點，滿足條件則加入新集合，新集合的點的依然按照y座標排序好
if (Y[i].x >= leftBorder && Y[i].x <= rightBorder)
newY[newYLen++] = Y[i];
}
for (i = 0; i<newYLen; i++) { //在新集合中，判斷每個點跟它後面的7個點的距離作比較，若小於則更新最小距離
for (j = 1; j <= 7; j++) {
if ((i + j)<newYLen) { //加入條件，防止越界
double dis = Distance(newY[i], newY[i + j]);
if (dis < minDis) {
minDis = dis;
P1 = newY[i];
P2 = newY[i + j];
}
}
}
}
delete YL;
delete YR;
delete newY;
return minDis;
}
void Initial(int N, myPoint* X, myPoint* Y) {
srand(N); //種子為資料的規模
int i;
for (i = 0; i < N; i++) { //產生隨機數，X集和Y集中點的內容是一樣的，只是排序的方法不同而已
X[i].x = rand();
X[i].y = rand();
Y[i].x = X[i].x;
Y[i].y = X[i].y;
}
}
值得注意的是，在對中間區域的點按照y座標進行排序的過程中，必須充分利用之前有序的陣列Y，在O(n)的複雜度下產生新的陣列。之前在網上參考了相關的過程，發現有一些演算法就沒有解決好這麼一個問題。他們對陣列Y中抽取出來的點進行快速排序，那麼合併過程的複雜度就變成O(nlogn)，再加上logn 層遞迴呼叫，使得演算法的總體複雜度變成了O（n（logn）2），並不能達到我們想要達到的O（nlogn）的複雜度。
最後給出main函式和執行結果
#include<iostream>
#include<math.h>
#include<stdlib.h>
#include<time.h>
const int Count = 100; //重複計算的次數，解決小規模時時間太短無法得出時間的問題
int main()
{
int N;
cout << "請輸入問題的規模：";
cin >> N;
myPoint* X = new myPoint[N]; //X,Y兩個點集的內容相同，區別在於X點集是按照x座標順序排序，而Y點集市安裝y座標順序排序
myPoint* Y = new myPoint[N];
int i;
clock_t start, end;
double ave = 0.0;
double minDis = 0.0;
myPoint P1, P2; //距離最小的兩個點的座標
Initial(N, X, Y); //初始化陣列
for (i = 0; i < Count; i++) { //重複count次
start = clock();
minDis += force(0, N - 1, X,P1,P2);
end = clock();
ave += (double)(end - start);
}
ave /= Count;
minDis /= Count;
cout << "暴力演算法: 最短距離為：" << minDis << " 時間為：" <<ave<<" ms"<< endl;
cout << "最小距離的兩個點的座標為：(" << P1.x << "," << P1.y << "),(" << P2.x << "," << P2.y << ")"<<endl;
ave = 0.0;
minDis = 0.0;
for (i = 0; i < Count; i++) { //重複count次
start = clock();
quickSort(X, 0, N - 1, 1); //將點集按照x座標的升序排序，並儲存到X陣列中
quickSort(Y, 0, N - 1, 2); //將點集按照y座標的升序排序，並儲存到Y陣列中
minDis+= divide_conquer(0, N - 1, X, Y,P1, P2);
end = clock();
Initial(N, X, Y); //重新初始化，供下一次使用
ave += (double)(end - start);
}
ave /= Count;
minDis /= Count;
cout << "分治演算法: 最短距離為：" <<minDis<< " 時間為：" << ave << " ms" << endl;
cout << "最小距離的兩個點的座標為：(" << P1.x << "," << P1.y << "),(" << P2.x << "," << P2.y << ")" << endl;
delete[]X;
delete[]Y;
return 0;
}執行結果：
