[题解]POJ3675 Telescope——求多边形与圆相交部分的面积
题意简述
多测。每次给定一个\(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\)上的任意一点的坐标了：
显然当\(t\in[0,1]\)时。\(P\)是在线段\(AB\)上的；\(t>1\)时在\(AB\)延长线上，\(t<0\)时在\(AB\)反向延长线上。
由此我们想到，如果求出让\(P\)落在圆上的\(t\)，就可以用此时\(P\)和\(A,B\)的位置关系进行讨论了。
所以可以列一个关于\(t\)的一元二次方程：
化简并整理可得：
解这个方程。
然后根据根的情况分类：
- 无解：直线\(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 {x*t,y*t};}
}fir,pre,cur;
double cross(vec a,vec b){return a.x*b.y-b.x*a.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-=2*Pi;
if(theta<-Pi) theta+=2*Pi;
return theta;
}//a在b的顺时针>0，逆时针<0，pi/-pi是反向，=0是同向
//扇形面积↓
double S_sect(double theta,double r){return theta*r*r/2;}
double S_tri_circ(vec OA,vec OB,double r){
double a=OA.x*OA.x+OB.x*OB.x-2*OA.x*OB.x+OA.y*OA.y+OB.y*OB.y-2*OA.y*OB.y;
double b=2*OA.x*OB.x+2*OA.y*OB.y-2*OA.x*OA.x-2*OA.y*OA.y;
double c=OA.x*OA.x+OA.y*OA.y-r*r;
double delta=b*b-4*a*c;
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;
}
