HDUOJ 1392凸包graham算法

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
const int Max = 1100;

struct Point
{
    int x;
    int y;
};
int num;
Point p[Max];   ///原始点集
Point ch[Max];  ///凸包点集

///p0p1叉乘p0p2
int xmult(const Point &p0, const Point &p1, const Point &p2)
{
    return (p1.x-p0.x)*(p2.y-p0.y) - (p2.x-p0.x)*(p1.y-p0.y);
}

double dis(const Point &a, const Point &b)
{
    return sqrt((double)(a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

///逆时针顺序把极角排序
bool cmp(const Point &a, const Point &b)
{
    int ret = xmult(p[0],a,b);
    if(ret>0)
        return true;
    if(ret==0&&dis(p[0],a)<dis(p[0],b))//保留凸包上的点，共线取最近
        return true;
    return false;
}
//如果我们只要求保留极点而去除在边上的点

//我们只需在取外侧的点的时候 碰到共线的点取最远的

//相反 如果我们要保留所有在边上的点我们只需要在共线的点中取最近的

/**
p: 原始点序列
n: 原始点的个数
ch: 栈,用来保存凸包上的点
top: 凸包上点的个数
*/
void Graham(Point *p, int n, Point *ch, int &top)
{
    int i,k=0;
    for(i=1; i < n; i++)//取纵坐标最小的点，纵坐标相同取横坐标最小，此点必在凸包上
    {
        if(p[k].y>p[i].y || ((p[k].y==p[i].y)&&p[k].x>p[i].x))
            k = i;
    }
    swap(p[k],p[0]);//p[0]为此点
    sort(p+1,p+n,cmp);//
    top = 0;
    ch[top++] = p[0];
    ch[top++] = p[1];
    ch[top++] = p[2];
    for(i = 3; i < n; ch[top++]=p[i++])
        while(top>=2 && xmult(ch[top-2],ch[top-1],p[i])<0)    ///是否向左
            top--;
}

int main()
{
    //freopen("C:\\Users\\Sky\\Desktop\\1.in","r",stdin);
    int i;
    int n,l;
    int t;
    while(scanf("%d",&n),n)
    {
        for(i = 0; i < n; i++)
            scanf("%d %d",&p[i].x,&p[i].y);
        if(n==1)
        {
            printf("0.00\n");
            continue;
        }
        else if(n==2)
        {
            printf("%.2lf\n",dis(p[0],p[1]));
            continue;
        }
        Graham(p,n,ch,num);
        double result=0.0;
        for(i = 0; i < num-1; i++)
            result += dis(ch[i],ch[i+1]);
        result += dis(ch[0],ch[num-1]);
        printf("%.2lf\n",result);
    }
    return 0;
}

View Code

关于凸包博文，

引用出处http://www.cnblogs.com/Booble/default.html?page=1作者Master_Chivu

====================================================================

一.卷包裹算法(Gift Wrapping Algorithm)的特性

前面提到过卷包裹算法的复杂度问题

由于卷包裹算法是两重循环实现的 因此很好分析它的复杂度

1 while true do 2 begin 3 k: = 0 ; 4 inc(m); ch[m]: = j; 5 for i: = 1 to n do 6 if (i <> j) and ((k = 0 ) or 7 cmp(p[j],p[k],p[i])) 8 then k: = i; 9 if k = temp then break; 10 j: = k; 11 end ;

内部循环N次 外循环的次数决定于凸包上的点数H

所以卷包裹算法复杂度为O(HN) 这个复杂度很有特点

考虑一个随机点集的凸包上的点数往往很少 但是最坏情况下是O(N)级别的

比如所有的点都在一个圆上的时候

这就决定了卷包裹算法很适合随机的点集 但是如果是刻意构造的数据或是比较特殊的数据

就会达到O(N^2)的最坏复杂度 这是我们不愿意看到的

下面介绍最坏情况下复杂度更好的Graham扫描算法(Graham Scan Algorithm)

====================================================================

二.Graham扫描算法(Graham Scan Algorithm)

Graham扫描算法维护一个凸壳 通过不断在凸壳中加入新的点和去除影响凸性的点 最后形成凸包

The Graham scan is a method of computing the convex hull of a finite set of points in the plane with time complexity O(n log n). It is named after Ronald Graham, who published the original algorithm in 1972.[1] The algorithm finds all vertices of the convex hull ordered along its boundary.

http://en.wikipedia.org/wiki/Graham_scan

算法主体由两部分组成 先是排序 后是扫描 分块讲解一下

----------------------------------------------------------------------------------------------------------

1.点集排序

为了得到加入新点的顺序 Graham扫描法的第一步是对点集排序

排序是对杂乱的点集进行了梳理 这也是这种算法能够得到更高效率的根本原因

排序的方法也有两种 极角坐标排序(极角序) 和 直角坐标排序(水平序)

前者好理解一些 但是在实现的时候 后者更方便

先说极角序 为了极角排序 我们先得得到一个参考点

一般的 我们取最左边(横坐标最小)的点作为参考点 如果有多个这样的点就取最下面的(纵坐标最小)

看这样一个例子 这是一个任意给出的平面点集:

参考点的定义:在横坐标最小的情况下取纵坐标最小的点

所以所有的点只能在这个黄色的半平面中 而且正上方为闭(可取得) 正下方为开(不可取)

这就决定了参考点的性质:点集中任意两点和参考点所成的到角为锐角

这样我们取得参考点 然后再考虑极角排序

极角排序以参考点为极角坐标系原点 各个点的极角为关键字

由于上面我们得到的参考点的性质 我们可以设所有点的极角均在(-90,90]之间

排序完成后应该是这样的:

比较极角我们仍然可以利用向量的叉积

叉积在这里已经介绍了 http://www.cnblogs.com/Booble/archive/2011/02/28/1967179.html

同样由于参考点的性质 所有向量之间的到角都是在180度以内 不会产生错误

----------------------------------------------------------------------------------------------------------

2.Graham的栈扫描

Graham的扫描是一个很优美的过程 用到的数据结构也很简单 仅仅是一个栈而已

核心的思想是按照排好的序 依次加入新点得到新的边

如果和上一条边成左转关系就压栈继续 如果右转就弹栈直到和栈顶两点的边成左转关系 压栈继续

实现的时候我们不用存边 只需要含顺序在栈里存点 相邻两点就是一条边

由于我们时时刻刻都保证栈内是一个凸壳 所以最后扫描完毕 就得到了一个凸包

下面还是继续上面的那个样例 演示一下栈扫描的过程

这样Graham扫描算法基本完成

复杂度是排序O(Nlog2N) 扫描O(N) {每个点仅仅出入栈一次}

合起来是一个O(Nlog2N)的算法 很优秀

 ----------------------------------------------------------------------------------------------------------

3.双重共线点难题

和卷包裹算法一样 我们同样还要考虑共线点问题

而且Graham扫描算法的共线问题更复杂 所以需要仔细考虑

i).排序时的共线问题

如果极角相同 我们应该怎么定先后呢?

我们得加上第二关键字距离 比如极角相同 距离参考点近的先

不过不管是近的先还是 远的先 开始和结束的两条边总是矛盾的

我们必须对其中一条特殊处理 除了结束边外距离近的先 结束边上距离远的先

这就是为什么极角排序不是很好实现的原因了 下面会介绍一下水平序

ii).扫描时的共线问题

这个和卷包裹算法的处理放法如出一辙

如果需要保留共线的点就在到角相同时取距离最近的

如果仅仅需要凸包极点就取距离最远的

----------------------------------------------------------------------------------------------------------

4.直角坐标排序(水平序)

直角坐标排序方法没有了极角排序的不足

以横坐标为第一关键字 纵坐标为第二关键字 排序点集

然后从第一个点开始 分别利用Graham扫描生成左链和右链

需要注意以下两点:

i).左链和右链的旋转方向是相反的

ii).注意生成第二条链要忽略第一条链上的点

由于水平序的CMP函数比较简单 代码也更短

还有一件有意思的事 Graham扫描算法是1972年提出的 卷包裹算法是1973年提出的

其实不奇怪 这两个算法本来也没有优劣 所用的思想不同 各自善于处理不同的情况

Graham扫描所用时间在点集随机时 还不如卷包裹算法快

正如第一节所说 卷包裹算法已经可以处理大部分情况 也是一个可取的算法

实际应用中点集更趋向于均匀 而不是集中在凸包上

----------------------------------------------------------------------------------------------------------

最后给一下比较难实现的极角排序代码

GrahamScan

{ $inline on } const zero = 1e-6 ; maxn = 100000 ; type point = record x,y:extended; end ; var p,ch: array [ 1 ..maxn] of point; i,j,n,m,temp:longint; function sgn(x:extended):longint; inline; begin if abs(x) < zero then exit( 0 ); if x < 0 then sgn: =- 1 else sgn: = 1 ; end ; function cross(a,b,c,d:point):extended; inline; begin cross: = (b.x - a.x) * (d.y - c.y) - (d.x - c.x) * (b.y - a.y); end ; function dist(a,b:point):extended; inline; begin dist: = sqr(a.x - b.x) + sqr(a.y - b.y); end ; function cmp1(a,b,c:point):boolean; inline; var temp:longint; begin temp: = sgn(cross(a,b,a,c)); if temp <> 0 then exit(temp > 0 ); { *B } cmp1: = dist(a,b) < dist(a,c); end ; function cmp2(a,b,c:point):boolean; inline; var temp:longint; begin temp: = sgn(cross(a,b,b,c)); if temp <> 0 then exit(temp < 0 ); { *B } cmp2: = dist(a,b) < dist(a,c); { *A } end ; procedure swap( var a,b:point); inline; var temp:point; begin temp: = a; a: = b; b: = temp; end ; procedure sort(l,r:longint); var i,j:longint; x:point; begin i: = l; j: = r; x: = p[(l + r) >> 1 ]; repeat while cmp1(p[ 1 ],p[i],x) do inc(i); while cmp1(p[ 1 ],x,p[j]) do dec(j); if not (i > j) then begin swap(p[i],p[j]); inc(i); dec(j); end ; until i > j; if i < r then sort(i,r); if l < j then sort(l,j); end ; begin assign(input, ' Hull.in ' ); reset(input); assign(output, ' Hull2.out ' ); rewrite(output); readln(n); for i: = 1 to n do begin readln(p[i].x,p[i].y); if (p[i].x < p[ 1 ].x) or (sgn(p[i].x - p[ 1 ].x) = 0 ) and (p[i].y < p[ 1 ].y) then swap(p[ 1 ],p[i]); end ; sort( 2 ,n); while i > 2 do begin if sgn(cross(p[ 1 ],p[i],p[ 1 ],p[i - 1 ])) <> 0 then break; dec(i); end ; temp: = (i + n + 1 ) >> 1 ; for j: = n downto temp do swap(p[j],p[i + n - j]); m: = 1 ; ch[ 1 ]: = p[ 1 ]; inc(n); p[n]: = p[ 1 ]; for i: = 2 to n do begin while m > 1 do begin if not cmp2(ch[m - 1 ],ch[m],p[i]) then break; dec(m); end ; inc(m); ch[m]: = p[i]; end ; for i: = 1 to m - 1 do writeln(ch[i].x: 0 : 2 , ' ' ,ch[i].y: 0 : 2 ); close(input); close(output); end .
* A Change Direction * B Remove Colinear Points

====================================================================

三.快速凸包算法(Quickhull Algorithm)

对比Graham扫描算法和卷包裹算法

我们发现 Graham扫描算法在凸包上的点很密集时仍然适用

卷包裹算法在凸包上点集随机分布时是很高效的

那么有没有两个优点都具备的算法呢?

是有的! 快速凸包算法(Quickhull Algorithm)就是这样的一个算法

快速凸包算法是一个和快速排序(Quicksort Algorithm)神似的算法

尽管快速排序的最坏复杂度可以达到O(N^2)

但是有着极小的常数 实现方便 思路优美 绝大多数情况特别高效的快速排序 还是赢得了更多人的青睐

快速凸包算法也是这样 尽管可以构造一个数据使之达到O(N^2)的复杂度

但这需要刻意针对程序经过分析 才能做到 是实际应用中根本不会碰到的情况

在点集均匀分布时 快速凸包的复杂度更是达到了O(N) 是上面两种算法难以企及的

在绝大多数情况下 平均复杂度是O(Nlog2N) 也很高效

快速凸包继承了快速排序分治的思想 这是一个递归的过程

------------------------------------------------------------------------------------

伪代码如下:

 

1 void 快速凸包(P:点集 , S:向量 /* S.p,S.q:点) */ ){ 2 　　 /* P 在 S 左侧 */ 3 　　选取 P 中距离 S 最远的 点 Y ; 4 　　向量 A <- { S.p , Y } ; 向量 B <- { Y , S.q } ; 5 　　点集 Q <- 在 P 中 且在 A 左侧的点 ; 6 　　点集 R <- 在 P 中 且在 B 左侧的点 ; /* 划分 */ 7 　　快速凸包 ( Q , A ) ; /* 分治 */ 8 　　输出 (点 Y) ; /* 按中序输出 保证顺序 */ 9 　　快速凸包 ( P , B ) ; /* 分治 */ 10 }

初始化就选取 最左下和最右上的点 划分好 然后调用两次快速凸包函数分别求上下凸包

其中 选取和划分都需要用到向量的叉乘 注意方向

另外 存储点集可以用数组里的连续一段 传参的时候就传左右下标

划分的时候就是给数组交换元素 使新的点集 变成连续的一段

这里有一个很好的演示

http://www.cs.princeton.edu/courses/archive/spr10/cos226/demo/ah/QuickHull.html

还要补充说明一下

快速凸包在所有点都在圆周上的时候还是O(Nlog2N) 不过会比Graham扫描算法慢一些

可以说 这种数据是Graham扫描法的最好情况 一遍走完就行了

构造快速凸包的最坏情况就是使划分不均等 和构造快速排序最坏情况一样

贴以下我很辛苦写出来的代码吧 为了好看些 还心血来潮用cpp写了...

QuickHull

#include < iostream > #include < math.h > #define maxn 100000 #define zero 1e-12 #define sgn(x) (fabs(x)<zero?0:(x>0?1:-1)) #define cross(a,b,c) ((b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x)) #define cmp(a,b) (a.x<b.x || sgn(a.x-b.x)==0 && a.y<b.y)
using namespace std; struct point{ double x,y; }p[maxn]; double s[maxn];
void hull( int l, int r,point a,point b){ int x = l,i = l - 1 ,j = r + 1 ,k; for (k = l;k <= r;k ++ ){ double temp = sgn(s[x] - s[k]); if (temp < 0 || temp == 0 && cmp(p[x],p[k])) x = k; } point y = p[x]; for (k = l;k <= r;k ++ ){ s[ ++ i] = cross(p[k],a,y); if (sgn(s[i]) > 0 ) swap(p[i],p[k]); else i -- ; } for (k = r;k >= l;k -- ){ s[ -- j] = cross(p[k],y,b); if (sgn(s[j]) > 0 ) swap(p[j],p[k]); else j ++ ; } if (l <= i) hull(l,i,a,y); printf( " %.4lf %.4lf\n " ,y.x,y.y); if (j <= r) hull(j,r,y,b); }
int main(){ freopen( " CH2D.in " , " r " ,stdin); freopen( " CH2D1.out " , " w " ,stdout); int n,i,x = 0 ; scanf( " %d " , & n); for (i = 1 ;i <= n;i ++ ){ scanf( " %lf %lf " , & p[i].x, & p[i].y); if (x == 0 || cmp(p[i],p[x])) x = i; } swap(p[ 1 ],p[x]); printf( " %.4lf %.4lf\n " ,p[ 1 ].x,p[ 1 ].y); hull( 2 ,n,p[ 1 ],p[ 1 ]); return 0 ; }

====================================================================

四.凸包算法复杂度下界

(引自<算法艺术与信息学竞赛>)

====================================================================

五.高效算法的应用

这里有一个平面最远点对的问题 可以利用凸包解决

Poj 2187 http://poj.org/problem?id=2187

由于最远点对必然在凸包上

我们先求凸包 然后枚举凸包上的点 不过这个复杂度最坏是O(N^2)的

不过凸包在很多情况下会改善问题的平均复杂度

凸包上的点通常很少 所以这个问题也可以过

篇幅关系 在下一节会介绍降低最坏复杂度的旋转卡壳算法

====================================================================

文中部分图片来源

http://westhoffswelt.de/blog/0040_quickhull_introduction_and_php_implementation.html

http://en.wikipedia.org/wiki/Graham_scan

==============================

一.简单枚举算法的不足

上一次介绍了一个基本的求平面最远点对的算法

即先求点集的凸包 然后枚举凸包上的点来求最远点集

这是利用了凸包上的点相比 点集中的点 一般是很少的 平均情况很好 并且我们也能AC这个问题

但是这是有局限性的 当凸包上的点达到O(N)的级别时 凸包的优化作用就不存在了

不过我们还要考虑到 凸包还起了对凸包上点集排序的作用

凸包有很多的优美的性质 我们可以加以利用 以得到更加高效的算法

旋转卡壳算法就是利用凸包特性的一类解决问题的方法

==============================

二.旋转卡壳算法

旋转卡(qiǎ)壳算法(Rotating Calipers Algorithm):

是解决一些与凸包有关问题的有效算法 就像一对卡壳卡住凸包旋转而得名

Every time one blade of the caliper lies flat against an edge of the polygon, it forms an antipodal pair with the point or edge touching the opposite blade. It turns out that  the complete "rotation" of the caliper around the polygon detects all  antipodal pairs and may be carried out in O(n) time.

http://en.wikipedia.org/wiki/Rotating_calipers

(图片来自:http://cgm.cs.mcgill.ca/~orm/rotcal.html)

被一对卡壳正好卡住的对应点对称为对踵点(Antipodal point)

http://en.wikipedia.org/wiki/Antipodal_point

可以证明对踵点的个数不超过3N/2个 也就是说对踵点的个数是O(N)的

对踵点的个数也是我们下面解决问题时间复杂度的保证

上第一个图是卡壳的一般情况 卡住两点 图二是卡住一条边和一个点

由于实现中 卡住两点的情况不好处理 我们通常关注第二种情况

在第二种情况中 我们可以看到 一个对踵点和对应边之间的距离比其他点要大

也就是一个对踵点和对应边所形成的三角形是最大的 下面我们会据此得到对踵点的简化求法

看一下官方的伪代码:

当时我看完了 就一个字 长... 我最讨厌冗长的程序了...

begin p0: = pn; q: = NEXT[p]; while (Area(p,NEXT[p],NEXT[q]) > Area(p,NEXT[p],q)) do q: = NEXT[q]; q0: = q; while (q ! = p0) do begin p: = NEXT[p]; Print(p,q); while (Area(p,NEXT[p],NEXT[q]) > Area(p,NEXT[p],q) do begin q: = NEXT[q]; if ((p,q) ! = (q0,p0)) then Print(p,q) else return end ; if (Area(p,NEXT[p],NEXT[q]) = Area(p,NEXT[p],q)) then if ((p,q) ! = (q0,p0)) then Print(p,NEXT[q]) else Print(NEXT[p],q) end end .

几经折腾 终于找到了一个不错的实现:http://www.cnblogs.com/DreamUp/archive/2010/09/16/1828131.html

不过不是很好理解 这里作一下说明

1 ch[m + 1 ]: = ch[ 1 ]; j: = 2 ; 2   for i: = 1 to m do 3 begin 4 while cross(ch[i],ch[j],ch[i + 1 ]) < cross(ch[i],ch[j + 1 ],ch[i + 1 ]) do 5 begin inc(j); if j > m then j: = 1 ; end ; 6 writeln(ch[i].x,' ',ch[i].y, ' ' ,ch[j].x,' ',ch[j].y); 7 end ;

上面就是旋转卡壳寻找对踵点的过程

其中叉积函数Cross(A,B,C:Point):Real 返回AB到AC的二维定义下的叉积

这里主要用到了叉积求三角形面积的功能

我们对于一条对应边<CH i,CH Next[i]>求出距离这条边最远的点CH j

则由上面第二种情况可知 CH i 和 CH j 为一对对踵点 这样让 CH i 绕行凸包一周即可得到所有的对踵点

下面面这个图 由于本人的gif图制作水平拙劣 所以不好看

需要的可以下载几何画板察看原版GSP文件 点击这里下载GSP文件

接下来考虑 如何得到距离每条对应边的的最远点呢?

稍加分析 我们可以发现 凸包上的点依次与对应边产生的距离成单峰函数

具体证明可以从凸包定义入手 用反证法解决

这样我们再找到一个点 使下一个点的距离小于当前的点时就可以停止了

而且随着对应边的旋转 最远点也只会顺着这个方向旋转 我们可以从上一次的对踵点开始继续寻找这一次的

由于内层while循环的执行次数取决于j增加次数 j最多增加O(N)次

所以求出所有对踵点的时间复杂度为O(N)

还有有两点需要注意:

1.上面这段代码及代码的分析都是需要凸包上没有三点共线的

2.Next[i] 不需要手动求 在原代码中有很好的处理

最后指出网上很多文章的一个错误 一个点的对踵点并不是离这个点最远的点!

这样子的点对是根本不满足对踵点的性质的 即最为重要的单峰分布性质

下图是一个反例:

==============================

三.旋转卡壳算法的简单应用

至此我们终于可以更高效的解决平面最远点对问题了

有一个很重要的结论是 最远点对必然属于对踵点对集合

那么我们先求出凸包 然后求出对踵点对集合 然后选出距离最大的即可

用这个算法可以47ms AC这个问题 算上凸包的时间 总复杂度为O(Nlog2N)

代码如下:

Maxd

{ $inline on } { $optimization on } const maxn = 50000 ; type point = record x,y:longint; end ; var n,i,x,m,ans,j:longint; ch,p: array [ 1 ..maxn + 1 ] of point; s: array [ 1 ..maxn] of longint; function cross(a,b,c:point):longint; inline; begin cross: = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x); end ; function dist(a,b:point):longint; inline; begin dist: = sqr(a.x - b.x) + sqr(a.y - b.y); end ; function cmp(a,b:point):boolean; inline; begin cmp: = (a.x < b.x) or (a.x = b.x) and (a.y < b.y); end ; function max(a,b:longint):longint; begin if a > b then max: = a else max: = b; end ; procedure swap(a,b:longint); inline; var x:point; begin x: = p[a]; p[a]: = p[b]; p[b]: = x; end ; procedure hull(l,r:longint; a,b:point); var x,i,j,k:longint; y:point; begin x: = l; y: = p[l]; for k: = l to r do if (s[x] < s[k]) or (s[x] = s[k]) and (cmp(y,p[k])) then begin x: = k; y: = p[k]; end ; i: = l - 1 ; j: = r + 1 ; for k: = l to r do begin inc(i); s[i]: = cross(p[k],a,y); if s[i] > 0 then swap(i,k) else dec(i); end ; for k: = r downto l do begin dec(j); s[j]: = cross(p[k],y,b); if s[j] > 0 then swap(j,k) else inc(j); end ; if l <= i then hull(l,i,a,y); inc(m); ch[m]: = y; if j <= r then hull(j,r,y,b); end ; begin assign(input, ' Maxd.in ' ); reset(input); assign(output, ' Maxd.out ' ); rewrite(output); readln(n); for i: = 1 to n do begin readln(p[i].x,p[i].y); if (x = 0 ) or cmp(p[i],p[x]) then x: = i; end ; swap( 1 ,x); m: = 1 ; ch[ 1 ]: = p[ 1 ]; hull( 2 ,n,p[ 1 ],p[ 1 ]); ch[m + 1 ]: = ch[ 1 ]; j: = 2 ; ans: = 0 ; for i: = 1 to m do begin while cross(ch[i],ch[j],ch[i + 1 ]) < cross(ch[i],ch[j + 1 ],ch[i + 1 ]) do begin inc(j); if j > m then j: = 1 ; end ; ans: = max(ans,dist(ch[i],ch[j])); end ; writeln(ans); close(input); close(output); end .

引用出处http://www.cnblogs.com/Booble/default.html?page=1作者Master_Chivu

转载于:https://www.cnblogs.com/Skyxj/p/3188740.html