平面最近點對問題(分治)

 

平面最近點對問題是指：在給出的同一個平面內的所有點的坐標，然后找出這些點中最近的兩個點的距離.

 

方法1：窮舉

 

1）算法描述：已知集合S中有n個點，一共可以組成n(n-1)/2對點對，蠻力法就是對這n(n-1)/2對點對逐對進行距離計算，通過循環求得點集中的最近點對
2）算法時間復雜度：算法一共要執行 n(n-1)/2次循環，因此算法復雜度為O(n²)

代碼實現：

利用兩個for循環可實現所有點的配對,每次配對算出距離然后更新最短距離.

 

for (i=0 ; i < n ;i ++)｛
     for(j= i+1 ; j<n ;j ++)｛
              點i與點j的配對
       ｝
｝

　　 

方法2：分治

 

1) 把它分成兩個或多個更小的問題；

 2) 分別解決每個小問題； 

3) 把各小問題的解答組合起來，即可得到原問題的解答。小問題通常與原問題相似，可以遞歸地使用分而治之策略來解決。

 

在這里介紹一種時間復雜度為O(nlognlogn)的算法。其實，這里用到了分治的思想。將所給平面上n個點的集合S分成兩個子集S1和S2，每個子集中約有n/2個點。然后在每個子集中遞歸地求最接近的點對。在這里，一個關鍵的問題是如何實現分治法中的合並步驟，即由S1和S2的最接近點對，如何求得原集合S中的最接近點對。如果這兩個點分別在S1和S2中，問題就變得復雜了。

為了使問題變得簡單，首先考慮一維的情形。此時，S中的n個點退化為x軸上的n個實數x1，x2，...，xn。最接近點對即為這n個實數中相差最小的兩個實數。顯然可以先將點排好序，然后線性掃描就可以了。但我們為了便於推廣到二維的情形，嘗試用分治法解決這個問題。

假設我們用m點將S分為S1和S2兩個集合，這樣一來，對於所有的p(S1中的點)和q(S2中的點)，有p<q。

遞歸地在S1和S2上找出其最接近點對{p1,p2}和{q1,q2}，並設

　　　　　　　　　　　　　　　　　　d = min{ |p1-p2| , |q1-q2| } 　　　　　　　　　　　　　　　　　　

由此易知，S中最接近點對或者是{p1,p2}，或者是{q1,q2}，或者是某個{q3,p3}，如下圖所示。

    如果最接近點對是{q3,p3}，即|p3-q3|<d，則p3和q3兩者與m的距離都不超過d，且在區間(m-d,d]和(d,m+d]各有且僅有一個點。這樣，就可以在線性時間內實現合並。

    此時，一維情形下的最近點對時間復雜度為O(nlogn)。

 

在二維的情況下:

我們仿照一維的情況先把所有點按照x(橫坐標)從左到右升序排列.

以X橫坐標中間的點作為分界線.將平面的點分成左邊和右邊,以上圖為例,分為左邊5個右邊5個.

然后找到左邊的點中最近點對的距離d1,和右邊最近點對的距離d2。

令d=min{d1,d2}.那么d是當前整個平面的最近點對的距離嗎?

答案不是!!

 

當前的d1,d2都是兩邊各自的點的最近距離,但是沒有考慮到兩邊的點相互配對的情況,即是點對一個在左邊區域,一個在右邊區域.

如果我們這個時候把兩邊相互配對的所有情況全部羅列出來那么時間復雜度就變為第一種方法的O(n²).

這個時候我們便可以用上面找到的d來限制配對.即是明顯超過d的兩點不需要配對.如下：

那么在范圍內的只有三個點,也就是說只有這個三個點相互配對的距離才可能比d小.那么再按照方法一一樣在這些點中找到最短距離再跟d去比較.小的就是當前整個平面中所考慮的所有點中最近點對的距離。

然而在以上問題處理上,有一個問題尚未解決就是如何找到兩邊區域中的最近點對的距離d1,d2 ?

我們可以發現在處理上面這個問題的時候，和最開始處理所有平面的點最近點距離的問題類似，只是點的數目減半了.

那么我們可以遞歸以上問題.直到划分的區域中只有一個或者兩個點.這樣,該區域的最近點對距離就是無窮大或者該兩點的距離。這也是遞歸終止的條件。

 

代碼實現：

#include <iostream>  
#include <cstdio>  
#include <cstring>  
#include <cmath>  
#include <algorithm>  
using namespace std;
const double inf = 1e20;
const int maxn = 100005;

struct Point{
    double x, y;
}point[maxn];

int n, mpt[maxn];


//以x為基准排序
bool cmpxy(const Point& a, const Point& b){
    if (a.x != b.x)
        return a.x < b.x;
    return a.y < b.y;
}

bool cmpy(const int& a, const int& b){
    return point[a].y < point[b].y;
}

double min(double a, double b){
    return a < b ? a : b;
}

double dis(int i, int j){
    return sqrt((point[i].x - point[j].x)*(point[i].x - point[j].x) + (point[i].y - point[j].y)*(point[i].y - point[j].y));
}

double Closest_Pair(int left, int right){
    double d = inf;
    if (left == right)
        return d;
    if (left + 1 == right)
        return dis(left, right);
    int mid = (left + right) >> 1;
    double d1 = Closest_Pair(left, mid);
    double d2 = Closest_Pair(mid + 1, right);
    d = min(d1, d2);
    int i, j, k = 0;
    //分離出寬度為d的區間  
    for (i = left; i <= right; i++){
        if (fabs(point[mid].x - point[i].x) <= d)
            mpt[k++] = i;
    }
    sort(mpt, mpt + k, cmpy);
    //線性掃描  
    for (i = 0; i < k; i++){
        for (j = i + 1; j < k && point[mpt[j]].y - point[mpt[i]].y<d; j++){
            double d3 = dis(mpt[i], mpt[j]);
            if (d > d3)    d = d3;
        }
    }
    return d;
}

int main(){
    while (~scanf("%d", &n) && n){
        for (int i = 0; i < n; i++)
            scanf("%lf %lf", &point[i].x, &point[i].y);
        sort(point, point + n, cmpxy);
        printf("%.2lf\n", Closest_Pair(0, n - 1) / 2);
    }
    return 0;
}