平面上最近点对

        在二维平面上的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拆分左右两部分为S_L和S_R，L一般取点集S中所有点的中间点的x坐标来划分，这样可以保证S_L和S_R中的点数目各为n/2，

（否则以其他方式划分S，有可能导致S_L和S_R中点数目一个为1，一个为n-1，不利于算法效率，要尽量保持树的平衡性）

依次找出这两部分中的最小点对距离：δ_L和δ_R，记S_L和S_R中最小点对距离δ = min（δ_L，δ_R），如图1：

 

   

     以L为中心，δ为半径划分一个长带，最小点对还有可能存在于S_L和S_R的交界处，如下图2左图中的虚线带，p点和q点分别位于S_L和S_R的虚线范围内，在这个范围内，p点和q点之间的距离才会小于δ，最小点对计算才有意义。

    

 

Figure 2

 

 

      对于S_L虚框范围内的p点，在S_R虚框中与p点距离小于δ的顶多只有六个点，就是图二右图中的2个正方形的6的顶点。这个可以反推证明，如果右边这2个正方形内有7个点与p点距离小于δ，例如q点，则q点与下面正方形的四个顶点距离小于δ，则和δ为S_L和S_R中的最小点对距离

相矛盾。因此对于S_L虚框中的p点，不需求出p点和右边虚线框内所有点距离，只需计算S_R中与p点y坐标距离最近的6个点，就可以求出最近点对，节省了比较次数。

（否则的话，最坏情形下，在S_R虚框中有可能会有n/2个点，对于S_L虚框中的p点，每次要比较n/2次，浪费了算法的效率）

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>/*核心是分治算法1. 分别根据点的 x,y 值进行排序2. 在 x 轴上划一道垂线, 将点均分成两半3. 假设最近点对都在左/右部分, 递归计算左/右半部分的最短距离并返回较小值 dis4. 假设最近点对分别在左右两个部分, 横跨中心的竖线. 中心线为中心, 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;
}