题意简述
多测。每次给定一个\(N\)边形(保证相邻输入的顶点在多边形上也是邻接的),再给定一个以\((0,0)\)为圆心,半径为\(r\)的圆。
请计算出多边形和圆相交部分的面积(保留\(2\)位小数)。
- \(3\le N\le 50\)
- \(0.1\le r\le 1000\)
- \(-1000\le x_i,y_i\le 1000\)。
输入格式见原题面。
思路分析。
首先简化一下问题。我们先不考虑圆,从原点出发,对这个多边形进行三角剖分。
如图,依次连接各顶点,形成\(\vec{a},\vec{b},\vec{c},\vec{d},\vec{e},\vec{f},\vec{g}\)。接下来重复下面的操作:
- 连接相邻的两个向量\(\vec{x},\vec{y}\)形成一个三角形,如果\(\vec{x}\)到\(\vec{y}\)是顺时针,面积增加该三角形面积,否则减小该三角形面积(如图,红色部分是减去的,蓝色部分是增加的)。
我们发现这种方法可以很容易地处理各种多边形。如果有凹陷,被加的次数和被减的次数一定相等,相当于被抵消了。
如果再加上圆呢?那就把计算“三角形的面积”改成“三角形在圆内的面积”不就好了嘛。
计算“一个端点包括原点的三角形,在圆内的面积”,我们呢想到分类讨论三角形与圆的位置关系。
怎样优雅而有条理的讨论呢?我们可以用向量运算+代数的方式解决。
对于一个由向量\(\vec{OA}=(x_1,y_1),\vec{OB}=(x_2,y_2)\)组成的三角形\(AOB\),我们设\(\vec{AP}=t\times \vec{AB}\)。这样我们就可以用\(t\)表示出直线\(AB\)上的任意一点的坐标了:
\[\begin{aligned} \vec{OP}&=\vec{OA}+\vec{AP}\\ &=\vec{OA}+t\times \vec{AB}\\ &=\vec{OA}+t\times (\vec{OB}-\vec{OA})\\ &=(1-t)\times\vec{OA}+t\times \vec{OB}\\ &=(1-t)\times (x_1,y_1)+t\times (x_2,y_2)\\ &=((1-t)x_1+tx_2,(1-t)y_1+ty_2) \end{aligned}\]显然当\(t\in[0,1]\)时。\(P\)是在线段\(AB\)上的;\(t>1\)时在\(AB\)延长线上,\(t<0\)时在\(AB\)反向延长线上。
由此我们想到,如果求出让\(P\)落在圆上的\(t\),就可以用此时\(P\)和\(A,B\)的位置关系进行讨论了。
所以可以列一个关于\(t\)的一元二次方程:
\[\begin{aligned} x_P^2+y_P^2-r^2&=0\\ ((1-t)x_1+tx_2)^2+((1-t)y_1+ty_2)^2-r^2&=0\\ \end{aligned}\]化简并整理可得:
\[\begin{aligned} &(x_1^2+x_2^2-2x_1x_2+y_1^2+y_2^2-2y_1y_2)t^2\\ &+(2x_1x_2+2y_1y_2-2x_1^2-2y_1^2)t\\ &+(x_1^2+y_1^2-r^2)=0 \end{aligned}\\\]解这个方程。
然后根据根的情况分类:
- 无解:直线\(AB\)与圆无交点。答案就是图中扇形的面积。
- 有\(2\)个相同的解:直线\(AB\)与圆相切。答案仍然是图中扇形的面积。
- 有\(2\)个不同的解:设\(t_1,t_2\)分别表示较小和较大的解,进一步分类讨论:
- \(\textbf{Case 1,2:}\quad t_1>1,t_2>1\)或\(t_1<0,t_2<0\):
这两种情况可以放在一起考虑,都是\(P_1,P_2\)在线段\(AB\)外且同侧,即\(A,B\)在圆外且同侧,此时面积仍然是扇形。另一种情况的图就不放了。
- \(\textbf{Case 3:}\quad t_1\in [0,1],t_2>1\)
这种情况表示\(A\)在圆外\(B\)在圆内,即\(P_1\)在线段\(AB\)上而\(P_2\)不在。此时的面积是一个扇形+一个三角形。
- \(\textbf{Case 4:}\quad t_1<0,t_2\in [0,1]\)
和上面反过来,\(A\)在圆内,\(B\)在圆外,即\(P_2\)在线段\(AB\)上而\(P_1\)不在,同样是一个扇形+一个三角形。
- \(\textbf{Case 5:}\quad t_1<0,t_2>1\)
这说明\(A,B\)都在圆内部,即\(P_1,P_2\)都在线段\(AB\)外且不同侧,此时答案就是\(AOB\)的面积。
- \(\textbf{Case 6:}\quad t_1\in[0,1],t_2\in[0,1]\)
说明\(A,B\)都在圆外部且不同侧,即\(P_1,P_2\)都在线段\(AB\)上,此时答案是\(2\)个扇形+\(1\)个三角形。
- \(\textbf{Case 1,2:}\quad t_1>1,t_2>1\)或\(t_1<0,t_2<0\):
就酱。
计算两向量面积需要用叉积,看到有些题解用海伦公式计算。但海伦公式更适用于已知三边长度,这种只给出两个向量的情况,还需要用勾股定理计算三边,再套公式求解。代码相较叉积求面积较为冗长,且要多次求平方&开根,所以速度和精度都有所丢失。因此计算向量形成的面积还是推荐用叉积。
计算两向量夹角、计算两向量叉积、计算端点包含原点的三角形在园内的面积,都需要用两个向量作为参数提供。所以需要斟酌一下答案的正负性。一般都和叉积统一,即:如果\(\vec{a}\)在\(\vec{b}\)的顺时针是正数,逆时针则是负数。这样统一之后就不容易混淆了。别忘了计算后求一下绝对值,因为节点遍历顺序的正反决定结果的正负。
还有,保留\(2\)位小数别忘了。
Code
注:POJ编译器较老,代码需要经过一番修改才能通过编译,这里就不放修改前的代码了。
点击查看代码
``#include<bits/stdc++.h>
define Pi 3.1415926535897932384626
using namespace std;
struct vec{
double x,y;
vec operator+ (const vec &b){return {x+b.x,y+b.y};}
vec operator- (const vec &b){return {x-b.x,y-b.y};}
vec operator* (const double t){return {xt,yt};}
}fir,pre,cur;
double cross(vec a,vec b){return a.xb.y-b.xa.y;}
double R,ans;
int n;
double angle(vec a,vec b){//a相对b的角度([-pi,pi])
double theta=atan2(a.x,a.y)-atan2(b.x,b.y);
if(theta>Pi) theta-=2Pi;
if(theta<-Pi) theta+=2Pi;
return theta;
}//a在b的顺时针>0,逆时针<0,pi/-pi是反向,=0是同向
//扇形面积↓
double S_sect(double theta,double r){return thetarr/2;}
double S_tri_circ(vec OA,vec OB,double r){
double a=OA.xOA.x+OB.xOB.x-2OA.xOB.x+OA.yOA.y+OB.yOB.y-2OA.yOB.y;
double b=2OA.xOB.x+2OA.yOB.y-2OA.xOA.x-2OA.yOA.y;
double c=OA.xOA.x+OA.yOA.y-rr;
double delta=bb-4ac;
if(delta<=0) return S_sect(angle(OA,OB),r);
double t1=(-b-sqrt(delta))/2/a,t2=(-b+sqrt(delta))/2/a;//解方程得到t
if((t1<0&&t2<0)||(t1>1&&t2>1)) return S_sect(angle(OA,OB),r);//Case 1,2
if(t2>1&&t1>=0&&t1<=1){//Case 3
vec OP1=OA+(OB-OA)t1;
return cross(OP1,OB)/2+S_sect(angle(OA,OP1),r);
}
if(t1<0&&t2>=0&&t2<=1){//Case 4
vec OP2=OA+(OB-OA)t2;
return cross(OA,OP2)/2+S_sect(angle(OP2,OB),r);
}
if(t1<0&&t2>1) return cross(OA,OB)/2;//Case 5
vec OP1=OA+(OB-OA)t1,OP2=OA+(OB-OA)t2;//Case 6
return S_sect(angle(OA,OP1),r)+S_sect(angle(OP2,OB),r)+cross(OP1,OP2)/2;
}
int main(){
while(cin>>R>>n){
ans=0;
cin>>fir.x>>fir.y;
pre=fir;
while(--n){
cin>>cur.x>>cur.y;
ans+=S_tri_circ(pre,cur,R);
pre=cur;
}
ans+=S_tri_circ(cur,fir,R);
cout<<fixed<<setprecision(2)<<fabs(ans)<<"\n";
}
return 0;
}`
</details>