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];
(* 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]}
  • $\begingroup$ are the 'n` spheres themselves potentially self intersecting? If not this can readily be done analytically. $\endgroup$ – george2079 Apr 19 '14 at 14:05
  • $\begingroup$ @george2079 Judging from the output of Graphics3D[Sphere@@@coords[50]], the spheres may intersect. $\endgroup$ – Mark McClure Apr 19 '14 at 14:20
  • $\begingroup$ @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 ... $\endgroup$ – dakta Apr 19 '14 at 14:40
  • 1
    $\begingroup$ 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. $\endgroup$ – george2079 Apr 20 '14 at 17:01
  • $\begingroup$ @george2079 Sounds like a good idea. I will look into this in more detail ... $\endgroup$ – 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; 
     {i, 1, Length[centers]}];
test[phi_?NumericQ, theta_?NumericQ, 
   ptsRadii : {{{_, _, _}, _} ..}] := 
  ctest[phi, theta, Sequence @@ Transpose[ptsRadii]];
  Sin[phi] test[phi, theta, coords[50]], {phi, 0, Pi}, {theta, 0, 2 Pi},
  Method -> "AdaptiveMonteCarlo", PrecisionGoal -> 2] // AbsoluteTiming

(* Out: {4.800880, 1.72108} *)
|improve this answer|||||
  • 1
    $\begingroup$ 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. $\endgroup$ – george2079 Apr 21 '14 at 16:39
  • $\begingroup$ @george2079 I agree. He said he wanted fast, though, so fast I went. :) Feel free to edit. $\endgroup$ – Mark McClure Apr 21 '14 at 16:44
  • $\begingroup$ @MarkMcClure Thank you for your answer. It is fit for the intentended purpose. $\endgroup$ – dakta Apr 26 '14 at 9:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.