算法竞赛笔记.二维计算几何.基本概念/凸包/半平面交

概念

基本定义

点/向量

struct Point{
    double x,y;
    Point operator+(const Point& b)const{return {x+b.x,y+b.y};}
    Point operator-(const Point& b)const{return {x-b.x,y-b.y};}
    Point operator*(double p)const{return {x*p,y*p};}
};
typedef Point Vector;

直线

struct Line{
    Point p;
    Vector v;
    double ang;
    Line(){}
    Line(Point p,Vector v):p(p),v(v){ang=atan2(v.y,v.x);}
};

平面向量的叉积

double cross(Point a1,Point a2,Point b1,Point b2){
    return (a2.x-a1.x)*(b2.y-b1.y)-(b2.x-b1.x)*(a2.y-a1.y);
}

判断点是否在线段上

bool isOnSeg(Point p,Point A,Point B){
    if(dcmp(cross(p,A,p,B)))return 0;
    return dcmp(min(A.x,B.x)-p.x)<=0&&
           dcmp(p.x-max(A.x,B.x))<=0&&
           dcmp(min(A.y,B.y)-p.y)<=0&&
           dcmp(p.y-max(A.y,B.y))<=0;
}

判断点是否在多边形里

int isInPolygon(Point p,Point* s,int m){
    int wn=0;
    for(int i=1;i<=m;i++){
        if(isOnSeg(p,s[i],s[i+1]))return -1;
        int k=dcmp(cross(s[i],s[i+1],s[i],p));
        int d1=dcmp(s[i].y-p.y);
        int d2=dcmp(s[i+1].y-p.y);
        if(k>0&&d1<=0&&d2>0)wn++;
        if(k<0&&d2<=0&&d1>0)wn--;
    }
    return wn!=0;
}

判断线段是否有公共点

bool isSegsIntersect(Point A,Point B,Point C,Point D){
    int d1=dcmp(cross(A,B,A,C));
    int d2=dcmp(cross(A,B,A,D));
    int d3=dcmp(cross(C,D,C,A));
    int d4=dcmp(cross(C,D,C,B));
    if(d1*d2<0&&d3*d4<0)return 1;
    if(d1==0&&isOnSeg(C,A,B))return 1;
    if(d2==0&&isOnSeg(D,A,B))return 1;
    if(d3==0&&isOnSeg(A,C,D))return 1;
    if(d4==0&&isOnSeg(B,C,D))return 1;
    return 0;
}

凸包

凸包是平面上能包含所有给定点的凸多边形。凸包的周长最小。

Andrew 算法求凸包

将所有点按横坐标排序，首先从横坐标最小的点向右（逆时针）出发，使用一个单调栈来维护下凸壳，然后同理从右向左维护上凸壳。

int Andrew(int n,Point *p,Point *s){
    int m=0;
    sort(p+1,p+n+1,[&](Point a,Point b){
        return a.x<b.x||(a.x==b.x&&a.y<b.y);
    });
    for(int i=1;i<=n;i++){
        while(m>=2&&cross(s[m-1],s[m],s[m],p[i])<=0)m--;
        s[++m]=p[i];
    }
    int k=m;
    for(int i=n-1;i>=1;i--){
        while(m>k&&cross(s[m-1],s[m],s[m],p[i])<=0)m--;
        s[++m]=p[i];
    }
    if(m>1)m--;
    s[m+1]=s[1];
    return m;
}

半平面交（S&I 算法）

一条直线把平面分为两个半平面，若干个半平面的交集即为半平面交。

S&I 算法求半平面交的思路可以概括为以下四步：
将所有凸多边形的边看作有向直线（规定左侧为半平面内部），按其方向向量的极角从小到大排序；若极角相同，仅保留最靠左的那条。

利用双端队列（deque）按序处理直线。每加入一条新线，若队列末尾（或队首）两条线的交点落在新线的右侧，说明旧线已失效，将其从队列中弹出。

所有直线入队后，还需用队首直线检查并剔除队尾多余的线，反之亦然，以保证最后剩下的直线能围成一个封闭的凸多边形。

求出队列中相邻直线的交点，这些交点即为半平面交的顶点，最后利用叉积（鞋带公式）计算该多边形的面积。

int dcmp(double x){
    if(fabs(x)<eps)return 0;
    return x<0?-1:1;
}
double cross(Vector a,Vector b){
    return a.x*b.y-a.y*b.x;
}
struct Line{
    Point p;
    Vector v;
    double ang;
    Line(){}
    Line(Point p,Vector v):p(p),v(v){ang=atan2(v.y,v.x);}
};
//极角排序,极角相同则靠左的(限制更严)排前面
bool cmp(const Line& a,const Line& b){
    if(dcmp(a.ang-b.ang)==0)return dcmp(cross(a.v,b.p-a.p))<0;
    return a.ang<b.ang;
}
//判断点p是否在线L的右侧(或线上)
bool onRight(Line L,Point p){
    return dcmp(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 n;
Line L[N],q[N];
Point p[N];
int main(){
    ios::sync_with_stdio(false);cin.tie(0);
    if(!(cin>>n))return 0;
    int cnt=0;
    for(int i=1;i<=n;i++){
        int m;cin>>m;
        vector<Point> poly(m);
        for(int j=0;j<m;j++)cin>>poly[j].x>>poly[j].y;
        for(int j=0;j<m;j++)L[++cnt]=Line(poly[j],poly[(j+1)%m]-poly[j]);
    }
    sort(L+1,L+cnt+1,cmp);
    int tot=0;
    //极角去重
    for(int i=1;i<=cnt;i++){
        if(i>1&&dcmp(L[i].ang-L[i-1].ang)==0)continue;
        L[++tot]=L[i];
    }
    int head=1,tail=0;
    for(int i=1;i<=tot;i++){
        //新加入的线把队尾的交点切掉了
        while(head<tail&&onRight(L[i],getIntersection(q[tail-1],q[tail])))tail--;
        //新加入的线把队首的交点切掉了
        while(head<tail&&onRight(L[i],getIntersection(q[head],q[head+1])))head++;
        q[++tail]=L[i];
        //判断平行反向导致闭合空间为空
        if(head<tail&&dcmp(cross(q[tail-1].v,q[tail].v))<=0){
            cout<<fixed<<setprecision(3)<<0.000<<"\n";
            return 0;
        }
    }
    //首尾互相更新去冗余
    while(head<tail&&onRight(q[head],getIntersection(q[tail-1],q[tail])))tail--;
    while(head<tail&&onRight(q[tail],getIntersection(q[head],q[head+1])))head++;
    //不足以构成多边形
    if(tail-head+1<3||dcmp(cross(q[tail].v,q[head].v))<=0){
        cout<<fixed<<setprecision(3)<<0.000<<"\n";
        return 0;
    }
    int m=0;
    for(int i=head;i<tail;i++)p[++m]=getIntersection(q[i],q[i+1]);
    p[++m]=getIntersection(q[tail],q[head]);//收尾相接的交点
    //计算多边形面积
    double ans=0;
    for(int i=1;i<=m;i++)ans+=cross(p[i],p[i%m+1]);
    ans=fabs(ans)/2.0;
    cout<<fixed<<setprecision(3)<<ans<<"\n";
    return 0;
}