算法导论求n个点的最小距离

算法导论求n个点的最小距离
2010-01-20 17:23
在中文算法导论649页
算法：
0：把所有的点按照横坐标排序
1：用一条竖直的线L将所有的点分成两等份
2：递归算出左半部分的最近两点距离d1，右半部分的最近两点距离d2，取
d=min(d1,d2)
3：算出“一个在左半部分，另一个在右半部分”这样的点对的最短距离d3。

4：结果＝min(d1,d2,d3)
关键就是这第3步。

貌似这需要n^2的时间，把左边每个点和右边每个点都对比一下。

其实不然。

秘密就在这里。

首先，两边的点，与分割线L的距离超过d的，都可以扔掉了。

其次，即使两个点P1,P2（不妨令P1在左边，P2在右边）与分割线L的距离（水平距离）都小于d，如果它们的纵坐标之差大于d，也没戏。

就是这两点使得搜索范围大大减小：
对于左半部分的，与L的距离在d之内的，每个P1来说：右半部分内，符合以上两个条件的点P2最多只有6个！
原因就是：
d是两个半平面各自内，任意两点的最小距离，因此在同一个半平面内，任何两点距离都不可能超过d。

我们又要求P1和P2的水平距离不能超过d，垂直距离也不能超过d，在这个d*2d 的小方块内，最多只能放下6个距离不小于d的点。

因此，第3步总的比较距离的次数不超过n*6。

第3步的具体做法是：
3.1 删除所有到L的距离大于d的点。

O(n)
3.2 把右半平面的点按照纵坐标y排序。

O(nlogn)
3.3 对于左半平面内的每个点P1，找出右半平面内纵坐标与P1的纵坐标的差在d以内的点P2，计算距离取最小值，算出d3。

O(n*6) = O(n)
因为3.2的排序需要O(nlogn)，
所以整个算法的复杂度就是O(n((logn)^2))。

改进：
我们对3.2这个排序的O(nlogn)不太满意。

既然整个算法是递归的，我们可以利用第2步的子递归中已经排好序的序列，在第3.2部归并这两个子列，这样3.2的复杂度变成了O(n)。

这样，整个算法就是O(nlogn)的。

代码如下: VC6.0下编译通过
#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX 10000
typedefstruct point{
intx,y;
}POINT;
double delta = MAX ;
inttotalnum;
POINT mem_p1,mem_p2;
int cmx(const void *a,const void *b) //比较x坐标
{
return ((POINT *)a)->x-((POINT *)b)->x;
}
intcmy(const void *a,const void *b) //比较y坐标
{
return ((POINT *)a)->y-((POINT *)b)->y;
}
double min(double p1,double p2)
{
return p1>p2?p2:p1;
}
double dist(POINT s1,POINT s2) //求两点的距离
{
double dx = s1.x-s2.x;
double dy = s1.y-s2.y;
return sqrt(dx*dx + dy*dy);
}
void rec_cl_pair(POINT a[],inti,int j)
{
double temp,xx;
int k;
if(j-i<3)//小于或等于三个点的时候可以直接求得最小距离
{
qsort(a+i,j-i+1,sizeof(a[0]),cmy);//按Y坐标从小到大排列 xx = dist(a[i],a[i+1]);//两个点的距离
if(j-i==1)//只有两个点的时候直接返回
{
if(xx<delta)
{
mem_p1 = a[i];
mem_p2 = a[i+1];
delta = xx;
}
return;
}
double t = dist(a[i+1],a[i+2]); //有三个点的情况
if(t<delta)
{
mem_p1 = a[i+1];
mem_p2 = a[i+2];
delta = t;
}
t = dist(a[i],a[i+2]);
if(t<delta)
{
mem_p1 = a[i];
mem_p2 = a[i+2];
delta = t;
}
return;
}
k=(i+j)/2;
double middle=a[k].x;//中线点
rec_cl_pair(a,i,k);//求左边点的最小距离
rec_cl_pair(a,k+1,j);//求右边点的最小距离
int t=0;
POINT v[100];
for(k=i;k<=j;k++)
{
if(a[k].x>middle-delta && a[k].x < middle+ delta)
{
t++;
v[t] = a[k];//存入离中线距离小于delta的点
}
} //t-1为剩余点的个数
qsort(v+1,t,sizeof(v[1]),cmy);//以Y坐标的大小进行排序
for(k=1;k<t;k++)
{
for(int s=k+1;s<min(double(t),double(k+7));s++)
{
temp = dist(v[k],v[s]);
if(temp<delta)
{
delta = temp;
mem_p1=v[k];
mem_p2 = v[s];
}
}
}
}
void close_pair(POINT a[])
{
int n = totalnum;
qsort(a+1,n,sizeof(a[1]),cmx);//a接收的是a+1地址，按X坐标从小到大排列
rec_cl_pair(a,1,n);
}
void main(intargc, char *argv[])
{
POINT *a;
scanf("%d",&totalnum); //输入点的总数
a = (POINT *)malloc(sizeof(point)*(totalnum+1));//多申请一个内存空间a[0]不用
for(int i=1;i<=totalnum;i++)
scanf("%d%d",&a[i].x,&a[i].y);//输入点的X和Y坐标
close_pair(a);
printf("最近点对的距离为:%.3f\n",delta);//地址从a+1开始
printf("最短距离的两个点
为:\n(%d,%d)\n(%d,%d)\n",mem_p1.x,mem_p1.y,mem_p2.x,mem_p2.y);
}。