计算几何模板集

 // 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)
", 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)
", 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==2居然要特判。。。什么傻逼
        else if(n==2) printf("%.2lf
", (node[0]-node[1]).len());
        else
            printf("%.2lf
",  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)
", 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
", 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)
", 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
", res);
    }
    return 0;
}

---

  （未完）