Volume and Surface Area of a Compound Geometric Object
Show older comments
I am having difficulty finding the non-overlapping volume and surface area of a solid comprised of three spheres. Specifically, I am considering three unit spheres whose centers lie on the corners of an equilateral triangle of side length sqrt(2). The tricky part is the triple overlap in the center. There are some papers in the literature that provide closed-form solutions, but I haven't found them to be consistent with one another (maybe I'm doing something wrong), so I don't know which is correct. So, I thought that I would do it numerically with masks and an alpha shape (hacky code below), but I am failing miserably. I have never used alpha shapes before, so I might just have a fundamental misunderstanding of the very basics. Any suggestions, please? Maybe there's a better way? Thanks, in advance.
% Create a "solid" object from three unit-radius spheres whose centers lie
% on the corners of an equilateral triangle of side length sqrt(2). Compute
% the volume and surface area of the resulting solid.
%
% Expected solution: V ~ 11.2 L^3, SA ~ 28.1 L^2
% Preliminaries
close all;
clear all;
clc;
% Sphere centroids and radii
c1 = [0 0 0];
c2 = [sqrt(2) 0 0];
c3 = [sqrt(2)/2 sqrt(3/2) 0];
radius = 1;
% Generate high-resolution sphere
[X1,Y1,Z1] = sphere(200);
X1 = X1(:);
Y1 = Y1(:);
Z1 = Z1(:);
% Create points for two other spheres
X2 = X1 + c2(1);
X3 = X1 + c3(1);
Y3 = Y1 + c3(2);
% Create masks to eliminate interior points
m12 = sum(([X1 Y1 Z1]-c2).^2,2)<radius^2;
m13 = sum(([X1 Y1 Z1]-c3).^2,2)<radius^2;
mm1 = logical(m12+m13);
m21 = sum(([X2 Y1 Z1]-c1).^2,2)<radius^2;
m23 = sum(([X2 Y1 Z1]-c3).^2,2)<radius^2;
mm2 = logical(m21+m23);
m31 = sum(([X3 Y3 Z1]-c1).^2,2)<radius^2;
m32 = sum(([X3 Y3 Z1]-c2).^2,2)<radius^2;
mm3 = logical(m31+m32);
% Apply masks
P1 = [X1 Y1 Z1];
P1 = P1(~mm1,:);
P2 = [X2 Y1 Z1];
P2 = P2(~mm2,:);
P3 = [X3 Y3 Z1];
P3 = P3(~mm3,:);
% Array holding all exterior points of the object
P = [P1; P2; P3];
P = unique(P,"rows");
% Plot and view, just to be sure
figure;
plot3(P(:,1),P(:,2),P(:,3),".")
axis equal
% Make alpha shape, have a look
figure;
shp = alphaShape(P(:,1),P(:,2),P(:,3));
pc = criticalAlpha(shp,"one-region");
shp.Alpha = pc;
plot(shp)
axis equal
% Compute volume and surface area of the solid
shp.volume
shp.surfaceArea
4 Comments
Moe Szyslak
on 27 May 2026
Torsten
on 28 May 2026
Determining the volume with Monto Carlo simulation should be no problem.
Concerning the surface area: The total surface area of the three spheres (with overlapping parts) is 3 * 4*pi*R^2. To determine the "outer" surface area of the three overlapping spheres, first choose randomly which of the three spheres you want to choose. Then generate a random point on the surface of this sphere. If this random point is in the interior of one of the other two spheres, reject it. Else add 1 to a variable "counter". The outer surface area A of the three intersecting spheres will then be approximately A = counter/total number of trials * 3 * 4*pi*R^2.
If you don't know, you can google on how to generate uniformly distributed random points on the surface of a sphere.
Moe Szyslak
on 28 May 2026
Accepted Answer
Code for volume and surface area. You might want to try to vectorize the for-loops.
c1 = [0 0 0];
c2 = [sqrt(2) 0 0];
c3 = [sqrt(2)/2 sqrt(3/2) 0];
R = 1;
N = 1e8;
% Volume
counter = 0;
xmin = min([c1(1),c2(1),c3(1)]) - R;
xmax = max([c1(1),c2(1),c3(1)]) + R;
ymin = min([c1(2),c2(2),c3(2)]) - R;
ymax = max([c1(2),c2(2),c3(2)]) + R;
zmin = min([c1(3),c2(3),c3(3)]) - R;
zmax = max([c1(3),c2(3),c3(3)]) + R;
Vref = (xmax-xmin)*(ymax-ymin)*(zmax-zmin);
for i = 1:N
x = xmin + rand*(xmax-xmin);
y = ymin + rand*(ymax-ymin);
z = zmin + rand*(zmax-zmin);
point = [x y z];
if norm(point-c1) <= R || norm(point-c2) <= R || norm(point-c3) <= R
counter = counter + 1;
end
end
V = counter/N * Vref
% Surface area
counter = 0;
for i = 1:N
D1 = randi(3);
point = randn(1,3);
D2 = point./sqrt(sum(point.^2));
if D1 == 1
D2 = D2 + c1;
if norm(D2-c2) < R || norm(D2-c3) < R
continue
end
counter = counter + 1;
elseif D1 == 2
D2 = D2 + c2;
if norm(D2-c1) < R || norm(D2-c3) < R
continue
end
counter = counter + 1;
else
D2 = D2 + c3;
if norm(D2-c1) < R || norm(D2-c2) < R
continue
end
counter = counter + 1;
end
end
A = counter/N * 3 * 4*pi*R^2
7 Comments
Moe Szyslak
on 28 May 2026
Torsten
on 28 May 2026
Can you show your analytical formula for the volume ? I think I would not be able to compute it by integration.
Moe Szyslak
on 28 May 2026
Edited: Moe Szyslak
on 28 May 2026
It's something of a kluge. I start with three spheres then subtract six lenses -- that part, of course, is straightforward. The trickier part (trigonometrically) is adding the triple overlap volume back in. For that, I use the equations presented in this paper. The equations simplify greatly since all spheres are of the same radius and separation distance, but the end result is still somewhat gnarly:
where W = 1/3 of the triple overlap, R = radius (1), D = separation distance (sqrt(2)), and \theta = angle between centroids (\pi/3). I think this should hold for any R<D<sqrt(3) as long as the symmetry isn't changed. There are other treatments of this problem in the literature (e.g., see here and here), but they aren't always obvious to find. The second of those references results in the following formula, which is somewhat simpler, but has the R=1 condition "baked in":
Here V is the total triple overlap volume and a = sqrt(2) is the separation distance (D in the first eqn). Ultimately, I want an approach that works on other compound sphere solids as well, so I am interested in finding a general approach to this class of problem, a task with which you have helped immensely.
BTW -- I did vectorize the volume calculation. On my laptop the for loop was taking ~26 s, parfor was at ~6 s, and the vectorized version is ~5 s. I did not do any robust speed-checking, just back-of-the-envelope type stuff. I haven't gotten around to vectorizing the SA calc yet. Thanks again.
John D'Errico
on 28 May 2026
Ah, reading your comment, I see your comment about the triply overlapped region.
A last comment about my approach to get the surface area:
If the three parts of the outer surface don't have the same area as in the above case, choosing them equally likely with D1 = randi(3) will give better results in smaller area parts than in larger area parts of the outer surface because the relationship "partial surface area / number of points used for resolution" is smaller. This could distort the total result. You could handle each part of the outer surface separately without randomization, but even then, you should give different values of N for different parts of the outer surface depending on their size. Maybe preestimating the area of the different parts of the outer surface by using a small value for N for each of them first and then doing the actual computation with N-values scaled by this result for the different area parts could be a way to go.
Paul
on 29 May 2026
Is the Monte Carlo approach necessary?
For computing the surface area would it be valid instead to just place N_i uniform points on the i'th sphere, then reject points that are inside the other spheres and then compute the proportion? Maybe the N_i are selected proportional to the surface areas of the indivdiual spheres. Is this what you meant by "You could handle each part of the outer surface separately without randomization, but even then, you should give different values of N for different parts of the outer surface depending on their size."
Similarly for the volume, would it be valid instead to just use a uniform grid of points within the bounding cube?
I didn't know about the randn approach to picking random points on the sphere, so thanks for showing that.
Is this what you meant by "You could handle each part of the outer surface separately without randomization, but even then, you should give different values of N for different parts of the outer surface depending on their size."
I meant to keep random generation of points on the surface of the spheres, but to cancel the random choice of one of the three spheres by the D1 = randi(3) command. Instead, I suggest looping over the three spheres in question one by one (which is not necessary in this simple case since the outer surfaces of the three spheres are equal by symmetry).
For computing the surface area would it be valid instead to just place N_i uniform points on the i'th sphere, then reject points that are inside the other spheres and then compute the proportion? Similarly for the volume, would it be valid instead to just use a uniform grid of points within the bounding cube?
I'm not an expert in this field, but intuitively I think that using uniform grids to compute surface area and volume needs more evaluation points to get the same accuracy compared to the Monte Carlo approach.
Here is a more qualified answer:
It stresses the advantages of Monte Carlo over uniform gridding especially for higher-dimensional problems.
More Answers (1)
John D'Errico
on 27 May 2026
Edited: John D'Errico
on 28 May 2026
First, alpha shapes tend to be poor at representing volumes with sharp internal corners, like those between the intersection of two adjacent spheres.
Second, there is no need to drop points that lie in the intersection. The alpha shape handles that directly. So much of what you did was wasted.
% Sphere centroids and radii
c1 = [0 0 0];
c2 = [sqrt(2) 0 0];
c3 = [sqrt(2)/2 sqrt(3/2) 0];
radius = 1;
But more importantly, we need to look at just one sphere. Make sure we can walk before we try to run. We need to understand what is happening.
[X,Y,Z] = sphere(200);
A unit sphere, the volume should be approximately 4/3*pi, surface area 4*pi. See how well it does.
XYZ = [X(:),Y(:),Z(:)];
S = alphaShape(XYZ,inf)
An alpha radius of 1 should be quite reasonable, certainly for a nice convex domain, but I can use inf. This will give us the delaunay triangulation.
volume(S)
That seems about right, as it should be a little of an underestimate of 4/3*pi. For this many points, we are only low by a tiny fraction, around 0.02%.
4/3*pi
volume(S)/(4/3*pi)
I'm happy.
surfaceArea(S)
surfaceArea(S)/(4*pi)
I was a little worried about the surface area, that possibly the convex domain of the alpha shape with inf as the alpha radius might be a little crinkly, sort of a fractal surface.
However, what happens when you did as you did?
a = criticalAlpha(S,'one-region')
S.Alpha = a
surfaceArea(S)
volume(S)
And what happened? It turns out that the resulting alpha shape is a bit holy, with essentially a zero volume, and twice the expected surface area. The result is one region, BUT it is garbage. You are seeing an area that is the sum of the outside of a sphere, plus the inside of that sphere. Anyway, I think we should understand what happened in the above.
Howabout the case of three spheres? Again, be very careful, as alpha shapes suck at approximating those sharp interior corners where the spheres will intersect.
XYZ3 = [XYZ + c1;XYZ + c2; XYZ + c3];
S3 = alphaShape(XYZ3,inf); % again, this is equivalent to the delaunay triangulation
Again, there is no need to delete the internal points. In fact, that helps you a little, as now the critical alpha ends up a little smaller. But even so, using the critical alpha is a poor choice. You will end up with a holy mess again.
criticalAlpha(S3,'one-region')
S3.Alpha = 0.75;
volume(S3)
surfaceArea(S3)
The volume was a little high, as an alpha radius of 0.75 is still large to approximate those intersecting corners, but it is sufficiently large that we don't get an internal hole.
The issue is, we will find it difficult to do tremendously better using an alpha shape. Again, alpha shapes are poor at handing those sharp edges between a pair of spheres. Monte Carlo is one solution. For example...
box = [min(XYZ3);max(XYZ3)];
n = 1e8;
MCsample = rand(n,3).*diff(box,[],1) + box(1,:);
rownorm = @(X) sqrt(sum(X.^2,2));
MCvol = prod(diff(box,[],1))/n*sum(rownorm(MCsample - c1) <= 1 | rownorm(MCsample - c2) <= 1 | rownorm(MCsample - c3) <= 1)
That was not very good, but Monte Carlo is known for being easy to do, yet not very good in terms of accuracy. The error goes down only with sqrt(n). How to get the surface area? That will be more difficult, but how would I solve the problem? I'd probably be lazy, and ressurect an old set of MATLAB tools for working with simplicial complexes, written many years ago for solving all sorts of such problems in multiple dimensions on color science problems. The idea would be to start with the tiling of one sphere in three dimensions. The result would be a 2-manifold of triangles, in a 3-dimensinoal cartesian space. Think of that manifold as the convex huill of the sphere. We can do this three times, one for each sphere. Call the three spheres S1, S2, and S3. Each of them represents a simple closed convex domain, and we can easily find if a point lies inside or outside of each of them, or where a line segment intersects one such domain.
Using this olde suite of codes (I say olde here because it represents 25-30 year old MATLAB code, code I have not used or maintained or literally even touched in that long of a time. A testament to MATLAB, in that it still works perfectly.) I then created 3 surface manifolds.
S1 = tilesphere(1,1001,1999,c1)
S1 =
struct with fields:
domain: [1997002×3 double]
tessellation: [3994000×3 double]
range: [1997002×0 double]
S2 = tilesphere(1,1001,1999,c2)
S2 =
struct with fields:
domain: [1997002×3 double]
tessellation: [3994000×3 double]
range: [1997002×0 double]
S3 = tilesphere(1,1001,1999,c3)
S3 =
struct with fields:
domain: [1997002×3 double]
tessellation: [3994000×3 double]
range: [1997002×0 double]
Each of them is just a fine surface tiling of a sphere, in fractions of a degree in both lattitude and longitude.
We can compute the surface area of each as just the sum of the areas of each triangle. In this case, since they are triangles, the area is computed using a volume tool, as if they were triangles in 2-dimensions. (Again, remember I am workign with a 2-manifold here.)
bodyvolume(S1)
ans =
12.5663253891995
We would expect to see an area of 4*pi of course, since this is the surface of a very finely tiled unit 3-sphere.
4*pi
ans =
12.5663706143592
Next, for each simplex in S1, I'll need to discard all triangles that lie entirely in either of the other spheres. And if a triangle lies partly in a sphere, I can dissect that triangle into either one or two triangles that lie external to the sphere. This of this as a setdiff operation. The result will be a new 2-manifold. It will not be closed anymore, like a convex hull. But just compute the area of each triangle in the three resulting manifolds, and sum the areas up.
Give me a minute to write that code...
3 Comments
Moe Szyslak
on 28 May 2026
John D'Errico
on 28 May 2026
Edited: John D'Errico
on 28 May 2026
You have not seen them ... because I've not given them out. Sigh. And the reason for that is, they take some time to learn. And they are now at least 18 years old, based on the time stamp from one of the codes. So to the extent I'd give them away, I won't maintain them.
A simplicial complex is just a set of nodes and simplexes, all stuffed into a structure. Remember, these codes go back roughly 20 years.
The tile sphere creates a surface tiling of a sphere in 3-d. The arguments are radius, number of lattitudinal nodes, and number of longitudinal nodes, then the center point as a vector. Essentially, I tiled the surface every 1 degree in both directions there.
c1 = [0 0 0];
c2 = [sqrt(2) 0 0];
c3 = [sqrt(2)/2 sqrt(3/2) 0];
S1 = tilesphere(1,181,361,c1)
S1 =
struct with fields:
domain: [64622×3 double]
tessellation: [129240×3 double]
range: [64622×0 double]
As you can see, the result has 64622 nodes in a 3-d domain, triangulated into a 2-manifold with 129240 triangles. When I would actually do this with the goal of higher resolution, perhaps at a 1 degree resolution.
The range field allowed me to use these complexes as mappings from one color space to another. It also allowed me to do surface fits, numerical integration, and a large variety of things I will not get into here.
plotsc is a code to plot the complex. The figure below is a smaller tiling, so you can see the basic shape as a tiling.
S1small = tilesphere(1,21,40,c1);
plotsc(S1small)
Now, two spheres will intersect in a plane, if they intersect, and one does not contain the other. The equation of the plane depends on the radii of the spheres, and the centers of those sphere, and how far apart are the spheres. If we have two spheres of the same radius, and centers at c1 and c2, the plane will be defined by
pip12 = (c1+c2)/2 % point in the plane
pip12 =
0.70711 0 0
normal12 = (c1-c2);normal12 = normal12/norm(normal12) % normal vector to the plane
normal12 =
-1 0 0
Similarly, there are two other planes of interest.
pip13 = (c1+c3)/2
pip13 =
0.35355 0.61237 0
pip23 = (c2+c3)/2
pip23 =
1.0607 0.61237 0
normal13 = (c1-c3);normal13 = normal13/norm(normal13)
normal13 =
-0.5 -0.86603 0
normal23 = (c2-c3);normal23 = normal23/norm(normal23)
normal23 =
0.5 -0.86603 0
Next, I would perform a "truncation" operation,
S1subsmall = planartruncate(S1small,pip12,normal12)
S1subsmall =
struct with fields:
domain: [772×3 double]
tessellation: [1452×3 double]
range: [772×0 double]
S1subsmall = planartruncate(S1subsmall,pip13,normal13)
S1subsmall =
struct with fields:
domain: [744×3 double]
tessellation: [1364×3 double]
range: [744×0 double]
Planartruncate uses a plane of interest, and discards the part of a simplicial complex that lies on one side of the plane, keeping everything on the other side. In this case, the result after two such truncations, is a sphere, with two spherical caps removed.
Perform similar operations on each sphere.
As it turns out though, since this problem is fully symmetrical, we need do it only once. Just compute the area of the manifold shown above, then multiply by 3.
bodyvolume(S1subsmall)*3
ans =
27.873
Next, do the same operation on a much more finely tiled sphere, as we want to see how well we can do.
S1 = tilesphere(1,1000,2000,c1);
S1sub = planartruncate(S1,pip12,normal12);
S1sub = planartruncate(S1sub,pip13,normal13);
bodyvolume(S1sub)*3
ans =
27.976
That suggests the surface area is slightly under 28. And I do very much trust my code, though it is pushing 20 years old. Interestingly, I see it generates a result virtually identical to that given by Torsten.
Finally, reading your comment about the region of triple intersection, if I add back in that area, I can get this:
27.976 + 0.146512896829364
ans =
28.1225128968294
So it looks like I can get the same result as you did using purely numerical tools.
Moe Szyslak
on 28 May 2026
Categories
Find more on Surface and Mesh Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
