`
Eastsun
  • 浏览: 303500 次
  • 性别: Icon_minigender_1
  • 来自: 天津
社区版块
存档分类
最新评论

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

阅读更多

前言 :因为前几天做了一个有关凸包的题,并答应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代码
  • 下载次数: 65
分享到:
评论
12 楼 crackerwang 2007-08-23  
弄了几天.硬着头皮看英文资料我终于把PKU上的2187给过了.
这次你真的不用写了.
非常谢谢.以后在请教你其他问题
11 楼 crackerwang 2007-08-17  
现在一搜索"凸包 直径"就到你这了.
你这个构造凸包的方法不错比那个算及角的要减少精度误差.
算直径的还不是很明白.
10 楼 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)    
求周长会不会有影响?
上次在做一个题目的时候发现加等号错不加就对.一直没有找到这个特殊的数据
9 楼 Eastsun 2007-08-16  
...
你说“终于明白了”。
我还以为我可以不用写第二篇了。

不过你还是先自己google一下吧,我现在不在学校不方便写。
8 楼 crackerwang 2007-08-16  
你要是有时间不如把凸包直径的求法也写一下吧 
7 楼 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 问题的一切秘奥,均展现在我们的眼前。
于是某人说:赞美政委大能!算法的光辉必将照亮一切黑暗前路。
我说:善。
6 楼 crackerwang 2007-06-26  
貌似RMQ是一个很典型的算法.不知道为什么找不到具体的.都是找到一些什么分类的东西.The LCA Problem Revisted 这个我找到了.
你的RMQ是在哪本书上看到的..
难道是传说中的算法导论??
5 楼 Elminster 2007-06-25  
不过扫了一眼这道题:

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

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

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

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

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



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

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

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

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

1 楼 crackerwang 2007-06-24  
哈哈..
你终于弄上来了...
学习一下先 

相关推荐

Global site tag (gtag.js) - Google Analytics