[POJ][2451][Uyuw's Concert][半平面交]

朱泽园专为他那篇nlogn算法解决半平面交问题的论文而出的题目。

 

/*
题目：Uyuw's Concert
题目来源：POJ 2451
题目难度：中等
题目内容或思路：
    半平面交的模板题
    不加读入优化300+ms，加了读入优化后94ms，刷到前几了
做题日期：2011.2.24
*/
#include <cstdio>
#include <cstdlib>
#include <climits>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <string>
#include <queue>
#include <map>
#include <vector>
#include <bitset>
#include <cmath>
#include <set>
using namespace std;

const double eps = 1e-8;
const double pi = acos(-1.0);
const int N = 20010;
const double maxl = 10000;

struct cpoint {
    double x, y;
    cpoint(double xx = 0, double yy = 0): x(xx), y(yy) {};
};

int dcmp(double x) {
    if (x < -eps) return -1; else return x > eps;
}

double cross(cpoint p0, cpoint p1, cpoint p2) { // p0p1 与 p0p2 叉积
    return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
}

bool EqualPoint(cpoint a, cpoint b) {
    return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0;
}

//==============半平面交==================//

struct cvector {
    cpoint s, e;
    double ang, d;
};

void setline(double x1, double y1, double x2, double y2, cvector &v) {
    v.s.x = x1; v.s.y = y1;
    v.e.x = x2; v.e.y = y2;
    v.ang = atan2(y2 - y1, x2 - x1);
    if (dcmp(x1 - x2))    v.d = (x1 * y2 - x2 * y1) / fabs(x1 - x2);
    else                v.d = (x1 * y2 - x2 * y1) / fabs(y1 - y2);
}

//判向量平行
bool parallel(const cvector &a, const cvector &b) {
    double u = (a.e.x - a.s.x) * (b.e.y - b.s.y) 
            - (a.e.y - a.s.y) * (b.e.x - b.s.x);
    return dcmp(u) == 0;
}

//求两向量(直线)交点 (两向量不能平行或重合)
cpoint CrossPoint(const cvector &a, const cvector &b) {
    cpoint res;
    double u = cross(a.s, a.e, b.s), v = cross(a.e, a.s, b.e);
    res.x = (b.s.x * v + b.e.x * u) / (u + v);
    res.y = (b.s.y * v + b.e.y * u) / (u + v);
    return res;
}

//半平面交排序函数[优先顺序: 1.极角 2.前面的直线在后面的左边]
static bool VecCmp(const cvector &l, const cvector &r) {
    if (dcmp(l.ang - r.ang)) return l.ang < r.ang;
    return l.d < r.d;
}

cvector deq[N]; //用于计算的双端队列

//        获取半平面交的多边形（多边形的核）
//        注意:1.半平面在向量左边, 2.函数会改变vec[]中的值
//函数运行后如果n[即返回多边形的点数量]为0则
//不存在半平面交的多边形（不存在区域或区域面积无穷大）
void HalfPanelCross(cvector vec[], int n, cpoint cp[], int &m) {
    int i, tn; m = 0;
    sort(vec, vec + n, VecCmp);
    for(i = tn = 1; i < n; ++i) //平面在向量左边的筛选
        if(dcmp(vec[i].ang - vec[i - 1].ang) != 0)
            vec[tn++] = vec[i];
    n = tn;
    int bot = 0, top = 1;
    deq[0] = vec[0];
    deq[1] = vec[1]; // vec[]大小不可小于2
    for (i = 2; i < n; ++i) {
        if (parallel(deq[top], deq[top - 1]) ||
            parallel(deq[bot], deq[bot + 1]) ) return ;
        while ( bot < top && dcmp( cross(vec[i].s, vec[i].e,
            CrossPoint(deq[top], deq[top - 1])) ) < 0 ) top--;
        while ( bot < top && dcmp( cross(vec[i].s, vec[i].e, 
            CrossPoint(deq[bot], deq[bot + 1])) ) < 0 ) bot++;
        deq[++top] = vec[i];
    }
    while ( bot < top && dcmp( cross(deq[bot].s, deq[bot].e, 
        CrossPoint(deq[top], deq[top - 1])) ) < 0 ) top--;
    while ( bot < top && dcmp( cross(deq[top].s, deq[top].e, 
        CrossPoint(deq[bot], deq[bot + 1])) ) < 0 ) bot++;
    if (top <= bot + 1) return ; // 两条或两条以下的直线，面积无穷大
    for (i = bot; i < top; i ++)
        cp[m++] = CrossPoint(deq[i], deq[i + 1]);
    if (bot < top + 1)
        cp[m++] = CrossPoint(deq[bot], deq[top]);
    m = unique(cp, cp + m, EqualPoint) - cp;
    for (i = 0; i < m; ++i) {
        if (dcmp(cp[i].x) == 0) cp[i].x = 0;
        if (dcmp(cp[i].y) == 0) cp[i].y = 0;
    }
}

double PolygonArea(cpoint p[], int n) {
    if (n < 3) return 0;
    double s = p[0].y * (p[n - 1].x - p[1].x);
    for (int i = 1; i < n; ++i)
        s += p[i].y * (p[i - 1].x - p[(i + 1) % n].x);
    return fabs(s / 2); // 顺时针方向s为负
}

int n, m;
cvector v[N];
cpoint cp[N];

double readdouble() {
    char c;
    while (c = getchar(), c != '-' && c != '.' && !isdigit(c));
    int f = 1;
    if (c == '-') {
        f = -1;
        c = getchar();
    }
    double p = 1, res = 0;
    if (c != '.') {
        res = c - '0';
        while (isdigit(c = getchar())) {
            res = res * 10 + c - '0';
        }
    }
    if (c != '.') return res * f;
    while (isdigit(c = getchar())) {
        res = res * 10 + c - '0';
        p *= 10;
    }
    return res * f / p;
}

void solve() {
    setline(0, 0, maxl, 0, v[0]);
    setline(maxl, 0, maxl, maxl, v[1]);
    setline(maxl, maxl, 0, maxl, v[2]);
    setline(0, maxl, 0, 0, v[3]);
    n += 4;
    double x1, x2, y1, y2;
    for (int i = 4; i < n; ++i) {
        x1 = readdouble();
        y1 = readdouble();
        x2 = readdouble();
        y2 = readdouble();
        setline(x1, y1, x2, y2, v[i]);
    }
    HalfPanelCross(v, n, cp, m);
    if (m < 3)
        printf("0.0\n");
    else
        printf("%.1f\n", PolygonArea(cp, m));
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("D:\\in.txt", "r", stdin);
#endif
    while (scanf("%d", &n) != EOF) {
        solve();
    }
    return 0;
}