# Surface area of intersecting spheres

Given a sphere of radius 1 centered at the origin and $n$ spheres with radii $r_i$ centered at predefined coordinates, $c_i$, in space, I am after the surface area of the unit sphere that is not intersected by any of the surrounding spheres. E.g. given the coordinates

c = {{{1.2,0,0},1}, {{0,-1.7,0},1.2}, {{0.7,1,0},0.9}, {{-0.5,-0.5,-0.5},1}}


I am interested in the (potentially) visible blue surface area in the graphic generated by:

{Blue, Sphere[{0,0,0},1], Red, Sphere[#,#2]&@@@c} // Graphics3D


In principle, I can obtain this area by evaluating the following integral

fun = Function[{t,p},
Evaluate[ Times @@
(UnitStep[Total[({Sin[t]*Cos[p],Sin[t]*Sin[p],Cos[t]}-#)^2]-#2^2]& @@@ c)
]
]
NIntegrate[Evaluate[fun[t,p]*Sin[t]], {t,0,Pi}, {p,0,2Pi}] // AbsoluteTiming


However, this approach is slow (in particular for a large numbers of spheres), and has convergence issues.

In an attempt to speed up the integral evaluation, I have devised the following code snippet, which is based on a naive (non-adaptive) application of the trapezoidal rule:

n = 1000;
theta = N@Pi/(n-1)*Range[0,n-1];
phi = 2*N@Pi/(2*n-1)*Range[0,2*n-1];
int1 = ConstantArray[1.,n];
int1[[{1,-1}]] = 0.5;
int2 = ConstantArray[1.,2*n];
int2[[{1,-1}]]=0.5;
(* Calculating ptsT is slow. However, it could be pre-calculated once ... *)
ptsT=Transpose[Outer[{Sin[#]*Cos[#2],Sin[#]*Sin[#2],Cos[#]}&, theta, phi], {3,2,1}];

(
lv = Transpose @
Fold[#1*UnitStep[Total[(ptsT - #2[[1]])^2] - #2[[2]]^2] &,
ConstantArray[1., {2 n, n}], c];
int1.((lv.int2)*Sin[theta]) *Pi/(n - 1) * 2*Pi/(2 n - 1)
) // AbsoluteTiming


How can this calculation be sped up? I am also interested in increasing the number of surrounding spheres to approximately 50.

The following code snippet will generate $n$ spheres, which might be useful in comparing the performance of different approaches:

coords[n_] := Transpose@{
1.1 * Table[
With[{y = (2*i + 1)/n - 1, phi = i*(Pi*(3 - Sqrt[5]))},
{Cos[phi]*#, y, Sin[phi]*#} &[Sqrt[1 - y*y]]
],
{i, 0, n - 1}],
ConstantArray[0.3, n]}

• are the 'n spheres themselves potentially self intersecting? If not this can readily be done analytically. – george2079 Apr 19 '14 at 14:05
• @george2079 Judging from the output of Graphics3D[Sphere@@@coords[50]], the spheres may intersect. – Mark McClure Apr 19 '14 at 14:20
• @george2079 Yes, the spheres are self-intersecting. I think that even for the case of intersecting spheres an analytical solution is possible using a power diagram. However, I am currently after a numerical solution ... – dakta Apr 19 '14 at 14:40
• the best approach most likely involves generating a triangular surface mesh (I don't think any of NIntegrates methods will do that unfortunately.) Thats a bit involved to work up as an answer, but if you do a search there are a few meshing questions on this site. – george2079 Apr 20 '14 at 17:01
• @george2079 Sounds like a good idea. I will look into this in more detail ... – dakta Apr 26 '14 at 9:32

I don't know your speed or precision requirements but here's an approach that yields a low precision estimate to your 50 sphere problem in a few seconds. It's based on the fact that the surface area of a sphere can be computed via $$\int_0^{2\pi}\int_0^{\pi} \sin(\varphi) \, d\varphi \, d\theta.$$ We'll simply write a test function to determine when a point is close to one of the spheres and use this to restrict the domain of integration.

ctest = Compile[{{phi, _Real}, {theta, _Real},
{centers, _Real, 2}, {radii, _Real, 1}},
Module[{p3d, result},
result = 1.0;
p3d = {Cos[theta] Sin[phi], Sin[theta] Sin[phi], Cos[phi]};
Do[If[Norm[p3d - centers[[i]]] < radii[[i]], result = 0.0;
Break[]],
{i, 1, Length[centers]}];
result]];
test[phi_?NumericQ, theta_?NumericQ,
ptsRadii : {{{_, _, _}, _} ..}] :=

• for comparison Method -> "LocalAdaptive" crunches away for 2-3 minutes (without errors) and gives 1.69185` (1.8million evaluations). I'd be inclined to thinks that's more accurate, of course without an analytic solution thats just a hunch. – george2079 Apr 21 '14 at 16:39