轉自：http://blog.csdn.net/liufeng_king/article/details/8484284

問題場景：在應用中，常用諸如點、圓等簡單的幾何對象代表現實世界中的實體。在涉及這些幾何對象的問題中，常需要了解其鄰域中其他幾何對象的信息。例如，在空中交通控制問題中，若將飛機作為空間中移動的一個點來看待，則具有最大碰撞危險的2架飛機，就是這個空間中最接近的一對點。這類問題是計算幾何學中研究的基本問題之一。

問題描述：給定平面上n個點，找其中的一對點，使得在n個點的所有點對中，該點對的距離最小。嚴格地說，最接近點對可能多於1對。為了簡單起見，這裏只限於找其中的一對。

1、一維最接近點對問題

算法思路：

這個問題很容易理解，似乎也不難解決。我們只要將每一點與其他n-1個點的距離算出，找出達到最小距離的兩個點即可。然而，這樣做效率太低，需要O(n^2)的計算時間。在問題的計算復雜性中我們可以看到，該問題的計算時間下界為Ω(nlogn)。這個下界引導我們去找問題的一個θ(nlogn)算法。采用分治法思想，考慮將所給的n個點的集合S分成2個子集S1和S2，每個子集中約有n/2個點，然後在每個子集中遞歸地求其最接近的點對。在這裏，一個關鍵的問題是如何實現分治法中的合並步驟，即由S1和S2的最接近點對，如何求得原集合S中的最接近點對，因為S1和S2的最接近點對未必就是S的最接近點對。如果組成S的最接近點對的2個點都在S1中或都在S2中，則問題很容易解決。但是，如果這2個點分別在S1和S2中，則對於S1中任一點p，S2中最多只有n/2個點與它構成最接近點對的候選者，仍需做n^2/4次計算和比較才能確定S的最接近點對。因此，依此思路，合並步驟耗時為O(n^2)。整個算法所需計算時間T(n)應滿足:　T(n)=2T(n/2)+O(n^2)。它的解為T(n)=O(n^2)，即與合並步驟的耗時同階，這不比用窮舉的方法好。從解遞歸方程的套用公式法，我們看到問題出在合並步驟耗時太多。這啟發我們把註意力放在合並步驟上。

設S中的n個點為x軸上的n個實數x1,x2,..,xn。最接近點對即為這n個實數中相差最小的2個實數。我們顯然可以先將x1,x2,..,xn排好序，然後，用一次線性掃描就可以找出最接近點對。這種方法主要計算時間花在排序上，在排序算法已經證明，時間復雜度為O(nlogn)。然而這種方法無法直接推廣到二維的情形。因此，對這種一維的簡單情形，我們還是嘗試用分治法來求解，並希望能推廣到二維的情形。假設我們用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中的一個點。同理，(m,m+d]中也至多包含S中的一個點。由圖可以看出，如果(m-d,m]中有S中的點，則此點就是S1中最大點。同理，如果(m,m+d]中有S中的點，則此點就是S2中最小點。因此，我們用線性時間就能找到區間(m-d,m]和(m,m+d]中所有點，即p3和q3。從而我們用線性時間就可以將S1的解和S2的解合並成為S的解。也就是說，按這種分治策略，合並步可在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}。然而，這樣選取分割點m，有可能造成劃分出的子集S1和S2的不平衡。例如在最壞情況下，|S1|=1，|S2|=n-1，由此產生的分治法在最壞情況下所需的計算時間T(n)應滿足遞歸方程:
T(n)=T(n-1)+O(n)
它的解是T(n)=O(n^2)。這種效率降低的現象可以通過分治法中“平衡子問題”的方法加以解決。即通過適當選擇分割點m，使S1和S2中有大致相等個數的點。自然地，我們會想到用S的n個點的坐標的中位數來作分割點。在選擇算法中介紹的選取中位數的線性時間算法使我們可以在O(n)時間內確定一個平衡的分割點m。

本程序確定平衡點采用m=[max(S)+min(S)]/2方法。如果需要利用中位數作分割點，看結合筆者博文《0005算法筆記——線性時間選擇》改寫。

一維最接近臨近點對問題程序清單如下：

//2d10-1 一維最鄰近點對問題  
#include "stdafx.h"  
#include <ctime>  
#include <iostream>   
using namespace std;   
  
const int L=100;  
//點對結構體  
struct Pair  
{  
    float d;//點對距離  
    float d1,d2;//點對坐標  
};  
float Random();  
int input(float s[]);//構造S  
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));  
    int m;  
    float s[L];  
    Pair d;  
    m=input(s);  
    d=Cpair(s,0,m-1);  
    cout<<endl<<"最近點對坐標為： (d1:"<<d.d1<<",d2:"<<d.d2<<")";  
    cout<<endl<<"這兩點距離為： "<<d.d<<endl;  
    return 0;  
}  
  
  
float Random()  
{  
    float result=rand()%10000;  
     return result*0.01;  
}  
  
int input(float s[])  
{  
    int length;  
    cout<<"輸入點的數目： ";  
    cin>>length;  
    cout<<"點集在X軸上坐標為：";  
    for(int i=0;i<length;i++)  
    {  
        s[i]=Random();  
        cout<<s[i]<<" ";  
    }  
      
    return length;  
}  
  
  
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、二維最接近點對問題

將以上過程推廣到二維最接近點對問題，設S中的點為平面上的點，它們都有2個坐標值x和y。為了將平面上點集S線性分割為大小大致相等的2個子集S1和S2，我們選取一垂直線l:x=m來作為分割直線。其中m為S中各點x坐標的中位數。由此將S分割為S1={p∈S|px≤m}和S2={p∈S|px>m}。從而使S1和S2分別位於直線l的左側和右側，且S=S1∪S2 。由於m是S中各點x坐標值的中位數，因此S1和S2中的點數大致相等。遞歸地在S1和S2上解最接近點對問題，我們分別得到S1和S2中的最小距離d1和d2。現設d=min(d1,d2)。若S的最接近點對(p,q)之間的距離d(p,q)<d則p和q必分屬於S1和S2。不妨設p∈S1，q∈S2。那麽p和q距直線l的距離均小於d。因此，我們若用P1和P2分別表示直線l的左邊和右邊的寬為d的2個垂直長條，則p∈S1，q∈S2，如圖所示:

距直線l的距離小於d的所有點

在一維的情形，距分割點距離為d的2個區間(m-d,m](m,m+d]中最多各有S中一個點。因而這2點成為唯一的末檢查過的最接近點對候選者。二維的情形則要復雜些，此時，P1中所有點與P2中所有點構成的點對均為最接近點對的候選者。在最壞情況下有n2/4對這樣的候選者。但是P1和P2中的點具有以下的稀疏性質，它使我們不必檢查所有這n^2/4對候選者。考慮P1中任意一點p,它若與P2中的點q構成最接近點對的候選者，則必有d(p,q)<d。滿足這個條件的P2中的點有多少個呢？容易看出這樣的點一定落在一個d×2d的矩形R中，如下圖所示:

包含點q的dX2d矩形R

由d的意義可知P2中任何2個S中的點的距離都不小於d。由此可以推出矩形R中最多只有6個S中的點。事實上，我們可以將矩形R的長為2d的邊3等分，將它的長為d的邊2等分，由此導出6個（d/2）×（2d/3）的矩形。如左圖所示:

矩陣R中點的稀疏性

若矩形R中有多於6個S中的點，則由鴿舍原理易知至少有一個δ×2δ的小矩形中有2個以上S中的點。設u,v是這樣2個點，它們位於同一小矩形中，則：

因此d(u,v)≤5d/6<d 。這與d的意義相矛盾。也就是說矩形R中最多只有6個S中的點。圖4(b)是矩形R中含有S中的6個點的極端情形。由於這種稀疏性質，對於P1中任一點p，P2中最多只有6個點與它構成最接近點對的候選者。因此，在分治法的合並步驟中，我們最多只需要檢查6×n/2=3n對候選者，而不是n^2/4對候選者。這是否就意味著我們可以在O(n)時間內完成分治法的合並步驟呢？現在還不能作出這個結論，因為我們只知道對於P1中每個S1中的點p最多只需要檢查P2中的6個點，但是我們並不確切地知道要檢查哪6個點。為了解決這個問題，我們可以將p和P2中所有S2的點投影到垂直線l上。由於能與p點一起構成最接近點對候選者的S2中點一定在矩形R中，所以它們在直線l上的投影點距p在l上投影點的距離小於d。由上面的分析可知，這種投影點最多只有6個。因此，若將P1和P2中所有S的點按其y坐標排好序，則對P1中所有點p，對排好序的點列作一次掃描，就可以找出所有最接近點對的候選者，對P1中每一點最多只要檢查P2中排好序的相繼6個點。

程序清單如下：

//2d10-2 二維最鄰近點對問題  
#include "stdafx.h"  
#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];   
}  

程序結果：

【分治法】最接近點對問題（轉）