POJ 2451 Uyuw‘s Concert（半平面交）

题目：

（0，0）-（10000，10000）的平面，给出n条直线切，求最后的面积。

思路：

半平面交裸题，注意给的点的顺序就好。

#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
using namespace std;
const double eps = 1e-6;
struct Point
{
    double x, y;
    Point(double x = 0, double y = 0):x(x),y(y) {}
};
typedef Point Vector;
Vector operator + (Vector A, Vector B){return Vector(A.x+B.x, A.y+B.y);}
Vector operator - (Point A, Point B){return Vector(A.x-B.x, A.y-B.y);}
Vector operator * (Vector A, double p){return Vector(A.x*p, A.y*p);}
Vector operator /(Vector A, double p){return Vector(A.x/p,A.y/p);}
int sgn(double x)
{
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
double Dot(Vector A, Vector B){return A.x*B.x + A.y*B.y;}
double Cross(Vector A, Vector B){return A.x*B.y-A.y*B.x;}
double Length(Vector A){return sqrt(Dot(A, A));}
Vector Normal(Vector A) //向量A左转90°的单位法向量
{
    double L = Length(A);
    return Vector(-A.y/L, A.x/L);
}
struct Line
{
    Point p;//直线上任意一点
    Vector v;//方向向量，它的左边就是对应的半平面
    double ang;//极角，即从x轴正半轴旋转到向量v所需要的角（弧度）
    Line() {}
    Line(Point p, Vector v) : p(p), v(v)
    {
        ang = atan2(v.y, v.x);
    }
    bool operator < (const Line& L) const  //排序用的比较运算符
    {
        return ang < L.ang;
    }
} line[1000];
//点p在有向直线L的左侧
bool OnLeft(Line L, Point p)
{
    return Cross(L.v, p - L.p) >= 0;//等于零的时候在线上
}
//两直线交点。假定交点唯一存在
Point GetIntersection(Line a, Line b)
{
    Vector u = a.p - b.p;
    double t = Cross(b.v, u)/Cross(a.v, b.v);
    return a.p + a.v*t;
}
//半平面交的主过程
int HalfplaneIntersection(Line* L, int n, Point* poly)
{
    sort(L, L + n);//按照极角排序
    int fst = 0, lst = 0;//双端队列的第一个元素和最后一个元素
    Point *P = new Point[n];//p[i] 为 q[i]与q[i + 1]的交点
    Line *q = new Line[n];//双端队列
    q[fst = lst = 0] = L[0];//初始化为只有一个半平面L[0]
    for(int i = 1; i < n; ++i)
    {
        while(fst < lst && !OnLeft(L[i], P[lst - 1])) --lst;
        while(fst < lst && !OnLeft(L[i], P[fst])) ++fst;
        q[++lst] = L[i];
        if(sgn(Cross(q[lst].v, q[lst - 1].v)) == 0)
        {
            //两向量平行且同向，取内侧一个
            --lst;
            if(OnLeft(q[lst], L[i].p)) q[lst] = L[i];
        }
        if(fst < lst)
            P[lst - 1] = GetIntersection(q[lst - 1], q[lst]);
    }
    while(fst < lst && !OnLeft(q[fst], P[lst - 1])) --lst;
    //删除无用平面
    if(lst - fst <= 1) return 0;//空集
    P[lst] = GetIntersection(q[lst], q[fst]);//计算首尾两个半平面的交点
    //从deque复制到输出中
    int m = 0;
    for(int i = fst; i <= lst; ++i) poly[m++] = P[i];
    return m;
}
Point poly[20500];
void init()
{
    Point a[4];
    a[0] = Point(0, 0);
    a[1] = Point(0,10000);
    a[2] = Point(10000, 10000);
    a[3] = Point(10000,0);
    for(int i =0; i<4; i++)
    {
        line[i] = Line(a[i],a[i]- a[(i+1)%4]);
    }
}

int main()
{
    Point a,b;
    int cnt = 4;
    char s[20];

    init();
    int n;
    scanf("%d", &n);
    for(int i = 0;i<n;i++)
    {
        scanf("%lf%lf%lf%lf",&a.x, &a.y, &b.x,&b.y);
        line[cnt++] = Line(a, b-a);
    }
    int m = HalfplaneIntersection(line, cnt, poly);
    double area = 0;
    for(int i = 0;i<m;i++)
    {
        area += Cross(poly[i], poly[(i+1)%m]);
    }
    area = fabs(area)/2;
    printf("%.1f\n",area);
    return 0;
}