# How I calculate the volume of multiple intersecting spheres?

I would like to know how to calculate the volume of the union of multiple possibly intersecting spheres (using Sphere[]). Please see the figure:

I was trying to get a mesh/delaunayMesh out of the spheres and apply "Volume" which didn't work.

Edit: Fixed spelling (Feb 6th).

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey Feb 5 '15 at 0:49
• Please include any code you have developed to try to solve this problem. Also, please reformat your code you have displayed according to the guidelines in meta1027. Doing so will encourage more people to consider your question. – bbgodfrey Feb 5 '15 at 0:50
• Do you want the volume of the intersection of the spheres, or the volume of the union of the spheres (which happen to be intersecting)? – Rahul Feb 5 '15 at 1:22
• The volume of the union of the spheres. – Derb Feb 5 '15 at 1:23
• Most of the solutions posted here work fine for 2 to 4 nodes (spheres). But any number of spheres which is higher brings up the following message:RegionMeasure::nmet: Unable to compute the measure of region RegionUnion[Ball[{100.,100.,100.},30.],Ball[{120.,120.,120.},30.],Ball[{130.,130.,130.},30.],Ball[{140.,140.,140.},30.],<<1>>,<<1>>,Ball[{190.,190.,190.},30.],Ball[{200.,200.,200.},30.],Ball[{210.,210.,210.},30.],Ball[{230.,230.,230.},30.]]. >> Any idea what it could be? – Derb Feb 6 '15 at 20:55

spheres = {
Sphere[{50, 50, 50}, 25],
Sphere[{70, 70, 70}, 25]};

rgn = RegionUnion @@
(spheres /. Sphere -> Ball);

RegionMeasure[rgn // N]


122585.

Volume[rgn // N]


122585.

EDIT: Example with more elements

rgn2 = RegionUnion[
Ball[{100., 100., 100.}, 30.], Ball[{120., 120., 120.}, 30.],
Ball[{130., 130., 130.}, 30.], Ball[{140., 140., 140.}, 30.],
Ball[{190., 190., 190.}, 30.], Ball[{200., 200., 200.}, 30.],
Ball[{210., 210., 210.}, 30.], Ball[{230., 230., 230.}, 30.]];

RegionMeasure[rgn2]


590768.

Volume[rgn2]


590768.

• This helped me! The calculation is fast and I can extend the solution to multiple spheres. Many thanks! – Derb Feb 5 '15 at 17:51
• Again, thank you, this helped me. It works nicely for up to 4 to 5 randomly distributed spheres. But any number of spheres which is higher brings up the following message:RegionMeasure::nmet: Unable to compute the measure of region RegionUnion[Ball[{100.,100.,100.},30.],Ball[{120.,120.,120.},30.],Ball[{130.,130.,130.},30.],Ball[{140.,140.,140.},30.],<<1>>,<<1>>,Ball[{190.,190.,190.},30.],Ball[{200.,200.,200.},30.],Ball[{210.,210.,210.},30.],Ball[{230.,230.,230.},30.]]. >> Any idea what it could be? I have this problem also for two more solution proposed here. – Derb Feb 6 '15 at 20:51
• Example above with eight balls worked without a problem. – Bob Hanlon Feb 7 '15 at 0:50
• Bob, yes, it worked for me too. I included more two balls and it didn't work anymore. Same error. I think it must be my system. Does it work for you with 10? rgn2 = RegionUnion[Ball[{100., 100., 100.}, 30.], Ball[{120., 120., 120.}, 30.], Ball[{130., 130., 130.}, 30.], Ball[{140., 140., 140.}, 30.], Ball[{190., 190., 190.}, 30.], Ball[{200., 200., 200.}, 30.], Ball[{210., 210., 210.}, 30.], Ball[{230., 230., 230.}, 30.], Ball[{135., 130., 130.}, 30.], Ball[{140., 140., 140.}, 30.], Ball[{195., 190., 190.}, 30.], Ball[{200., 200., 200.}, 30.]]; RegionMeasure[rgn2] – Derb Feb 7 '15 at 2:38
• Neither RegionMeasure nor Volume worked for nine or more Balls on my system. Perhaps the level of complexity causes the calculation to exceed some time limit. You might try subdividing the region, sum the volume of the subregions, and subtract the intersections of the subregions. – Bob Hanlon Feb 7 '15 at 3:52

Let the spheres have radius $r1$ and $r2$ and their centers be separated by distance $d$. There are four cases:

• $r1+r2 < d$ (separate spheres): $V = {4 \pi \over 3} (r1^3 + r2^3)$
• $r1 > r2 \wedge d + r2 < r1$ (sphere 2 within sphere 1): $V = {4 \pi \over 3} r1^3$
• $r2 > r1 \wedge d + r1 < r2$ (sphere 1 within sphere 2): $V = {4 \pi \over 3} r2^3$
• $r1 + r2 <d \wedge (d + r2 > r1 \vee d + r1 > r2)$ (partially intersecting spheres): $V = {4 \pi \over 3}(r1^3 + r2^3) - V_{cap1} - V_{cap2}$ (see below).

This last result comes from realizing that when two spheres partially intersect, we can define a plane through the circle defined by the spheres' intersecting surfaces. Then there are two "caps" that are "overcounted": the "cap" of sphere 1 within sphere 2, and the "cap" of sphere 2 within sphere 1. The cap of sphere 1 within sphere 2 has volume

$V_{cap1} = \int_{h1}^r \pi r^2(x) dx$ ,

where

$r^2(x) = r1^2 - x^2$.

Thus $V_{cap1} = \frac{1}{3} \pi (\text{h1}-\text{r1})^2 (\text{h1}+2 \text{r1})$.

Likewise, we have

$V_{cap2} = \frac{\pi (d-\text{h1}+2 \text{r2}) \left(d^2-2 d \text{r2}+\text{r1}^2-\text{r2}^2\right)^2}{12 d^2}$.

Here $h1$ is the distance from sphere 1's center to the plane, and likewise for sphere 2. Note that here $h1 + h2 = d$. We solve for $h1$ by $\sqrt{r1^2 - h1^2} = \sqrt{r2^2 - (d - h1)^2}$ then the total volume is the sum of the volumes of the individual spheres (${4 \pi \over 3}ri^3$) minus the two overcounted "caps" given by $h1$ and $h2 = d - h1$.

Of course this computes lightning fast.

• the analytic approach should be preferred ( for only two spheres in any case). You ought to post the mathematica code though.. – george2079 Feb 5 '15 at 20:20

Unfortunately as of 10.0.2, RegionIntersection is not implemented for MeshRegion nor BoundaryMeshRegion objects embedded in 3D. But you could use ImplicitRegion[] as follows:

r1 = ImplicitRegion[(x - 50)^2 + (y - 50)^2 + (z - 50)^2 <= 25^2, {x,
y, z}];
r2 = ImplicitRegion[(x - 70)^2 + (y - 70)^2 + (z - 70)^2 <= 25^2, {x,
y, z}];
Volume[RegionIntersection[r1, r2]]
RegionPlot3D[{r1, r2}, PlotPoints -> 50]

• Note that calls to Volume[] can sometimes take a long time! – M.R. Feb 5 '15 at 0:40
• I run the code. It really took long. :-( At the end, I got the following message: Unable to compute the volume of region. Any idea? – Derb Feb 5 '15 at 1:17
ball[{x0_, y0_, z0_}, r_] := (x - x0)^2 + (y - y0)^2 + (z - z0)^2 <= r^2.;
region = ImplicitRegion[ball[{50, 50, 50}, 25] && ball[{75, 75, 75}, 25], {x, y, z}];
Volume[region]

(*1683.46*)

• ball[{x0_, y0_, z0_}, r_] := (x - x0)^2 + (y - y0)^2 + (z - z0)^2 <= r^2.; region = ImplicitRegion[ ball[{50, 50, 50}, 25] && ball[{75, 75, 75}, 10], {x, y, z}]; Volume[region] (* 0 *) which is incorrect. You probably want NAND. – David G. Stork Feb 5 '15 at 2:41
• There is no intersection in your case. OP asked for the intersection not the union. – Ivan Feb 5 '15 at 3:49
 NIntegrate[
Boole[
Norm[ {x, y, z}  - {50, 50, 50}] < 25 ||
Norm[ {x, y, z}  - {70, 70, 70}] < 25  ] ,
{x, 0, 100}, {y, 0, 100}, {z, 0, 100}] // Chop


122585.