72 views (last 30 days)

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

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.

Opportunities for recent engineering grads.

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

Start Hunting!
## 0 Comments

Sign in to comment.