0007算法笔记——【分治法】最接近点对问题

转载链接:http://blog.csdn.net/liufeng_king/article/details/8484284   

问题场景:在应用中,常用诸如点、圆等简单的几何对象代表现实世界中的实体。在涉及这些几何对象的问题中,常需要了解其邻域中其他几何对象的信息。例如,在空中交通控制问题中,若将飞机作为空间中移动的一个点来看待,则具有最大碰撞危险的2架飞机,就是这个空间中最接近的一对点。这类问题是计算几何学中研究的基本问题之一。

      问题描述给定平面上n个点,找其中的一对点,使得在n个点的所有点对中,该点对的距离最小。严格地说,最接近点对可能多于1对。为了简单起见,这里只限于找其中的一对。

      1、一维最接近点对问题

       算法思路

       这个问题很容易理解,似乎也不难解决。我们只要将每一点与其他n-1个点的距离算出,找出达到最小距离的两个点即可。然而,这样做效率太低,需要O(n^2)的计算时间。在问题的计算复杂性中我们可以看到,该问题的计算时间下界为Ω(nlogn)。这个下界引导我们去找问题的一个θ(nlogn)算法。采用分治法思想,考虑将所给的n个点的集合S分成2个子集S1和S2,每个子集中约有n/2个点,然后在每个子集中递归地求其最接近的点对。在这里,一个关键的问题是如何实现分治法中的合并步骤,即由S1和S2的最接近点对,如何求得原集合S中的最接近点对,因为S1和S2的最接近点对未必就是S的最接近点对。如果组成S的最接近点对的2个点都在S1中或都在S2中,则问题很容易解决。但是,如果这2个点分别在S1和S2中,则对于S1中任一点p,S2中最多只有n/2个点与它构成最接近点对的候选者,仍需做n^2/4次计算和比较才能确定S的最接近点对。因此,依此思路,合并步骤耗时为O(n^2)。整个算法所需计算时间T(n)应满足: T(n)=2T(n/2)+O(n^2)。它的解为T(n)=O(n^2),即与合并步骤的耗时同阶,这不比用穷举的方法好。从解递归方程的套用公式法,我们看到问题出在合并步骤耗时太多。这启发我们把注意力放在合并步骤上。      

      设S中的n个点为x轴上的n个实数x1,x2,..,xn。最接近点对即为这n个实数中相差最小的2个实数。我们显然可以先将x1,x2,..,xn排好序,然后,用一次线性扫描就可以找出最接近点对。这种方法主要计算时间花在排序上,在排序算法已经证明,时间复杂度为O(nlogn)。然而这种方法无法直接推广到二维的情形。因此,对这种一维的简单情形,我们还是尝试用分治法来求解,并希望能推广到二维的情形。假设我们用x轴上某个点m将S划分为2个子集S1和S2,使得S1={x∈S|x≤m};S2={x∈S|x>m}。这样一来,对于所有p∈S1和q∈S2有p<q。递归地在S1和S2上找出其最接近点对{p1,p2}和{q1,q2},并设d=min{|p1-p2|,|q1-q2|},S中的最接近点对或者是{p1,p2},或者是{q1,q2},或者是某个{p3,q3},其中p3∈S1且q3∈S2。如图所示。

     如果S的最接近点对是{p3,q3},即|p3-q3|<d,则p3和q3两者与m的距离不超过d,即|p3-m|<d,|q3-m|<d,也就是说,p3∈(m-d,m],q3∈(m,m+d]。由于在S1中,每个长度为d的半闭区间至多包含一个点(否则必有两点距离小于d),并且m是S1和S2的分割点,因此(m-d,m]中至多包含S中的一个点。同理,(m,m+d]中也至多包含S中的一个点。由图可以看出,如果(m-d,m]中有S中的点,则此点就是S1中最大点。同理,如果(m,m+d]中有S中的点,则此点就是S2中最小点。因此,我们用线性时间就能找到区间(m-d,m]和(m,m+d]中所有点,即p3和q3。从而我们用线性时间就可以将S1的解和S2的解合并成为S的解。也就是说,按这种分治策略,合并步可在O(n)时间内完成。这样是否就可以得到一个有效的算法了呢?还有一个问题需要认真考虑,即分割点m的选取,及S1和S2的划分。选取分割点m的一个基本要求是由此导出集合S的一个线性分割,即S=S1∪S2 ,S1∩S2=Φ,且S1={x|x≤m};S2={x|x>m}。容易看出,如果选取m=[max(S)+min(S)]/2,可以满足线性分割的要求。选取分割点后,再用O(n)时间即可将S划分成S1={x∈S|x≤m}和S2={x∈S|x>m}。然而,这样选取分割点m,有可能造成划分出的子集S1和S2的不平衡。例如在最坏情况下,|S1|=1,|S2|=n-1,由此产生的分治法在最坏情况下所需的计算时间T(n)应满足递归方程:
   T(n)=T(n-1)+O(n)
    它的解是T(n)=O(n^2)。这种效率降低的现象可以通过分治法中“平衡子问题”的方法加以解决。即通过适当选择分割点m,使S1和S2中有大致相等个数的点。自然地,我们会想到用S的n个点的坐标的中位数来作分割点。在选择算法中介绍的选取中位数的线性时间算法使我们可以在O(n)时间内确定一个平衡的分割点m。

       本程序确定平衡点采用m=[max(S)+min(S)]/2方法。如果需要利用中位数作分割点,看结合笔者博文《0005算法笔记——线性时间选择》改写。

     一维最接近临近点对问题程序清单如下:

[cpp]  view plain  copy
  1. //2d10-1 一维最邻近点对问题  
  2. #include "stdafx.h"  
  3. #include <ctime>  
  4. #include <iostream>   
  5. using namespace std;   
  6.   
  7. const int L=100;  
  8. //点对结构体  
  9. struct Pair  
  10. {  
  11.     float d;//点对距离  
  12.     float d1,d2;//点对坐标  
  13. };  
  14. float Random();  
  15. int input(float s[]);//构造S  
  16. float Max(float s[],int p,int q);  
  17. float Min(float s[],int p,int q);  
  18. template <class Type>  
  19. void Swap(Type &x,Type &y);  
  20. template <class Type>  
  21. int Partition(Type s[],Type x,int l,int r);  
  22. Pair Cpair(float s[],int l,int r);  
  23.   
  24. int main()  
  25. {  
  26.     srand((unsigned)time(NULL));  
  27.     int m;  
  28.     float s[L];  
  29.     Pair d;  
  30.     m=input(s);  
  31.     d=Cpair(s,0,m-1);  
  32.     cout<<endl<<"最近点对坐标为: (d1:"<<d.d1<<",d2:"<<d.d2<<")";  
  33.     cout<<endl<<"这两点距离为: "<<d.d<<endl;  
  34.     return 0;  
  35. }  
  36.   
  37.   
  38. float Random()  
  39. {  
  40.     float result=rand()%10000;  
  41.      return result*0.01;  
  42. }  
  43.   
  44. int input(float s[])  
  45. {  
  46.     int length;  
  47.     cout<<"输入点的数目: ";  
  48.     cin>>length;  
  49.     cout<<"点集在X轴上坐标为:";  
  50.     for(int i=0;i<length;i++)  
  51.     {  
  52.         s[i]=Random();  
  53.         cout<<s[i]<<" ";  
  54.     }  
  55.       
  56.     return length;  
  57. }  
  58.   
  59.   
  60. float Max(float s[],int l,int r)//返回s[]中的最大值  
  61. {  
  62.     float s_max=s[l];  
  63.     for(int i=l+1;i<=r;i++)  
  64.         if(s_max<s[i])  
  65.             s_max=s[i];  
  66.     return s_max;  
  67. }  
  68.   
  69. float Min(float s[],int l,int r)//返回s[]中的最小值  
  70. {  
  71.     float s_min=s[l];  
  72.     for(int i=l+1;i<=r;i++)   
  73.         if(s_min>s[i])  
  74.             s_min=s[i];  
  75.     return s_min;  
  76. }  
  77.   
  78. template <class Type>  
  79. void Swap(Type &x,Type &y)  
  80. {  
  81.     Type temp = x;  
  82.     x = y;  
  83.     y = temp;  
  84. }  
  85.   
  86. template <class Type>  
  87. int Partition(Type s[],Type x,int l,int r)  
  88. {  
  89.     int i = l - 1,j = r + 1;  
  90.   
  91.     while(true)  
  92.     {  
  93.         while(s[++i]<x && i<r);  
  94.         while(s[--j]>x);  
  95.         if(i>=j)  
  96.         {  
  97.             break;  
  98.         }  
  99.         Swap(s[i],s[j]);  
  100.     }  
  101.     return j;  
  102. }  
  103.   
  104. //返回s[]中的具有最近距离的点对及其距离  
  105. Pair Cpair(float s[],int l,int r)  
  106. {  
  107.     Pair min_d={99999,0,0};//最短距离  
  108.   
  109.     if(r-l<1) return min_d;  
  110.     float m1=Max(s,l,r),m2=Min(s,l,r);  
  111.   
  112.     float m=(m1+m2)/2;//找出点集中的中位数  
  113.   
  114.     //将点集中的各元素按与m的大小关系分组  
  115.     int j = Partition(s,m,l,r);  
  116.   
  117.     Pair d1=Cpair(s,l,j),d2=Cpair(s,j+1,r);//递归  
  118.     float p=Max(s,l,j),q=Min(s,j+1,r);  
  119.   
  120.     //返回s[]中的具有最近距离的点对及其距离  
  121.     if(d1.d<d2.d)  
  122.     {  
  123.         if((q-p)<d1.d)  
  124.         {  
  125.             min_d.d=(q-p);  
  126.             min_d.d1=q;  
  127.             min_d.d2=p;  
  128.             return min_d;  
  129.         }  
  130.         else return d1;  
  131.     }  
  132.     else  
  133.     {  
  134.         if((q-p)<d2.d)  
  135.         {  
  136.             min_d.d=(q-p);  
  137.             min_d.d1=q;  
  138.             min_d.d2=p;  
  139.             return min_d;  
  140.         }  
  141.         else return d2;  
  142.     }  
  143. }  
程序运行结果如下:


      该算法的分割步骤和合并步骤总共耗时O(n)。因此,算法耗费的计算时间T(n)满足递归方程:

      

       解此递归方程可得T(n)=O(nlogn)。

         2、二维最接近点对问题

     将以上过程推广到二维最接近点对问题,设S中的点为平面上的点,它们都有2个坐标值x和y。为了将平面上点集S线性分割为大小大致相等的2个子集S1和S2,我们选取一垂直线l:x=m来作为分割直线。其中m为S中各点x坐标的中位数。由此将S分割为S1={p∈S|px≤m}和S2={p∈S|px>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)<d则p和q必分属于S1和S2。不妨设p∈S1,q∈S2。那么p和q距直线l的距离均小于d。因此,我们若用P1和P2分别表示直线l的左边和右边的宽为d的2个垂直长条,则p∈S1,q∈S2,如图所示:


     距直线l的距离小于d的所有点

       在一维的情形,距分割点距离为d的2个区间(m-d,m](m,m+d]中最多各有S中一个点。因而这2点成为唯一的末检查过的最接近点对候选者。二维的情形则要复杂些,此时,P1中所有点与P2中所有点构成的点对均为最接近点对的候选者。在最坏情况下有n2/4对这样的候选者。但是P1和P2中的点具有以下的稀疏性质,它使我们不必检查所有这n^2/4对候选者。考虑P1中任意一点p,它若与P2中的点q构成最接近点对的候选者,则必有d(p,q)<d。满足这个条件的P2中的点有多少个呢?容易看出这样的点一定落在一个d×2d的矩形R中,如下图所示:


包含点q的dX2d矩形R

     由d的意义可知P2中任何2个S中的点的距离都不小于d。由此可以推出矩形R中最多只有6个S中的点。事实上,我们可以将矩形R的长为2d的边3等分,将它的长为d的边2等分,由此导出6个(d/2)×(2d/3)的矩形。如左图所示:

矩阵R中点的稀疏性

     若矩形R中有多于6个S中的点,则由鸽舍原理易知至少有一个δ×2δ的小矩形中有2个以上S中的点。设u,v是这样2个点,它们位于同一小矩形中,则:

因此d(u,v)≤5d/6<d 。这与d的意义相矛盾。也就是说矩形R中最多只有6个S中的点。图4(b)是矩形R中含有S中的6个点的极端情形。由于这种稀疏性质,对于P1中任一点p,P2中最多只有6个点与它构成最接近点对的候选者。因此,在分治法的合并步骤中,我们最多只需要检查6×n/2=3n对候选者,而不是n^2/4对候选者。这是否就意味着我们可以在O(n)时间内完成分治法的合并步骤呢?现在还不能作出这个结论,因为我们只知道对于P1中每个S1中的点p最多只需要检查P2中的6个点,但是我们并不确切地知道要检查哪6个点。为了解决这个问题,我们可以将p和P2中所有S2的点投影到垂直线l上。由于能与p点一起构成最接近点对候选者的S2中点一定在矩形R中,所以它们在直线l上的投影点距p在l上投影点的距离小于d。由上面的分析可知,这种投影点最多只有6个。因此,若将P1和P2中所有S的点按其y坐标排好序,则对P1中所有点p,对排好序的点列作一次扫描,就可以找出所有最接近点对的候选者,对P1中每一点最多只要检查P2中排好序的相继6个点。

程序清单如下:

[cpp]  view plain  copy
  1. //2d10-2 二维最邻近点对问题  
  2. #include "stdafx.h"  
  3. #include<time.h>  
  4. #include<iostream>   
  5. #include<cmath>  
  6.   
  7. using namespace std;  
  8. const int M=50;  
  9.   
  10. //用类PointX和PointY表示依x坐标和y坐标排好序的点  
  11. class PointX {  
  12.     public:   
  13.         int operator<=(PointX a)const  
  14.         { return (x<=a.x); }  
  15.         int ID; //点编号  
  16.         float x,y; //点坐标   
  17. };  
  18.   
  19. class PointY {   
  20.     public:   
  21.         int operator<=(PointY a)const  
  22.         { return(y<=a.y); }  
  23.         int p; //同一点在数组x中的坐标   
  24.         float x,y; //点坐标  
  25. };  
  26.   
  27. float Random();  
  28. template <class Type>  
  29. float dis(const Type&u,const Type&v);   
  30.   
  31. bool Cpair2(PointX X[], int n,PointX& a,PointX& b, float& d);  
  32. void closest(PointX X[],PointY Y[],PointY Z[], int l, int r,PointX& a,PointX& b,float& d);  
  33.   
  34. template <typename Type>   
  35. void Copy(Type a[],Type b[], int left,int right);  
  36.   
  37. template <class Type>  
  38. void Merge(Type c[],Type d[],int l,int m,int r);  
  39.   
  40. template <class Type>  
  41. void MergeSort(Type a[],Type b[],int left,int right);  
  42.   
  43. int main()  
  44. {   
  45.     srand((unsigned)time(NULL));  
  46.     int length;  
  47.   
  48.     cout<<"请输入点对数:";  
  49.     cin>>length;  
  50.   
  51.     PointX X[M];  
  52.     cout<<"随机生成的二维点对为:"<<endl;  
  53.   
  54.     for(int i=0;i<length;i++)  
  55.     {  
  56.         X[i].ID=i;  
  57.         X[i].x=Random();  
  58.         X[i].y=Random();  
  59.         cout<<"("<<X[i].x<<","<<X[i].y<<") ";  
  60.     }  
  61.   
  62.     PointX a;   
  63.     PointX b;   
  64.     float d;  
  65.   
  66.     Cpair2(X,length,a,b,d);   
  67.   
  68.     cout<<endl;  
  69.     cout<<"最邻近点对为:("<<a.x<<","<<a.y<<")和("<<b.x<<","<<b.y<<") "<<endl;  
  70.     cout<<"最邻近距离为: "<<d<<endl;  
  71.   
  72.     return 0;  
  73. }  
  74.   
  75. float Random()  
  76. {  
  77.     float result=rand()%10000;  
  78.     return result*0.01;  
  79. }  
  80.   
  81. //平面上任意两点u和v之间的距离可计算如下  
  82. template <class Type>  
  83. inline float dis(const Type& u,const Type& v)  
  84. {  
  85.     float dx=u.x-v.x;  
  86.     float dy=u.y-v.y;   
  87.     return sqrt(dx*dx+dy*dy);   
  88. }  
  89.   
  90. bool Cpair2(PointX X[], int n,PointX& a,PointX& b,float& d)  
  91. {  
  92.     if(n<2) return false;  
  93.   
  94.     PointX* tmpX = new PointX[n];  
  95.     MergeSort(X,tmpX,0,n-1);  
  96.   
  97.     PointY* Y=new PointY[n];   
  98.     for(int i=0;i<n;i++) //将数组X中的点复制到数组Y中  
  99.     {   
  100.         Y[i].p=i;  
  101.         Y[i].x=X[i].x;  
  102.         Y[i].y=X[i].y;  
  103.     }   
  104.   
  105.     PointY* tmpY = new PointY[n];  
  106.     MergeSort(Y,tmpY,0,n-1);  
  107.   
  108.     PointY* Z=new PointY[n];  
  109.     closest(X,Y,Z,0,n-1,a,b,d);   
  110.   
  111.     delete []Y;   
  112.     delete []Z;  
  113.     delete []tmpX;  
  114.     delete []tmpY;  
  115.     return true;   
  116. }  
  117. void closest(PointX X[],PointY Y[],PointY Z[], int l, int r,PointX& a,PointX& b,float& d)   
  118. {   
  119.     if(r-l==1) //两点的情形   
  120.     {  
  121.         a=X[l];  
  122.         b=X[r];  
  123.         d=dis(X[l],X[r]);  
  124.         return;   
  125.     }   
  126.   
  127.     if(r-l==2) //3点的情形   
  128.     {  
  129.         float d1=dis(X[l],X[l+1]);  
  130.         float d2=dis(X[l+1],X[r]);   
  131.         float d3=dis(X[l],X[r]);   
  132.   
  133.         if(d1<=d2 && d1<=d3)   
  134.         {   
  135.             a=X[l];  
  136.             b=X[l+1];  
  137.             d=d1;  
  138.             return;  
  139.         }   
  140.   
  141.         if(d2<=d3)  
  142.         {   
  143.             a=X[l+1];  
  144.             b=X[r];  
  145.             d=d2;  
  146.         }   
  147.         else {   
  148.             a=X[l];   
  149.             b=X[r];   
  150.             d=d3;   
  151.         }  
  152.         return;   
  153.     }   
  154.   
  155.     //多于3点的情形,用分治法   
  156.     int m=(l+r)/2;   
  157.     int f=l,g=m+1;   
  158.   
  159.     //在算法预处理阶段,将数组X中的点依x坐标排序,将数组Y中的点依y坐标排序  
  160.     //算法分割阶段,将子数组X[l:r]均匀划分成两个不想交的子集,取m=(l+r)/2  
  161.     //X[l:m]和X[m+1:r]就是满足要求的分割。  
  162.     for(int i=l;i<=r;i++)  
  163.     {  
  164.         if(Y[i].p>m) Z[g++]=Y[i];   
  165.         else Z[f++]=Y[i];  
  166.     }  
  167.   
  168.     closest(X,Z,Y,l,m,a,b,d);  
  169.     float dr;  
  170.   
  171.     PointX ar,br;  
  172.     closest(X,Z,Y,m+1,r,ar,br,dr);   
  173.   
  174.     if(dr<d)  
  175.     {  
  176.         a=ar;   
  177.         b=br;   
  178.         d=dr;   
  179.     }   
  180.   
  181.     Merge(Z,Y,l,m,r);//重构数组Y  
  182.   
  183.     //d矩形条内的点置于Z中  
  184.     int k=l;   
  185.     for(int i=l;i<=r;i++)  
  186.     {  
  187.         if(fabs(X[m].x-Y[i].x)<d)  
  188.         {   
  189.             Z[k++]=Y[i];   
  190.         }  
  191.     }  
  192.   
  193.     //搜索Z[l:k-1]  
  194.     for(int i=l;i<k;i++)   
  195.     {   
  196.         for(int j=i+1;j<k && Z[j].y-Z[i].y<d;j++)   
  197.         {   
  198.             float dp=dis(Z[i],Z[j]);  
  199.             if(dp<d)   
  200.             {   
  201.                 d=dp;   
  202.                 a=X[Z[i].p];  
  203.                 b=X[Z[j].p];   
  204.             }  
  205.         }  
  206.     }   
  207. }  
  208.   
  209. template <class Type>  
  210. void Merge(Type c[],Type d[],int l,int m,int r)  
  211. {  
  212.     int i = l,j = m + 1,k = l;  
  213.     while((i<=m)&&(j<=r))  
  214.     {  
  215.         if(c[i]<=c[j])  
  216.         {  
  217.             d[k++] = c[i++];  
  218.         }  
  219.         else  
  220.         {  
  221.             d[k++] = c[j++];  
  222.         }  
  223.     }  
  224.   
  225.     if(i>m)  
  226.     {  
  227.         for(int q=j; q<=r; q++)  
  228.         {  
  229.             d[k++] = c[q];  
  230.         }     
  231.     }  
  232.     else  
  233.     {  
  234.         for(int q=i; q<=m; q++)  
  235.         {  
  236.             d[k++] = c[q];  
  237.         }  
  238.     }  
  239. }  
  240.   
  241. template <class Type>  
  242. void MergeSort(Type a[],Type b[],int left,int right)  
  243. {  
  244.     if(left<right)  
  245.     {  
  246.         int i = (left + right)/2;  
  247.         MergeSort(a,b,left,i);  
  248.         MergeSort(a,b,i+1,right);  
  249.         Merge(a,b,left,i,right);//合并到数组b  
  250.         Copy(a,b,left,right);//复制回数组a         
  251.     }  
  252. }  
  253.   
  254. template <typename Type>   
  255. void Copy(Type a[],Type b[], int left,int right)  
  256. {  
  257.     for(int i=left;i<=right;i++)   
  258.         a[i]=b[i];   
  259. }  

程序结果:

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值