bzoj2618: [Cqoi2006]凸多边形 计算几何：半平面交

bzoj2618: [Cqoi2006]凸多边形

Description

逆时针给出n个凸多边形的顶点坐标，求它们交的面积。例如n=2时，两个凸多边形如下图：

则相交部分的面积为5.233。

Input

第一行有一个整数n，表示凸多边形的个数，以下依次描述各个多边形。第i个多边形的第一行包含一个整数mi，表示多边形的边数，以下mi行每行两个整数，逆时针给出各个顶点的坐标。

Output

输出文件仅包含一个实数，表示相交部分的面积，保留三位小数。

Sample Input

2
6
-2 0
-1 -2
1 -2
2 0
1 2
-1 2
4
0 -3
1 -1
2 2
-1 0

Sample Output

5.233

HINT

100%的数据满足：2<=n<=10，3<=mi<=50，每维坐标为[-1000,1000]内的整数

题解

裸半平面交
用的是zzy大神的增量法
每次扔一条直线交点判叉积
好吧我承认这波是我水的~

代码

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<cstring>
using namespace std;
const double eps = 1e-8;
const int N = 1100;
int read() {
    char ch = getchar(); int x = 0, f = 1;
    while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9') {x = x * 10 + ch - '0'; ch = getchar();}
    return x * f;
}
double sgn(double x) {return x > eps ? 1 : (x < -eps ? -1 : 0);}
struct P {
    double x, y;
    P(double a = 0, double b = 0) : x(a), y(b) {}
    P operator - (P a) const {return P(x - a.x, y - a.y);}
    P operator + (P a) const {return P(x + a.x, y + a.y);}
    P operator * (double a) {return P(x * a, y * a);}
    P operator / (double a) {return P(x / a, y / a);}
    double operator * (P a) const {return x * a.y - y * a.x;}
}ret[N], p[N];
double xmult(P a, P b, P c) {return (b - a) * (c - a);}
double angle(P a) {return atan2(a.y, a.x);}

struct line {
    P a, b; double ang;
    line(P c = P(0, 0), P d = P(0, 0), double e = 0) : a(c), b(d), ang(e) {}
    P operator + (line c) {
        double w1 = xmult(c.a, a, b), w2 = xmult(c.b, b, a);
        return (c.a * w2 + c.b * w1) / (w1 + w2);
    }
}l[N];
int top, id[N], q[N], cnt;
double ans;
bool cmp(int a, int b) {
    int t = sgn(l[a].ang - l[b].ang);
    if(!t) return sgn(xmult(l[a].a, l[b].a, l[b].b)) > 0;
    return t < 0;
}

void HPI() {
    int L = 1, R = 2, tot = 1;
    for(int i = 1;i <= top; ++i) id[i] = i;
    sort(id + 1, id + top + 1, cmp);
    for(int i = 2;i <= top; ++i)
        if(sgn(l[id[i]].ang - l[id[i - 1]].ang) > 0)
            id[++tot] = id[i];
    q[1] = id[1]; q[2] = id[2];
    for(int i = 3;i <= tot; ++i) {
        while(L < R && sgn(xmult(l[q[R - 1]] + l[q[R]], l[id[i]].a, l[id[i]].b)) < 0) --R;
        while(L < R && sgn(xmult(l[q[L + 1]] + l[q[L]], l[id[i]].a, l[id[i]].b)) < 0) ++L;
        q[++R] = id[i];
    }
    while(L < R && sgn(xmult(l[q[R - 1]] + l[q[R]], l[q[L]].a, l[q[L]].b)) < 0) --R;
    while(L < R && sgn(xmult(l[q[L + 1]] + l[q[L]], l[q[R]].a, l[q[R]].b)) < 0) ++L;
    q[R + 1] = q[L];
    for(int i = L;i <= R; ++i) ret[++cnt] = l[q[i]] + l[q[i + 1]];
}

int main()
{
    int n = read();
    for(int i = 1;i <= n; ++i) {
        int x = read();
        for(int j = 1; j <= x; ++j) p[j].x = read(), p[j].y = read(); 
        p[x + 1] = p[1];
        for(int j = 1;j <= x; ++j) l[++top] = line(p[j], p[j + 1], angle(p[j] - p[j + 1]));
    }
    HPI();
    if(cnt >= 3) {
        ret[cnt+1] = ret[1];
        for(int i = 1; i <= cnt; ++i) ans += ret[i] * ret[i + 1];
    }
    ans = fabs(ans) / 2;
    printf("%.3lf\n", ans);
    return 0;
}