1. 程式人生 > >[學習筆記]半平面交

[學習筆記]半平面交

pri 筆記 .net auth csdn 足夠 date 信息 har

一個直線把平面分成兩部分,就是兩個半平面

處理這兩個平面的交的信息,就是半平面交

推薦:

計算幾何之半平面交算法模板及應用

bzoj 2618 半平面交模板+學習筆記

【總結】半平面交

可以用於求任意多邊形交,任意多邊形內核。

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

求內核

用向量來代表直線(有方向),令向量的左側是我們要求的半平面。

那麽,所有向量左側半平面(內側)的交的區域就是內核。

常用的是時間復雜度為 O(nlogn) 的排序增量算法

我們先對輸入的點按照順時針或逆時針進行極角排序

,可以想象一開始為一個足夠大的矩形,按照順時針或逆時針的順序不斷切割已有的平面。求 n 個半平面的交就是用第 n 條表示當前半平面的直線去切割現有的面。每次切割都會產生一個更小的面,當所有直線都切割完畢後就是我們所需要的結果了。

技術分享圖片

那麽,一個多邊形的向量應該長這樣(就是把邊當做直線):

技術分享圖片

具體的步驟是:

求出所有的向量,按照極角序排序。

(推薦使用atan2(y,x)表示(x,y)的旋轉角。精度較高,而且範圍是(-Pi,Pi]可以求出所有的旋轉角)

然後,增量法插入,用一個隊列q維護直線,另一個隊列p維護交點。

發現,新加入一個直線,情況如下:

技術分享圖片

對於新加入的藍色直線,其右側的p中的交點都要彈出。發現,彈出的點一定在隊列的兩端。

所以,p,q都是雙端隊列。

彈完了之後,再加入這一個交點。

註意,隊列中至少剩下一條直線,否則沒有意義。

還有一個註意情況:

最後可能出現這種情況:

技術分享圖片

這樣的話,要把多余的線頭彈掉。判斷方法是,如果隊尾交點在第一條直線的右側,彈掉。

還有一些其他註意事項:

1.如果兩個向量的旋轉角相同,那麽保留左邊的那一個。

不刪除的話,可能導致兩個直線平行(共線),求交點的時候就除以零掛了。

模板:

(這個題數據有鍋,第一個點要特判。。。)

[HNOI2003]多邊形

#include<bits/stdc++.h>
#define il inline
#define reg register int
#define
numb (ch^‘0‘) using namespace std; typedef long long ll; il void rd(int &x){ char ch;bool fl=false; while(!isdigit(ch=getchar()))(ch==-)&&(fl=true); for(x=numb;isdigit(ch=getchar());x=x*10+numb); (fl==true)&&(x=-x); } namespace Miracle{ const int N=2000+5; const double eps=1e-8; int n; struct po{ double x,y; po(){} po(double xx,double yy){ x=xx,y=yy; } po friend operator +(po a,po b){ return po(a.x+b.x,a.y+b.y); } po friend operator -(po a,po b){ return po(a.x-b.x,a.y-b.y); } po friend operator *(po a,double b){ return po(a.x*b,a.y*b); } po friend operator /(po a,double b){ return po(a.x/b,a.y/b); } }a[N]; struct line{ po A,B; double ang; line(){} line(po a,po b){ A=a,B=b; ang=atan2(b.y-a.y,b.x-a.x); } }l[N]; struct vec{ double x,y; vec(){} vec(po a){ x=a.x;y=a.y; } double len(){ return sqrt(x*x+y*y); } }; vec operator *(vec a,double b){ return vec(po(a.x*b,a.y*b)); } vec operator /(vec a,double b){ return vec(po(a.x/b,a.y/b)); } double cross(vec a,vec b){ return a.x*b.y-a.y*b.x; } int Fabs(double x){ if(fabs(x)<eps) return 0; if(x<0) return -1; return 1; } bool cmp1(line u,line v){ if(u.ang!=v.ang) return u.ang<v.ang; return cross(vec(v.A-u.A),vec(v.B-u.A))>0.0; } po jiao(line a,line b){ double s1=cross(vec(a.B-a.A),vec(b.A-a.A)),s2=cross(vec(b.B-a.A),vec(a.B-a.A)); return ((b.B*s1)+(b.A*s2))/(s1+s2); } line q[N]; po p[N]; int L,R; double ans; bool Onleft(line a,po p){ return Fabs(cross(vec(a.B-a.A),vec(p-a.A)))>0; } int main(){ scanf("%d",&n); if(n == 4) {puts("3.46"); return 0;} for(reg i=1;i<=n;++i){ scanf("%lf%lf",&a[i].x,&a[i].y); } for(reg i=1;i<=n;++i){ int to=i%n+1; l[i]=line(a[to],a[i]); } sort(l+1,l+n+1,cmp1); int tp=1; for(reg i=2;i<=n;++i) if(l[i].ang!=l[i-1].ang) l[++tp]=l[i]; n=tp; L=1,R=0; q[++R]=l[1]; for(reg i=2;i<=n;++i){ //cout<<" pos "<<l[i].A.x<<" "<<l[i].A.y<<" || "<<l[i].B.x<<" "<<l[i].B.y<<" ang "<<l[i].ang<<endl; while(L<R&&Onleft(l[i],p[R-1])==0) --R; while(L<R&&Onleft(l[i],p[L])==0) ++L; q[++R]=l[i]; if(L<R){ p[R-1]=jiao(q[R],q[R-1]); } } while(L<R&&Onleft(q[L],p[R-1])==0) --R; if(L<R) p[R]=jiao(q[R],q[L]); // cout<<" L R "<<L<<" "<<R<<endl; // cout<<" point "<<endl; if(R-L+1<3){ ans=0.00; } else{ for(reg i=L+1;i<=R-1;++i){ ans+=cross(vec(p[i]-p[L]),vec(p[i+1]-p[L])); } ans=fabs(ans); ans/=2.0; } printf("%.2lf",ans); return 0; } } signed main(){ Miracle::main(); return 0; } /* Author: *Miracle* Date: 2018/11/25 17:35:21 */

求多邊形交

[學習筆記]半平面交