计算几何及其应用——凸包问题
其实与计算几何中的最小圆覆盖问题很类似,凸包问题探究的是如何构造可以覆盖给定点集最小的凸多边形。
我们先从人脑的思维来分析一下这个问题,所谓凸包,起名字包含了两个关键的信息。
1.凸:这里所求作的是凸多边形,这是很关键的一点。因为在构造的时候可能会有下图的疑问。
右边的图的面积岂不是更小?但是我们还是要构造出左边的图,因为题设给出了——是凸多边形。
2.包(这里谈Jarvis算法):这个字的含义在于这个凸多边形必须要把给定的点集都给包住,这个条件使我们思考如何构造这个凸多边形的开始。既然想包住,我们显然要从点集中最靠外围的点开始构造(这显而易见),这里我们从最靠左下的点A开始。因为这个点已经是最下面的了,所以我们做出一条以A为端点的水平线射线,在这条射线下是不可能有点的, 我们开始在[0,180)旋转这个射线来“扫描”符合要求的点,最先出现的点(设为B)便是我们想要的,于是我们再以这个点为端点再做一条射线,开始循环操作。这里在扫描点的时候,无非会出现下面两种情况:
可以看到,在定量描述A、B、C的相对位置的时候,我们引入的向量及其叉积。因为在一般情况下我们都会知道点的坐标。由此我们就有了“选点”的依据。在这种情况下,我们可以保证选出来的点构成的凸多边形一定能包住所有的点。
有了上面对凸包概念更深层次的理解,我们来结合一个问题更深入的探讨凸包问题的应用。
阅读完此题后发现,这里其实是一个间接应用凸包的一个问题,关键就在于题设中说围墙必须离城堡不少于L。我们进行简单地作图会得到下面的图形。
我们以每个节点为圆心,L为半径做出了n个圆,然后将凸包的各个边进行平移,便会得到外围带有圆弧的图形。上图形成了一个类环图形。我们这里采用高等数学中极限的思想进行理解:内环和外环都有无数个点,并且他们一一对应,而且对应点之间的距离都是L。如果某组对应点之间的距离大于L,那么在外环中那个点的对应位置会形成一个小小的"凸起",显然周长就增大了。因此这种方式得到的围墙周长一定是最小的。以上是对“最小周长”的一个简单的证明。
有了这层证明,我们会发现,外围的周长其实是凸包的周长加上以L为半径的周长,这就把整个问题转化到了凸包问题上来。
构造凸包的算法实现:
1.上文提及构造凸包需要从最外围的点开始选,因此在输入点之后进行排序是必不可少的。
2.有了上述关于选点的叉积判别式,不难从点集中筛选出凸包的各个顶点。
代码如下。
#include<cmath> #include<algorithm> #include<stdio.h> #define size 1000 using namespace std; struct pint { int x,y; }x[size]; int n,l,ans[size],cnt,sta[size],tail; bool cmp(pint a,pint b) { return (a.y < b.y || (a.y==b.y && a.x<b.x)); //对点集进行排序 } bool CrossLeft(pint p1,pint p2,pint p3) { //比较是否是左转,如果共线不算是左转 return ((p3.x-p1.x)*(p2.y-p1.y)-(p2.x-p1.x)*(p3.y-p1.y))<0; //选点 } void make_Convex_hull() { tail = cnt = 0; sort(x,x + n,cmp); sta[tail++] = 0; sta[tail++] = 1; for(int i = 2;i < n;i++) //对排序后的点集进行遍历 { while(tail > 1 && !CrossLeft(x[sta[tail-1]] , x[sta[tail-2]] , x[i])) tail--; sta[tail++]=i; //这一层的循环会在点集的最高点处结束,此时凸包构造了一半 } for(int i = 0;i < tail;i++) //记录筛选出来的点 ans[cnt++] = sta[i]; tail = 0; sta[tail++] = n - 1; sta[tail++] = n - 2; for(int i = n-3 ;i >= 0;i--) //从最高点开始继续进行构造 { while(tail>1 && !CrossLeft(x[sta[tail-1]],x[sta[tail-2]],x[i])) tail--; sta[tail++]=i; //记录筛选出来的点 } for(int i=0;i<tail;i++) ans[cnt++]=sta[i]; } int main() { int t; while(scanf("%d",&t)!=EOF) { while(t--) { scanf("%d%d",&n,&l); for(int i = 0;i < n;i++) scanf("%d%d",&x[i].x , &x[i].y); make_Convex_hull(); double re=4 * acos(0.0) * l; for(int i = 0;i < cnt - 1;i++) re += sqrt((x[ans[i]].x - x[ans[i + 1]].x) *( x[ans[i]].x - x[ans[i + 1]].x) +(x[ans[i]].y - x[ans[i + 1]].y) * (x[ans[i]].y - x[ans[i + 1]].y)); printf("%.0lf\n",re); if(t) printf("\n"); } } return 0; }
然我们再看一道用Jarvis算法实现的凸包问题。
参考代码。
#include<iostream> #include<algorithm> #include<stdio.h> #include<cmath> using namespace std; const int Max = 105; struct Point { double x , y; }p[Max]; int n ,res[Max],top; bool cmp(Point a,Point b) { if(a.y == b.y) return a.x < b.x; return a.y < b.y; } bool mult(Point sp , Point ep , Point op) { return(sp.x - op.x)*(ep.y-op.y) >= (ep.x-op.x)*(sp.y-op.y); } void Graham() { int i , len; top = 1; sort(p,p+n,cmp); if(n == 0) return; res[0] = 0; if(n == 1) return; res[1] = 1; if(n == 2) return; res[2] = 2; for(i = 2;i < n;i++) { while(top && mult(p[i],p[res[top]],p[res[top-1]])) top--; res[++top] = i; } len = top; res[++top] = n - 2; for(i = n - 3;i >= 0;i--) { while(top != len && mult(p[i],p[res[top]],p[res[top-1]])) top--; res[++top] = i; } }
今天来讨论构建凸包的另一个算法——Graham算法。
相对Jarvis算法选取先选最下面的点,这里Graham算法先选取最左端的点,然后对点集进行极角排序(Jarvis算法按照纵坐标大小来排序,而Graham算法是按照该点与极点连线的斜率进行排序。),然后模拟栈的机理,起始的时候在栈中放入两个点,然后加入下一个点,判断是否满足左转(否则无法形成凸角),不满足的话剔除掉栈顶端的点,然后加入下一个点。遍历点集之后,此时栈中的顶部的点一定是点集中最靠上的点,然后再设置一次遍历,完成对凸包的构造。
参考下图可以更好地理解这个算法的过程。
我们来通过一个具体的凸包问题来体会一下这种算法的编程实现。(Problem source : hdu1392)
题目大意:给你一些树的点,现在需要用一条管道把所有的数都圈起来,这里让你计算需要多长的管道。
编程实现:我们可以看到,这是很明显的凸包问题。而这里我们考虑有Graham算法来实现它,其编程实现的一个重点就是极角排序和模拟栈选凸包点的过程,完成了这两个重要的过程,所求的凸包就记录在我们开的栈空间中,然后就不难求出周长了。
极角排序:根据其定义,我们可以轻松的将其转化成叉积的形式。
参考代码如下。
#include<stdio.h> #include<algorithm> #include<cmath> using namespace std; const int maxn = 110; const double eps = 1e-8; int n , stack[maxn]; struct point { double x , y; } p[maxn]; double dis(point a , point b) { return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)); } double Cross_left(point a , point b , point c) { return (a.x - c.x)*(b.y - c.y) - (a.y - c.y)*(b.x - c.x); } bool cmp(point a , point b) //极角排序 { double x = Cross_left(a , b , p[0]); if(x > 0) return true; else if(fabs(x) <= eps && (dis(b,p[0]) >= dis(a,p[0]))) return true; else return false; } double Graham() { int i , k = 0 , top = 2; point pmin=p[0] ; for(i = 1;i < n;i++) { if(p[i].y < pmin.y || (p[i].y - pmin.y <= eps && p[i].x < pmin.x)) { pmin = p[i] , k = i; //找极点 } } swap(p[0] , p[k]); sort(p+1,p+n,cmp); for(i = 0;i < 3;i++) stack[i] = i; for(i = 3;i < n;i++) { while(Cross_left(p[stack[top]] , p[i] , p[stack[top-1]]) <= 0)//共线情况不入栈 { top--; if(top == 0) break; } stack[++top] = i; } double cir = 0;top++; for(i = 0;i < top;i++) cir += dis(p[stack[i]] , p[stack[(i+1)%top]]); return cir; } int main() { while(scanf("%d",&n)!= EOF && n) { for(int i = 0;i < n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); if(n == 1) printf("0.00\n"); else if(n == 2) printf("%.2lf\n",dis(p[0],p[1])); else printf("%.2lf\n",Graham()); } }
基于我们对凸包问题的探讨,我们再来看一道与凸包问题有关的题目。(Problem source : hdu 4978)
题目大意:平面内有无数条间隔为D的平行线,现在有一枚直径为D的硬币,内有n个点,并且满足任意三个点不在一条直线,任意两点之间有连线的条件,现在问你抛掷这枚硬币,硬币内的线段与平行线相交的概率是多少。
数理分析:我们这里首先给出一个概率论的模型。
蒲丰抛针实验:将长度为L的针抛入间距为D(L<D)的平行线中,针与平行线相交的概率P = 2L/πD。
这是概率论中求π的近似值的一个经典几何概型,在以后讨论概率论的文章中我们会详细的介绍,这里不展开这个模型的证明。
在题目描述中给出的“直径为D内”的限定,其实就呼应了蒲丰抛针实验中L < D的限制条件,并且我们需要把n个点构成的图形看做一个整体,显然我们需要构造一个凸包(因为它呈现出这个整体的边界,我们只需讨论边界与平行线相交的概率即可)。
这样,我们就把蒲丰抛针和凸包两个模型结合在了一起,但是在具体求解的过程中,我们应该怎样操作呢?
假设我们已经做出了包含n个顶点的凸包,第i条边的边长为Li,那么单独来看这条边,它与平行线相交的概率Pi = (2 Li) / (πD),考虑到每条直线与改边相交的同时,一定会和该凸包的第j条边相交,因此我们遍历每一条边的时候,应该乘上1/2,即P = 1/2∑Pi = L/πD,这里的L,即是凸包的周长。
有了这层数理逻辑,我们只需通过我们学习过的Graham-scan算法构造出凸包并求出周长,然后带入上面的式子即可。
参考代码如下。
#include<stdio.h> #include<algorithm> #include<cmath> using namespace std; const double pi = acos(-1.0); const int maxn = 110; const double eps = 1e-8; int n , stack[maxn]; struct point { double x , y; } p[maxn]; double dis(point a , point b) { return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y)); } double Cross_left(point a , point b , point c) { return (a.x - c.x)*(b.y - c.y) - (a.y - c.y)*(b.x - c.x); } bool cmp(point a , point b) //极角排序 { double x = Cross_left(a , b , p[0]); if(x > 0) return true; else if(fabs(x) <= eps && (dis(b,p[0]) >= dis(a,p[0]))) return true; else return false; } double Graham() { int i , k = 0 , top = 2; point pmin=p[0] ; for(i = 1;i < n;i++) { if(p[i].y < pmin.y || (p[i].y - pmin.y <= eps && p[i].x < pmin.x)) { pmin = p[i] , k = i; //找极点 } } swap(p[0] , p[k]); sort(p+1,p+n,cmp); for(i = 0;i < 3;i++) stack[i] = i; for(i = 3;i < n;i++) { while(Cross_left(p[stack[top]] , p[i] , p[stack[top-1]]) <= 0)//共线情况不入栈 { top--; if(top == 0) break; } stack[++top] = i; } double cir = 0;top++; for(i = 0;i < top;i++) cir += dis(p[stack[i]] , p[stack[(i+1)%top]]); return cir; } int main() { int ncases; int t = 1; double d; scanf("%d",&ncases); while(ncases--) { scanf("%d%lf",&n,&d); for(int i = 0;i < n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); printf("Case #%d: ",t++); if(n == 1) printf("0.0000\n"); else if(n == 2) printf("%.4lf\n",2*dis(p[0],p[1])/(pi*d)); else printf("%.4lf\n",Graham()/(pi*d)); } }
基于我们对凸包问题的探讨,我们现在思考这样一个军事问题,在领土内,如何设置一个放哨点,使其可以看到领土的每一个角落?
抽象化得看这个模型的话,就是求一个平面几何图形的内核。而内核是什么呢?即该点与多边形上的任意一个点的连线都在该几何图形的内部,想想看,内核是否就是上面那个军事问题中符合要求的放哨点?
而如何用计算机设计算法来找到这个内核呢,或者说内核区域呢?这里我们就要用到半平面求交。
关于半平面求交的定义很简单,即给出一条直线ax+by+c = 0,和一个平面a。我们取ax + by + c > 0 (当然也可以小于零)并且属于a的平面区域a'。
那么我们回到这个求解内核的实际问题当中来,我们得到顺指针排列的多边形的点集V。
基于内核的定义,需要点与边上的任意点的连线都在多边形内部,那么我们考虑依次遍历所有的边。然后以改边所在直线为准,进行半平面求交,这里计算ax+by+c大于0还是小于0是根据第一个边来确定的。经过半平面求交计算求得的平面区域,表示该区域的点与当前遍历的边上任意点的连线是在多边形内部的(凹角区域除外,但遍历到凹角所在边时可排除),然后遍历多边形的边,再连续半平面求交,因此,当我们遍历了所有的边的时候,剩下的平面区域满足与所有边上的任意点的连线都在多边形内部,即是所谓的内核区域。
因此我们可以概括出如下的算法步骤:
step1:得到多边形顶点集V,元素按照顺时针或者逆时针排列。
step2:i = 1,得到边e1 = v1v2,计算e1所在直线l1:a1x + b1y + c1 = 0 ,计算l[1]与多边形半平面的交(并确定直线取大于0还是小于0),得到平面flat[2].
step3:i∈[2,n],遍历i,得到边ei = viv(i+1),计算其所在直线l[i]:ai*x + bi*y + c = 0,计算li与平面flat[i]半平面的交。
最终flat[n]即是内核区域。
我们通过一个具体的题目来实现以下上述的算法过程。(pku 3335)
Description
This year, ACM/ICPC World finals will be held in a hall in form of a simple polygon. The coaches and spectators are seated along the edges of the polygon. We want to place a rotating scoreboard somewhere in the hall such that a spectator sitting anywhere on the boundary of the hall can view the scoreboard (i.e., his line of sight is not blocked by a wall). Note that if the line of sight of a spectator is tangent to the polygon boundary (either in a vertex or in an edge), he can still view the scoreboard. You may view spectator's seats as points along the boundary of the simple polygon, and consider the scoreboard as a point as well. Your program is given the corners of the hall (the vertices of the polygon), and must check if there is a location for the scoreboard (a point inside the polygon) such that the scoreboard can be viewed from any point on the edges of the polygon.
Input
The first number in the input line, T is the number of test cases. Each test case is specified on a single line of input in the form n x1 y1 x2 y2 ... xn yn where n (3 ≤ n ≤ 100) is the number of vertices in the polygon, and the pair of integers xi yi sequence specify the vertices of the polygon sorted in order.
Output
The output contains T lines, each corresponding to an input test case in that order. The output line contains either YES or NO depending on whether the scoreboard can be placed inside the hall conforming to the problem conditions.
参考代码如下。
#include <iostream> #include <cstdio> #include <algorithm> #include <cmath> using namespace std; const double eps = 1e-8; const int maxn = 105; int dq[maxn], top, bot, pn, order[maxn], ln; struct Point { double x, y; } p[maxn]; struct Line { Point a, b; double angle; } l[maxn]; int dblcmp(double k) { if (fabs(k) < eps) return 0; return k > 0 ? 1 : -1; } double multi(Point p0, Point p1, Point p2) { return (p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x); } bool cmp(int u, int v) { int d = dblcmp(l[u].angle-l[v].angle); if (!d) return dblcmp(multi(l[u].a, l[v].a, l[v].b)) < 0; return d < 0; } void getIntersect(Line l1, Line l2, Point& p) { double dot1,dot2; dot1 = multi(l2.a, l1.b, l1.a); dot2 = multi(l1.b, l2.b, l1.a); p.x = (l2.a.x * dot2 + l2.b.x * dot1) / (dot2 + dot1); p.y = (l2.a.y * dot2 + l2.b.y * dot1) / (dot2 + dot1); } bool judge(Line l0, Line l1, Line l2) { Point p; getIntersect(l1, l2, p); return dblcmp(multi(p, l0.a, l0.b)) > 0; } void addLine(double x1, double y1, double x2, double y2) { l[ln].a.x = x1; l[ln].a.y = y1; l[ln].b.x = x2; l[ln].b.y = y2; l[ln].angle = atan2(y2-y1, x2-x1); order[ln] = ln; ln++; } void halfPlaneIntersection() { int i, j; sort(order, order+ln, cmp); for (i = 1, j = 0; i < ln; i++) if (dblcmp(l[order[i]].angle-l[order[j]].angle) > 0) order[++j] = order[i]; ln = j + 1; dq[0] = order[0]; dq[1] = order[1]; bot = 0; top = 1; for (i = 2; i < ln; i++) { while (bot < top && judge(l[order[i]], l[dq[top-1]], l[dq[top]])) top--; while (bot < top && judge(l[order[i]], l[dq[bot+1]], l[dq[bot]])) bot++; dq[++top] = order[i]; } while (bot < top && judge(l[dq[bot]], l[dq[top-1]], l[dq[top]])) top--; while (bot < top && judge(l[dq[top]], l[dq[bot+1]], l[dq[bot]])) bot++; } bool isThereACore() { if (top-bot > 1) return true; return false; } int main() { int t, i; scanf ("%d", &t); while (t--) { scanf ("%d", &pn); for (i = 0; i < pn; i++) scanf ("%lf%lf", &p[i].x, &p[i].y); for (ln = i = 0; i < pn-1; i++) addLine(p[i].x, p[i].y, p[i+1].x, p[i+1].y); addLine(p[i].x, p[i].y, p[0].x, p[0].y); halfPlaneIntersection(); if (isThereACore()) printf ("YES\n"); else printf ("NO\n"); } return 0; }
基于上文对构造凸包的探讨,下面我们来讨论如何判断一个多边形是否是凸包。(Problem source : pku 1584)
Description
Input
Output
题目大意:给出一个圆的半径和圆心,并顺次(顺时针或者逆时针都可)给出n个点。首先判断该n边形是否是凸包,如果是,判断该圆是否被抱在这个凸包当中。
数理分析:通过问题描述,我们不难看出我们需要解决如下两个问题。
①判断一个多边形是否是凸包。
②判断圆是否包含在其内部。
针对第一个问题,我们依旧利用在计算几何中强有力的工具——叉积。我们遍历凸包相邻三点,叉积的值均为非正或非负即表明是一个凸包。这里值得一提的是叉积等于0这种情况,它表示三点共线,如果求解的凸包要求顶点不能共线,只需把上面的结论改成均为正或负即可。
针对第二个问题,我们再分解其判断的流程。
步骤一:要判断这个圆是否在凸包的内部,首先得判断该点是否在凸包的内部,同样可以利用上一段所描述的叉积来实现。
步骤二:然后我们想到,圆在凸包内部的临界情况是内切,因此我们可以通过圆心到各边长的距离和半径的比较来判断这个圆是否在凸包的内部。
这里我们可以看到,如果圆心不在凸包内部,在步骤二的判断过程中依然会出现符合要求的情况,但是圆却不在凸包的内部,这佐证了步骤一存在的必要性。
编程实现:有了上述的数理分析,我们在编程实现的时候需要一些比较巧的编程技巧以简化代码。在判断是否是凸包和判断圆心是否在凸包内部的时候我么你都用到了遍历点求叉积的过程,由于这里没有告诉究竟是顺时针给出点集还是逆时针给出,因此我们设置一个标记数组s(初始化为1),叉乘为0、正、负分别记为0、1、2,得到结果后,记s[x] = 0,(x = 0 、1、 2),我们再利用位运算符“|”,当s[1] | s[2] = 0的时候,表明同时出现了叉乘为正和为负的情况,此时即可判断不符合要求。这就很好解决了点集给出时的方向不确定性的问题。
参考代码如下。
#include<iostream> #include<cstdio> #include<math.h> #define eps 1e-8 #define _sign(x) (x>eps ? 1 : x < -eps ? 2 : 0) struct Point{double x , y;}; struct Line{double a , b , c;}; struct Circle{double r;Point center;}; double xmult(Point p1,Point p2,Point p0) { return (p1.x-p0.x)*(p2.y-p0.y) - (p2.x-p0.x)*(p1.y-p0.y); } int isConvex(int n,Point *p) { int i , s[3] = {1,1,1}; for(i = 0;i < n && s[1]|s[2];i++) s[_sign(xmult(p[(i+1)%n],p[(i+2)%n],p[i]))] = 0; return s[1]|s[2]; } int insideConvex(Point q,int n,Point *p) { int i , s[3] = {1,1,1}; for(i = 0;i < n && s[1]|s[2];i++) s[_sign(xmult(p[(i+1)%n],q,p[i]))] = 0; return s[1]|s[2]; } double poToLine(Point p , Line l) { return fabs(l.a*p.x + l.b*p.y +l.c)/sqrt(l.a*l.a + l.b*l.b); } Line twoPoLine(Point p1, Point p2) { Line l; l.a = p1.y - p2.y; l.b = p2.x - p1.x; l.c = p1.x*p2.y - p2.x*p1.y; return l; } Point p[10000]; int main() { int n; Circle c; while(scanf("%d",&n) != EOF,n > 2) { int i; scanf("%lf%lf%lf",&c.r,&c.center.x,&c.center.y); for(i = 0;i < n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); if(isConvex(n,p)) { if(insideConvex(c.center,n,p)) { for(i = 0;i < n;i++) { Line l = twoPoLine(p[i],p[(i+1)%n]); if(poToLine(c.center,l)-c.r< 0.0) {printf("PEG WILL NOT FIT\n");break;} } if(i == n) {printf("PEG WILL FIT\n");} } else printf("PEG WILL NOT FIT\n"); } else printf("HOLE IS ILL-FORMED\n"); } return 0; }
我们再来看一个有关凸包的问题。(Problem source : pku 2007)
Description
A closed polygon is called convex if the line segment joining any two points of the polygon lies in the polygon. Figure 1 shows a closed polygon which is convex and one which is not convex. (Informally, a closed polygon is convex if its border doesn't have any "dents".)
![](https://images2015.cnblogs.com/blog/894523/201603/894523-20160316170129631-416476419.jpg)
The first property is that the vertices of the polygon will be confined to three or fewer of the four quadrants of the coordinate plane. In the example shown in Figure 2, none of the vertices are in the second quadrant (where x < 0, y > 0).
To describe the second property, suppose you "take a trip" around the polygon: start at (0, 0), visit all other vertices exactly once, and arrive at (0, 0). As you visit each vertex (other than (0, 0)), draw the diagonal that connects the current vertex with (0, 0), and calculate the slope of this diagonal. Then, within each quadrant, the slopes of these diagonals will form a decreasing or increasing sequence of numbers, i.e., they will be sorted. Figure 3 illustrates this point.
![](https://images2015.cnblogs.com/blog/894523/201603/894523-20160316170218599-2017990963.jpg)
![](https://images2015.cnblogs.com/blog/894523/201603/894523-20160316170258365-1624404114.jpg)
Input
Output
题目大意:无序给出凸包的顶点集,其中一个顶点有(0,0),现在要求你按照逆时针输出该凸包的顶点,从(0,0)开始。
数理分析:在前文讨论构造凸包的高效算法Graham-Scan算法的时候,我们蜻蜓点水般地点了一下极角排序这个概念,这里通过这个问题更加详细的介绍一下这个概念的含义。 极角排序是基于极坐标来说的。对于极坐标:选取平面内任意一点o做坐标原点,那么平面内的任意点A都可以用(x,y)来表示,其中x = ρcosθ,y = ρsinθ。这里ρ是A到o的距离,θ则是oA与过o的水平线的正方向的夹角。某两点的θ相同,那么我们选取ρ较小的排在前面。 那么我们根据点集V中各元素极角θ的大小进行排序,然后再次访问该点集,会发现我们访问的点总体上是呈逆时针顺序的分布了。
那么我们再回这个问题,由于这里说明了给出的一定是能够组成凸包的点集,且从(0,0)开始逆时针输出,那么我们就很自然与上面极角排序的模型进行关联,其实本质上就是以笛卡尔坐标系下的原点为极坐标系下的原点,将给出的点集进行极角排序,然后依次输出,便一定是按照逆时针方向输出凸包顶点。 知道了要用极角排序这一点了,接下来的问题就是如何实现极角排序了。方法有很多,而根据极角排序的定义,在这里给出一种利用计算几何中最常用的工具——叉积,刚好能够将其实现。 编程实现:我们通过设置一个bool函数cmp,来给出o、v1、v2的相对位置,然后借用c++中方面的sort排序即可快速完成对整个点集的极角排序。
参考代码如下(值得注意的是从题目描述来看,这道题目从屏幕输出来看是不受参数限制的,因此在编码的时候用EOF控制结束输入即可):
#include<iostream> #include<cmath> #include<stdio.h> #include<algorithm> #define EPS 1e-8 using namespace std; struct point{double x ,y;}; point convex[50]; double cross(point p1, point p2, point q1,point q2) { return (q2.y - q1.y)*(p2.x - p1.x) - (q2.x - q1.x)*(p2.y - p1.y); } bool cmp(point a , point b) { point o; o.x = o.y = 0; return cross(o,b,o,a) < 0; } int main() { int cnt = 0; while(scanf("%lf%lf",&convex[cnt].x,&convex[cnt].y) != EOF) { cnt++; } sort(convex + 1 , convex + cnt , cmp); for(int i = 0;i < cnt;i++) cout<<"("<<convex[i].x<<","<<convex[i].y<<")"<<endl; return 0; }
通过上文对凸包问题的初步了解,这里我们将介绍旋转卡壳的相关概念。
其概念非常简单、很好理解。
首先介绍支撑线,给定一个凸多边形,如果直线L过凸多边形的顶点,并且该凸多边形的其余顶点都分布在该直线L的一侧,则我们称这条直线是这条凸多边形的一条支撑线。
在支撑线的基础上,如果凸多边形存在两条平行的支撑线,那么这两条支撑线分别经过的两个点叫做一对对踵点。
那么,一个凸多边形的所有成对的对踵点当中,存在一组两点相距距离最大的,那么这个长度我们成为凸包的直径。
而这些和旋转卡壳这个名字有什么关系呢?通过其字面意思,联系我们上文中给出的凸包直径的概念,如果我们找到表征凸包直径的对踵点,并画出两个点所在的两条支撑线,此时该凸包其实是可以绕着直径,在两条支撑线当中任意旋转的,这其实就是“旋转卡壳”想表达的意思,明白算法名称的含义,将更加得有利于我们理解算法在干什么。
我们再抽象化的概括一下问题,即给定一个凸包,求其直径。或者也可以表述成,给定点集,求该点集的任意两元素的最大距离。
知道了算法的内涵了,我们就要思考算法是怎么具体实现的了。
从问题描述上来看,似乎可以用穷举来实现,但是随着点集元素的增大,这种暴力穷举的方法的时间复杂度也会激增。
有没有更优化的求解方法呢?
我们逆向的去思考这个问题,假设我们已经给出了凸包的一对对踵点,但两点的连线不一定是凸包的直径,我们根据这两点也是可以做出两条互相平行的支撑线,随后我们让这两条互相平行沿着凸包的边旋转,显然,这样下去会遍历出所有成对的对踵点所在支撑线平行的情况,而又由于其实状态下两条互相平行线之间的距离是确定的,因此当支撑线越靠近凸包的边,此时夹在两条平行线之间的对踵点连线也是越“歪斜”的(可以通过作图实验或者思维想象一下),也就越长,那么我们显然能够看到它的终态:平行线中的一条与凸包的一条边v(i)v(i+1)重合,那么凸包的直径就是max(v(q)v(i) , v(q)v(j))。(v(q)表示另一条支撑线经过的那个顶点)。
这样在按照某个方向旋转互相平行的支撑线的时候,一个点的对踵点也是按照这个方向旋转的,这在维护最大值的时候,相对于穷举算法,便避免了大量的重复计算。
我们通过一个题目来具体的实现一下这个算法。(Problem source : pku 2187)
Total Submissions: 32797 | Accepted: 10177 |
Description
Even though Bessie travels directly in a straight line between pairs of farms, the distance between some farms can be quite large, so she wants to bring a suitcase full of hay with her so she has enough food to eat on each leg of her journey. Since Bessie refills her suitcase at every farm she visits, she wants to determine the maximum possible distance she might need to travel so she knows the size of suitcase she must bring.Help Bessie by computing the maximum distance among all pairs of farms.
Input
* Lines 2..N+1: Two space-separated integers x and y specifying coordinate of each farm
Output
题目大意:给定一个点集,让你其中两点之间最大的距离。
数理分析:基于上面的分析,这里就可以恰到好处的用到凸包问题的旋转卡壳算法了。
编程实现:由于时间原因,笔者在这里暂且只贴出代码而没有给出注释和编程实现的分析,在以后修正文章的时候会详细给出分析的。
#include <cmath> #include<stdio.h> #include <algorithm> #include <iostream> using namespace std; #define MAXN 50005 struct Point { int x, y; bool operator < (const Point& _P) const { return y<_P.y||(y==_P.y&&x<_P.x); }; }pset[MAXN],ch[MAXN]; double cross(Point a,Point b,Point o) { return (a.x - o.x) * (b.y - o.y) - (b.x - o.x) * (a.y - o.y); } void convex_hull(Point *p,Point *ch,int n,int &len) { sort(p, p+n); ch[0]=p[0]; ch[1]=p[1]; int top=1; for(int i=2;i<n;i++) { while(top>0&&cross(ch[top],p[i],ch[top-1])<=0) top--; ch[++top]=p[i]; } int tmp=top; for(int i=n-2;i>=0;i--) { while(top>tmp&&cross(ch[top],p[i],ch[top-1])<=0) top--; ch[++top]=p[i]; } len=top; } int dist2(Point a,Point b) { return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y); } int rotating_calipers(Point *ch,int n) { int q=1,ans=0; ch[n]=ch[0]; for(int p=0;p<n;p++) { while(cross(ch[p+1],ch[q+1],ch[p])>cross(ch[p+1],ch[q],ch[p])) q=(q+1)%n; ans=max(ans,max(dist2(ch[p],ch[q]),dist2(ch[p+1],ch[q+1]))); } return ans; } int main() { int n, len; while(scanf("%d", &n)!=EOF) { for(int i = 0;i < n;i++) { scanf("%d %d",&pset[i].x,&pset[i].y); } convex_hull(pset,ch,n,len); printf("%d\n",rotating_calipers(ch,len)); } return 0; }