一個直線把平面分成兩部分,就是兩個半平面
處理這兩個平面的交的信息,就是半平面交
推薦:
計算幾何之半平面交算法模板及應用
bzoj 2618 半平面交模板+學習筆記
【總結】半平面交
可以用於求任意多邊形交,任意多邊形內核。
(內核:如果多邊形中存在一個區域使得在區域中可以看到多邊形中任意位置(反之亦然),則這個區域就是多邊形的核。可以用半平面交來求解。)
求內核
用向量來代表直線(有方向),令向量的左側是我們要求的半平面。
那么,所有向量左側半平面(內側)的交的區域就是內核。
常用的是時間復雜度為 O(nlogn) 的排序增量算法。
我們先對輸入的點按照順時針或逆時針進行極角排序,可以想象一開始為一個足夠大的矩形,按照順時針或逆時針的順序不斷切割已有的平面。求 n 個半平面的交就是用第 n 條表示當前半平面的直線去切割現有的面。每次切割都會產生一個更小的面,當所有直線都切割完畢后就是我們所需要的結果了。
那么,一個多邊形的向量應該長這樣(就是把邊當做直線):
具體的步驟是:
求出所有的向量,按照極角序排序。
(推薦使用atan2(y,x)表示(x,y)的旋轉角。精度較高,而且范圍是(-Pi,Pi]可以求出所有的旋轉角)
然后,增量法插入,用一個隊列q維護直線,另一個隊列p維護交點。
發現,新加入一個直線,情況如下:
對於新加入的藍色直線,其右側的p中的交點都要彈出。發現,彈出的點一定在隊列的兩端。
所以,p,q都是雙端隊列。
彈完了之后,再加入這一個交點。
注意,隊列中至少剩下一條直線,否則沒有意義。
還有一個注意情況:
最后可能出現這種情況:
這樣的話,要把多余的線頭彈掉。判斷方法是,如果隊尾交點在第一條直線的右側,彈掉。
還有一些其他注意事項:
1.如果兩個向量的旋轉角相同,那么保留左邊的那一個。
不刪除的話,可能導致兩個直線平行(共線),求交點的時候就除以零掛了。
模板:
(這個題數據有鍋,第一個點要特判。。。)
#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 */
求多邊形交
由於最終的區域,必須在所有多邊形的內部,即所有向量的內側。
所以直接求半平面交即可。
(當然多虧這是些凸多邊形,否則暴力半平面交顯然就錯了。而且我並不知道怎么做。。。)
#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=1000+5; const double eps=1e-8; struct po{ double x,y; po(){} po(double xx,double yy){ x=xx,y=yy; } }a[N]; po operator +(po a,po b){ return po(a.x+b.x,a.y+b.y); } po operator -(po a,po b){ return po(a.x-b.x,a.y-b.y); } po operator *(po a,double b){ return po(a.x*b,a.y*b); } po operator /(po a,double b){ return po(a.x/b,a.y/b); } double cross(po a,po b){ return a.x*b.y-a.y*b.x; } struct line{ po A,B; po V; double ang; line(){} line(po a,po b){ A=a,B=b; V=b-a; ang=atan2(b.y-a.y,b.x-a.x); } bool friend operator <(line a,line b){ if(a.ang!=b.ang) return a.ang<b.ang; return cross(b.A-a.A,a.V)>0; } }l[N]; bool Onleft(line a,po b){ return cross(a.V,b-a.A)>0; } po jiao(line a,line b){ po A=a.A,B=a.B,C=b.A,D=b.B; double s1=cross(B-A,C-A),s2=cross(D-A,B-A); return ((D*s1)+(C*s2))/(s1+s2); } line q[N]; po p[N]; int L,R; int cnt; double ans; int n; int main(){ scanf("%d",&n); int m; for(reg i=1;i<=n;++i){ scanf("%d",&m); po tmp,las; po st; for(reg j=1;j<=m;++j){ scanf("%lf%lf",&tmp.x,&tmp.y); if(j!=1) { //cout<<las.x<<" "<<las.y<<" || "<<tmp.x<<" "<<tmp.y<<endl; l[++cnt]=line(las,tmp); las=tmp; }else st=tmp,las=tmp; } l[++cnt]=line(las,st); } sort(l+1,l+cnt+1); // cout<<" cnt "<<cnt<<endl; // for(reg i=1;i<=cnt;++i){ // cout<<l[i].ang<<endl; // } n=1; for(reg i=2;i<=cnt;++i) if(l[i].ang!=l[i-1].ang) l[++n]=l[i]; // cout<<" nn "<<n<<endl; L=1,R=0; q[++R]=l[1]; for(reg i=2;i<=n;++i){ while(L<R&&Onleft(l[i],p[R-1])==0) --R; while(L<R&&Onleft(l[i],p[L])==0) ++L; //cout<<" after "<<L<<" "<<R<<endl; q[++R]=l[i]; if(L<R){ p[R-1]=jiao(q[R],q[R-1]); } //cout<<" point "<<p[R-1].x<<" "<<p[R-1].y<<endl; } while(L<R&&Onleft(q[L],p[R-1])==0) --R; if(L<R) p[R]=jiao(q[R],q[L]); if(R-L+1<3){ ans=0.00; }else{ for(reg i=L+1;i<=R-1;++i){ ans+=cross(p[i]-p[L],p[i+1]-p[L]); } ans=fabs(ans); ans/=2.0; } printf("%.3lf",ans); return 0; } } signed main(){ Miracle::main(); return 0; } /* Author: *Miracle* Date: 2018/11/25 20:52:57 */