半平面的核

【题意】给出很多个半平面，这里每个半平面由线段组成，都是指向线段方向的左边表示有
(x1 - x) * (y2 - y) - (x2 - x) * (y1 - y) >=0 ( >=0 表示左边，<=0 表示右边)
要你求个半平面的核，就是所有半平面所围成的面积
【算法】O(nlogn)的半平面交算法，最后统计出得到的多边形的点，然后利用叉积公式求出面积就行了

#include<cstdio>
#include<vector>
#include<cmath>
#include<algorithm>
using namespace std;
const double eps=1e-10,big=10000.0;
const int maxn = 20010;
struct point
{
    double x,y;
};
struct polygon   //存放最后半平面交中相邻边的交点，就是一个多边形的所有点
{
    int n;
    point p[maxn];
};
struct line   //半平面，这里是线段
{
    point a,b;
};
double at2[maxn];
int ord[maxn],dq[maxn+1],lnum;

int n;
polygon pg;

line ls[maxn]; //半平面集合
inline int sig(double k)   //判是不是等于0，返回-1，0，1，分别是小于，等于，大于
{
    return (k < -eps)? -1: k > eps;
}
//叉积>0代表在左边，<0代表在右边，=0代表共线
//e是否在o->s的左边onleft(sig(multi))>=0
inline double multi(point o, point s, point e)
{
    //构造向量，然后返回叉积
    return (s.x-o.x)*(e.y-o.y)-(e.x-o.x)*(s.y-o.y);
}
//直线求交点
point isIntersected(point s1, point e1, point s2, point e2)
{
    double dot1,dot2;
    point pp;
    dot1 = multi(s2,e1,s1);
    dot2 = multi(e1,e2,s1);
    pp.x = (s2.x * dot2 + e2.x * dot1) / (dot2 + dot1);
    pp.y = (s2.y * dot2 + e2.y * dot1) / (dot2 + dot1);
    return pp;
}
//象限排序
inline bool cmp(int u,int v)
{
    if(sig(at2[u]-at2[v])==0)
        return sig(multi(ls[v].a,ls[v].b,ls[u].b))>=0;
    return at2[u]<at2[v];
}
//判断半平面的交点在当前半平面外
bool judgein(int x,int y,int z)
{
    point pnt = isIntersected(ls[x].a, ls[x].b, ls[y].a, ls[y].b); //求交点
    return  sig(multi(ls[z].a,ls[z].b,pnt)) < 0;  
    //判断交点位置，如果在右面，返回true，如果要排除三点共线，改成<=
}
//半平面交
void HalfPlaneIntersection(polygon &pg)   //预处理
{
    int n = lnum , tmpn , i;
    /*  对于atan2(y,x)
        结果为正表示从 X 轴逆时针旋转的角度，结果为负表示从 X 轴顺时针旋转的角度。
        　　 atan2(a, b) 与 atan(a/b)稍有不同，atan2(a,b)的取值范围介于 -pi 到 pi 之间（不包括 -pi），
        　　 而atan(a/b)的取值范围介于-pi/2到pi/2之间（不包括±pi)*/
    for(i = 0 ; i < n ; i ++)
    {
        //atan2(y,x)求出每条线段对应坐标系的角度
        at2[i] = atan2( ls[i].b.y - ls[i].a.y, ls[i].b.x - ls[i].a.x);
        ord[i] = i;
    }
    sort(ord , ord + n , cmp);
    for (i = 1 , tmpn = 1 ; i < n ; i++) //处理重线的情况
        if( sig(at2[ord[i-1]] - at2[ord[i]]) != 0 ) ord[tmpn++] = ord[i];
    n = tmpn;
    //圈地
    int bot = 1,top = bot + 1; //双端栈，bot为栈底，top为栈顶
    dq[bot] = ord[0];
    dq[top] = ord[1]; //先压两根线进栈
    for(i = 2 ; i < n ; i ++)
    {
        //bot < top 表示要保证栈里至少有2条线段，如果剩下1条，就不继续退栈
        //judgein，判断如果栈中两条线的交点如果在当前插入先的右边，就退栈
        while( bot < top && judgein(dq[top-1] , dq[top] , ord[i]) ) top--;
        //对栈顶要同样的操作
        while( bot < top && judgein(dq[bot+1] , dq[bot] , ord[i]) ) bot++;
        dq[++top] = ord[i];
    }
    //最后还要处理一下栈里面存在的栈顶的线在栈底交点末尾位置，或者栈顶在栈尾两条线的右边
    while( bot < top && judgein(dq[top-1] , dq[top] , dq[bot]) ) top--;
    while( bot < top && judgein(dq[bot+1] , dq[bot] , dq[top]) ) bot++;
    //最后一条线是重合的
    dq[--bot] = dq[top];
    //求多边形
    pg.n = 0;
    for(i = bot + 1; i <= top ; i++) //求相邻两条线的交点
        pg.p[pg.n++] = isIntersected(ls[dq[i-1]].a, ls[dq[i-1]].b, ls[dq[i]].a,ls[dq[i]].b);
}

inline void add(double a,double b,double c,double d)
{
    //添加线段
    ls[lnum].a.x = a;
    ls[lnum].a.y = b;
    ls[lnum].b.x = c;
    ls[lnum].b.y = d;
    lnum++;
}
int main()
{
    int n,i;
    scanf("%d",&n);
    double a,b,c,d;
    for(i = 0 ; i < n ; i ++)
    {
        //输入代表一条向量(x = (c - a),y = (d - b));
        scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
        add(a,b,c,d);
    }
    //下面是构造一个大矩形边界
    add(0,0,big,0);//down
    add(big,0,big,big);//right
    add(big,big,0,big);//up
    add(0,big,0,0);//left
    HalfPlaneIntersection(pg); //求半平面交
    double area = 0;
    n = pg.n;
    ///最后多边形的各个点保存在pg里面
    for(i = 0 ; i < n ; i ++)
        area += pg.p[i].x * pg.p[(i+1)%n].y - pg.p[(i+1)%n].x * pg.p[i].y;
        //x1 * y2 - x2 * y1用叉积求多边形面积
    area=fabs(area)/2.0; 
    //所有面积应该是三角形面积之和，而叉积求出来的是四边形的面积和，所以要除2

    printf("%.1f\n",area);
    return 0;
}