numerical integration - 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]}

Answer

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 : {{{_, _, _}, _} ..}] := 
  ctest[phi, theta, Sequence @@ Transpose[ptsRadii]];
NIntegrate[
  Sin[phi] test[phi, theta, coords[50]], {phi, 0, Pi}, {theta, 0, 2 Pi},

  Method -> "AdaptiveMonteCarlo", PrecisionGoal -> 2] // AbsoluteTiming

(* Out: {4.800880, 1.72108} *)