在二维平面上的n个点中,如何快速的找出最近的一对点,就是最近点对问题。
一种简单的想法是暴力枚举每两个点,记录最小距离,显然,时间复杂度为O(n^2)。
在这里介绍一种时间复杂度为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)。
在二维情形下,类似的,利用分治法,但是难点在于如何实现线性的合并?
由上图可见,形成的宽为2d的带状区间,最多可能有n个点,合并时间最坏情况下为n^2,。但是,P1和P2中的点具有以下稀疏的性质,对于P1中的任意一点,P2中的点必定落在一个d X 2d的矩形中,且最多只需检查六个点(鸽巢原理)。
这样,先将带状区间的点按y坐标排序,然后线性扫描,这样合并的时间复杂度为O(nlogn),几乎为线性了。
算法描述:已知集合S中有n个点,分治法的思想就是将S进行拆分,分为2部分求最近点对。算法每次选择一条垂线L,将S拆分左右两部分为SL和SR,L一般取点集S中所有点的中间点的x坐标来划分,这样可以保证SL和SR中的点数目各为n/2,
(否则以其他方式划分S,有可能导致SL和SR中点数目一个为1,一个为n-1,不利于算法效率,要尽量保持树的平衡性)
依次找出这两部分中的最小点对距离:δL和δR,记SL和SR中最小点对距离δ = min(δL,δR),如图1:
以L为中心,δ为半径划分一个长带,最小点对还有可能存在于SL和SR的交界处,如下图2左图中的虚线带,p点和q点分别位于SL和SR的虚线范围内,在这个范围内,p点和q点之间的距离才会小于δ,最小点对计算才有意义。
Figure 2
对于SL虚框范围内的p点,在SR虚框中与p点距离小于δ的顶多只有六个点,就是图二右图中的2个正方形的6的顶点。这个可以反推证明,如果右边这2个正方形内有7个点与p点距离小于δ,例如q点,则q点与下面正方形的四个顶点距离小于δ,则和δ为SL和SR中的最小点对距离
相矛盾。因此对于SL虚框中的p点,不需求出p点和右边虚线框内所有点距离,只需计算SR中与p点y坐标距离最近的6个点,就可以求出最近点对,节省了比较次数。
(否则的话,最坏情形下,在SR虚框中有可能会有n/2个点,对于SL虚框中的p点,每次要比较n/2次,浪费了算法的效率)
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*
核心是分治算法
1. 分别根据点的 x,y 值进行排序
2. 在 x 轴上划一道垂线, 将点均分成两半
3. 假设最近点对都在左/右部分, 递归计算左/右半部分的最短距离
并返回较小值 dis
4. 假设最近点对分别在左右两个部分, 横跨中心的竖线. 中心线为中心,
2*dis 为宽度画一个矩形, 横跨中心线的最近点对 candidate 都在这个矩形内.
将这些点按照 y 值的大小加入到数组中. 遍历数组中的点, 将该点与其后的
7 个点计算距离, 返回最小距离
5. 为什么仅和 7 个点作对比呢. 因为已经假设 dis 是左右不分最近点对的最小值,
这就说明在一个长(宽)为 dis 的正方形内, 至多有 4 个点. 长为 dis*2,
宽为 dis 的长方形至多 8 个.
*/
const double INF = 1e20;//10^20,科学计数方法
const int N = 100005;
struct Point
{
double x;
double y;
}point[100005]; //point用来存储输入的坐班点
int n;
int tmpt[100005];//tmpt[]用来保存,中间区域点的顺序
int cmpxy(const void *a, const void *b) //对所有的点进行横坐标升序排序
{ //横坐标相同的按纵坐标升序
Point *c=(Point *)a;
Point *d=(Point *)b;
if(c->x!= d->x)
return c->x-d->x;
return c->y-d->y;
}
int cmpy(const void *a, const void * b) //中间区域排序
{ //对距离mid点横向距离少于d的点,进行纵坐标升序排序
int c=*(int*)a, d=*(int*)b;
return point[c].y - point[d].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)/2; //取中间位置
double d1 = Closest_Pair(left,mid);//求左边最小距离
double d2 = Closest_Pair(mid+1,right);//求右边最小距离
d = min(d1,d2);
// printf("d1=%.4f d2=%.4f d01=%.4f\n",d1,d2,d);
int i,j,k=0;
//分离出距离中心横向宽度为d的点区间
for(i=left;i<=right;i++)
{
if(fabs(point[mid].x-point[i].x)<=d)
tmpt[k++]=i;
}
qsort(tmpt,k,sizeof(int),cmpy);//分离出来的区间点,纵坐标进升序排序
// for(i=0;i<k;i++) printf("tmpt[%d]=%d ",i,tmpt[i]);
// printf("分界\n");
//线性扫描,求解中间位置两侧的最小两点距离
for(i = 0;i<k;i++)
{
for(j=i+1;j<k;j++)//理解!!
{
if(point[tmpt[j]].y-point[tmpt[i]].y<d)
{
double d3 = dis(tmpt[i],tmpt[j]);
// printf("d3=%.4f\n",d3);
if(d > d3)
d = d3;
}
}
}
// printf("d02=%.4f\n",d);
return d;
}
int main()
{
freopen("nearest.in","r",stdin);
freopen("nearest.out","w",stdout);
scanf("%d",&n);
for(int i = 0; i < n; i++)
scanf("%lf%lf",&point[i].x,&point[i].y);
qsort(point,n,sizeof(Point),cmpxy);
/* for(int i = 0; i < n; i++)
printf("%lf %lf\n",point[i].x,point[i].y);
*/
printf("%.4lf\n",Closest_Pair(0,n-1));
return 0;
}