计算几何模板集
// 2019.8.13
// 昨天的内容是计算几何部分点、线、多边形、向量、凸包等基本概念的介绍和模板的引入，之前的博客已有涉及，不再重复记录。
// 今天同样是ddy带来了计算几何进阶知识点的讲解。
半平面交
见我的博客 F - Feng Shui一题。
凸包 && 旋转卡壳
#include<iostream> #include<cstdio> #include<cmath> #include<algorithm> using namespace std; const int maxn = 5000; typedef double db; struct P { db x, y; P(db x=0, db y=0):x(x), y(y) {} P operator-(const P& a) { return P(x-a.x, y-a.y); } P operator+(const P& a) { return P(x+a.x, y+a.y); } P operator*(const db a) { return P(a*x, a*y); } db len() { return sqrt(x*x+y*y); } void print() { printf("(%.2lf %.2lf)\n", x, y); } }node[maxn]; int n; // 向量a,b叉积 axb db cross(const P& a, const P& b) { return a.x*b.y - b.x*a.y; } // 按极角排序比较函数 bool cmp(P a, P b) { db x = cross(a-node[0], b-node[0]); if(x>0) return 1; if(x==0) return (a-node[0]).len()<(b-node[0]).len(); return 0; } // 凸包 P stack[maxn]; int top; // 0~cnt-1 // 返回凸包上点的个数 // 点按逆时针放入stack int graham() { // 1. 找到最左下方的点作为起始点 int k=0; for(int i=1;i<n;i++) { if(node[i].y<node[k].y || node[i].y==node[k].y && node[i].x<node[k].x) k = i; } swap(node[0], node[k]); // 2. 将全部点按照极角排序 sort(node+1, node+n, cmp); // 3. 保存凸包上的点 stack[0] = node[0]; stack[1] = node[1]; int top = 1; for(int i=2;i<n;i++) { while(top>=1 && cross(stack[top]-stack[top-1], node[i]-stack[top-1])<=0) --top; stack[++top] = node[i]; } for(int i=0;i<=top;i++) { stack[i].print(); } return top+1; } // 旋转卡壳 // 返回凸包的直径 // c0: 圆心位置 P c0; db rotateCalipers() { db ans = 0; int cnt = graham(); stack[cnt] = stack[0]; int q = 1; for(int p=0;p<cnt;p++) { while(cross(stack[p+1]-stack[p], stack[q]-stack[p])<cross(stack[p+1]-stack[p], stack[q+1]-stack[p])) q = (q+1)%cnt; db dis1 = (stack[p]-stack[q]).len(); db dis2 = (stack[p+1]-stack[q+1]).len(); // 求直径的圆心 /* if(dis1<dis2) { if(ans<dis2) { c0 = (stack[p+1]+stack[q+1])*0.5; ans = dis2; } } else { if(ans<dis1) { c0 = (stack[p]+stack[q])*0.5; ans = dis1; } } */ ans = max(ans, max(dis1, dis2)); } return ans; }
时间复杂度
对于Graham算法求凸包：第一步O(n)，第二步排序O(nlogn)，第三步O(n)。
总体时间复杂度O(nlogn)。
旋转卡壳时间复杂度O(n)。
应用
求平面上点集的直径，或最远点对，即点之间的最大距离。
例：HDU 1392 求凸包周长
AC代码
#include<iostream> #include<cstdio> #include<cmath> #include<algorithm> using namespace std; const int maxn = 110; typedef double db; struct P { db x, y; P(db x=0, db y=0):x(x), y(y) {} P operator-(const P& a) { return P(x-a.x, y-a.y); } P operator+(const P& a) { return P(x+a.x, y+a.y); } P operator*(const db a) { return P(a*x, a*y); } double len() { return sqrt(x*x+y*y); } void print() { printf("(%.2lf %.2lf)\n", x, y); } }node[maxn]; int n; // 向量a,b叉积 axb db cross(const P& a, const P& b) { return a.x*b.y - b.x*a.y; } // 凸包 P stack[maxn]; int top; // 0~cnt-1 bool cmp(P a, P b) { db x = cross(a-node[0], b-node[0]); if(x>0) return 1; if(x==0) return (a-node[0]).len()<(b-node[0]).len(); return 0; } // 返回凸包上点的个数 // 点按逆时针放入stack double graham() { // 1. 找到最左下方的点 int k=0; for(int i=1;i<n;i++) { if(node[i].y<node[k].y || node[i].y==node[k].y && node[i].x<node[k].x) k = i; } swap(node[0], node[k]); // 2. 按照极角排序 sort(node+1, node+n, cmp); /* for(int i=0;i<n;i++) { node[i].print(); } */ // 3. 保存凸包上的点 stack[0] = node[0]; stack[1] = node[1]; int top = 1; for(int i=2;i<n;i++) { while(top>=1 && cross(stack[top]-stack[top-1], node[i]-stack[top-1])<=0) --top; stack[++top] = node[i]; } stack[++top] = stack[0]; double ans = 0; for(int i=1;i<=top;i++) { ans += (stack[i]-stack[i-1]).len(); } return ans; } int main() { while(scanf("%d", &n)!=EOF && n) { for(int i=0;i<n;i++) { scanf("%lf %lf", &node[i].x, &node[i].y); } if(n==1) printf("0.00\n"); // n==2居然要特判。。。什么傻逼 else if(n==2) printf("%.2lf\n", (node[0]-node[1]).len()); else printf("%.2lf\n", graham()); } return 0; }
最小圆覆盖
题目背景：HDU 3007 Buried memory
题意：求最小半径的圆覆盖平面上全部点。
题解：有一种神奇的复杂度为O(n)的随机增量法解决这个问题。
算法流程
- 先对点集random_shuffle处理
- 假设已经有了前i-1个点的最小圆覆盖Ci−1，检验点pi是否在Ci−1内。如果在，则Ci=Ci−1；否则，pi一定在圆Ci上，进行下一步
- 包括点pi在内的前j-1个点(j<i)的最小圆覆盖为Cj−1，检验点pj是否在Cj−1内。如果在，则Cj=Cj−1；否则，pj一定在圆Cj上，进行下一步
- 包括点pi和点pj在内的前k-1个点(k<j)的最小圆覆盖为Ck−1，检验点pk是否在Ck−1内。如果在，则Ck=Ck−1，否则，pk一定在圆Ck上。返回
AC代码（模板）：
#include<iostream> #include<cstdio> #include<cmath> #include<algorithm> using namespace std; const int maxn = 510; const double eps = 1e-8; typedef double db; struct P { db x, y; P(db x=0, db y=0):x(x), y(y) {} P operator-(const P& a) { return P(x-a.x, y-a.y); } P operator+(const P& a) { return P(x+a.x, y+a.y); } P operator*(const db a) { return P(x*a, y*a); } P operator/(const db a) { return P(x/a, y/a); } double len() { return sqrt(x*x+y*y); } void print() { printf("(%.2lf %.2lf)\n", x, y); } }node[maxn]; struct Circle { P center; double r; }; int n; // 向量a,b叉积 axb db cross(const P& a, const P& b) { return a.x*b.y - b.x*a.y; } // 求三角形外接圆 Circle CircleOfTriangle(P p1, P p2, P p3) { double Bx = p2.x - p1.x, By = p2.y - p1.y; double Cx = p3.x - p1.x, Cy = p3.y - p1.y; double D = 2*(Bx*Cy - By*Cx); double cx = (Cy*(Bx*Bx + By*By) - By*(Cx*Cx + Cy*Cy))/D + p1.x; double cy = (Bx*(Cx*Cx + Cy*Cy) - Cx*(Bx*Bx + By*By))/D + p1.y; Circle c; c.center = P(cx, cy); c.r = (p1-c.center).len(); return c; } // 最小圆覆盖 void min_cover_circle() { random_shuffle(node, node+n); Circle c; c.center = node[0]; c.r = 0; for(int i=1;i<n;i++) { if((node[i]-c.center).len()>c.r+eps) { c.center = node[i]; c.r = 0; for(int j=0;j<i;j++) { if((node[j]-c.center).len()>c.r+eps) { c.center = (node[i]+node[j])/2; c.r = (node[i]-node[j]).len()/2; for(int k=0;k<j;k++) { if((node[k]-c.center).len()>c.r+eps) { c = CircleOfTriangle(node[i], node[j], node[k]); } } } } } } printf("%.2lf %.2lf %.2lf\n", c.center.x, c.center.y, c.r); } int main() { while(scanf("%d", &n)!=EOF && n) { for(int i=0;i<n;i++) { scanf("%lf %lf", &node[i].x, &node[i].y); } min_cover_circle(); } return 0; }
三分套三分
题目背景：HDU 3400 Line belt
题意：平面上有两条线段，AB，CD，在线段上移动速度为v1，v2，在平面上移动速度为v3，求 A 移动到 D 的最短时间。
题解：设P，Q分别在线段AB，CD上，移动路径一定为AP-PQ-QD。本题应该满足凸函数性质（why？），于是每段上三分求最优。
注意点A，B，C，D等为全局变量，main函数里用P A(ax, ay)赋值（产生了局部变量A）让我debug半天。。。
AC代码：
#include<iostream> #include<cstdio> #include<cmath> #include<algorithm> using namespace std; const double eps = 1e-6; typedef double db; struct P { db x, y; P(db x=0, db y=0):x(x), y(y) {} P operator-(const P& a) { return P(x-a.x, y-a.y); } P operator+(const P& a) { return P(x+a.x, y+a.y); } P operator*(const db a) { return P(a*x, a*y); } double len() { return sqrt(x*x+y*y); } void print() { printf("(%.2lf %.2lf)\n", x, y); } }A, B, C, D, AB, CD; int ax, ay, bx, by, cx, cy, dx, dy; int v1, v2, v3; // 计算通过AP-PQ-QD需要的时间 double cal(double t1, double t2) { double res = 0; res += AB.len()*t1/v1; res += CD.len()*(1-t2)/v2; P PQ = (CD*t2 + C) - (AB*t1 + A); res += PQ.len()/v3; return res; } // AB上取点A+t*AB，三分CD求最优 double trisearch(double t) { double l = 0, r = 1; double res = 0; while(r-l>eps) { double mid1 = l + (r-l)/3; double mid2 = l + (r-l)/3*2; double res1 = cal(t, mid1); double res2 = cal(t, mid2); if(res1<res2) { r = mid2; res = res1; } else { l = mid1; res = res2; } } return res; } int main() { int t; cin>>t; while(t--) { scanf("%d %d %d %d", &ax, &ay, &bx, &by); scanf("%d %d %d %d", &cx, &cy, &dx, &dy); scanf("%d %d %d", &v1, &v2, &v3); A = P(ax, ay); B = P(bx, by); C = P(cx, cy); D = P(dx, dy); AB = B-A; CD = D-C; // 三分AB求最优解 double l = 0, r = 1; double res = 0; while(r-l>eps) { double mid1 = l + (r-l)/3; double mid2 = l + (r-l)/3*2; double res1 = trisearch(mid1); double res2 = trisearch(mid2); if(res1<res2) { r = mid2; res = res1; } else { l = mid1; res = res2; } } printf("%.2lf\n", res); } return 0; }
（未完）
