計算幾何-半平面交

阿新 • • 發佈：2022-05-29

計算幾何-半平面交

半平面

平面內的一條直線把這個平面分成兩部分，每一部分對這個平面來說，都叫做半平面。包括這條直線的半平面叫做閉半平面，否則叫做開半平面。

解析式為 $Ax + By +C >=0$或$Ax + By +C <=0$。

在計算幾何中用向量表示，整個題統一以向量的左側或右側為半平面。

半平面交

半平面交就是多個半平面的交集。半平面交是一個點集。

它可以理解為向量集中每一個向量的右側的交，或者是下面方程組的解。

\[ \begin{equation} \begin{cases} A1x+B1y+C1\ge0\newline A2x+B2y+C2\ge0\newline ~~~~~~~~~~~~~~\dots \end{cases} \end{equation} \]

多邊形的核

如果一個點集中的點與多邊形上任意一點的連線與多邊形沒有其他交點，那麼這個點集被稱為多邊形的核。

把多邊形的每條邊看成是首尾相連的向量，那麼這些向量在多邊形內部方向的半平面交就是多邊形的核。

求法

D&C演算法

該演算法是基於分治思想的：

將$n$個半平面分成兩個$n/2$的集合;
對兩個子集和遞迴求解半平面交;
將前一步算出來的兩個交利用平面掃描法求解。

時間複雜度$(n \log n)$這個演算法並不常用，主要介紹的是下面這個。

S&I演算法

該演算法是在2006年有中國隊隊員朱澤園提出來的“排序增量法”。

假設給出$n$條直線，求這$n$條直線的左方半平面的交集：

首先對這$n$條直線按極角排序；
用一個佇列去維護半平面的交集，和相鄰兩條直線的交點；
每次加入新的直線時判斷是否有交點在該直線的右面，如果是則彈出直線，先判隊尾再判隊首，注意判斷平行情況；
最後佇列中的交集即為半平面交。

【模板】半平面交

Code:

#include<bits/stdc++.h>
#define eps 1e-8
using namespace std;
const int maxn=1010;
struct geometric{
    double x,y;
    geometric(double X=0,double Y=0):x(X),y(Y) {}
    friend geometric operator + (const geometric a,const geometric b){return geometric(a.x+b.x,a.y+b.y);} 
    friend geometric operator - (const geometric a,const geometric b){return geometric(a.x-b.x,a.y-b.y);} 
    friend geometric operator * (const geometric a,double p){return geometric(a.x*p,a.y*p);}
    friend geometric operator / (const geometric a,double p){return geometric(a.x/p,a.y/p);}
    double dis(geometric a,geometric b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
    double dot(geometric a1,geometric a2,geometric b1,geometric b2){return (a2.x-a1.x)*(b2.x-b1.x)+(a2.y-a1.y)*(b2.y-b1.y);}
    double cross(geometric a1,geometric a2,geometric b1,geometric b2){return (a2.x-a1.x)*(b2.y-b1.y)-(a2.y-a1.y)*(b2.x-b1.x);}
    double corner(geometric a1,geometric a2,geometric b1,geometric b2){return dot(a1,a1,b1,b2)/(dis(a1,a2)*dis(b1,b2));}
    double area(geometric a1,geometric a2,geometric b1,geometric b2){return fabs(cross(a1,a2,b1,b2));}
    double angle(geometric a){return atan2(a.y,a.x);}
}opt;
int n,m,tot,head=1,tail=1;double ans;
geometric data[maxn],origin,T[maxn];
struct line{
    geometric A,B;double An;
    line(geometric a,geometric b):A(a),B(b) {An=opt.angle(B);}
    line(){}
    bool operator < (const line &a)const{return An<a.An;}
    geometric sdot(line a,line b){
        geometric c=a.A-b.A;
        double k=opt.cross(origin,b.B,origin,c)/opt.cross(origin,a.B,origin,b.B);
        return a.A+a.B*k;
    }
}q[maxn],p[maxn],take;
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        scanf("%d",&m);
        for(int j=1;j<=m;j++)
        scanf("%lf%lf",&data[j].x,&data[j].y);
        for(int j=1;j<=m;j++)
        {
            if(j==m)p[++tot]=line(data[j],data[1]-data[j]);
            else p[++tot]=line(data[j],data[j+1]-data[j]);
        }
    }
    sort(p+1,p+tot+1);
    q[head]=p[head];
    for(int i=2;i<=tot;i++)
    {
        while(head<tail&&opt.cross(origin,p[i].B,p[i].A,T[tail-1])<=eps)tail--;
        while(head<tail&&opt.cross(origin,p[i].B,p[i].A,T[head])<=eps)head++;
        q[++tail]=p[i];
        if(fabs(opt.cross(origin,q[tail].B,origin,q[tail-1].B))<=eps)
        {
            tail--;
            if(opt.cross(origin,q[tail].B,q[tail].A,p[i].A)>eps)q[tail]=p[i];
        }
        if(head<tail)T[tail-1]=take.sdot(q[tail-1],q[tail]);
    }   
    while(head<tail&&opt.cross(origin,q[head].B,q[head].A,T[tail-1])<=eps)tail--;
    if(tail-head>1)
    T[tail]=take.sdot(q[head],q[tail]);
    for(int i=head;i<=tail;i++)
    {
        if(i==tail)ans+=opt.cross(origin,T[i],origin,T[head]);
        else ans+=opt.cross(origin,T[i],origin,T[i+1]);
    }
    printf("%.3lf",ans/2);
    return 0;
}