圆的面积并

《圆的并》解题报告

1、题目描述①

给定n(1≤n≤1000)个圆，求n个圆的并的面积。圆的坐标和半径的范围是-10000到10000，答案精确到小数点后面6位。

2、算法分析

关于求圆的并的方法有很多。这中间有很多是求近似值的算法，例如随机算法，割平面算法，这些算法都很简单，当精度要求不高时，使用这些算法是比较好的选择。但是，如果精度要求高一点，这些算法的复杂度也会大大的提高。例如这里的数据范围，用随机算法时最坏情况下至少需要随机1010个点，割平面算法的复杂度也不止1000*1010的级别。

这道题，我们只有用理论上能求出准确值的算法（不考虑计算时的精度误差）。由于圆的并是一个很复杂的形状，直接计算它的面积不是很容易，想法当然是把它们分割成若干部分，每一个部分求出它们的面积，而每一部分的面积都不是很难算。这样就不难想到离散交点的方法，把所有的交点按照x坐标排序后，对任意两个交点，用一条条的竖直线去分割图形：

图1

①题目来源：经典问题

这样，平面被分成若干条形区域，怎样计算条形区域的面积呢？由于没有交点，所以这个变得简单。

看图中绿色方框框住的区域里面的圆，就不难发现，它可以分成若干梯形和弓形的面积（下面的左图）：

图2

考虑一般的情况，一个竖直区域的左右两条分界线，一个圆与它们相交的情况有3种：不与任意一条分界线相交、与其中一条相交，与两条都相交（相切不算相交）。其实如果把每个圆的两条竖直切线都拿来离散，那么就只需要考虑圆与两条竖直线都相交的情况了。

两条分界线切割圆时，在两条分割线中会形成两条弧，上面和下面各一条。这些弧之间不可能有交点，于是把它们的左的端点从上到下排序，左端点相同的按照右端点排序。记一个层次，然后再从上往下扫描，遇到上边界就将层次加1，遇到下边界就将层次减1。到了层次为0的时候，就形成了一个独立区域。这个区域的面积可以通过两个弓形的面积和一个梯形的面积相接计算出来。如图2的右，形成了两个独立区域，粉红色表示层次为1，绿色表示层次为2，蓝色表示层次为3。

再来分析算法的时间复杂度，最坏情况下，交点个数的级别是O(n 2)的，共分成了O(n 2)个区域，然后每一个区域都与n 个圆相交，加上排序，复杂度是O(n 3log 2n)的。

总感觉到题目描述很简单，应该复杂度降低一些。上面的算法很直观，而下面说的算法也比较直观，可以把复杂度降为O(n 2log 2n)。

试想，如果人来做此题，而不使用计算机计算，那么，人会采取什么样的方法呢？看下面的图：

两个圆相交的情况，可以把面积分成两个弓形的面积之和。（上面两个图）

下面是三个圆两两相交的情况，可以把面积看成是3个弓形与一个三角形面积之和。（下面两个图）

图3

看图4，4个圆，如上图所示，可以看成是8个弓形（外面4个，里面4个）的面积，再加上一个4边形面积（外面的4边形），再减去一个4边形面积（里面的4边形）。

上面的图形虽然简单，人能够很快的看出它是什么样的图形拼接而成，但是，如果圆的个数多了一点，还能不能变得如此简单？

从图3的3个圆的情况可以看出，有些交点是没有用到的！虽然3个圆有6个交点，但是，只有3个交点在计算中起到了作用。观察这3个点与另外3个点的区别。这3个圆的周围没有完全被蓝色部分包围！推广到一般的情况，也可以这样做。

首先要做一些预处理，如果一个圆完全被另一个圆包围，那么这一个圆可以删除。删除后，如果一个圆是孤立的圆，不与其它任何圆相交，就可以把这个圆的面积现算出来，也将它删去。剩下的工作就是求交点了。

对于一个圆来说，其他的某些圆覆盖了他的圆弧上的某一段，即若干个区间。那么所有的圆覆盖的部分也是若干个区间。如果有k 个区间被覆盖了，也就会有k 个区间没有被覆盖。这2k 个区间被2k 个点分开。

图中黑色部分是被覆盖的部分，即“看不见的”，而绿色部分是为覆盖的部分，即“露在外面的”，所以在计算弓形的面积时，绿色部分所围成的弓形一定会计算在总面积中间。所以我们用若干有向弦把这些部分分开，有向是指：沿圆的逆时针方向，如图5的右边所示。

对所有的圆都这样处理了以后，只看这些连线，有什么发现？不难发现任意

图5

一个分割点一定会有两条连线！一条连线连进来，一条连出去。并且也只会有两条连线。这个试着画就可以了，例如如果3个圆经过通过同一个点，那么其中一定有一个圆，这个点的两边的弧都被覆盖了，即这个点不是这个圆的分割点！

如果任意一个分割点都会有两条连线，那么，这些连线之间形成了若干多边形。这些多边形就是我们要求的多边形的面积。但是，这些多边形中有的面积是负的，这个只需要看这个多边形的连线是顺时针的还是逆时针的，顺的为正，逆的为负。

用一个例子来说明这个算法。

然后再把一个被完全包含的圆和一个孤立的圆特殊处理删掉后，就变成了右边的图。注意还有一个圆它没有被任意一个圆完全包围，但是它却没有圆弧露在外面。

算法就是这样，虽然没有严谨的证明来写，但是这样，感性的认识多于理性的思考，更容易使人理解。有些东西，千言万语也表达不出来，而有些东西，无法用言语表达，也不需要用言语表达。

该是分析时间复杂度的时候了。求每一个圆被覆盖的区间的复杂度是O(nlogn)，因为需要排序，所以求所有的圆的覆盖区间也只需要O(n 2logn)了。而算法其余的部分，例如前面的预处理，构造这些有向边，

以及计算多边形的面积。

都不会高于这个复杂度。不过似乎有个猜想，分割点的个数的级别是多少？上界是O(n2)，但似乎很难达到这个上界。

[参考文献]

《算法艺术与信息学竞赛》——刘汝佳黄亮著[讨论]

似乎讨论的结果都是前面一种方法。

[感谢]

周源、刘汝佳

3、程序

const

chash=12343;

inputfile='area.in';

outputfile='area.out';

zero=1e-8;

var

n,nodes:integer;

x,y,r:array[1..1000]of extended;

ans:extended;

inter:array[0..1000,1..2]of extended;

link,next:array[1..10000]of integer;

node:array[1..10000,1..2]of extended;

first:array[0..chash-1]of integer;

procedure init;

var

i:integer;

begin

assign(input,inputfile);reset(input);

readln(n);

for i:=1 to n do readln(x[i],y[i],r[i]);

close(input);

end;

function sqrdist(x1,y1,x2,y2:extended):extended;

begin

sqrdist:=sqr(x1-x2)+sqr(y1-y2);

end;

function dist(x1,y1,x2,y2:extended):extended;

begin

dist:=sqrt(sqr(x1-x2)+sqr(y1-y2))

end;

procedure prepare;

var

i,j:integer;

b:array[1..1000]of boolean;

begin

fillchar(b,sizeof(b),1);

for i:=1 to n do

for j:=i+1 to n do

if (x[i]=x[j])and(y[i]=y[j])and(r[i]=r[j]) then begin

b[i]:=false;break

end;

for i:=1 to n do

for j:=1 to n do if (i<>j)and(r[i]+zero

if dist(x[i],y[i],x[j],y[j])<=r[j]-r[i] then begin

b[i]:=false;break

end;

j:=0;

for i:=1 to n do if b[i] then begin

inc(j);x[j]:=x[i];y[j]:=y[i];r[j]:=r[i]

end;

n:=j

end;

function getangle(x,y:extended):extended;

begin

if x<-zero then getangle:=arctan(y/x)+pi

else if x>zero then

if y>0 then getangle:=arctan(y/x) else getangle:=arctan(y/x)+pi*2

else if y>0 then getangle:=pi/2 else getangle:=pi*3/2

end;

procedure getcross(i,j:integer;var t1,t2:extended);

var

a,b,c,a1,b1,c1,x1,y1,x2,y2,t,l:extended;

begin

a:=(x[i]-x[j])*2;b:=(y[i]-y[j])*2;

c:=sqr(r[j])-sqr(r[i])+sqr(x[i])-sqr(x[j])+sqr(y[i])-sqr(y[j]);

a1:=b;b1:=-a;

c1:=a1*x[i]+b1*y[i];

t:=a*b1-b*a1;

x1:=(c*b1-b*c1)/t;y1:=(a*c1-c*a1)/t;

l:=sqrt(sqr(r[i])-sqr(x1-x[i])-sqr(y1-y[i]));

t:=sqrt(sqr(a1)+sqr(b1));

x2:=x1+l*a1/t;y2:=y1+l*b1/t;

x1:=x1*2-x2;y1:=y1*2-y2;

t1:=getangle(x1-x[i],y1-y[i]);

t2:=getangle(x2-x[i],y2-y[i]);

if t2

t:=(t1+t2)/2;

if dist(x[j],y[j],x[i]+r[i]*cos(t),y[i]+r[i]*sin(t))>r[j] then begin t:=t1;t1:=t2;t2:=t;

if t2

end;

end;

procedure sort(l,r:integer);

var

i,j:integer;

k1,k2:extended;

begin

i:=l;j:=r;

k1:=inter[(l+r) shr 1,1];k2:=inter[(l+r) shr 1,2];

while i<=j do begin

while (inter[i,1]+zero

(abs(inter[i,1]-k1)k2+zero) do inc(i);

while (inter[j,1]>k1+zero)or

(abs(inter[j,1]-k1)

if i<=j then begin

inter[0]:=inter[i];inter[i]:=inter[j];inter[j]:=inter[0];

inc(i);dec(j)

end;

end;

if l

if i

end;

function getwhere(x,y:extended):integer;

var

i,t:integer;

begin

t:=trunc(abs(x+y+zero)*100000) mod chash;

i:=first[t];

while i<>0 do begin

if (abs(node[i,1]-x)

getwhere:=i;exit

end;

i:=next[i]

end;

inc(nodes);

node[nodes,1]:=x;node[nodes,2]:=y;

next[nodes]:=first[t];

first[t]:=nodes;

getwhere:=nodes

end;

function getchord(r,a:extended):extended;

begin

getchord:=sqr(r)/2*(a-sin(a))

end;

procedure getnode;

var

i,j,k,top,t1,t2:integer;

begin

nodes:=0;

for i:=1 to n do begin

top:=0;

for j:=1 to n do

if (i<>j)and (dist(x[i],y[i],x[j],y[j])+zero

inc(top);

getcross(i,j,inter[top,1],inter[top,2]);

end;

if top>0 then begin

sort(1,top);

k:=0;

for j:=1 to top do

if (k=0)or(inter[j,1]>inter[k,2]) then begin

inc(k);inter[k]:=inter[j]

end

else

if inter[j,2]>inter[k,2] then inter[k,2]:=inter[j,2];

top:=k;

while (top>0)and(inter[top,2]+zero>inter[1,1]+pi*2) do begin

if inter[top,1]-pi*2

inter[1,1]:=inter[top,1]-pi*2;

dec(top)

end;

if top>0 then begin

for j:=1 to top-1 do begin

ans:=ans+getchord(r[i],inter[j+1,1]-inter[j,2]);

t1:=getwhere(x[i]+r[i]*cos(inter[j+1,1]),y[i]+r[i]*sin(inter[j+1,1]));

t2:=getwhere(x[i]+r[i]*cos(inter[j,2]),y[i]+r[i]*sin(inter[j,2]));

link[t1]:=t2;

end;

ans:=ans+getchord(r[i],inter[1,1]+pi*2-inter[top,2]);

t1:=getwhere(x[i]+r[i]*cos(inter[1,1]),y[i]+r[i]*sin(inter[1,1]));

t2:=getwhere(x[i]+r[i]*cos(inter[top,2]),y[i]+r[i]*sin(inter[top,2]));

link[t1]:=t2;

end

end

else ans:=ans+pi*r[i]

end;

end;

procedure work;

var

i,j:integer;

visited:array[1..10000]of boolean;

begin

ans:=0;

getnode;

fillchar(visited,sizeof(visited),0);

for i:=1 to nodes do if not visited[i] then begin

j:=i;

repeat

visited[j]:=true;

ans:=ans+(node[link[j],1]*node[j,2]-node[j,1]*node[link[j],2])/2;

j:=link[j];

until j=i

end;

assign(output,outputfile);rewrite(output);

writeln(ans:0:6);

close(output)

end;

begin

init;

prepare;

work

end.