MATLAB: How to find volume of intersected part when two spheres are intersected using Matlab

image analysisimage processingmathematicsMATLAB

using following code, two spheres are intersected. i want to find the intersected volume which is circle. The graph shows the intersected spheres. Thanks for all cooperation in advance.
radius = 10; 
numPoints = 1000;	% Use a large value.

% Get a 3-by-numPoints list of (x, y, z) coordinates.

r = randn(3, numPoints);
r = bsxfun(@rdivide, r, sqrt(sum(r.^2,1)));
r = radius * r;
x = r(1,:) +3 ; 
y = r(2,:)+2; 
z = r(3,:)+4; 
figure(1)
scatter3(x, y, z);
axis square; % Make sure the aspect ratio is maintained as it's displayed and rotated.

xlabel('X', 'FontSize', 20);
ylabel('Y', 'FontSize', 20);
zlabel('Z', 'FontSize', 20);
% Enlarge figure to full screen.

set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
%
%msgbox('Now use the circular arrow icon on the toolbar to rotate the sphere.');  
%% 2nd sphere
radius = 13; 
numPoints = 1000;	% Use a large value.
% Get a 3-by-numPoints list of (x, y, z) coordinates.
r = randn(3, numPoints);
r = bsxfun(@rdivide, r, sqrt(sum(r.^2,1)));
r = radius * r;
% Extract the x, y, and z coordinates from the array.
x = r(1,:)+13; 
y = r(2,:)+12; 
z = r(3,:)+11; 
% Display the shell of points
hold on
scatter3(x, y, z);
axis square; % Make sure the aspect ratio is maintained as it's displayed and rotated.
xlabel('X', 'FontSize', 20);
ylabel('Y', 'FontSize', 20);
zlabel('Z', 'FontSize', 20);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);

Best Answer

You can use my spheresegmentvolume utility to solve this. It is on the file exchange.
https://www.mathworks.com/matlabcentral/fileexchange/18641-spheresegmentvolume
The idea is, if two spheres intersect, then you can compute the volume of TWO spherical caps. The intersected volume is the sum of those two caps.
For example, consider two spheres, of different radii and centers, but such that the spheres do intersect, with some finite volume.
I'll be arbitrary here for my example. Suppose sphere1 has center [0 0 0], with a radius of 2. I'll create sphere 2 with a center at [4,0,0]. and a radius of 3.
The distance between centers is thus 4 units, and the sum of the radii is 5, which is greater than the distance between centers. So there is a clear intersection between the two spheres. Had the sum of the radii been smaller than the distance between centers, then no intersection could exist.
As you can also appreciate, as long as one of the spheres is not entirely inside the other, then the surfaces of the two spheres will always intersect in a circle. You can test to see if either sphere lives entirely inside the other easily enough too. In that case, the volume of intersection will be just the volume of the smaller sphere.
A non-trivial intersection volume can now be visualized as the sum of two spherical caps.
We can reduce the problem simply to this case, where there are two spheres, separated by a distance D. The distance is all that matters. So I will arbitrarily put one sphere centered at the origin, thus [0,0,0], and the second sphere centered at the point [D,0,0].
The equations of those two spheres are
x^2 + y^2 + z^2 = R1^2
(x-D)^2 + y^2 + z^2 = R2^2
Now we wish to find the plane of intersection of the spheres, IF there is some non-trivial intersection at all. So just subtract the sphere equations. This yields a simple equation, where many terms drop out.
-2*D*x = R2^2 - R1^2 - D^2
Solving for x, we can then find the location in x of the plane of intersection.
x = (D^2 - R2^2 + R1^2)/(2*D)
This is implemented in the code below. As I said, it requires ONLY knowledge of the radii and the distance between centers.
function V = sphereintersectvolume(D,R1,R2)
% V = sphereintersectvolume(D,R1,R2)
% arguments: (input) 
% D - distance between centers of the two spheres
% R1, R2 - radii of the pair of spheres 
if any([D,R1,R2] < 0)
  error('all radii and distances should in general be non-negative numbers.')
end
% we need to know which radius is smaller of the two
Rmin = min(R1,R2);
Rmax = max(R1,R2);
if D > (Rmin + Rmax)
  V = 0;
elseif Rmax >= (D + Rmin)
  % this triggers if one sphere lies entirely inside the other.
  % then the volume returned is just the volume of the smaller sphere.
  V = 4/3*pi*Rmin^3;
else
  % There is a non-trivial intersection
  X = (D^2 - Rmax^2 + Rmin^2)/(2*D);
  V = spheresegmentvolume([X,Rmin],3,Rmin) + spheresegmentvolume([D - X,Rmax],3,Rmax);
end
Does this work? Of course. Lets try it out. I'll consider two spheres here, with radii 2 and 3. Now, assume that we start with the circles centers very close together.
If the circles have centers that are no larger than 1 unit apart, then the intersection volume is just the volume of the smaller sphere. That is because the difference in radii is 1 unit here.
Circles that are farther apart than 5 units in this case will have a null intersection. So, as the distance moves from 0 to 5 units and above here, we should see the volume of intersection decrease from 33.5103216382911 (the volume of the smaller sphere) to zero.
Dlist = [0.1 0.99 1.01 1.5 1.99 2.01 2.5 2.99 3 3.01 4 4.99 5 6];
for D = Dlist
  disp([D,sphereintersectvolume(D,R1,R2)])
end
                       0.1          33.5103216382911
                      0.99          33.5103216382911
                      1.01           33.508456385047
                       1.5           30.466903755126
                      1.99          24.8636525053106
                      2.01          24.6162550618698
                       2.5          18.4895817633149
                      2.99          12.6782340788667
                         3          12.5663706143592
                      3.01          12.4548329430293
                         4          3.46884188833873
                      4.99      0.000376697840158655
                         5                         0
                         6                         0
As you can see in the test, that is exactly what happens. In this case, if the distance was no larger than 1, the volume is that of the smaller sphere. If the distance is larger than 5=2+3, then the volume is zero.
Related Question