如何計算n個圓的聯集面積
前言
一般人第一次遇到這個問題，可能會想要想辦法用排容原理，找圓之間交疊的凸包之類的…。
然而我只要舉一個例子，你就會發現我們就算把凸包找出來了，我們也非常難知道找到的凸包到底是要減掉還是要加上這個面積。
來源
我們會發現包含中空的那個四邊形，我們不容易知道到底那塊要不要加上。
想法
我們可以利用線積分，對圓的聯集的邊界逆時鐘積分，就可以得到面積。
首先我們需要得到邊界，也就是沒有被別的圓包含的圓弧。要得到這個我們可以對於某個圓$c_i$去遍歷所有別的圓$c_j,i\neq j$，把交疊得到的圓弧記下來，那麼沒被記下來的圓弧就是我們要積分的部分了。
對於圓$a,b$，我們只需要考慮三種狀況($d=|a.center-b.center|$):
- same(a.radius+b.radius,d)
- d<=abs(a.radius-b.radius)+eps
- d<abs(a.radius+b.radius)-eps
接著我們可以利用$d,a.radius,b.radius$和餘弦定理來計算圓弧的角度範圍(以平行$x$軸為角度$0$)。
接著我們來看對於一個圓$c_i$，圓心在$(x_0,y_0)$，對於圓弧$[0,\theta]$積分所得到的值:
$\int_\Gamma 𝑥𝑑𝑦−𝑦𝑑𝑥=\int_\Gamma (𝑥_0+𝑟\cos\theta)\frac{dy}{d\theta}−(𝑦_0+𝑟\sin\theta)\frac{dx}{d\theta} d\theta$
$=\int_\Gamma (𝑥_0+𝑟\cos\theta)(r\cos\theta)−(𝑦_0+𝑟\sin\theta)(-r\sin\theta) d\theta$
$=𝑟[𝑥_0\sin\theta−𝑦_0\cos\theta+r\theta]_0^\theta=𝑟[𝑥_0\sin\theta−𝑦_0\cos\theta+𝑦_0+r\theta]$
最後我們直接上code，可以直接看CoverSegment和CircleUnionArea這兩個function。
程式碼:
const db eps=1e-8,pi=acos(-1);
bool same(db a,db b){return abs(a-b)<eps;}
db sq(db x){return x*x;}
struct P{
db x,y;
P():x(0),y(0){}
P(db x,db y):x(x),y(y){}
P operator+(P b){return P(x+b.x,y+b.y);}
P operator-(P b){return P(x-b.x,y-b.y);}
P operator*(db b){return P(x*b,y*b);}
P operator/(db b){return P(x/b,y/b);}
db operator*(P b){return x*b.x+y*b.y;}
db operator^(P b){return x*b.y-y*b.x;}
db abs(){return hypot(x,y);}
P unit(){return *this/abs();}
P spin(db o){
db c=cos(o),s=sin(o);
return P(c*x-s*y,s*x+c*y);
}
db angle(){return atan2(y,x);}
};
struct C{
P c;
db r;
C(P c=P(0,0),db r=0):c(c),r(r){}
};
vector<pair<db,db>> CoverSegment(C &a,C &b){
db d=(a.c-b.c).abs();
vector<pair<db,db>> res;
if(same(a.r+b.r,d));
else if(d<=abs(a.r-b.r)+eps){
if(a.r<b.r)res.pb({0,2*pi});
}else if(d<abs(a.r+b.r)-eps){
db o=acos((sq(a.r)+sq(d)-sq(b.r))/(2*a.r*d)),z=(b.c-a.c).angle();
if(z<0)z+=2*pi;
db l=z-o,r=z+o;
if(l<0)l+=2*pi;if(r>2*pi)r-=2*pi;
if(l>r)res.pb({l,2*pi}),res.pb({0,r});
else res.pb({l,r});
}return res;
}
db CircleUnionArea(vector<C> c){//circles should be distinct
db res=0;
rep(i,0,SZ(c)){
db w=0;vector<pair<db,db>> s={{2*pi,9}},z;
rep(j,0,SZ(c))if(i!=j){
z=CoverSegment(c[i],c[j]);
for(auto &e:z)s.pb(e);
}sort(all(s));
auto f=[&](db t){return c[i].r*(c[i].c.x*sin(t)-c[i].c.y*cos(t)+c[i].c.y+c[i].r*t);};
for(auto &e:s){
if(e.fi>w)res+=f(e.fi)-f(w);
w=max(w,e.se);
}
}return res/2;
}
const int _n=110;
db t,n,r,x,y;
vector<C> a;
vector<PII> ps;
main(void) {ios_base::sync_with_stdio(0);cin.tie(0);cout.tie(0);
cin>>n>>r;rep(i,0,n){cin>>x>>y;ps.pb({x,y});}
int nn=unique(all(ps))-ps.begin();rep(i,0,nn)a.pb(C({ps[i].fi,ps[i].se},r));
cout<<fixed<<setprecision(2)<<round(CircleUnionArea(a)*100)/100<<'\n';
return 0;
}
標頭、模板請點Codepad看
Codepad