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

조회 수: 29(최근 30일)
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]);

채택된 답변

John D'Errico
John D'Errico 2020년 3월 7일
편집: John D'Errico 2020년 3월 8일
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
M.S. Khan
M.S. Khan 2020년 3월 9일
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.

댓글을 달려면 로그인하십시오.

추가 답변(1개)

darova
darova 2020년 3월 7일
편집: darova 2020년 3월 7일
I used cosine theorem to calculate a1 and a2
Volume of sphere sector for first sphere can be found as:
See attached script
  댓글 수: 3

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by