BZOJ2178 圆的面积并

前言

由于笔者被这道题目虐了很久，感觉心生不爽，所以写篇题解造福一下大众。希望别起到反效果就好了。

题解

这里的做法是计算直接算圆弧的积分。

首先比较坑的两个点（现在想想一点都不坑，只是自己菜）：

• 被包含的圆是直接不计算贡献的。
• 如果两圆重合，在排除被包含圆的时候可能会互相影响导致被直接删掉。

然后这里主要探讨的是实现的方法。

首先我们要得到两圆的关系，这个就直接比较 (r_1,r_2,dis(O_1,O_2)) 就行了。

int Circle_Circle_Check(Circle a,Circle b){
	double dis=len(a.O-b.O);
	if(b.r>=dis+a.r) return -1;//ain b
	if(a.r>=dis+b.r) return 1;//bin a
	if(dis<a.r+b.r&&a.r<dis+b.r&&b.r<dis+a.r) return 2;
	return 0;
}

然后关于计算两圆的交点，这里有一个只需要一次开方操作的做法，精度还是蛮高的。

大体思想就是利用余弦定理算出角度，然后利用角度找到垂足，然后求一下垂直的向量就好了。

但是按照上面直接做是要开多次方的，然后我们整理一下式子，发现有几个开方是多余的。

pair<Point,Point> Circle_Circle_Get(Circle a,Circle b){
	pair<Point,Point> res;
	Vector v=b.O-a.O;double dis_2=len_2(v);
	double Cos=(a.r*a.r+dis_2-b.r*b.r)/(2*a.r*dis_2);
	Vector v1=v*a.r*Cos,v2=Rotate_90(v)*sqrt((a.r*a.r-len_2(v1))/dis_2);
	return res.first=a.O+v1-v2,res.second=a.O+v1+v2,res;
}

根据求到的交点我们再将点转换为以这个圆的圆心为出发点的向量，这样我们就可以利用 atan2 函数直接极角排序，然后一个扫描就可以求出哪些弧是并集的表面 。

for(int j=0;j<(int)a.size();++j){
    if(i==j) continue;
    int tmp=Circle_Circle_Check(a[i],a[j]);
    if(tmp==2){
        pair<Point,Point> res=Circle_Circle_Get(a[i],a[j]);
        Vector v1=res.first-a[i].O,v2=res.second-a[i].O;
        if(atan2(v1.y,v1.x)<=atan2(v2.y,v2.x)){
            bag.push_back(make_pair(v1,1));
            bag.push_back(make_pair(v2,-1));
        }
        else{
            bag.push_back(make_pair((Vector){-a[i].r,-1e-8},1));
            bag.push_back(make_pair(v2,-1)),bag.push_back(make_pair(v1,1));
            bag.push_back(make_pair((Vector){-a[i].r,0},-1));
        }
    }
}

然后我们就求弧的积分，转换为弓形加上有向梯形面积，然后弓形的面积又可以转变为扇形加上有向三角形的面积，然后利用叉积就可以计算了。

int tmp=0;bool flag=true;
Vector lst=(Vector){-a[i].r,-1e-8};
for(int j=0;j<(int)bag.size();++j){
    // printf("%lf %lf %d
",bag[j].first.x,bag[j].first.y,bag[j].second);
    if(!tmp&&flag){
        Point p1=a[i].O+lst,p2=a[i].O+bag[j].first;
        res+=(p1.x-p2.x)*(p1.y+p2.y)/2;
        double tmp1=atan2(lst.y,lst.x),tmp2=atan2(bag[j].first.y,bag[j].first.x);
        res+=((tmp2-tmp1)*a[i].r*a[i].r+(bag[j].first^lst))/2;
        flag=false;
    }
    tmp+=bag[j].second;
    if(j+1<(int)bag.size()){
        double tmp1=atan2(bag[j].first.y,bag[j].first.x);
        double tmp2=atan2(bag[j+1].first.y,bag[j+1].first.x);
        if(tmp1==tmp2) continue;
    }
    if(!tmp) lst=bag[j].first,flag=true;
}

完整代码

#include<bits/stdc++.h>
using namespace std;
const int N=1e3+5;
const double Pi=acos(-1);
int n;double res=0;
struct Point{double x,y;};
struct Vector{double x,y;};
struct Circle{Point O;double r;}b[N];vector<Circle> a;
Vector operator + (Vector a,Vector b){return (Vector){a.x+b.x,a.y+b.y};}
Vector operator - (Point a,Point b){return (Vector){a.x-b.x,a.y-b.y};}
Vector operator * (Vector a,double b){return (Vector){a.x*b,a.y*b};}
Vector operator / (Vector a,double b){return (Vector){a.x/b,a.y/b};}
Point operator + (Point a,Vector b){return (Point){a.x+b.x,a.y+b.y};}
Point operator - (Point a,Vector b){return (Point){a.x-b.x,a.y-b.y};}
double operator ^ (Vector a,Vector b){return a.x*b.y-a.y*b.x;}
Vector Rotate_90(Vector a){return (Vector){-a.y,a.x};}
double len_2(Vector a){return a.x*a.x+a.y*a.y;}
double len(Vector a){return sqrt(len_2(a));}
bool cmp(pair<Vector,int> a,pair<Vector,int> b){
	return atan2(a.first.y,a.first.x)<atan2(b.first.y,b.first.x);
}
int Circle_Circle_Check(Circle a,Circle b){
	double dis=len(a.O-b.O);
	if(b.r>=dis+a.r) return -1;
	if(a.r>=dis+b.r) return 1;
	if(dis<a.r+b.r&&a.r<dis+b.r&&b.r<dis+a.r) return 2;
	return 0;
}
pair<Point,Point> Circle_Circle_Get(Circle a,Circle b){
	pair<Point,Point> res;
	Vector v=b.O-a.O;double dis_2=len_2(v);
	double Cos=(a.r*a.r+dis_2-b.r*b.r)/(2*a.r*dis_2);
	Vector v1=v*a.r*Cos,v2=Rotate_90(v)*sqrt((a.r*a.r-len_2(v1))/dis_2);
	return res.first=a.O+v1-v2,res.second=a.O+v1+v2,res;
}
int main(){
	cin>>n;
	for(int i=1;i<=n;++i){
		scanf("%lf%lf%lf",&b[i].O.x,&b[i].O.y,&b[i].r);
	}
	for(int i=1;i<=n;++i){
		bool flag=true;
		for(int j=1;j<=n;++j){
			if(i==j) continue;
			if(Circle_Circle_Check(b[i],b[j])<0){
				if(Circle_Circle_Check(b[j],b[i])<0&&i<j) continue;
				flag=false;break;
			}
		}
		if(flag) a.push_back(b[i]);
	}
	for(int i=0;i<(int)a.size();++i){
		// printf("---------------
");
		// printf("%lf %lf %lf
",a[i].O.x,a[i].O.y,a[i].r);
		vector<pair<Vector,int> > bag;
		for(int j=0;j<(int)a.size();++j){
			if(i==j) continue;
			int tmp=Circle_Circle_Check(a[i],a[j]);
			if(tmp==2){
				pair<Point,Point> res=Circle_Circle_Get(a[i],a[j]);
				Vector v1=res.first-a[i].O,v2=res.second-a[i].O;
				if(atan2(v1.y,v1.x)<=atan2(v2.y,v2.x)){
					bag.push_back(make_pair(v1,1));
					bag.push_back(make_pair(v2,-1));
				}
				else{
					bag.push_back(make_pair((Vector){-a[i].r,-1e-8},1));
					bag.push_back(make_pair(v2,-1)),bag.push_back(make_pair(v1,1));
					bag.push_back(make_pair((Vector){-a[i].r,0},-1));
				}
			}
		}
		sort(bag.begin(),bag.end(),cmp);
		int tmp=0;bool flag=true;
		Vector lst=(Vector){-a[i].r,-1e-8};
		for(int j=0;j<(int)bag.size();++j){
			// printf("%lf %lf %d
",bag[j].first.x,bag[j].first.y,bag[j].second);
			if(!tmp&&flag){
				Point p1=a[i].O+lst,p2=a[i].O+bag[j].first;
				res+=(p1.x-p2.x)*(p1.y+p2.y)/2;
				double tmp1=atan2(lst.y,lst.x),tmp2=atan2(bag[j].first.y,bag[j].first.x);
				res+=((tmp2-tmp1)*a[i].r*a[i].r+(bag[j].first^lst))/2;
				flag=false;
			}
			tmp+=bag[j].second;
			if(j+1<(int)bag.size()){
				double tmp1=atan2(bag[j].first.y,bag[j].first.x);
				double tmp2=atan2(bag[j+1].first.y,bag[j+1].first.x);
				if(tmp1==tmp2) continue;
			}
			if(!tmp) lst=bag[j].first,flag=true;
		}
		Point p1=a[i].O+lst,p2=a[i].O+(Vector){-a[i].r,0};
		res+=(p1.x-p2.x)*(p1.y+p2.y)/2;
		double tmp1=atan2(lst.y,lst.x),tmp2=atan2(0,-a[i].r);
		res+=((tmp2-tmp1)*a[i].r*a[i].r+((Vector){-a[i].r,0}^lst))/2;
	}
	return printf("%.3lf
",res),0;
}
/*
2
0 0 10
20 10 20
*/