[SPOJ CIRU]The area of the union of circles(自适应Simpson积分求圆并面积)

时间:2022-09-05 14:47:00

题目链接

http://acm.hust.edu.cn/vjudge/problem/viewProblem.action?id=16754

题目大意

求圆的面积并

思路

首先膜拜一发福大核武的神做法
http://hi.baidu.com/aekdycoin/item/b8ff6adc73c0e71dd78ed0d6
这种做法在applepi大爷的计算几何讲义中也有提到过,非常神奇优美,但是很不好写,毕竟不是所有人都像核武大爷那样神犇,没有很强的可推广性。

再来膜拜一发七爷的自适应Simpson积分题解
http://hi.baidu.com/sevenkplus/item/8cf47779e1e83b376cc37caa

自适应Simpson积分来求圆并的做法,就是用 f(x) 来表示一根穿过了横坐标为 x 的扫描线,被圆形覆盖的长度,那么 f(x) 的图像是在一段一段的区间中是连续的,并且函数的曲线肯定是在 x 轴的上方。而且显然最终的圆并面积就是 f(x) 的图像与 x 轴圈起来的面积,我们可以用数值积分来求这个图像的面积,数值积分分三种:矩形、梯形、Simpson,其中Simpson积分的做法就是把函数图像用一根二次函数曲线来拟合,精度最高。
而自适应Simpson积分,简单地说就是在当前区间直接拟合的精度足够的情况下,就不会继续递归到左右两半区间再做Simpson积分,提高了算法的效率。

看起来这个玩意就是乱搞嘛是不是= =,似乎看上去很不靠谱,但是实际上非常好用,除非一些特殊的反例(http://wapapp.baidu.com/shuxk/item/bd1210284b5d2e869c63d1b7),不光是圆的面积并,只要你会求出某种图形盖在扫描线上的区间,几乎任何玩意的面积并都是可以求的,实为杀人越货之利器23333。

那么我们就是要把所有圆分成一坨一坨圆,每一坨圆是连在一块的,对于每一坨圆,它们的 f(x) 图像肯定是连续的,就能用自适应Simpson积分把这一坨圆的面积求出来,累加答案即可。
当然还要注意怎么求一根扫描线被圆覆盖的长度,这和之前我做的下落的圆盘那题差不多,也是个比较简单的区间贪心问题,而找一坨一坨圆也是一样的贪心,只不过是在 x 轴上做贪心而已,怎么做相信大家都会,渣渣我就不说了。

嘴巴上讲讲做法当然很简单,但是我还是花了一个白天的时间才A掉这题,因此还要注意一些细节:
1、EPS开1e-6足矣,1e-7也可以,不过速度慢了4倍
2、此题很卡精度和细节,判 < EPS神马的必须取等号,比如说写成 < EPS就会WA,写成 <= EPS就没问题,不过在1e-7精度下没问题,坑爹啊
3、找区间并的长度的贪心很简单,但是还是要注意,不能掉以轻心(看我注释中的感叹号就是我WA的地方)
4、要做一个预处理,删掉所有被其他大圆完全覆盖的小圆,否则会WA,不要问我为什么,我也不知道,亲测WA,扒了别人代码才发现问题的。

代码


//自适应Simpson积分
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <cmath>

#define MAXN 10000
#define EPS 1e-6

using namespace std;

int n;

struct Point //点
{
double x,y;
Point(){}
Point(double _x,double _y):x(_x),y(_y){}
bool operator<(const Point &b)const { return x==b.x?y<b.y:x<b.x; }
};

struct Circle //圆
{
Point o;
double r;
Circle(){}
Circle(Point _o,double _r):o(_o),r(_r){}
}circles[MAXN],tmp[MAXN];

bool cmp2(Circle a,Circle b)
{
return a.r>b.r;
}

bool cmp(Circle a,Circle b)
{
if(fabs(a.o.x-a.r-b.o.x+b.r)<=EPS) return a.o.x+a.r<b.o.x+b.r;
return a.o.x-a.r<b.o.x-b.r;
}

pair<double,double>intervals[MAXN];
int tot=0;
int st,ed;

bool check(Circle a,Circle b) //检查b是否被套在a里面
{
return (a.o.x-b.o.x)*(a.o.x-b.o.x)+(a.o.y-b.o.y)*(a.o.y-b.o.y)<=(a.r-b.r)*(a.r-b.r);
}

void prework()
{
sort(circles+1,circles+n+1,cmp2);
int cnt=0;
for(int i=1;i<=n;i++)
{
bool flag=true;
for(int j=1;j<=cnt;j++)
if(check(tmp[j],circles[i]))
{
flag=false;
break;
}
if(flag) tmp[++cnt]=circles[i];
}
n=cnt;
for(int i=1;i<=n;i++)
circles[i]=tmp[i];
}

void getCircleIntersec(Circle a,double x) //求x=x这条平行于y轴的直线穿过圆a的长度
{
intervals[++tot]=make_pair(a.o.y-sqrt(a.r*a.r-(x-a.o.x)*(x-a.o.x)),a.o.y+sqrt(a.r*a.r-(x-a.o.x)*(x-a.o.x)));
}

double f(double x) //求扫描线x被圆覆盖部分的长度
{
tot=0;
for(int i=st;i<=ed;i++)
if(x<circles[i].o.x+circles[i].r&&x>circles[i].o.x-circles[i].r)
getCircleIntersec(circles[i],x);
sort(intervals+1,intervals+tot+1);
double ans=0,start=-1e12,end=-1e12; //!!!!!!!!!!!
for(int i=1;i<=tot;i++)
{
if(intervals[i].first>=end)
{
ans+=end-start;
start=intervals[i].first;
end=intervals[i].second;
}
else end=max(end,intervals[i].second);
}
ans+=end-start;
return ans;
}

double calc(double len,double fL,double fM,double fR) //求长度为len的[L,R]区间,中点为M的Simpson近似面积
{
return (fL+4*fM+fR)*len/6;
}

double Simpson(double L,double M,double R,double fL,double fM,double fR,double sqr) //Simpson积分求区间[L,R]的面积并,f(L)=L,f(R)=R,f(M)=M,把[L,R]当成整体来拟合得到的面积是sqr
{
double M1=(L+M)/2,M2=(M+R)/2;
double fM1=f(M1),fM2=f(M2);
double g1=calc(M-L,fL,fM1,fM),g2=calc(R-M,fM,fM2,fR);
if(fabs(sqr-g1-g2)<=EPS) //把当前区间分成2半再拟合得到的答案差别很小,就不再递归下去了
return g1+g2;
return Simpson(L,M1,M,fL,fM1,fM,g1)+Simpson(M,M2,R,fM,fM2,fR,g2);
}

int main()
{
double ans=0;
scanf("%d",&n);
for(int i=1;i<=n;i++)
scanf("%lf%lf%lf",&circles[i].o.x,&circles[i].o.y,&circles[i].r);
prework();
sort(circles+1,circles+n+1,cmp);
for(int i=1,j;i<=n;i++)
{
double L=circles[i].o.x-circles[i].r,R=circles[i].o.x+circles[i].r;
for(j=i+1;j<=n;j++)
{
if(circles[j].o.x-circles[j].r>R) break;
else R=max(R,circles[j].o.x+circles[j].r); //!!!!!!!!!!!!!!
}
double M=(L+R)/2;
st=i,ed=j-1;
i=j-1;
double fL=f(L),fM=f(M),fR=f(R);
ans+=Simpson(L,M,R,fL,fM,fR,calc(R-L,fL,fM,fR));
}
printf("%.3lf\n",ans);
return 0;
}