查找平面见最近点对

题目：

给定平面上n个点，找其中的一对点，使得在n个点组成的所有点对中，该点对间的距离最小。

算法核心思想：

严格的讲，最接近点对可能多于1对，为简单起见，只找其中的1对作为问题的解。简单的说，只要将每一点与其它n-1个点的距离算出，找出达到最小距离的2点即可。但这样效率太低，故想到分治法来解决这个问题。也就是说，将所给的平面上n个点的集合S分成2个子集S1和S2，每个子集中约有n/2个点。然后在每个子集中递归的求其最接近的点对。这里，关键问题是如何实现分治法中的合并步骤，即由S1和S2的最接近点对，如何求得原集合S中的最接近点对。如果组成S的最接近点对的2个点都在S1中或都在S2中，则问题很容易解决，但如果这2个点分别在S1和S2中，问题就不那么简单了。

S中的点为平面上的点，它们都有两个坐标值x和y。为了将平面上点集S线形分割为大小大致相等的两个子集S1和S2，选取一垂直线l：x=m来作为分割直线。其中m为S中各点x坐标的中位数。由此将S分割为S1={p∈S|x（p）<=m}和S2={p∈S|x（p）>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必分属于S1和S2。不妨设p∈S1，q∈S2。p和q距离直线l的距离均小于d。因此，若用P1和P2分别表示直线l的左边和右边的宽为d的两个垂直长条。则p∈P1，q∈P2，如下图所示。



在二维条件下，P1中所有点与P2中所有点构成的点对均为最接近点对的候选者。由d的意义可知，P2中任何两个S中的点的距离都不小于d。由此可推出矩阵R中最多只有6个S中的点。由此稀疏性质，对于P1中任一点p，P2中最多只有6个点与它构成最接近点对的候选者。因此，在分治法的合并步骤中，最多只需要检查6*n/2=3*n个候选者。但并不确切知道要检查哪6个点。为解决这问题，可以将p和P2中所有S2的点投影到垂直线l上。由于能与p点一起构成最接近点对候选者的S2中的点一定在d*2d的矩形中，所以它们在直线l上的投影点距p在l上投影点的距离小于d。由上述分析可知，这种投影点最多有6个。因此，若将P1和P2中所有S中点按其y坐标排好序，则对P1中所有点，对排好序的点列作一次扫描，就可以找出所有最接近点对的候选者。对P1中每一点最多只要检查P2中排好序的相继6个点。

算法步骤如下：

Step1:将点对按x坐标排序；

Step2:如果|S|=1,d=-1;

Step3:如果|S|=2,d=这两个点的距离；

Step4:如果|S|=3,d=min{这三个点两两之间的距离};

Step5:mid=S中各点x坐标的中位数；

Step6:ans=min{递归的按中位数划分的左右两个子区间求出的最段距离}；

Step7:Q用来记录所有在S[mid]左右各ans范围内的坐标点，并且按y轴排序，由前面的分析可知这样的点最多有7个；

Step8:求出Q中距离最短的两个点d1,ans=min{d1,ans};

代码如下：

// cloest.cpp : Defines the entry point for the console application.
//

//#include "stdafx.h"
#include <cstdio>
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <limits>
#include<algorithm>
using namespace std;
#define N 5000
struct point
{
double x , y;
}p
;
point q
;//保存筛选的坐标点
point close[2];//保存最近点对

bool cmpx(const point &c , const point &b)
{
if(c.x==b.x) return c.y<=b.y;
else return c.x < b.x;
}
bool cmpy(const point &c , const point &b)
{
return c.y <= b.y;
}

double square(double x)
{
return x*x;
}
double dis(point &a , point &b)//求两个点之间的距离
{
return sqrt( (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}
double min(double a , double b)
{
return a < b ? a : b;
}
double GetThreePoint(int low,int high)
{
static double help = numeric_limits<double>::max()/10;//该语句是关键,核心在Static！
double calcu;
int i,j;
for (i = low;i<high;i++)
{
for (j = i + 1;j<=high;j++)
{
calcu = dis(p[i],p[j]);
if (calcu<help)
{
help = calcu;
close[0] = p[i];
close[1] = p[j];
//cout<<"help = "<<help<<endl;
}
}
}
return help;
}
double closest(int low , int high)//分治法求最近点对
{
double d2,d1,d3;
static double ans = numeric_limits<double>::max()/10;//关键在Static！
if(high-low<=2)//只有小于等于3个点的情况
{
return GetThreePoint(low,high);
}
int mid = (low + high)>>1;
ans = min(ans,min( closest(low , mid) , closest(mid + 1 , high) ));    //分治法进行递归求解
if(ans == -1) return ans;
int i , j , cnt = 0;
for(i = low ; i <= high ; ++i)   //把x坐标在p[mid].x-ans~p[mid].x+ans范围内的点取出来
{
if(p[i].x >= p[mid].x - ans && p[i].x <= p[mid].x + ans)
{
q[cnt++] = p[i];       //保存在mid点左右各ans距离的点
}
}
sort(q,q+cnt,cmpy);//按y坐标进行升序排序
for(i = 0 ; i < cnt-1 ; ++i)//求在mid点左右各ans距离的点中的最近点对
{
for (j = i + 1; j < cnt && q[j].y - q[i].y < ans; ++j) {
d3 = dis(q[i], q[j]);
//cout<<"ans = "<<ans<<endl;
//cout<<"d3 = "<<d3<<endl;
if (d3 <=ans) {
ans = d3;
close[0] = q[i];
close[1] = q[j];
//cout<<"q[j] = "<<q[j].x<<" "<<q[j].y<<endl;
}
}
}
return ans;
}
void qongju(point p[],int n)//穷举求出最短点对，用于比较分治法的结果是否正确
{
double len,num;
num = 100000;
int n1,n2;
int i,j;
for(i=0;i<n;i++)
for(j=i+1;j<n;j++)
{
len = dis(p[i],p[j]);
if(len<num)
{
num = len;
n1 = i;
n2 = j;
}
}
printf("穷举法求最近点对距离为：%.1lf\n",num);
printf("最近点对为：[%.0lf,%.0lf],[%.0lf,%.0lf]\n",p[n1].x,p[n1].y,p[n2].x,p[n2].y);
}
int main(void)
{
int i,n;
printf("请输入点的对数：");
scanf("%d",&n);
for(i=0;i<n;i++)
{
p[i].x = rand()%(n+20);
p[i].y = rand()%(n+20);
}
sort(p ,p+n, cmpx);//按x坐标排序
printf("after sorting array is:-----------------------------------------\n");
for(i = 0 ; i < n ; ++i)
{
printf("第%d个点：[%.0lf,%.0lf]\n",i,p[i].x,p[i].y);
}
double distance = closest(0 , n - 1);
if(distance != -1)
{
printf("分治法最近点对距离为：%.1lf\n",closest(0 , n - 1));
printf("最近点对为：[%.0lf,%.0lf],[%.0lf,%.0lf]\n",close[0].x,close[0].y,close[1].x,close[1].y);
}
else printf("只有少于两个点！\n");
qongju(p,n);
printf("请输入点的对数：");
return 0;
}