MATLAB Answers

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

72 views (last 30 days)
M.S. Khan
M.S. Khan on 7 Mar 2020
Commented: darova on 10 Mar 2020
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]);

  0 Comments

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 7 Mar 2020
Edited: John D'Errico on 8 Mar 2020
You can use my spheresegmentvolume utility to solve this. It is on the file exchange.
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.

  8 Comments

Show 5 older comments
M.S. Khan
M.S. Khan on 9 Mar 2020
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
John D'Errico on 9 Mar 2020
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
M.S. Khan on 9 Mar 2020
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.

Sign in to comment.

More Answers (1)

darova
darova on 7 Mar 2020
Edited: darova on 7 Mar 2020
I used cosine theorem to calculate a1 and a2
Volume of sphere sector for first sphere can be found as:
See attached script

  3 Comments

M.S. Khan
M.S. Khan on 8 Mar 2020
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
M.S. Khan on 10 Mar 2020
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!