分治法解决最近点对问题

问题

给定平面上n个点，找其中的一对点，使得在n个点的所有点对中，该点对的距离最小。严格地说，最接近点对可能多于1对。为了简单起见，这里只限于找其中的一对。

原理（这段为抄袭https://blog.csdn.net/liufeng_king/article/details/8484284）

设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个点。（抄袭结束）

伪代码

输入：按x坐标排列的n(n>=2)个点的集合S={(x1,y1),(x2,y2),...,(xn,yn)}

输出：最近点的距离

1.如果n==2，则返回（x1,y1）和（x2,y2）之间的距离，算法结束；

2.如果n==3，则返回（x1,y1）、（x2,y2）和（x3,y3）之间的最小距离，算法结束；（此步必要在于，若n==3划分之后必然有一半为n==1，导致无法正确执行递归）

3.划分：m==S中各点x坐标的中位数；

4.d1 = 计算{(x1,y1),...,(xm,ym)}的最近对距离；

5.d2 = 计算{(xm,ym),...,(xn,yn)}的最近对距离；

6.d = min(d1,d2);

7.依次考察集合S中的点p（x，y），如果（x<=xm 并且x>=xm-d），则将点p放入集合P1中；如果（x>xm 并且x<=xm+d），则将点p放入集合P2中;

8.将集合P1和P2按y坐标升序排列；

9.对集合P1和P2中的每个点p（x，y），在y坐标区间[y,y+d]内最对取出6个候选点，计算与点p的最近距离d3；

10.返回min{d，d3}；

注：对于第7，8步，由于任何排序算法需要O（nlgn）的复杂度，但这里需要O（n）的复杂度才可使整个算法的复杂度为O（nlgn），因此采用将归并算法融入求解过程中。对第9步，实际情况是如下图所示

c++实现

/*程序： 最近对的距离作者：Moyu
*/
#include<iostream>
#include<vector>
#include<cmath>
#include<algorithm>
#include<fstream>
#include<sstream>using namespace std;struct Point{double x;double y;
};
inline bool Compx(const Point &p1, const Point &p2)
{return p1.x < p2.x;
}
inline bool Compy(const Point &p1, const Point &p2)
{return p1.y < p2.y;
}
inline double Distance(const Point &a, const Point &b)
{return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
void Merge(vector<Point> &v, int lo, int m, int hi)
{vector<Point> vl(v.begin()+lo,v.begin()+m);int i = lo;int j = 0;int k = m;while(i < hi){if(j < vl.size() && (k == hi || vl[j].y <= v[k].y))v[i++] = vl[j++];if(k < hi && (j == vl.size() || vl[j].y > v[k].y))v[i++] = v[k++];}
}
/*函数：点集最近对的距离参数：vx：以x排序的点集
*/
double Closest(vector<Point> &vx, int lo, int hi)
{if(hi - lo == 2){if(vx[lo].y > vx[hi-1].y){swap(vx[lo],vx[hi-1]);}return Distance(vx[lo],vx[hi-1]);}if(hi - lo == 3){sort(vx.begin()+lo,vx.begin()+hi,Compy); double d1 = Distance(vx[lo],vx[lo+1]);double d2 = Distance(vx[lo],vx[hi-1]);double d3 = Distance(vx[lo+1],vx[hi-1]);return min({d1,d2,d3});}int m = (lo + hi) / 2; double mx = vx[m].x;double dl = Closest(vx,lo,m);double dr = Closest(vx,m,hi);double d = min(dl,dr);Merge(vx,lo,m,hi);vector<Point> vp;for(int i = lo; i < hi; ++i){if(abs(vx[i].x - mx) < d)vp.push_back(vx[i]);}for(int i = 0; i < vp.size(); ++i){for(int j = i + 1; j < vp.size(); ++j){if(vp[j].y - vp[i].y >= d)break;else{double dm = Distance(vp[i],vp[j]);if(dm < d)d = dm;}}}return d;
}
int main()
{vector<Point> v;ifstream file("point.txt",ifstream::in);string line;while(getline(file,line)){Point p;stringstream liness(line);liness >> p.x >> p.y;v.push_back(p);}sort(v.begin(),v.end(),Compx);cout << Closest(v,0,v.size()) << endl;return 0;
}

箴言录：

堂堂正正做人，踏踏实实做事。