SPOJ CIRU SPOJ VCIRCLE 圓的面積並問題

SPOJ VCIRCLE SPOJ CIRU

兩道題都是給出若干圓 就面積並，數據規模和精度要求不同。

求圓面積並有兩種常見的方法，一種是Simpson積分，另一種是幾何法。

在這里給出幾何方法。

PS.以下算法基於正方向為逆時針

 

考慮上圖中的藍色圓,綠色的圓和藍色的圓交於 A,B 2個交點 ,我們在逆時針系下考慮，那么 可以知道 對於藍色的圓，它對應的某個 角度區間被覆蓋了

假設 區間為 [A, B], 並且角度是按照 圓心到交點的 向量的 極角來定義 (為了方便，我一般都把角度轉化到 [0,2pi]區間內) 那么可以知道在這種 標識情況下，可能存在以下情況

這種區間的跨度如何解決呢?實際上十分簡單，只需要把[A,B] 拆成 [A, 2PI], [0,B]即可,也就是所謂的添加虛擬點

下面介紹一下 對於我們當前所求任務的實際運用( 利用上述做法)

首先 對於所給的N個圓，我們可以進行去冗雜，表現為:
(1) 去掉被包含(內含 or 內切)的小圓 ()
(2) 去掉相同的圓

枚舉一個圓，並對於剩下的圓和它求交點，對於所求的的交點，可以得到一個角度區間 [A,B], 當然區間如果跨越了(例如 [1.5PI, 0.5PI],注意這里是有方向的) 2PI那么需要拆 區間 

可以知道，最后區間的並集必然是最后 所有圓和當前圓的交集的一個邊界！

於是我們得到互補區間(所謂互補區間就是[0,2PI]內除去區間並的區間,可能有多個)

 

假設我們先枚舉了橙色的圓，那么得到了許多角度區間，可以知道綠色的和藍色的角度區間是“未被覆蓋的”，對於未被覆蓋的

圓弧染色！

而對於其他圓，我們也做同樣的步驟， 同時把他們未被覆蓋的角度區間的圓弧標記為黑色陰影

於是最終的結果就如下圖 (染色只考慮圓弧)

 

通過觀察不難發現,圓的並是許多的圓弧+ 許多多邊形的面積之和(注意這里為簡單多邊形,並且面積有正負之別！)

於是我們累加互補區間的圓弧面積到答案，並把互補區間確定的弦的有向面積累加到答案

(例如在上面的例子中，我們在掃描到橙色圓的這一步只需要累加藍色和綠色的圓弧面積 以及 藍色和綠色的有向面積，注意這里藍色和綠色的邊必然是最后那個多邊形的邊！)

這里涉及到一個問題，就是：
圓弧可能大於一半的圓，例如上圖中最大的圓，當然如果你推出了公式，那么實際上很容易發現無論圓弧長啥樣都能算出正確的答案！

 

代碼如下 

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<vector>
#include<queue>
#include<stack>
#include<map>
#include<algorithm>
#include<complex>
using namespace std;

const double EPS=1e-9,PI=acos(-1.0);

int cmp(double k)
{
 return k<-EPS ? -1:k>EPS ? 1:0;
}



inline double sqr(double x)
{
 return x*x;
}






struct point
{
 double x,y;
 point (){}
 point (double a,double b):x(a),y(b){}
 bool input()
 	{
 	 return scanf("%lf%lf",&x,&y)!=EOF;
	}
 friend point operator +(const point &a,const point &b)
 	{
 	 return point(a.x+b.x,a.y+b.y);
	}	
 friend point operator -(const point &a,const point &b)
 	{
 	 return point(a.x-b.x,a.y-b.y);
	}
 friend bool operator ==(const point &a,const point &b)
 	{
 	 return cmp(a.x-b.x)==0&&cmp(a.y-b.y)==0;
	}
 friend point operator *(const point &a,const double &b)
 	{
 	 return point(a.x*b,a.y*b);
	}
 friend point operator*(const double &a,const point &b)
 	{
 	 return point(a*b.x,a*b.y);
	}
 friend point operator /(const point &a,const double &b)
 	{
 	 return point(a.x/b,a.y/b);
	}
 double norm()
 	{
 	 return sqrt(sqr(x)+sqr(y));
	}
};


double cross(const point &a,const point &b)
{
 return a.x*b.y-a.y*b.x;
}

struct Circle
{
 point p;
 double r;
 bool operator <(const Circle &o)const
 	{
      if(cmp(r-o.r)!=0)return cmp(r-o.r)==-1;
      if(cmp(p.x-o.p.x)!=0)return cmp(p.x-o.p.x)==-1;
      return cmp(p.y-o.p.y)==-1;
 	}
 bool operator ==(const Circle &o)const
 	{
 	 return cmp(r-o.r)==0&&cmp(p.x-o.p.x)==0&&cmp(p.y-o.p.y)==0;
	}
};

point rotate(const point &p,double cost,double sint)
{
 double x=p.x,y=p.y;
 return point(x*cost-y*sint,x*sint+y*cost);
}

pair<point,point> crosspoint(point ap,double ar,point bp,double br)
{
 double d=(ap-bp).norm();
 double cost=(ar*ar+d*d-br*br)/(2*ar*d);
 double sint=sqrt(1.-cost*cost);
 point v=(bp-ap)/(bp-ap).norm()*ar;
 return make_pair(ap+rotate(v,cost,-sint),ap+rotate(v,cost,sint));
}

inline pair<point,point> crosspoint(const Circle &a,const Circle &b)
{
 return crosspoint(a.p,a.r,b.p,b.r);
}

const int maxn=2000;
struct Node
{
 point p;
 double a;
 int d;
 Node(const point &p,double a,int d):p(p),a(a),d(d){}
 bool operator <(const Node &o)const{
  return a<o.a;
 }
};

double arg(point p)
{
 return arg(complex<double>(p.x,p.y));
}

double solve(int m,Circle tc[],Circle c[])
{
 int n=0;
 sort(tc,tc+m);
 m=unique(tc,tc+m)-tc;//unique返回去重后的尾地址 
 for(int i=m-1;i>=0;i--)
 	{
     bool ok=true;
     for(int j=i+1;j<m;++j)
     	{
     	 double d=(tc[i].p-tc[j].p).norm();
     	 if(cmp(d-abs(tc[i].r-tc[j].r))<=0)
     	 	{
     	 	 ok=false;
     	 	 break;
			}
		}
	 if(ok)c[n++]=tc[i];
	}
 double ans=0;
 for(int i=0;i<n;++i)
 	{
 	 vector<Node> event;
 	 point boundary=c[i].p+point(-c[i].r,0);
 	 event.push_back(Node(boundary,-PI,0));
 	 event.push_back(Node(boundary,PI,0));
 	 for(int j=0;j<n;++j)
 	 	{
 	 	 if(i==j)continue;
 	 	 double d=(c[i].p-c[j].p).norm();
 	 	 if(cmp(d-(c[i].r+c[j].r))<0)
 	 	 	{
 	 	 	 pair<point,point> ret=crosspoint(c[i],c[j]);
 	 	 	 double x=arg(ret.first-c[i].p);
 	 	 	 double y=arg(ret.second-c[i].p);
 	 	 	 if(cmp(x-y)>0){
 	 	 	     event.push_back(Node(ret.first,x,1));
 	 	 	     event.push_back(Node(boundary,PI,-1));
 	 	 	     event.push_back(Node(boundary,-PI,1));
 	 	 	     event.push_back(Node(ret.second,y,-1));
				}else{
					event.push_back(Node(ret.first,x,1));
					event.push_back(Node(ret.second,y,-1));
				}
			}
		}
	 sort(event.begin(),event.end());
	 int sum=event[0].d;
	 for(int j=1;j<(int)event.size();++j)
	 	{
	 	 if(sum==0)
	 	 	{
	 	 	 ans+=cross(event[j-1].p,event[j].p)/2.;
	 	 	 double x=event[j-1].a;
	 	 	 double y=event[j].a;
	 	 	 double area=c[i].r*c[i].r*(y-x)/2;
	 	 	 point v1=event[j-1].p-c[i].p;
	 	 	 point v2=event[j].p-c[i].p;
	 	 	 area-=cross(v1,v2)/2.;
	 	 	 ans+=area;
			}
		 sum+=event[j].d;
		}	
	}
 return ans;
}

Circle c[maxn],tc[maxn];
int m;
int main()
{freopen("t.txt","r",stdin);
 scanf("%d",&m);
 for(int i=0;i<m;i++)
 	tc[i].p.input(),scanf("%lf",&tc[i].r);
 printf("%.5lf\n",solve(m,tc,c)+0.00000005);
 return 0;
}