2007-06-21

二维点集的凸包及其直径(1)

关键字: 凸包 Graham 水平排序

前言:因为前几天做了一个有关凸包的题,并答应crackerwang写个blog解释一下我的算法.因为我比较懒的原因,一直拖到现在才写.预计一共有两篇,第一篇介绍求二维点集凸包的O(N*logN)时间复杂度的算法.第二篇介绍求凸包直径的O(N)时间复杂度的算法.
下面首先给出http://acm.tju.edu.cn/toj/showp2847.html该题的C++代码,本文将使用Java代码来描述. 

cpp 代码
  1. /**  
  2. * TJU 2847 的C++解法  
  3. * Written By Eastsun   
  4. */  
  5. #include＜stdio.h＞   
  6. #include＜math.h＞   
  7. #include＜algorithm＞   
  8. #include＜functional＞   
  9. #define  S(arr,a,b,c) ((arr[b].x-arr[a].x)*(arr[c].y-arr[a].y)-(arr[c].x-arr[a].x)*(arr[b].y-arr[a].y))    
  10. #define  J(arr,a,b,c,d)  ((arr[b].x-arr[a].x)*(arr[d].y-arr[c].y)-(arr[d].x-arr[c].x)*(arr[b].y-arr[a].y))   
  11. #define  Q(x) ((x)*(x))   
  12. #define  D(arr,a,b) (Q(arr[a].x-arr[b].x)+Q(arr[a].y-arr[b].y))   
  13. using namespace std;   
  14.   
  15. struct Point{   
  16.     double x,y;   
  17. };   
  18.   
  19. struct myCmp:public binary_function＜Point,Point,bool＞{   
  20.     bool operator()(const Point& p1,const Point& p2)const{   
  21.         if(p1.x==p2.x) return p1.y＜=p2.y;   
  22.         else           return p1.x＜p2.x;   
  23.     }   
  24. };   
  25.   
  26. int main(){   
  27.     int n;   
  28.     while(scanf("%d",&n),n){   
  29.         Point *ps =new Point[n];   
  30.         for(int i=0;i＜n;i++) scanf("%lf%lf",&ps[i].x,&ps[i].y);   
  31.         sort(ps,ps+n,myCmp());   
  32.         int *v =new int[2*n];   
  33.         int p,q;   
  34.         p =q =n;   
  35.         for(int i=0;i＜n;i++){   
  36.             v[p] =v[q] =i;   
  37.             while(p-n＞=2&&S(ps,v[p],v[p-1],v[p-2])＞=0){v[p-1] =v[p]; p--;}   
  38.             while(n-q＞=2&&S(ps,v[q+2],v[q+1],v[q])＞=0){v[q+1] =v[q]; q++;}   
  39.             p++;   
  40.             q--;   
  41.         }   
  42.         int len =p -q -2;   
  43.         Point *pps =new Point[len+2];   
  44.         for(int i=q+1;i＜p;i++) pps[i-q] =ps[v[i]];    
  45.         pps[0] =pps[len];   
  46.         int i=0,j=1;   
  47.         while(J(pps,i,i+1,j,j+1)＞0) j++;   
  48.         double max =0;   
  49.         while(i＜=len){   
  50.             double det =J(pps,i,i+1,j,j+1);   
  51.             if(det＜0) i++;   
  52.             else if(det＞0) j =(j+1＞len?j+1-len:j+1);   
  53.             else{   
  54.                 i++;   
  55.                 continue;   
  56.             }   
  57.             double tmp =D(pps,i,j);   
  58.             if(tmp＞max) max =tmp;   
  59.         }   
  60.         delete[] ps;   
  61.         delete[] pps;   
  62.         delete[] v;   
  63.         printf("%.2lf\n",sqrt(max));   
  64.     }   
  65.     return 0;   
  66. }  

一:凸包的定义及其应用
  对于二维点集,其凸包(convex hull)是指包含所有点的最小的凸多边形.直观上看,如果用一条橡皮筋将所有的点圈住,当橡皮筋拉紧后的形状就是这些点的凸包.下面就是凸包的示意图:

convex hull 

       那么,我们为什么需要凸包这个概念呢?它又能解决什么问题?
       首先,凸包上的点相对原有的点集,我们可以想象,其数量将大大减少.研究表明,对于二维情形,凸包顶点数m(k) =O(n^1/3).更一般的,对于k维球体中均匀分布的n个点,其凸包顶点数m(k) =O(n^(k-1)/(k+1)),可见凸包可大大降低平均意义下的时空复杂度.
      另一方面,凸包相对原有点集增加了一个"序",原来是一个杂乱无章的点集,而现在是一个性质优美的凸多边形,研究起来方便很多.
 
      关于其应用,这里我们只针对TJU2847来说,原题是求一个点集中任意两点的距离的最大值.显然,如果我们直接考虑这个点集中任意两点的距离,时间复杂度是O(N^2),下面我们可以看到,当求出其凸包后,我们可以在O(k)时间内求出这个值(这里的k是指凸包顶点个数k<=N).再加上构造凸包的时间复杂度O(N*logN),我们在O(N*logN)时间复杂度内解决了这个问题.
 
二:构造凸包的算法
      二维点集构造凸包的算法有:卷包裹法(Gift-Wrapping),Graham-Scan扫描法以及分治算法,增量算法等等.这里我采用的是Graham-Scan算法的一种变形--x-y排序.(至于原本的极角排序的Graham-Scan算法,有兴趣的可以参阅《Introduction to Algorithms》第VII章的Computational Geometry一节,里面有详细的说明及正确性证明):
      首先将给定点集对x的值进行排序,x值相同的按y的值排序.简单来说就是从左到右,从下到上排序.
      排序后,记这些点为ps[0,1,..,k],显然点列的第一个点与最后一个点都在凸包上(想一想那个橡皮筋),然后我们从第一点开始到最后一个点进行如下处理:
      构造一个堆栈(数组)stack[],stack[0] =ps[0],然后维持栈中点都是按右手系旋转(或者说从stack[0]到stack[p],所有的点都是按逆时针排列),对于每一个新增加的点,都首先检查栈顶的两个点与其是不是保持右手系,如果是,把这个点加入栈;否则,将栈顶点去掉,继续检查...一直到符合要求为止.

      这样处理后的结果是:stack[]中得到一条连接ps[0]与ps[k]的,并且维持逆时针顺转的链.这个链就是我们要求的凸包的下半部分.用伪代码来描叙就是:

java 代码
  1. //ps中保存了已经排好序的点   
  2. Point[] stack =new Point[ps.length];   
  3. int index =0, p =0;   
  4. while(index
  5.     stack[p] =ps[index++];   
  6.     while(p>=2&& stack[p-2],stack[p-1],stack[p] 不是右手系){   
  7.          stack[p-1] =stack[p];   
  8.          p --;   
  9.     }   
  10.     p ++;
  11. }  

      显然,我们再从ps[k]到ps[0]做类似的处理就可以得到凸包的上半部分.当然,事实上我们可以把这两个工作一起做了:维持一个双头栈,使其头部的三个点为右手系,其尾部的三个点为左手系.这样经过一次扫描就可以得到整个凸包.伪代码如下:

java 代码
  1. //ps中保存了已经排好序的点   
  2. Point[] stack =new Point[2*ps.length];   
  3. int index =0, head =ps.length,tail =ps.length;   
  4. while(index
  5.     stack[head] =stack[tail] =ps[index++];   
  6.     while(head-ps.length>=2&& stack[head-2],stack[head-1],stack[head] 不是右手系){   
  7.          stack[head-1] =stack[head];   
  8.          head --;   
  9.     }   
  10.     while(ps.length-tail>=2&& stack[tail+2],stack[tail+1],stack[tail] 不是左手系){   
  11.          stack[tail+1] =stack[tail];   
  12.          tail ++;   
  13.     }   
  14.     head ++;   
  15.     tail --;   
  16. }  

            可以看出,整个算法并不复杂,或者可以说很优美.哦,等等,伪代码里面还有个判断"左手系","右手系"是怎么来的?这就涉及到一个数学概念:叉积

三:叉积
    叉积,又叫外积,向量积.这里我不对叉积做太深入的探讨,只做一些概念性的描叙(且没有考虑数学上的准确性,只是为了解决问题方便),有兴趣的可以找本大学的<解析几何>之类的数学书看看.
    定义:对于二维平面的两个向量a =[xa,ya],b =[xb,yb],其叉积 a×b =xa*yb -xb*ya (其实叉积的定义是在三维空间中的..)
    另一方面,对于点X1 =[x1,y1],X2 =[x2,y2],对应的向量 X1X2 =[x2-x1,y2-y1]
    性质:
    对于平面上的三个点A,B,C所组成的三角形ABC的面积是向量AB与AC叉积的绝对值的一半,即 2*S(A,B,C) =|AB×AC|,并且当A,B,C三点按逆时间顺转时,AB×CD为正;当A,B,C三点按顺时间顺转时,AB×AC为负;A,B,C三点共线时,其叉积为0.

 

四:具体代码
    有了上面的知识,下面给出具体的求凸包的Java代码.
    首先给出用于表示Point的类:

java 代码
  1. package convex;   
  2. public class Point implements Comparable＜Point＞{
  3.     public double x,y;   
  4.     public static double distanceSq(Point p1,Point p2){   
  5.         return (p2.x-p1.x)*(p2.x-p1.x)+(p2.y-p1.y)*(p2.y-p1.y);   
  6.     }   
  7.     public Point(){   
  8.         this(0.0,0.0);   
  9.     }   
  10.     public Point(double x,double y){   
  11.         this.x =x;   
  12.         this.y =y;   
  13.     }   
  14.     /**  
  15.     * 实现Comparable的compareTo方法,将Point按从左到右,从下到上排序  
  16.     */  
  17.     public int compareTo(Point p){   
  18.         if(x ==p.x) return y==p.y?0:(y>p.y?1:-1);   
  19.         else        return x>p.x?1:-1;   
  20.     }   
  21. }  

 

 下面给出求凸包的代码:
java 代码
  1. package convex;   
  2. import static java.util.Arrays.sort;   
  3. /**  
  4. * 用于求点集凸包以及求凸包直径的算法类  
  5. */  
  6. public class Algorithm{   
  7.     //避免生成该类实例   
  8.     private Algorithm(){}   
  9.     /**  
  10.     * 构造Point数组ps中从fromIndex到toIndex的Point的凸包,时间复杂度O(N*logN)  
  11.     * 注意: 该方法可能改变ps数组中从fromIndex到toIndex元素的顺序  
  12.     * 前置条件: 一个任意的二维点集,不同点的个数>=2  
  13.     *@param ps 保存二维点集的Point数组  
  14.     *@param fromIndex 点集在数组中的起始位置  
  15.     *@param toIndex   点集在数组中的结束位置  
  16.     *@param convex    用来保存凸多边形的顶点的Point数组,其顶点将按逆时针排列  
  17.     *@param offset    保存数据在convex数组中的起始位置  
  18.     *@return length   凸多边形的顶点数目  
  19.     */  
  20.     public static int getPointsConvexClosure(Point[] ps,int fromIndex,int toIndex,Point[] convex,int offset){   
  21.         sort(ps,fromIndex,toIndex);   
  22.         int len =toIndex -fromIndex;   
  23.         Point[] tmp =new Point[2*len];   
  24.         int up =len, down =len;   
  25.         for(int index =fromIndex;index
  26.             tmp[up] =tmp[down] =ps[index];   
  27.             while(len-up>=2&&multiply(tmp[up+2],tmp[up+1],tmp[up])>=0){   
  28.                 tmp[up+1] =tmp[up];    
  29.                 up++;   
  30.             }   
  31.             while(down-len>=2&&multiply(tmp[down-2],tmp[down-1],tmp[down])<=0){   
  32.                 tmp[down-1] =tmp[down];    
  33.                 down--;   
  34.             }   
  35.             up --;   
  36.             down ++;   
  37.         }   
  38.         System.arraycopy(tmp,up+1,convex,offset,down-up-2);   
  39.         return down-up-2;   
  40.     }   
  41.     /**
        * 计算向量ab与ac的叉积
        */
        private static double multiply(Point a,Point b,Point c){
            return multiply(a,b,a,c);
        }
        /**
        * 计算向量ab与cd的叉积
        */
        private static double multiply(Point a,Point b,Point c,Point d){
            return (b.x-a.x)*(d.y-c.y)-(d.x-c.x)*(b.y-a.y);
        }
  42. }  
  • convex.rar (3 KB)
  • 描述: 本文中使用的Java代码,并包含TJU2847的 Java代码
  • 下载次数: 25
评论
crackerwang 2007-08-23
弄了几天.硬着头皮看英文资料我终于把PKU上的2187给过了.
这次你真的不用写了.
非常谢谢.以后在请教你其他问题
crackerwang 2007-08-17
现在一搜索"凸包 直径"就到你这了.
你这个构造凸包的方法不错比那个算及角的要减少精度误差.
算直径的还不是很明白.
crackerwang 2007-08-17
你高估我了.
凸包那个算法还好.
不过还是有点小小的问题.
这个地方while(down-len>=2&&multiply(tmp[down-2],tmp[down-1],tmp[down])<=0){
tmp[down-1] =tmp[down];
down--;
}
改成while(down-len>=2&&multiply(tmp[down-2],tmp[down-1],tmp[down])<0)
求周长会不会有影响?
上次在做一个题目的时候发现加等号错不加就对.一直没有找到这个特殊的数据
Eastsun 2007-08-16
...
你说“终于明白了”。
我还以为我可以不用写第二篇了。

不过你还是先自己google一下吧，我现在不在学校不方便写。
crackerwang 2007-08-16
你要是有时间不如把凸包直径的求法也写一下吧
Elminster 2007-06-26
crackerwang 写道
貌似RMQ是一个很典型的算法.不知道为什么找不到具体的.都是找到一些什么分类的东西.The LCA Problem Revisted 这个我找到了.
你的RMQ是在哪本书上看到的..
难道是传说中的算法导论??


某年月日，某人问我：你可知道 RMQ 问题？
我说：不曾听过这个名字。
某人于是描述了一遍 RMQ 问题，并说：似乎有一个算法，可以做到预处理时间 O(n)，查询时间 O(1)。
我想了一回，便说：那定然是个很精巧的算法，我不能及。
某人说：请替我查访一下罢。
我说：好，我便查访一下罢。
于是我拜请 google 大神，照见一切以 RMQ 为名的先贤足迹，于 15 分钟后寻得 R.E.Tarjan、M.A.Bender 等先知留下的文献卷轴。于是 RMQ-LCA 问题的一切秘奥，均展现在我们的眼前。
于是某人说：赞美政委大能！算法的光辉必将照亮一切黑暗前路。
我说：善。
crackerwang 2007-06-26
貌似RMQ是一个很典型的算法.不知道为什么找不到具体的.都是找到一些什么分类的东西.The LCA Problem Revisted 这个我找到了.
你的RMQ是在哪本书上看到的..
难道是传说中的算法导论??
Elminster 2007-06-25
不过扫了一眼这道题：

http://acm.pku.edu.cn/JudgeOnline/problem?id=2823

这道题用 RMQ 应该是违规的。从题意来看，任意时刻你应该只能访问其中 k 个元素。而且事实上，这题用 RMQ 也是牛刀割鸡了。
Elminster 2007-06-24
crackerwang 写道
终于明白了..
但是我在做那个题目的时候是想构造一个最小外接圆,然后在求上面所有点的两两距离.最后因为没有看明白书上的dephi代码而做罢了..
看来eastsun大哥对算法和数据结构的功底不在一般计算机学生之下了..对你的佩服已经由如滔滔...(将来也要和你一样,嘿嘿)

我再问个典型算法的问题:RMQ-LCA,你会吗?这个好像是个很典型的算法,但是在网上却找不到很好的答案...你要是会就稍微在底下给我留言一下...(当然要是再写个BLOG就更加好了..嘿嘿)

具体的问题你可以参照这个:http://acm.pku.edu.cn/JudgeOnline/problem?id=2823

我用优先队列过是能过,但是效率就....



直接看 The LCA Problem Revisted 这篇论文好了，可以直接下载，很容易懂。
crackerwang 2007-06-24
当时(打拼音老是打错)我在做那个题目...(打拼音老是打错)
crackerwang 2007-06-24
终于明白了..
但是我在做那个题目的时候是想构造一个最小外接圆,然后在求上面所有点的两两距离.最后因为没有看明白书上的dephi代码而做罢了..
看来eastsun大哥对算法和数据结构的功底不在一般计算机学生之下了..对你的佩服已经由如滔滔...(将来也要和你一样,嘿嘿)

我再问个典型算法的问题:RMQ-LCA,你会吗?这个好像是个很典型的算法,但是在网上却找不到很好的答案...你要是会就稍微在底下给我留言一下...(当然要是再写个BLOG就更加好了..嘿嘿)

具体的问题你可以参照这个:http://acm.pku.edu.cn/JudgeOnline/problem?id=2823

我用优先队列过是能过,但是效率就....
Eastsun
搜索本博客
我的相册
Ba2baf55-8a44-300c-973d-5f0a8705818e-thumb
businessblacksteel1.png
共 60 张
最近加入圈子
存档
最新评论