首先说明,此篇博客是我在初学凸包时将看到过的所有很好很经典博客里简单易懂或者精彩的部分摘取组合而成的。每一篇博客侧重点不同,有些讲的比较粗略但是另一些博客很详细,所以我觉得将这些精华聚合在一起还是很重要的,每篇博客我都会标明出处。另外有些地方我也会增加一些我自己的见解,如果有不对的地方,希望各位大佬指出。如果您觉得好的话欢迎转载,但请说明出处:
http://blog.csdn.net/qq_34374664?viewmode=contents
转自:http://blog.csdn.net/bone_ace/article/details/46239187
首先,什么是凸包?
假设平面上有p0~p12共13个点,过某些点作一个多边形,使这个多边形能把所有点都“包”起来。当这个多边形是凸多边形的时候,我们就叫它“凸包”。如下图:
然后,什么是凸包问题?
我们把这些点放在二维坐标系里面,那么每个点都能用 (x,y) 来表示。
现给出点的数目13,和各个点的坐标。求构成凸包的点?
解一:穷举法(蛮力法)
时间复杂度:O(n³)。
思路:两点确定一条直线,如果剩余的其它点都在这条直线的同一侧,则这两个点是凸包上的点,否则就不是。
步骤:
- 将点集里面的所有点两两配对,组成 n(n-1)/2 条直线。
- 对于每条直线,再检查剩余的 (n-2) 个点是否在直线的同一侧。
如何判断一个点 p3 是在直线 p1p2 的左边还是右边呢?(坐标:p1(x1,y1),p2(x2,y2),p3(x3,y3))
当上式结果为正时,p3在直线 p1p2 的左侧;当结果为负时,p3在直线 p1p2 的右边。
解二:分治法
时间复杂度:O(n㏒n)。
思路:应用分治法思想,把一个大问题分成几个结构相同的子问题,把子问题再分成几个更小的子问题……。然后我们就能用递归的方法,分别求这些子问题的解。最后把每个子问题的解“组装”成原来大问题的解。
步骤:
- 把所有的点都放在二维坐标系里面。那么横坐标最小和最大的两个点 P1 和 Pn 一定是凸包上的点(为什么呢?用反证法很容易证明,这里不详讲)。直线 P1Pn 把点集分成了两部分,即 X 轴上面和下面两部分,分别叫做上包和下包。
- 对上包:求距离直线 P1Pn 最远的点,即下图中的点 Pmax 。
- 作直线 P1Pmax 、PnPmax,把直线 P1Pmax 左侧的点当成是上包,把直线 PnPmax 右侧的点也当成是上包。
- 重复步骤 2、3。
- 对下包也作类似操作。
然而怎么求距离某直线最远的点呢?我们还是用到解一中的公式:
设有一个点 P3 和直线 P1P2 。(坐标:p1(x1,y1),p2(x2,y2),p3(x3,y3))
对上式的结果取绝对值,绝对值越大,则距离直线越远。
注意:在步骤一,如果横坐标最小的点不止一个,那么这几个点都是凸包上的点,此时上包和下包的划分就有点不同了,需要注意。
解三:Jarvis步进法
时间复杂度:O(nH)。(其中 n 是点的总个数,H 是凸包上的点的个数)
思路:
- 纵坐标最小的那个点一定是凸包上的点,例如图上的 P0。
- 从 P0 开始,按逆时针的方向,逐个找凸包上的点,每前进一步找到一个点,所以叫作步进法。
- 怎么找下一个点呢?利用夹角。假设现在已经找到 {P0,P1,P2} 了,要找下一个点:剩下的点分别和 P2 组成向量,设这个向量与向量P1P2的夹角为 β 。当 β 最小时就是所要求的下一个点了,此处为 P3 。
注意:
- 找第二个点 P1 时,因为已经找到的只有 P0 一个点,所以向量只能和水平线作夹角 α,当 α 最小时求得第二个点。
- 共线情况:如果直线 P2P3 上还有一个点 P4,即三个点共线,此时由向量P2P3 和向量P2P4 产生的两个 β 是相同的。我们应该把 P3、P4 都当做凸包上的点,并且把距离 P2 最远的那个点(即图中的P4)作为最后搜索到的点,继续找它的下一个连接点。
解四:Graham扫描法
时间复杂度:O(n㏒n)
思路:Graham扫描的思想和Jarris步进法类似,也是先找到凸包上的一个点,然后从那个点开始按逆时针方向逐个找凸包上的点,但它不是利用夹角。
步骤:
- 把所有点放在二维坐标系中,则纵坐标最小的点一定是凸包上的点,如图中的P0。
- 把所有点的坐标平移一下,使 P0 作为原点,如上图。
- 计算各个点相对于 P0 的幅角 α ,按从小到大的顺序对各个点排序。当 α 相同时,距离 P0 比较近的排在前面。例如上图得到的结果为 P1,P2,P3,P4,P5,P6,P7,P8。我们由几何知识可以知道,结果中第一个点 P1 和最后一个点 P8 一定是凸包上的点。
(以上是准备步骤,以下开始求凸包)
以上,我们已经知道了凸包上的第一个点 P0 和第二个点 P1,我们把它们放在栈里面。现在从步骤3求得的那个结果里,把 P1 后面的那个点拿出来做当前点,即 P2 。接下来开始找第三个点: - 连接P0和栈顶的那个点,得到直线 L 。看当前点是在直线 L 的右边还是左边。如果在直线的右边就执行步骤5;如果在直线上,或者在直线的左边就执行步骤6。
- 如果在右边,则栈顶的那个元素不是凸包上的点,把栈顶元素出栈。执行步骤4。
- 当前点是凸包上的点,把它压入栈,执行步骤7。
- 检查当前的点 P2 是不是步骤3那个结果的最后一个元素。是最后一个元素的话就结束。如果不是的话就把 P2 后面那个点做当前点,返回步骤4。
最后,栈中的元素就是凸包上的点了。
以下为用Graham扫描法动态求解的过程:
解五:Melkman算法
说真的,这个算法我也还没有看清。网上的资料也少的可怜,我暂且把网上的解释截个图在这里,往后搞懂以后再回来补上。
或者有人看懂了的,希望不吝指教,不甚感激!
扩展:
以上讨论的只是二维的凸包,如果延生为三维、多维的凸包问题呢?如何求解?
不过首先,二维凸包可以用来解决围栏问题、城市规划问题、聚类分析等等。但是三维、多维的凸包可能的使用范畴有哪些?
附:快包算法代码(C语言):
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
比较常用的:
转自: http://blog.csdn.net/tmljs1988/article/details/7259331凸包问题 —— Graham扫描法:
(1)找出点集p[]中最左下的点p1,把p1同点集中其他各点用线段连接,并计算这些线段与水平线的夹角,然后按夹角从小到大和按到p1的距离从近到远排序(夹角范围为 [0, 180)度,而且可以删除相同夹角且距离p1较近的点,保留最远点,这样可减少计算量。因为直线上的非端点不是凸包的极点,即如果p1,p2,p3在一条直线上,则只取凸点p1,p3。p2不在端点,故可以去掉),得到新的节点序列p1,p2,...pn.依次连接这些点,得到一个多边形(已经逆时针,有所进展,但还需去掉不在凸包上的点)。此时p1是凸包的边界起点,p2和pn也是最终凸包的顶点,p[n+1]=p1(看成循环的)
(2)删除p3,p4,...p[n-1]中不在凸包上的点:
先把p1,p2,p3入栈S中,再依次循环(i = 3 -> n-1),若栈顶的两个点和当前的点p[i]这三点连线的方向向顺时针方向偏转,表明是凹的,应删除,则栈顶元素出栈(要循环判断,即可能前面的仍是凹的,还需再出栈,举例如下图),直到向逆时针方向偏转或者栈内只有2个元素了(p1p2),就把当前点p[i]入栈。
最后栈中的元素就是最终凸包上的点。
分析:一般会从最左下点p1开始,根据所有点斜率中最小的求下一个凸包点pa,再根据pa的所有点斜率中最小的来求下一个凸包点pb,依此类推,但这样就是三重循环(p1,pa,pb是一次,p1,pa,pb内部的排序有2次,共三重循环),这样效率不高(其实这就是卷包裹法,复杂度为O(NH),其中N是全部的点数 H是最终在凸包上的点数)。Graham扫描法只取所有点对p1的斜率,后面的点充分利用该斜率的信息,并作某些处理,进行改进,以提高效率。这是由多到少,由多个到1个的方法,并充分利用已知条件。
下面图例也可说明栈扫描的过程(摘自http://www.cnblogs.com/Booble/archive/2011/03/10/1980089.html#2065991):
----------------
---------------
--------------------
以下截图代码摘自《ACM程序设计培训教程 吴昊 中国铁道出版社》 :
―――――――――――
下面为网上转载:
对于一个有三个或以上点的点集Q,Graham扫描法的过程如下:
令p0为Q中Y-X坐标排序下最小的点
设<p1,p2,...pm>为对其余点按以p0为中心的极角逆时针排序所得的点集(如果有多个点有相同的极角,除了距p0最远的点外全部移除
压p0进栈S
压p1进栈S
压p2进栈S
for i ← 3 to m
do while 由S的栈顶元素的下一个元素、S的栈顶元素以及pi构成的折线段不拐向左侧
对S弹栈
压pi进栈S
return S;
若计算周长。将边相加 zoj1453
若计算面积。穷举三点构成的三角形。。找最大zoj2419
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MaxNode 100015
int stack[MaxNode];
int top;
typedef struct POINT
{
int x;
int y;
}POINT;
POINT point[MaxNode];
void swap(POINT point[],int i,int j)
{
POINT tmp;
tmp=point[i];
point[i]=point[j];
point[j]=tmp;
}
double multi(POINT p1,POINT p2,POINT p0) //叉乘
{
return ((p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x));
}
double distence(POINT p1,POINT p2) //p1,p2的距离
{
return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}
int cmp(const void *a,const void *b)
{
POINT *c=(POINT *)a;
POINT *d=(POINT *)b;
double k=multi(*c,*d,point[0]);
if(k<0) return 1;
else if(k==0&&distence(*c,point[0])>=distence(*d,point[0])) return 1; //极角相同,比距离
else return -1;
}
void grahamscan(int n)
{
int i,u;
u=0;
for(i = 1;i<= n-1;i++) //找到最左下的点p0
if((point[i].y < point[u].y)||(point[i].y==point[u].y&&point[i].x < point[u].x))
u=i;
swap(point,0,u);
qsort(point+1,n-1,sizeof(point[0]),cmp); //按极角排序
for(i = 0;i <= 2;i++) stack[i] = i; //p0,p1,p2入栈
top=2;
for(i = 3;i <= n-1;i++) //最终凸包的各顶点的编号依次存在stack[]中。
{
while(multi(point[i],point[stack[top]],point[stack[top-1]])>=0) //弹出非左转的点
{
if(top==0)break;
top--;
}
top++;
stack[top] = i;
}
}
//求凸包的面积
double polygonArea(int n,POINT p[])
{
double area;
int i;
area = 0;
for(i = 1;i <= n;i++){
area += (p[stack[i - 1]].x * p[stack[i % n]].y - p[stack[i % n]].x * p[stack[i - 1]].y);
}
return fabs(area) / 2;
}
int main()
{
int i;
for(i=0;i<10;i++)
scanf("%d%d",&(point[i].x),&(point[i].y));
grahamscan(10);
printf("%lf\n",polygonArea(10,point));
}
快速凸包:转自: http://www.cnblogs.com/Booble/archive/2011/03/10/1980089.html
三.快速凸包算法(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写了...
#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)的
不过凸包在很多情况下会改善问题的平均复杂度
凸包上的点通常很少 所以这个问题也可以过。