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

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]);  

Answer 1

[John D'Errico]

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.

[M.S. Khan]

Thanks Mr. John. God bless u for your guidance.

i am trying to run the script but its cannot be run giving some errors. could you please guide.

Also if i want to take different centre like (3,4,5) or (7, 8, 9) along with different Radii, it can be applied

[John D'Errico]

See that I changed my code, making it considerably simpler.

As far as any errors you got, I cannot possibly know what errors you get, since I cannot read your mind nor can I see inside your computer. Perhaps you just have no clue how to save an m-file function in MATLAB, or how to use one. Or, perhaps you never bothered to download spheresegmentvolume from the file exchange. Or you put it in some strange place. How can I delve into your mind?

Finally, I said that ALL that matters is the distance between the sphere centers. Can you compute the distance between two points?

Hint: A simple way to compute the Euclidean distance between two points is to compute the norm of the difference of the vectors. So if the centers are as you said, then what is the distance between them?

D = norm([3,4,5]-[7,8,9])
D =
          6.92820323027551

D is ALL that matters. Why? Suppose you translated both spheres over, each by the same amount, in ANY direction? Does the volume of intersection change? NOOOOO!!!!!! All that matters is the distance between centers. So you can center the pair of spheres anywhere in the universe, if you wish, perhaps on the surface of Mars, or pick some distant galaxy. All that matters is how far apart are the centers, and their radii.

[M.S. Khan]

Mr. John, if we take D = norm([3,4,5]-[7,8,9]), then equation for x = (D^2 - R2^2 + R1^2)/(2*D) will need to be changed because we will have equations for "y" and "z" as well.

[John D'Errico]

ARGH. NO. THINK. What is the distance between the points you gave as an example?

D = norm([3,4,5]-[7,8,9])
D =
          6.92820323027551

It is a scalar. Does the code i provide run? (and run properly?) Yes. I've now said (how many times?) that it needs to know only the distance between the centers.

Next, READ WHAT I WROTE. THINK.

"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]."

Effectively, given the distance between the centers, my code creates the problem internally as two spheres centered at the points [0 0 0], and [D,0,0]. Can you see that those two points lie at a distance of D from each other? It does not care if in the real world, those spheres live at two different points. ALL THAT MATTERS IS THE DISTANCE BETWEEN CENTERS. I've said it a few times more now. I offset one sphere by a distance D internally in the X coordinate only. You should see, if you actually look at what I wrote, that Y and Z drop out of the problem completely.

Then I do two lines of algebra, and I write:

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)

I don't care where the true centers of the sphere are. I ONLY CARE WHAT IS THE DISTANCE BETWEEN THEM, because that is the only pertinent piece of information needed.

I'm not sure why you are fixating on the question of y and z, but not thinking about what I have written, when I tell you that all that matters is the distance. Do you not understand that I am effectively changing the coordinate system of the problem in the mathematics as I wrote the equations? The volume inside a sphere does not change if I rotate the coordinate system. The volume of a rotated fragment of a sphere is still the same before and after any simple rotation of the coordinate system. A translation of the coordinate system changes nothing either.

I suppose I could have written the code to allow the user to provide the centers of the spheres? Then the very first thing I would have done in the code is to compute the distance between them, and then discard the true centers as irrelevant information anyway.

[M.S. Khan]

Thanks Mr. John. Code is running perfectly. Please dont mind for asking questions because your comments are of great importance. Do you mean if we translate or rotate by x or y or z, the distance of centers of two sphere will remain same. i will be very thankful if you have some relavent stuff for updating my knowlege, please send to me by my email: [e-mail]

Warm regards.

[M.S. Khan]

Mr. Jones, could you please check my coding? am i correct to input centeres and then using norm, i found the distance and showed the distance other than caps. Please guide me.

function [V, V_req] = sphereintersectvolume(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 
C1 = input('Centre sphere 1= ')
C2 = input('centre sphere 2 = ')
D = norm(C1 - C2)
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 
    V_max = 4/3* pi* Rmax^3; 
    V_min = 4/3* pi* Rmin^3;
    V_req = V_max + V_min - V;

end

[John D'Errico]

My name is John, not Jones.

Yes. If you move both spheres, translating them in ANY direction, or if you rotate the line between the pair of spheres in ANY way, how can that possibly change the volume of the intersection? Nothing changes, unless the distance between the spheres changes. So all that is pertinent to this question is what is the distance.

I'm not at all sure why you need to use input for the centers of the spheres, but do as you please. Input is perhaps the worst user interface tool ever created in MATLAB or in any language. Learn to avoid it.

Instead, learn to use the arguments to functions, passing in any necessary information as arguments. But that is your choice.

I assume you also want to compute the combined volume of the spheres less the intersection volume. That is indeed what V_req would represent.

[M.S. Khan]

Thanks Mr.john for all your guidance and cooperation. God bless u. I will avoid input next time and will try to use it as arguments of function. Please explain to me this function i.e. spheresegmentvolume([X,Rmin],3,Rmin). The sum of these two in the code gives volume of intersected part of spheres. What u think, the volume that i have calculated i.e the total volume minus intersected volume is correct.

Answer 2

[darova]

I used cosine theorem to calculate a1 and a2

Volume of sphere sector for first sphere can be found as:

See attached script

[M.S. Khan]

Thanks darova for guidance. the same procedrure for V will be repeated for 2nd sphere as well to find the required volume. i am trying on it.

[M.S. Khan]

Mr. Darova, please guide me, how you found volume of sector = 2/3 * pi *R1^3( 1 - Cos(a1)). and volume of cone. Waiting for your kind cooperation please.

[darova]

For sphere sector see this:

Cone volume: LINK