半平面交

github傳送門

簡介

Emmmm學完旋轉卡殼感覺自己已經是個廢人了..
修整了一個周末, 回來接著跟計算幾何勢力硬幹... (這個周末是不是有點長?)
今天就講講半平面交吧.
請自己回顧必修五 線性規劃相關知識...
什麽是半平面?
就是一條直線一側的點構成的集合..
用解析幾何的觀點來看就是滿足\(Ax+By+C<0\)這個不等式的點的集合.
那麽半平面交自然就是許多這樣的集合的交集咯~

最後就很像線性規劃做出來的可行域一樣...

半平面交可以長這樣

這樣

甚至這樣

也就是說, 半平面交的結果可能是無界的一大片, 可能是一個有界的凸多邊形, 也可能是點、線段、直線(沒圖╮(╯_╰)╭), 甚至空集

引例

題目

先來道不那麽標準的板子題:戳這裏..

分析

這個很顯然就是要求朝上的無界的一片的那種.. md就是個下凸殼啊...

所以我們按照凸殼的方式試著思考一下.
我們先將直線按斜率排序(都按斜截式給你了顯然不會不存在....) 如果斜率一樣截距小的可以直接不要了(肯定會被擋住..)
維護一個棧, 先將第一條直線壓入, 這條直線肯定能露出來(誰讓斜率最大了呢)
然後我們考慮新加入的直線對前面的直線的影響.

我們可以看到,當新加入一條斜率比較小的直線(綠色或黃色, 按理說應該是平行的, 如果看著不大平行就湊合看吧)

• 如果比前面兩條直線的交點要高(比如綠色) 那麽前一條直線就會被遮住(紅色左邊被綠遮, 右邊被黑遮...);
• 而如果比前面兩條直線的交點低(比如黃色) 那麽前一條直線就會有一部分露出來.

所以出現了綠色的情況我們就把紅色退棧. 然後繼續判斷.

是不跟Graham掃描法一個樣紙?

時間復雜度\(O(nlogn)\) 復雜度仍然卡在排序上. 這是整數而且斜率極差不超過1e6所以可以很方便的用基數排序降到\(O(n)\)虐場.. 然而窩還是用sort 因為省事...

現在唯一的問題就是判斷直線在點上面了(根據樣例可以得知交點一個點是不算子線段的, 也要算在上面).. 但這就是解析幾何的問題了.. (其實可以用類似於方向向量的方法解, 只是我不想再寫叉積了(想起被叉積支配的恐懼..))

我們令前兩條直線分別為\(y=k_1x+b_1, y=k_2x+b_2\)

代碼

下面貼上1A代碼:

#include <cmath>
#include <cstdio>
#include <algorithm>
const double eps=1e-9;
inline int dcmp(const double &a){
    if(fabs(a)<eps) return 0;
    return a<0?-1:1;
}
struct line{
    double k,b;
    int id;
}p[50050],stk[50050]; int tp,n;
inline bool cmp(const line &a,const line &b){
    if(!dcmp(a.k-b.k)) return a.b>b.b;
    return a.k>b.k;
}
inline bool cmp2(const line &a,const line &b){
    return a.id<b.id;
}
void halfPlanesIntersection(){
    stk[++tp]=p[0];
    for(int i=1;i<n;++i){
        if(!dcmp(p[i].k-stk[tp].k)) continue;
        while(tp>1&&dcmp((stk[tp].b-stk[tp-1].b)*(stk[tp-1].k-p[i].k)+(stk[tp-1].k-stk[tp].k)*(stk[tp-1].b-p[i].b))<1) --tp;
        stk[++tp]=p[i];
    }
}
int main(){
    scanf("%d",&n);
    for(int i=0;i<n;++i){
        scanf("%lf%lf",&p[i].k,&p[i].b);
        p[i].id=i+1;
    }
    std::sort(p, p+n, cmp);
    halfPlanesIntersection();
    std::sort(stk+1, stk+tp+1, cmp2);
    for(int i=1;i<=tp;++i)
        printf("%d ",stk[i].id);
}

Emmmm這個題還有一種SAO操作就是Trinkle大神提出的半平面交對偶轉凸包做法... orz
大體上就是把直線\(y=kx+b\)視為一個點\((k,b)\), 然後求一波上凸殼.
等下個題看下能不能用上吧.

正題

剛才那個題太解析幾何了, 拉個平時數學課不睡覺的人就能做..
我們來看更"計算幾何"的東西..

比如半平面除了能用線性規劃裏面的\(Ax+By+C<0\)以外, 還可以用向量的左邊來表示...
向量的左邊嘛, 當然就是叉積<0咯(被叉積支配的恐懼...)

方法一

有一種非常簡單的\(O(n^2)\)做法, 我們假設這個東西是有界的, 於是就從無窮遠處套個框..
然後用每一條邊去切割這個多邊形就行了.

我們枚舉一條邊的時候, 分別判斷其他的點是否在左邊.. 如果在左邊顯然就不在半平面交了.
然後有一個端點在左邊的邊會與正在枚舉的邊有交點(8), 我們就把這個新的交點加入我們的程序即可.
一共要枚舉\(n\)條邊, 每條邊要判斷\(n\)次, 所以時間復雜度顯然是\(O(n^2)\)的.
這種做法有個非常大的好處就是好理解....

代碼(無)

等用到了再寫..

方法二

還有一種方法叫分治.. 算法的流程大約是這樣的:

• 在最外面套個框

• 遞歸搞出前一半的半平面交
• 遞歸高出後一半的半平面交

• 利用凸多邊形\(O(n)\)求交(好像是CGI什麽的) 來處理整個的半平面交

時間復雜度也好分析: \(T(n)=2T(\frac n2)+O(n)=>O(nlogn)\)...
但是好像比較麻煩(而且凸多邊形求交真的不是搞一個半平面交???)
有更好用的\(O(nlogn)\)做法所以這個就不常見了...

代碼(無)

代碼應該也不太好寫 等用到了再寫...

方法三

隆重請出我們的方法三...
好像是某ZZY大神自己yy出來的?? orz
然後被某大神優化了一下就出現了現在的版本...
真的是高端大氣上檔次...
不過這樣對於我這種蒟蒻來講就不大好理解了..

求半平面交的方式和凸包很像很像..
先將直線的方向向量起點都平移到原點, 然後按極角排序..

但是直接用叉積就會出現一些問題(就是不喜歡atan2..), 比如紅色不一定能排在紫色的左邊了...(為了方便向量的先後關系應該如彩虹的顏色順序..)
所以我們要分類討論(幾何常見套路), 將向量分成x軸上方(包括x軸正半軸)和x軸下方(包括x軸負半軸)兩類. 很顯然第一類<第二類(因為要寫sort裏的cmp函數就用大於小於了).
然後對於每一類分別求叉積就好了.. 比如橙色在黃色的順時針方向所以橙色<黃色, 比如藍色在紫色的順時針方向所以藍色<紫色, 以此類推.
那麽極角相同呢?

由於我們要求的是交集, 所以要把右邊的排除. 我們可以通過\(\vec a\)的終點-\(\vec b\)的起點與\(\vec b\)做叉積, 來判斷直線的左右關系.

排完序之後就進入了正題.
我們還是看添加一條直線之後會造成什麽影響.

(P.S. 圖裏沒有畫箭頭, 方向都是↗)

根據上面不那麽標準的題目可以分析得到, 當我們新加入一條邊的時候, 便檢測一下前面兩條邊的交點在不在這條新加入的邊所對應的半平面上. 所以對於紅點來說, 它在粉色對應的半平面但不在紫色對應的半平面, 所以如果加入的是粉色那就加入了, 而如果加入的是紫色那麽淺綠是要退棧的..

而我們的半平面交不只是求一個凸殼(所以說引例不標準嘛), 所以圖形可能是封閉的, 那麽新加入的邊也就可能影響到最開始加入的邊...

所以我們要維護一個雙端隊列.. 每次分別判斷隊首隊尾兩條邊的交點在不在對應的半平面內. 如果不在就彈出去.

然後就出現了下一個問題, 就是求完以後, 可能會有很多冗余的直線.

比如我們求完黑色的半平面交以後, 紅色的也能成功的入隊, 但是對答案並沒有什麽貢獻. 隊首也是一樣, 所以我們要再掃, 把這些冗余的出隊.

由於每條直也最多只會進隊一次, 出隊一次, 所以復雜度是\(O(n)\)的,

總復雜度就自然是排序的\(O(nlogn)\)(所以也可以用基數排序做到\(O(n)\), 但我還是用sort..)

然後細節上來說還是比較麻煩的... 用到了好多前面的基礎知識, 然後發現自己全tm不記得了..

代碼(這次有了XD)

貼一道算是模板題的題目吧 bzoj2618
把多邊形的每一條邊拆成直線, 最後求一下凸多邊形的面積就行了(其實有點綜合的意思~)
WA竟然是因為求面積的時候\(pcnt\)打成了\(n\)... 是不該去打星際了..

代碼:

#include <cmath>
#include <cstdio>
#include <algorithm>
const int N=101010;
const double eps=1e-9;
inline int dcmp(const double &a){
    if(fabs(a)<eps) return 0;
    return a<0?-1:1;
}
struct vec{
    double x,y;
    vec(double X=0,double Y=0):x(X),y(Y){}
}pt[N],yM; int pcnt;
inline vec operator-(const vec &a,const vec &b){
    return vec(a.x-b.x,a.y-b.y);
}
inline double operator*(const vec &a,const vec &b){
    return a.x*b.y-a.y*b.x;
}
struct line{
    vec a,b;
    line(){}
    line(vec X,vec Y):a(X),b(Y){}
}p[N],dq[N]; int tp,bt,n; //雙端隊列
inline bool aboveX(const line &a) //是否在x軸上方
{
    if(!dcmp(a.b.y-a.a.y)) return dcmp(a.b.x-a.a.x)>0;
    return dcmp(a.b.y-a.a.y)>0;
}
inline bool cmpa(const line &a,const line &b){
    if(aboveX(a)!=aboveX(b)) return aboveX(a);
    if(!dcmp((a.b-a.a)*(b.b-b.a))) 
        return dcmp((a.b-a.a)*(b.b-a.a))<0;
    return dcmp((a.b-a.a)*(b.b-b.a))>0;
}
inline vec getLinesIntersect(const line &a,const line &b){
    double a1,a2,b1,b2,c1,c2,d;
    a1=a.b.y-a.a.y; b1=a.a.x-a.b.x; c1=a.a*a.b;
    a2=b.b.y-b.a.y; b2=b.a.x-b.b.x; c2=b.a*b.b;
    d=a1*b2-a2*b1;
    return vec((b2*c1-b1*c2)/d,(a1*c2-a2*c1)/d);
}
inline bool pd(const line &a,const line &b,const line &c){
    //判斷直線a和直線b的交點是否在c左邊(在半平面內)
    vec p=getLinesIntersect(a,b);
    return dcmp((p-c.a)*(c.b-c.a))>-1; //true表示不在半平面內
}
void halfPlaneIntersection(){
    int i,j; vec t;
    for(i=0,j=0;i<n;++i)
        if(dcmp((p[i].b-p[i].a)*(p[j].b-p[j].a)))
            p[++j]=p[i];
    n=j+1; //去掉極角相同
    dq[bt]=p[0];
    dq[++tp]=p[1]; //開始兩條直線入隊
    for(int i=2;i<n;++i){
        while(tp>bt&&pd(dq[tp],dq[tp-1],p[i])) --tp;
        while(bt<tp&&pd(dq[bt],dq[bt+1],p[i])) ++bt;
        dq[++tp]=p[i];
    }
    while(tp>bt&&pd(dq[tp],dq[tp-1],dq[bt])) --tp;
    while(bt<tp&&pd(dq[bt],dq[bt+1],dq[tp])) ++bt;
    dq[++tp]=dq[bt]; //已經保存好繞一圈的直線了
  
    //確定半平面交部分的頂點
    for(i=bt;i<tp;++i,++pcnt)
        pt[pcnt]=getLinesIntersect(dq[i],dq[i+1]);
}
double polyArea(double s=0){
    if(pcnt<3) return 0;
    pt[pcnt]=pt[0];
    for(int i=0;i<pcnt;++i)
        s=s+pt[i]*pt[i+1];
    return 0.5*fabs(s);
}
int main(){
    int T,w=0; 
    scanf("%d",&T);
    while(T--){
        int ps; scanf("%d",&ps);
        for(int i=1;i<=ps;++i) scanf("%lf%lf",&pt[i].x,&pt[i].y);
        pt[0]=pt[ps];
        for(int i=0;i<ps;++i) p[n].a=pt[i],p[n++].b=pt[i+1];
    }
    std::sort(p,p+n,cmpa);
    halfPlaneIntersection();
    double area=polyArea();
    printf("%.3lf",area);
}

其實想清楚的話也不是很難, 但是真的細啊OvO

大約就這樣吧?

【筆記篇】最良心的計算幾何學習筆記（六）