(二）计算几何基础之凸包与半平面交

凸包

形象解释：有很多头奶牛在田野上，现在希望你用一条绳子将这些奶牛给围起来，这样所形成的多边形称为凸包。

Andrew算法求凸包步骤：

1. 将给出的所有点排序。（横坐标为第一关键字，纵坐标为第二关键字）
2. 将整个凸包分成两个部分来求（上半部分和下半部分），先从左至右做上半部分，再从右至左做下半部分。（用栈来维护）
   上凸包看左侧点，下凸包看右侧点。
3. 最后用栈底更新一遍。（边界点）
   
   时间复杂度：O(nlogn)

code:

#include<bits/stdc++.h>
using namespace std;
#define x first
#define y second


const int N = 10010;
int n;
typedef pair<double, double> PDD;


PDD q[N];
int stk[N];
bool used[N];

PDD operator-(PDD a,PDD b)
{
    return {a.x-b.x,a.y-b.y};
}
double cross(PDD a,PDD b)
{
    return a.x*b.y-a.y*b.x;
}
double area(PDD a,PDD b,PDD c)
{
    return cross(b-a,c-a);
}
double get_dist(PDD a,PDD b)
{
    return sqrt(pow((b.x-a.x),2)+pow(b.y-a.y,2));
}
double andrew()
{
    sort(q,q+n);
    int top = 0;
    for(int i = 0;i<n;i++)
    {//从左往右求下凸包（第一遍）
        while(top>=2&&area(q[stk[top-1]],q[stk[top]],q[i])<=0)//这里如果共线算凸包就是<=，如果不算就是<，详情可见洛谷和acwing的算法进阶课
        {
            //说明当前判断到的点在栈顶元素所在边的逆时针方向
            //所以当前栈顶元素所在边应该删掉,当前栈顶元素也应该被删掉
            if (area(q[stk[top - 1]], q[stk[top]], q[i]) < 0)
                used[stk[top]] = false;
            top--;//栈顶指针需要--
        }
        stk[++top] = i;//把新的元素作为栈顶元素
        used[i] = true;//标记当前点
    }
    used[0] = false;
    for(int i = n-1;i>=0;i--)
    {//从右往左
        if(used[i])
            continue;
        while(top>=2&&area(q[stk[top-1]],q[stk[top]],q[i])<=0)
        {
            top--;
        }
        stk[++top] = i;
    }
    double res = 0;
    for(int i = 2;i<=top;i++)
    {
        res+=get_dist(q[stk[i-1]],q[stk[i]]);
    }
    return res;
}
int main()
{
    cin>>n;
    for(int i = 0;i<n;i++)
    {
        cin>>q[i].x>>q[i].y;
    }
    double res = andrew();
    printf("%.2lf\n",res);
    return 0;
}

例题一：信用卡凸包

原题链接

分析：其实就是求圆心围成的凸包再加上一个圆的周长。

#include<bits/stdc++.h>
using namespace std;

#define x first
#define y second

typedef pair<double, double> PDD;
const int N = 40010;
int n,cnt;
PDD q[N];
int stk[N],top;
bool used[N];
const double PI = acos(-1);

PDD rotate(PDD a,double b)
{
    return {a.x*cos(b)+a.y*sin(b),-a.x*sin(b)+a.y*cos(b)};
}
PDD operator-(PDD a,PDD b)
{
    return {a.x-b.x,a.y-b.y};
}
double cross(PDD a,PDD b)
{
    return a.x*b.y-a.y*b.x;
}
double area(PDD a,PDD b,PDD c)
{
    return cross(b-a,c-a);
}
double get_dist(PDD a,PDD b)
{
    double dx = b.x-a.x;
    double dy = b.y-a.y;
    return sqrt(dx*dx+dy*dy);
}
double andrew()
{
    sort(q, q + cnt);
    for (int i = 0; i < cnt; i ++ )
    {
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
        {
            // 凸包边界上的点即使被从栈中删掉，也不能删掉used上的标记
            if (area(q[stk[top - 1]], q[stk[top]], q[i]) < 0)
                used[stk[top -- ]] = false;
            else top -- ;
        }
        stk[ ++ top] = i;
        used[i] = true;
    }
    used[0] = 0;
    for (int i = cnt - 1; i >= 0; i -- )
    {
        if (used[i]) continue;
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
            top -- ;
        stk[ ++ top] = i;
    }

    double res = 0;
    for (int i = 2; i <= top; i ++ )
        res += get_dist(q[stk[i - 1]], q[stk[i]]);
    return res;
}
int main()
{
    cin>>n;
    double a,b,r;
    cin>>a>>b>>r;
    a = a/2-r,b = b/2-r;
    int dx[] = {1,1,-1,-1},dy[] = {1,-1,-1,1};
    while (n -- )
    {
        double x,y,z;
        cin>>x>>y>>z;
        for(int i = 0;i<4;i++)
        {
            auto t = rotate({dx[i]*b,dy[i]*a},-z);//
            q[cnt++]= {x+t.x,y+t.y};
        }
    }
    double res = andrew();
    printf("%.2lf\n",res+2*PI*r);
    return 0;
}

code:

半平面交

形象解释：
给定许多直线，求出这些直线左半边部分围成的一个凸多边形。
步骤：

1. 将向量按角度排序 (atan2(y,x)：计算角度,[-180,180])
2. 按顺序从左往右扫描所有的向量 。（用双端队列维护半平面交的轮廓） （如何维护？）

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>

#define x first
#define y second

using namespace std;

typedef pair<double,double> PDD;
const int N = 510;
const double eps = 1e-8;

int cnt;
struct Line
{
    PDD st, ed;
}line[N];
PDD pg[N], ans[N];
int q[N];

int sign(double x)
{
    if (fabs(x) < eps) return 0;
    if (x < 0) return -1;
    return 1;
}

int dcmp(double x, double y)
{
    if (fabs(x - y) < eps) return 0;
    if (x < y) return -1;
    return 1;
}

double get_angle(const Line& a)
{
    return atan2(a.ed.y - a.st.y, a.ed.x - a.st.x);
}

PDD operator-(PDD a, PDD b)
{
    return {a.x - b.x, a.y - b.y};
}

double cross(PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}

double area(PDD a, PDD b, PDD c)
{
    return cross(b - a, c - a);
}

bool cmp(const Line& a, const Line& b)
{
    double A = get_angle(a), B = get_angle(b);
    if (!dcmp(A, B)) return area(a.st, a.ed, b.ed) < 0;
    return A < B;
}

PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w)
{
    auto u = p - q;
    double t = cross(w, u) / cross(v, w);
    return {p.x + v.x * t, p.y + v.y * t};
}

PDD get_line_intersection(Line a, Line b)
{
    return get_line_intersection(a.st, a.ed - a.st, b.st, b.ed - b.st);
}

// bc的交点是否在a的右侧
bool on_right(Line& a, Line& b, Line& c)
{
    auto o = get_line_intersection(b, c);
    return sign(area(a.st, a.ed, o)) <= 0;
}

double half_plane_intersection()
{
    sort(line, line + cnt, cmp);
    int hh = 0, tt = -1;
    for (int i = 0; i < cnt; i ++ )
    {
        if (i && !dcmp(get_angle(line[i]), get_angle(line[i - 1]))) continue;
        while (hh + 1 <= tt && on_right(line[i], line[q[tt - 1]], line[q[tt]])) tt -- ;
        while (hh + 1 <= tt && on_right(line[i], line[q[hh]], line[q[hh + 1]])) hh ++ ;
        q[ ++ tt] = i;
    }
    while (hh + 1 <= tt && on_right(line[q[hh]], line[q[tt - 1]], line[q[tt]])) tt -- ;
    while (hh + 1 <= tt && on_right(line[q[tt]], line[q[hh]], line[q[hh + 1]])) hh ++ ;

    q[ ++ tt] = q[hh];
    int k = 0;
    for (int i = hh; i < tt; i ++ )
        ans[k ++ ] = get_line_intersection(line[q[i]], line[q[i + 1]]);
    double res = 0;
    for (int i = 1; i + 1 < k; i ++ )
        res += area(ans[0], ans[i], ans[i + 1]);
    return res / 2;
}

int main()
{
    int n, m;
    scanf("%d", &n);
    while (n -- )
    {
        scanf("%d", &m);
        for (int i = 0; i < m; i ++ ) scanf("%lf%lf", &pg[i].x, &pg[i].y);
        for (int i = 0; i < m; i ++ )
            line[cnt ++ ] = {pg[i], pg[(i + 1) % m]};
    }
    double res = half_plane_intersection();
    printf("%.3lf\n", res);

    return 0;
}