1. 程式人生 > 實用技巧 >【洛谷5467】[PKUSC2018] PKUSC(計算幾何)

【洛谷5467】[PKUSC2018] PKUSC(計算幾何)

點此看題面

  • 給定\(n\)個點和一個\(m\)個點構成的多邊形。
  • 隨機將多邊形旋轉一個角度,求落在多邊形內部的點數的期望。
  • \(n\le200,m\le500\)

能自己把這道題的做法口胡出來,但具體實現中一些地方還是參考了一下題解。

不過計算幾何應該都是口胡容易程式碼難吧。。。

解題思路

考慮求多邊形內點數的期望,由於點與點之間是獨立的,所以這就相當於是求每個點落在多邊形內的概率之和。

而一個點落在多邊形內的概率,就是\([0,2\pi)\)內所有能使它落在多邊形內的弧度總數除以\(2\pi\)

轉多邊形實在是太麻煩了,因此我們改為考慮點的旋轉,而這實際上就是一個圓。

那麼我們的問題就變成了以原點為圓心、以該點到原點的距離為半徑畫一個圓,圓上在多邊形內部的弧對應的角度總和。

這隻要先求出圓與多邊形上每一條邊的交點,然後分別考慮每相鄰兩個交點之間的弧是否在多邊形內部即可,而這也等價於這條弧的中點是否在多邊形內部。

口胡到這裡其實已經結束了,接下來具體介紹幾種函式的實現。

求線段與圓的交點

注意這道題中的圓必然以原點為圓心。

我們先求出原點到線段的垂足,如果原點到該垂足的距離\(l\)已經超過半徑\(r\)說明二者相離,無交點。

否則,我們分別按該線段的兩種方向,給這個垂足加上一個長度為\(\sqrt{r^2-l^2}\)的向量,然後看看得到的這個新點是不是線上段上就好了。

判斷點是否在多邊形內

這倒是一個我之前就寫過的東西,但現在發現有更簡單的寫法。

考慮從一個點隨便向一個方向引出一條線,那麼與多邊形的邊相交次數的奇偶性就意味著它是否在多邊形中。

但是,一個交點同時存在於兩條邊中,所以當這條線經過交點的時候就可能會出錯。

以前我傻乎乎地直接向右引一條線,每次要討論是否與端點相交,有些麻煩。

現在發現,其實只要隨便向一個刁鑽的角度引一條線,就不用考慮端點的問題了。

程式碼:\(O(nm^2)\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 200
#define M 500
#define DB double
#define eps 1e-8
using namespace std;
int n,m,cnt,x[N+5],y[N+5];DB Pi=acos(-1);
I int dcmp(Con DB& x) {return fabs(x)<eps?0:(x<0?-1:1);}
struct P
{
	#define CP Con P&
	DB x,y,ang;I P() {x=y=ang=0;}I P(Con DB& a,Con DB& b):x(a),y(b){ang=0;}
	I P operator + (CP o) Con {return P(x+o.x,y+o.y);}
	I P operator - (CP o) Con {return P(x-o.x,y-o.y);}
	I P operator * (Con DB& o) Con {return P(x*o,y*o);}
	I P operator / (Con DB& o) Con {return P(x/o,y/o);}
	I bool operator == (CP o) Con {return !dcmp(ang-o.ang);}
	I bool operator < (CP o) Con {return dcmp(ang-o.ang)<0;}
	I DB operator * (CP o) Con {return x*o.x+y*o.y;}
	I DB operator ^ (CP o) Con {return x*o.y-y*o.x;}
	I DB Len() {return sqrt(x*x+y*y);}
}p[M+5],f[2*M+5];
struct S
{
	#define CS Con S&
	P A,B;I S() {A=B=P();}I S(CP a,CP b):A(a),B(b){}
};
I P Foot(CS s) {return P(-(s.B-s.A).y,(s.B-s.A).x);}//垂線上的一點
I bool P_on_S(CP p,CS s)//判斷這條直線上的一點是否線上段上
{
	if(dcmp(p.x-min(s.A.x,s.B.x))<0||dcmp(p.x-max(s.A.x,s.B.x))>0) return 0;
	if(dcmp(p.y-min(s.A.y,s.B.y))<0||dcmp(p.y-max(s.A.y,s.B.y))>0) return 0;return 1;
}
I P L_with_L(CS X,CS Y)//求兩條直線交點
{
	DB s1=(X.B-X.A)^(Y.A-X.A),s2=(X.B-X.A)^(Y.B-X.A);
	return Y.A+(Y.B-Y.A)*s1/(s1-s2);//利用面積比得出長度比
}
I void S_with_C(CS s,Con DB& r)//求線段與圓的交點
{
	P t=L_with_L(S(P(),Foot(s)),s);if(t.Len()>r) return;//求出垂足,如果距離超過r說明相離
	P g=(s.B-s.A)/(s.B-s.A).Len()*sqrt(r*r-t.Len()*t.Len());//求出用於變化的向量
	P_on_S(t+g,s)&&(f[++cnt]=t+g,0),P_on_S(t-g,s)&&(f[++cnt]=t-g,0);//兩種方向
}
I bool P_in_G(CP x)//判斷點是否在多邊形內
{
	RI i,res=0;P t;S s(x,P(x.x+0.413,x.y+0.521));//向一個刁鑽的角度引一條線
	for(i=1;i<=m;++i) t=L_with_L(s,S(p[i],p[i+1])),//求出與邊所在直線的交點
		res^=P_on_S(t,S(p[i],p[i+1]))&&dcmp(t.x-x.x)>=0;return res;//判斷是否與邊有交
}
I DB Calc(CI x,CI y)//求解點(x,y)的答案
{
	RI i;DB r=P(x,y).Len();for(cnt=0,i=1;i<=m;++i) S_with_C(S(p[i],p[i+1]),r);//求出圓與多邊形每條邊的交點
	if(!cnt) return P_in_G(P(r,0))?2*Pi:0;for(i=1;i<=cnt;++i) f[i].ang=atan2(f[i].x,f[i].y);//特判沒有交點,判斷是否完全在多邊形內部
	sort(f+1,f+cnt+1),cnt=unique(f+1,f+cnt+1)-f-1,(f[cnt+1]=f[1]).ang+=2*Pi;
	DB res=0;P t;for(i=1;i<=cnt;++i) t=Foot(S(f[i],f[i+1])),//檢驗相鄰兩交點之間的弧
		P_in_G(t/t.Len()*r)&&(res+=f[i+1].ang-f[i].ang);return res;//取中點判斷是否在多邊形內部
}
int main()
{
	RI i;for(scanf("%d%d",&n,&m),i=1;i<=n;++i) scanf("%d%d",x+i,y+i);
	for(i=1;i<=m;++i) scanf("%lf%lf",&p[i].x,&p[i].y);p[m+1]=p[1];
	DB ans=0;for(i=1;i<=n;++i) ans+=Calc(x[i],y[i]);return printf("%.5lf\n",ans/(2*Pi)),0;//求所有點概率和
}