我的半平面交模板

昨天研究了一下午的半平面交算法，发现要看懂网上的模板真的是不容易啊。。。。orz，然后就从最基本的开始了，先去打印了一份zzy的n（logn）求半平面交的论文

研究了一发，大概思路懂了，然后看到网上有一份别人推荐的模板，就接着研究模板去了，后来看着模板，自己再涂涂改改，就整理出了我的第一份半平面交模板，可喜可贺~~~，下面就是我的模板还有注释一些我的理解

其实就是poj 的2451

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const double eps=1e-8;
int dq[200005],top,bot,pn;
double marea;
struct P
{
    double x,y;
    P (){}
    P (double x,double y):x(x),y(y){
    }
    P operator -(P p){
        return P(x-p.x,y-p.y);
    }
    P operator +(P p){
        return P(x+p.x,y+p.y);
    }
    P operator *(double d){
        return P(x*d,y*d);
    }
    double det(P p){
        return x*p.y-y*p.x;
    }
}p[200005];
struct Line
{
    P a,b;
    double angle;
}l[1000000];
//判断大小函数。
int dblcmp(double k)
{
    if(fabs(k)<eps)
    return 0;
    return k>0?1:-1;
}
//极角排序，相同的则按靠里面的放前
bool cmp(Line a,Line b)
{
    int d=dblcmp(a.angle-b.angle);
    if(!d)
    return dblcmp((b.a-a.a).det(b.b-a.a))>0;
    return d<0;
}
//求直线a,b的交点函数
P in(Line a,Line b)
{
    return a.a+(a.b-a.a)*((b.b-b.a).det(b.a-a.a) / (b.b-b.a).det(a.b-a.a));
}
//b,c为半平面的顶点的两条直线，若它们的交点在新加入的半平面外，则返回1；
bool judge(Line a,Line b,Line c)
{
    P r=in(b,c);
    return dblcmp((a.a-r).det(a.b-r))<0;
}
int n;
void HalfPlaneIntersect()
{
    int i,j;
    //按极角排序
    sort(l,l+n,cmp);
    for(i=0,j=0;i<n;i++)
    {
        if(dblcmp(l[i].angle-l[j].angle)>0)  //存在极角相同的则取第一个
        l[++j]=l[i];
    }
    n=j+1;
    //dq为双端队列，存线段的编号信息
    dq[0]=0;
    dq[1]=1;
    top=1;
    bot=0;
    for(i=2;i<n;i++)
    {
        //若顶端的两条直线的交点在新半平面l[i]外，那么这个顶点就弹出
        while(top>bot && judge(l[i],l[dq[top]],l[dq[top-1]])) top--;
        //旋转一圈左右，如果底部的交点在新加入的半平面外，底部弹出
        while(top>bot && judge(l[i],l[dq[bot]],l[dq[bot+1]])) bot++;
        dq[++top]=i;
    }
    //若顶部的交点不在底部的半平面内，则顶部弹出
    while(top>bot && judge(l[dq[bot]],l[dq[top]],l[dq[top-1]])) top--;
    //同理若底部的交点不在顶部的半平面内，则底部弹出
    while(top>bot && judge(l[dq[top]],l[dq[bot]],l[dq[bot+1]])) bot++;
    //底部加入进来，方便计算半平面交成的凸多边形顶点
    dq[++top]=dq[bot];
    for(pn=0,i=bot;i<top;i++,pn++)
    {
        p[pn]=in(l[dq[i+1]],l[dq[i]]);
    }
}
void getarea()
{
    marea=0;
    if(pn<3)
    marea=0;
    //面积计算公式,也可以用叉乘直接算
    for(int i=0;i<pn;i++)
    {
        marea+=(p[i].x*p[i+1].y-p[i].y*p[i+1].x);
    }
    marea+=(p[pn-1].x*p[0].y-p[pn-1].y*p[0].x);
    marea=fabs(marea);
    marea/=2.0;
}
int main()
{
    scanf("%d",&n);
    for(int i=0;i<n;i++)
    {
        scanf("%lf%lf%lf%lf",&l[i].a.x,&l[i].a.y,&l[i].b.x,&l[i].b.y);
        l[i].angle=atan2(l[i].b.y-l[i].a.y,l[i].b.x-l[i].a.x);
    }
    //添加题意中的边界
    l[n].a.x=0;l[n].a.y=0;l[n].b.x=10000;l[n].b.y=0;
    l[n].angle=atan2(l[n].b.y-l[n].a.y,l[n].b.x-l[n].a.x);
    n++;
    l[n].a.x=10000;l[n].a.y=0;l[n].b.x=10000;l[n].b.y=10000;
    l[n].angle=atan2(l[n].b.y-l[n].a.y,l[n].b.x-l[n].a.x);
    n++;
    l[n].a.x=10000;l[n].a.y=10000;l[n].b.x=0;l[n].b.y=10000;
    l[n].angle=atan2(l[n].b.y-l[n].a.y,l[n].b.x-l[n].a.x);
    n++;
    l[n].a.x=0;l[n].a.y=10000;l[n].b.x=0;l[n].b.y=0;
    l[n].angle=atan2(l[n].b.y-l[n].a.y,l[n].b.x-l[n].a.x);
    n++;
    //半平面交
    HalfPlaneIntersect();
    //求面积
    getarea();
    printf("%.1lf\n",marea);
    return 0;
}