计算几何_quick hull算法-CSDN博客

一.凸集&凸包

(下文中所有的集合若不作特殊说明都是指欧氏空间上的集合)

凸集(Convex Set):任意两点的连线都在这个集合内的集合就是一个凸集.

A set in Euclidean space is convex set if it contains all the line segments connecting any pair of its points.

http://mathworld.wolfram.com/Convex.html

凸包(Convex Hull):包含集合S的所有凸集的交集就是集合S的凸包.

The convex hull of a set of points S in N dimensions is the intersection of all convex sets containing S.

我们经常关注一个点集的凸包这也是计算几何学的一个基本问题

我们现在已经有成熟的算法可以求出平面点集的凸包和空间点集的凸包

甚至有的算法可以方便的求出任意维度欧氏空间内的一个点集的凸包

凸包有着优美而实用的性质我们可以利用凸包把一个杂乱的点集所包含的信息进行有效的概括梳理

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

二.平面点集的凸包的算法

(下文中所有凸包若不作特殊说明都是指平面点集的凸包)

有两种直观的理解凸包的方式

在木板上钉钉子用一个有弹性的橡皮筋去框住所有钉子

橡皮筋形成的图形就是这个钉子所构成的点集的凸包

还有一种理解我们用一根麻绳绑住一个外面的钉子然后拉着麻绳绕所有钉子一圈

这个麻绳最后也构成了点集的凸包

其中第二种理解是我们一个经典算法 卷包裹法(Gift Wrapping)的思路

卷包裹算法从一个必然在凸包上的点开始向着一个方向依次选择最外侧的点

当回到最初的点时所选出的点集就是所要求的凸包

这里还有两个问题不是很清楚:

1.怎么确定一个肯定在凸包上的点?

这个问题很好解决取一个最左边的也就是横坐标最小的点

如果有多个这样的点就取这些点里纵坐标最小的

这样可以很好的处理共线的情况

2.如何确定下一个点(即最外侧的点)?

我们需要利用向量的叉积来解决这个问题

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

向量的叉积(Cross Product)原本是三维空间中的问题在二维中也有巧妙的应用

http://mathworld.wolfram.com/CrossProduct.html

(下文中所有的叉积若不作特殊说明都是指二维中新定义的叉积

下文中所有的向量乘法若不作特殊说明都是指向量的叉积)

我们定义二维向量<x1,y1>和<x2,y2>的叉积为一个实数Cp=x1*y2-x2*y1

叉积有两条性质很常用也很好用

1.叉积的一半是一个三角形的有向面积

这个公式可以避免面积计算的误差如果点是整点那么所有运算都是整数

2.向量的叉积的符号代表着向量旋转的方向

向量的叉积是不满足交换律的

向量A乘以向量B 如果为正则为A逆时针旋转向B 否则为顺时针

当然这里A转向B的角总是考虑一个小于180度以内的角否则就会出错

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

有了向量我们就可以选取一个最外侧的点了

比如现在我们卷包裹卷到J点我们要选取一个最外侧的点

当然比较红色的到角可以直接得到最外侧的点不过不方便

我们可以考虑那个蓝色的到角

利用向量我们可以比较哪个点"更外侧"

比如点K和点I 我们利用向量JK乘以向量JI得到一个数这个数应该是负数说明I比K更外侧

两个向量的比较具有传递性 所以我们可以像N个数里取最大的数一样取出最外侧的

遍历所有点每个点都和现有最外侧的点比较得到新的最外侧的点

至此两个问题都得以解决我们可以写出满足一般要求的卷包裹算法了

不过还遗留有一个问题就是处理共线的问题

有时候我们需要凸包边上的点也考虑到有时候却需要去掉这些点

我们通常称在凸包顶点处的点为极点

如果我们只要求保留极点而去除在边上的点

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

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

这样整个卷包裹法终于完成了

还有两点要补充说明一下:

1.卷包裹算法的复杂度是O(NH)

N是全部的点数 H是最终在凸包上的点数

所以卷包裹算法很适合凸包上的点很少的时候通常随机数据很快

但是构造出的凸包上的点很多的数据这个算法就会很慢

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

2.卷包裹算法输出的点是有序的

这也是对二维凸包算法的一个基本要求

通常只有保证有序才能进行进一步的计算

通过改变CMP函数可以改变上文中提到的共线(取/不取)以及这里的序(顺时针/逆时针)的问题

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

三.凸包的面积和周长

凸包的一个简单算法介绍完了

来看一个具体的问题 Poj 1113 http://poj.org/problem?id=1113

问题是求到一个简单多边形距离最少为L的点的轨迹的周长

我们可以先求凸包然后得到计算公式

实际上圆的部分可以拼接起来变成一个整圆

这样就好算多了注意把共线的都去掉不然误差会很大的

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

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

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

内部循环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扫描所用时间在点集随机时还不如卷包裹算法快

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

伪代码如下:

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

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

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

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

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

这里有一个很好的演示

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

还要补充说明一下

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

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

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

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;
}

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

四.凸包算法复杂度下界

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