标签：return res 面积 BZOJ2178 Vector double Circle dis
前言
由于笔者被这道题目虐了很久，感觉心生不爽，所以写篇题解造福一下大众。希望别起到反效果就好了。
题解
这里的做法是计算直接算圆弧的积分。
首先比较坑的两个点（现在想想一点都不坑，只是自己菜）：
- 被包含的圆是直接不计算贡献的。
- 如果两圆重合，在排除被包含圆的时候可能会互相影响导致被直接删掉。
然后这里主要探讨的是实现的方法。
首先我们要得到两圆的关系，这个就直接比较 \(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;//a\in b
if(a.r>=dis+b.r) return 1;//b\in 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\n",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("---------------\n");
// printf("%lf %lf %lf\n",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\n",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\n",res),0;
}
/*
2
0 0 10
20 10 20
*/
标签：return,res,面积,BZOJ2178,Vector,double,Circle,dis 来源： https://www.cnblogs.com/Point-King/p/15194919.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。