如何計算n個圓的聯集面積

如何計算n個圓的聯集面積

前言

一般人第一次遇到這個問題，可能會想要想辦法用排容原理，找圓之間交疊的凸包之類的...。
然而我只要舉一個例子，你就會發現我們就算把凸包找出來了，我們也非常難知道找到的凸包到底是要減掉還是要加上這個面積。

來源
我們會發現包含中空的那個四邊形，我們不容易知道到底那塊要不要加上。

@copyright petjelinux 版權所有
觀看更多正版原始文章請至petjelinux的blog

想法

我們可以利用線積分，對圓的聯集的邊界逆時鐘積分，就可以得到面積。
首先我們需要得到邊界，也就是沒有被別的圓包含的圓弧。要得到這個我們可以對於某個圓\(c_i\)去遍歷所有別的圓\(c_j,i\neq j\)，把交疊得到的圓弧記下來，那麼沒被記下來的圓弧就是我們要積分的部分了。
對於圓\(a,b\)，我們只需要考慮三種狀況(\(d=|a.center-b.center|\)):

1. same(a.radius+b.radius,d)
2. d<=abs(a.radius-b.radius)+eps
3. 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