SPOJ 8073 CIRU - The area of the union of circles（圆面积并）

Description
给出n个圆的圆心坐标和半径，求这n个圆的面积并
Input
第一行一整数n表示圆的个数，之后n行每行三个整数x[i],y[i],r[i]表示第i个圆的圆心坐标和半径(1<=n<=1000,|x[i]|,|y[i]|,r[i]<=1000)
Output
输出这n个圆的面积并，结果保留小数点后三位
Sample Input
3
0 0 1
0 0 1
100 100 1
Sample Output
6.283
Solution
首先把被包含的小圆去掉，然后把独立的圆单独计算面积，对于剩下的圆，每个圆都和某些其他圆有相交部分，把一个相交部分用一个角度区间表示（假设逆时针从0~2PI，如果一个区间[a,b]跨过了2PI，就把该区间拆成[a,2PI]和[0,b]），这些角度区间可能还有交，所以把一个圆和其他圆相交部分排序然后合并有交的区间，然后就得到了若干不相交的区间，这些区间端点都是这个圆与其他圆的交点，在圆上找到这些交点并把相邻交点连起来后，一个圆就被分成了一个凸多边形和若干圆弧，凸多边形的面积显然要被统计到答案中，而对于圆弧部分，如果两个圆相交，那么这个相交部分那块圆弧其实是在另一个圆的凸多边形部分中的，所以这种圆弧不用统计面积，而对于非相交分布的圆弧是要统计到答案中的
Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
#define INF 0x3f3f3f3f
#define maxn 1111
#define eps 1e-8
#define PI acos(-1.0)
int sign(double x)
{
    if(fabs(x)<eps)return 0;
    return x>0?1:-1;
}
struct Point
{
    double x,y;
    Point(){};
    Point(double _x,double _y){x=_x,y=_y;}
    Point operator+(const Point &b)const{return Point(x+b.x,y+b.y);}
    Point operator-(const Point &b)const{return Point(x-b.x,y-b.y);}
    bool operator==(const Point &b)const{return sign(x-b.x)==0&&sign(y-b.y)==0;}
    double norm(){return sqrt(x*x+y*y);}
};
double det(Point a,Point b)
{
    return a.x*b.y-a.y*b.x;
}
double get_angle(Point a)
{
    return atan2(a.y,a.x);
}
struct Reg
{
    double st,ed;
    Reg(){};
    Reg(double _st,double _ed){st=_st,ed=_ed;}
    bool operator<(const Reg&b)const
    {
        return sign(st-b.st)<0||sign(st-b.st)==0&&sign(ed-b.ed)<0;
    }
};
struct Circle
{
    double r;//半径
    Point o;//圆心坐标
    Circle(){};
    Circle(double _r,Point _o){r=_r,o=_o;}
    vector<Reg>reg;//相交部分的区间表示
    Point get_point(double a){return Point(o.x+r*cos(a),o.y+r*sin(a));}
}C[maxn];
int Circle_Insert(Point o1,double r1,Point o2,double r2,Point *p)
{
    double d=(o1-o2).norm();
    if(sign(d)==0)
    {
        if(sign(r1-r2)==0)return -1;//重合
        return 0;//同心圆
    }
    if(sign(r1+r2-d)<0)return 0;//相离
    if(sign(fabs(r1-r2)-d)>0)return 0;//包含
    double ang1=atan2(o2.y-o1.y,o2.x-o1.x);
    double ang2=acos((r1*r1+d*d-r2*r2)/(2*r1*d));
    p[0]=Point(o1.x+r1*cos(ang1+ang2),o1.y+r1*sin(ang1+ang2));
    p[1]=Point(o1.x+r1*cos(ang1-ang2),o1.y+r1*sin(ang1-ang2));
    if(p[0]==p[1])return 1;//相切
    return 2;
}
bool mark[maxn];
double Solve(Circle *C,int n)
{
    double ans=0;
    memset(mark,0,sizeof(mark));
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            if(i!=j&&!mark[j])
            {
                double d=(C[i].o-C[j].o).norm();
                if(sign(d+C[i].r-C[j].r)<=0)
                {
                    mark[i]=1;
                    break;
                }
            }
    for(int i=1;i<=n;i++)
        if(!mark[i])
        {
            Point p[2];
            int flag=0;
            C[i].reg.clear();
            for(int j=1;j<=n;j++)
                if(i!=j&&!mark[j])
                {
                    int num=Circle_Insert(C[i].o,C[i].r,C[j].o,C[j].r,p);
                    if(num!=2)continue;
                    flag=1;
                    double st=get_angle(p[1]-C[i].o),ed=get_angle(p[0]-C[i].o);
                    if(sign(st)<0)st+=2*PI;
                    if(sign(ed)<0)ed+=2*PI;
                    if(sign(st-ed)>0)
                        C[i].reg.push_back(Reg(st,2*PI)),C[i].reg.push_back(Reg(0,ed));
                    else C[i].reg.push_back(Reg(st,ed));
                }
            if(!flag)//与其它圆不交的圆
            {
                ans+=PI*C[i].r*C[i].r;
                continue;
            }
            sort(C[i].reg.begin(),C[i].reg.end());
            int cnt=1;
            for(int j=1;j<C[i].reg.size();j++)
                if(sign(C[i].reg[cnt-1].ed-C[i].reg[j].st)>=0)
                    C[i].reg[cnt-1].ed=max(C[i].reg[cnt-1].ed,C[i].reg[j].ed);
                else C[i].reg[cnt++]=C[i].reg[j];
            C[i].reg.push_back(C[i].reg[0]);
            C[i].reg[cnt]=C[i].reg[0];
            for(int j=0;j<cnt;j++)
            {
                p[0]=C[i].get_point(C[i].reg[j].ed);
                p[1]=C[i].get_point(C[i].reg[j+1].st);
                ans+=0.5*det(p[0],p[1]);
                double ang=C[i].reg[j+1].st-C[i].reg[j].ed;
                if(sign(ang)<0)ang+=2*PI;
                ans+=0.5*C[i].r*C[i].r*(ang-sin(ang));
            }
        }
    return ans;
}
int main()
{
    int n;
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%lf%lf%lf",&C[i].o.x,&C[i].o.y,&C[i].r);
    printf("%.3f\n",Solve(C,n)+eps);
    return 0;
}