1. 程式人生 > >半平面交詳解

半平面交詳解

定義:

半平面:

顧名思義,就是平面的一半。一條直線會把平面分成兩部分,就是兩個半平面。對於半平面,我們可以用直線方程式如:\(ax + by >= c\) 表示,更常用的是用直線表示。

半平面交:

顧名思義,就是多個半平面求交集。其結果可能是一個凸多邊形、無窮平面、直線、線段、點等。

多邊形的核:

如果多邊形中存在一個區域使得在區域中可以看到多邊形中任意位置(反之亦然),則這個區域就是多邊形的核。可以用半平面交來求解。

極點:

\((x,y)\) 與原點的連線與 \(x\) 軸的夾角,其範圍為 [0,360].

為了方便求解我們假設所有的直線的左側為我們所需要的半平面。

一般來說求解半平面交有兩種方法

① 分治法 \(O(n\log_2 n)\)

② 增量法 \(O(n\log_2 n)\)

但是在這裡由於分治法常數較大,程式碼實現較第二種複雜,所以我們著重介紹第二種方法。

演算法流程

① 將所有直線極角排序,角度相同的保留下需要的一個

② 用一個雙端佇列儲存當前半平面交,每次通過判斷隊首與隊尾第一個交點是否滿足當前直線來更新

③ 先用隊尾判定隊首交點是否合法,再用隊首判斷隊尾交點是否合法

④ 最後求出來的半平面交是一個凸多邊形

極角排序

如下圖排序後

\((a,b),(c,d),(e,f),(g,h),(i,j)\)

所以我們可以更方便的逆時針依次構造,半平面交

由於我們規定了所有的直線的左側為我們所需要的半平面

所以極角相同的直線,我們保留最靠右的。

png

構造半平面交

我們可以依照以下流程來構造一個半平面交,並且構造完成的半平面交有多種情況

① 直線、線段、點不合法

② 凸多邊形,無窮平面(可以增加4個用於限制的半平面,使得平面變得有限)

我們維護兩個雙端佇列

一個儲存當前有用的直線(半平面),一個儲存半平面交的點。

我們依次加入每一條直線(半平面),在加入之前先將之前儲存了的點,但不是最終半平面交中的點彈出佇列

如下圖我們首先讓 直線\(AB\) 進入佇列,再加入直線\(FG\),並且求出 \(AB\)\(FG\)

的交點 \(H\),並且把它加入第二個佇列裡

那麼加入 直線\(CD\) 時我們會發現, \(H\) 這個點在半平面之外,那麼就將它彈出佇列,同時將這條邊也彈出佇列。

所以只要一個點在加入的這條直線的右邊我們就將它彈出。

png

為什麼我們要使用雙端佇列?

可以明顯的發現,雙端佇列的中的點是逆時針排列的

且滿足一定的單調性

當前加入一條直線,如果佇列中的 \(I\) 號點不會彈出佇列,那麼 \(J\) 也不會 (\(I>J\)

且最開始的點,也可能在半平面交之外

假設有一條直線在 \(AB\) 上方,但是它的極角比 \(AB\) 大,那麼就會讓 \(E\) 點,彈出

所以每次彈出,我們都從隊首和隊尾分別彈出。

刪去不合法的點

最後隊首和隊尾都會剩下一些點,它們在半平面交之外,分別使用第一個佇列的隊尾隊首判斷一下

按照之前的方法彈出即可

png

~論文PPT上的圖~

POJ2451模板題,求半平面交面積

AC程式碼,有註釋

#include<cmath>
#include<cstring>
#include<vector>
#include<cstdio>
#include<iostream>
#include<algorithm>
using namespace std;
const double eps=1e-6;
const int maxn=2e5+10;
const double Pi=acos(-1.00);
inline int dcmp(double x)
{
    if(x>eps)return 1;
    return x<-eps?-1:0;
}
struct Vector
{
    double x,y;
    Vector(double X=0,double Y=0)
    {
        x=X,y=Y;
    }
    bool operator == (const Vector &b)const 
    {
        return dcmp(x-b.x)==0&&dcmp(y-b.y)==0;
    }
    double angle()
    {
        return atan2(y,x);//求出極角
    }
};
typedef Vector Point;
Vector operator + (Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator - (Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator * (Vector a,double b){return Vector(a.x*b,a.y*b);}
Vector operator / (Vector a,double b){return Vector(a.x/b,a.y/b);}
struct Line
{
    Point s,t;
    double ang;
    Line(Point X=Vector(),Point Y=Vector())
    {
        s=X,t=Y,ang=(Y-X).angle();
    }
};
typedef Line Segment;
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;
}
bool is_parallel(Line a,Line b)//判斷a,b直線是否平行
{
    return dcmp(cross(a.t-a.s,b.t-b.s))==0;
}
Point intersection(Line a,Line b)//求出a,b的交點
{
    return a.s+(a.t-a.s)*(cross(b.t-b.s,a.s-b.s)/cross(a.t-a.s,b.t-b.s));
}
double area(Point *p,int n)//求出多邊形的面積
{
    double res=0;
    p[n+1]=p[1];
    for(int i=1;i<=n;i++)res+=cross(p[i],p[i+1]);
    return fabs(res/2);
}
bool operator < (const Line &a,const Line &b)//極角排序,如果極角相同則,選擇最靠左的直線
{
    double r=a.ang-b.ang;
    if(dcmp(r)!=0)return dcmp(r)==-1;
    return dcmp(cross(a.t-a.s,b.t-a.s))==-1;
}
bool OnRight(Line a,Point b)//檢查b是否在a直線的右邊
{
    return dcmp(cross(a.t-a.s,b-a.s))<0;
}
bool SI(Line *l,int n,Point *s,int &m)//增量法求半平面交
{
    static Line que[maxn];
    static Point que2[maxn];//兩個雙端佇列
    int head=0,tail=0;
    sort(l+1,l+1+n);
    que[0]=l[1];
    for(int i=2;i<=n;i++)
    if(dcmp(l[i].ang-l[i-1].ang)!=0)//極角相等的直線,取一個
    {
        if(head<tail&&(is_parallel(que[head],que[head+1])||is_parallel(que[tail],que[tail-1])))return false;//如果兩個直線共線,但是極角不同,則沒有半平面交
        while(head<tail&&OnRight(l[i],que2[tail-1]))tail--;//如果在直線右邊,刪除點
        while(head<tail&&OnRight(l[i],que2[head]))head++;
        que[++tail]=l[i];
        if(head<tail)que2[tail-1]=intersection(que[tail],que[tail-1]);//加入新點
    }
    while(head<tail&&OnRight(que[head],que2[tail-1]))tail--;//刪去多餘點
    while(head<tail&&OnRight(que[tail],que2[head]))head++;
    if(tail-head<=1)return false;//只有一個點或零個點,沒有半平面交
    que2[tail]=intersection(que[head],que[tail]);//加入最後一條邊,和第一條邊的交點
    m=0;
    for(int i=head;i<=tail;i++)s[++m]=que2[i];
    return true;
}
const double lim=10000;
int n,m;
Point p[maxn];
Line l[maxn];
double solve()
{
    Point a=Point(0,0);//加入最大限制,防止半平面交無限大
    Point b=Point(lim,0);
    Point c=Point(lim,lim);
    Point d=Point(0,lim);
    l[++n]=Line(a,b);
    l[++n]=Line(b,c);
    l[++n]=Line(c,d);
    l[++n]=Line(d,a);
    if(!SI(l,n,p,m))return 0;
    return area(p,m);
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        Point a,b;
        scanf("%lf%lf%lf%lf",&a.x,&a.y,&b.x,&b.y);
        l[i]=Line(a,b);
    }
    printf("%.1f\n",solve());
}