寻找最近点对问题(FCPP)
问题描述:
给定平面上n个点,找出距离最近的两个点。
思考过程:
1)对于这种问题,我们首先想到的求解方法就是求出所有点对的距离,并找出最近的那个,当然这个是个显而
易见的方法,具体过程大体可以警醒如下描述。
① 定义变量
N (点的数目),X[N] (点的x坐标值),Y[N](点的y坐标值),closest(最近的值),P(最近点对的一个点),Q(最近点对的另一个点)。
② 依次扫描并比较出最近点对
初始化closest=get_len(X[1],Y[1],X[2],Y[2]);
P=(X[1],Y[1]);Q=(X[2],Y[2]);
For i=1:N
For j=i+1:N
Tmp=get_len(X[i],Y[i],X[j],Y[j]);
If Tmp < closest
Closest=Tmp;
P=(X[i],Y[i]);
Q= (X[j],Y[j]);
End If
End For
End For
这样运行后就可以求出所有的点中的最近的点对值,然而这是一种Brute-force(蛮力法),当然时间复杂度也是非常不能令人满意的Q(n2)。
2)对于1)中的算法基本上很简单而且很容易想到,的确,对于一些问题,人们往往想的简单做的难,想的难就做的简单啦(呵呵)。在此,可以利用分治法来对问题求解。具体算法步骤如下:
① 将所有的点按X坐标值排序,然后平均等分成两个部分PL与PR。在左右分开处画个分割线part_line。(如下图)
② 找出PL中的最近值与PR中的最近值分别为dL,dR。
③ 令d=min(dL,dR);现在,问题来了,d一定最小么?
回答是相当确定的,就是肯定不是最小的。原因也显而易见,就是最近点对可能出现在PL与PR的分割处。
④ 根据③中的问题,可以以part_line为中心左右|d|的距离里的点考虑(为什么?超过这个距离显然不需要 考虑啦,因为距离一定会大于d啦~~~)。现在,问题又来了,如何取得这个范围内的最近点呢?蛮力?全搜索计算?对!不行的这是,可能所有的点都落在这个范围里,就此一项我的时间复杂度就又会回到Q(n2)。怎么办?
⑤ 对于④中的问题,解决方法是这样的,将左右|d|范围里的点按Y值排序,然后依次从每个点出发做水平线,并在其之上d的距离做水平线,这样得到了一个2d*d的矩形,显然,在此矩形之外的点无需开率啦,可以证明,矩形内最多容下8个点(包括自己本身)。所以意味着在排好序的Y坐标点中,只需要考虑其后的7个点就行了,此外已有人证明了事实上只需要考虑其后的4个点就行啦。可见,这些的时间复杂度最多Q(N)。
⑥ 纵观整个过程,首先要对X坐标排序时间为O(NlogN),这是分治之前的预处理。然后问题分成了N/2规模,然后用Q(N)的复杂度得到中间的最近点对,然后可以在常数时间内得到最终的点对与最近值。So。。。T(N)=2T(N/2)+O(N),可以得到该算法的时间为O(NlogN)(即使加上一开始的排序预处理的时间复杂度)。
算法具体实现过程:
①定义函数 FCPP(st,ed,X,Y,IND,d,p,q);
其中:
st,ed分别是X数组的最小和最大下标
X[]和Y[] 这个不用说,就是存放点的X和Y坐标的数组
IND[]是存放Y坐标对应X坐标的对应关系,比如(X[IND[i]],Y[i])表示一个点。(以上是输入的变量)
δ表示最近的值,p,q分别是最近点对的点 (以上是输出变量)
②具体运行过程为:
FCPP(j,k,X,Ypres,INDpres,d,p,q)
{
令m= k-j+1;//即m是X[j]到X[k]这一段中的数据元素的个数。
if(m≤3){
根据Ypres,INDpres和X数组找出全部3个(或2个)点,
求出其中的最近点对(p,q)及最小距离d,然后返回d,p,q。
}
if(m≥4){
进行分治前的准备工作:
引入4个数组YL,INDL(长度均为ém/2ù),
YR,INDR(长度均为ëm/2û),
YL存放左边点集中点的y坐标(按从小到大的次序),
INDL存放与之对应的x坐标在X数组中的下标,
YR和INDR的作用与YL和INDL类似。
LengthL←0;LengthR←0;/* LengthL记录放入左边数组YL中的点数*/
For i=1 to m do /*对m个点逐一进行检查,Q(m)时间*/
{
if INDpres[i]≤j+m/2-1 /*当前检查的点应分在左边的点集中*/
{ LengthL增1; /*左边的点数增1*/
YL[LengthL]←Ypres [i];/*当前检查点的y坐标放在左边Y数组中
INDL[LengthL]←INDpres [i];/*在左边索引数组中保持y坐标与x坐标的对应关系*/
} else { /*当前被检查的点应分在右边的点集中*/
LengthR增1; /*右边的点数增1*/
YR[LengthR]←Ypres[i];/*当前检查点的y坐标放在右边Y数组中*/
INDR[LengthR]←INDpres [i];/*在右边索引数组中保持y坐标与x坐标的对应关系*/
} /*当前检查点已放在左边或右边的Y数组中*/
} /*m个点检查完毕时,*/
/*YL和YR仍然保持了从小到大的次序*/
/*然后分治:*/
FCPP(j, j+m/2-1, X, YL, INDL,dL,pL,qL);
FCPP(j+m/2, k, X, YR, INDR,dR,pR,qR);
/*以上两个调用执行时间为2T(m/2)*/
令d=min{dL,dR }并记录下相应的(p,q); /*O(1)时间*/
/*求分割线的x坐标值:*/
xdiv=(X[j+m/2-1]+ X[j+m/2])/2; /*O(1)时间*/
/*引入数组Yd和INDd,长度最大为m,*/
/*Yd放带状域内点的y坐标,INDd放对应的x的下标*/
LengthYd←0;
For i=1 to m do /*对m个点逐一进行检查,Q(m)时间*/
{ if (xdiv-d < X[ INDpres[i] ] < xdiv+d)
/*当前点落到2d带状区域内*/
{ LengthYd增1;Yd [LengthYd]←Ypres [i];
/*记录下该点的y坐标*/
INDd [LengthYd]←INDpres [i];}
/*保持y坐标与x坐标的对应关系*/
} /*m个点检查完毕时,*/
/*Yd中的y坐标仍然保持了从小到大的次序*/
if LengthYd≤1 then return;/*只有一个或没有点落到带状区域内,无须算*/
For i=1 to LengthYd-1 do /* 至少有两个点落到*/
/*带状区域内时对每个点要逐一检查,*/
/*最后一个点因其后无点,故该点无需检查*/
{ ypres←Yd [i];xpres←X[INDd[i]];
u←min{i+4, LengthYd}; /*u是本次检查的最大下标*/
v←i+1;/* v是本次检查的游动变量,从i+1到u为止*/
while (v≤u & (Yd[v] - ypres) < d) /*y坐标值之差大于等于d就停止检查,因其后不可能有*/
{
计算点(xpres,ypres)与点(X[INDd[v]], Yd[v])之间的距离d’;if d’< d then {d←d’};
点p,q更新为(xpres,ypres)和(X[INDd[v]], Yd[v]);}
v增1; /*v所指的点检查完毕*/
}/*该while循环最多执行4次,故为O(1)时间*/
} /* LengthYd最大为m,故For循环执行时间最多为O(m)*/
}//end for
}//end if m>=4
}//end FCPP
实现的代码为(C语言):
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAXLEN 200
int init(float *x,float *y,int *ind)
{
int n,i;
while(1==scanf("%d",&n))
{
if(n > 1) break;
printf("n invalid!\n");
}
for(i=0;i<n;i++)
{
scanf("%f%f",x+i,y+i);
ind[i]=i;
}
return n;
}
void pre_sort(float *x,float *y,int *ind,int n)
{//must finish in O(nlogn) but here using normorl sort method
//bubble sort
int i,j;
int tmp;
float tmp_y,tmp_x;
for(i=0;i<n;i++)
for(j=n-1;j > i;j--)
if(x[j-1] > x[j])
{
//swap x
tmp_x=x[j-1];
x[j-1]=x[j];
x[j]=tmp_x;
//swap y
tmp_y=y[j-1];
y[j-1]=y[j];
y[j]=tmp_y;
//swap ind no need now!
/*tmp=ind[j-1];
ind[j-1]=ind[j];
ind[j]=tmp;*/
}
for(i=0;i<n;i++)
for(j=n-1;j > i;j--)
if(y[j-1] > y[j])
{
//swap y
tmp_y=y[j-1];
y[j-1]=y[j];
y[j]=tmp_y;
//swap ind
tmp=ind[j-1];
ind[j-1]=ind[j];
ind[j]=tmp;
}
}
float get_pair_len(float x1,float y1,float x2,float y2)
{
return sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
}
void find_closest_pair(int st,int ed,float *x,float *y,int *ind,float *closest,float *x1,float *y1,float *x2,float *y2)
{//key code
if(ed-st == 1){//only 2 point
*closest=get_pair_len(x[ind[0]],y[0],x[ind[1]],y[1]);
*x1=x[ind[0]];
*y1=y[0];
*x2=x[ind[1]];
*y2=y[1];
}
else if(ed-st ==2)
{//only 3 point
float len;
*closest=get_pair_len(x[ind[0]],y[0],x[ind[1]],y[1]);
*x1=x[ind[0]];
*y1=y[0];
*x2=x[ind[1]];
*y2=y[1];
len=get_pair_len(x[ind[0]],y[0],x[ind[2]],y[2]);
if(len < *closest){
*closest=len;
*x1=x[ind[0]];
*y1=y[0];
*x2=x[ind[2]];
*y2=y[2];
}
len=get_pair_len(x[ind[1]],y[1],x[ind[2]],y[2]);
if(len < *closest){
*closest=len;
*x1=x[ind[1]];
*y1=y[1];
*x2=x[ind[2]];
*y2=y[2];
}
}else{//at least 4 points
int i,cl,cr,ct;
int mid;
float t_c1,t_c2;
float t_c1_x1,t_c1_y1,t_c1_x2,t_c1_y2;
float t_c2_x1,t_c2_y1,t_c2_x2,t_c2_y2;
float *yl,*yr,*yct;
int *indl,*indr,*indct;
mid=(st+ed)/2;
yl=(float*)malloc((mid-st+1)*sizeof(float));
indl=(int*)malloc((mid-st+1)*sizeof(int));
yr=(float*)malloc((ed-mid)*sizeof(float));
indr=(int*)malloc((ed-mid)*sizeof(int));
if(!yl || !yr || !indl || !indr) exit(-2);//non enough mm
for(i=0,cl=0,cr=0;i<=ed-st;i++)
{//store new order y&ind
if(ind[i] <= mid){
yl[cl]=y[i];
indl[cl]=ind[i];
cl++;
}
else{
yr[cr]=y[i];
indr[cr]=ind[i];
cr++;
}
}
//divided
find_closest_pair(st,mid,x,yl,indl,&t_c1,&t_c1_x1,&t_c1_y1,&t_c1_x2,&t_c1_y2);
find_closest_pair(mid+1,ed,x,yr,indr,&t_c2,&t_c2_x1,&t_c2_y1,&t_c2_x2,&t_c2_y2);
if(t_c1 < t_c2){
*closest=t_c1;
*x1=t_c1_x1;
*y1=t_c1_y1;
*x2=t_c1_x2;
*y2=t_c1_y2;
}
else{
*closest=t_c2;
*x1=t_c2_x1;
*y1=t_c2_y1;
*x2=t_c2_x2;
*y2=t_c2_y2;
}//get the closest pair between the two part(L&R)
int v,ve;
float closest_tmp;
float part_line=(x[mid]+x[mid+1])/2.0;
yct=(float*)malloc((ed-st+1)*sizeof(float));
indct=(int*)malloc((ed-st+1)*sizeof(int));
//get all the y in the rect
for(i=0,ct=0;i<=ed-st;i++)
if(x[ind[i]]>part_line-*closest && x[ind[i]]<part_line+*closest ){
yct[ct]=y[i];
indct[ct]=ind[i];
ct++;
}
//check followed the 4 points
if(ct > 1){
for(i=0;i< ct-1;i++){
ve=(i+4 < ct-1 ? i+4:ct-1);
v=i+1;
while(v<=ve && yct[v]-yct[i] <=*closest){
closest_tmp=get_pair_len(x[indct[i]],yct[i],x[indct[v]],yct[v]);
if(closest_tmp < *closest){
*closest=closest_tmp;
*x1=x[indct[i]];
*y1=yct[i];
*x2=x[indct[v]];
*y2=yct[v];
}
v++;
}
}
}
free(yl);
free(yr);
free(indl);
free(indr);
free(yct);
free(indct);
}
}
int main()
{
float X[MAXLEN],Y[MAXLEN];
int IND[MAXLEN],n;
n=init(X,Y,IND);
float closest,x1,y1,x2,y2;
pre_sort(X,Y,IND,n);
find_closest_pair(0,n-1,X,Y,IND,&closest,&x1,&y1,&x2,&y2);
printf("the closest value is:%-4.2f (%-4.2f,%-4.2f) (%-4.2f,%-4.2f)\n",closest,x1,y1,x2,y2);
return 0;
}
对与这个问题还有个更好的解法,可以在O(N)的时间内搞定,用的随机算法,具体的在以后贴出来的,that's all~~~