How to find out the intersection points between two 3D ellipsoids?

37 views (last 30 days)
I am trying to find out the intersection points between two ellipsoids. There are four intersection points in total as can be seen in the figure below.
I plotted these ellipsoids using the ellipsoid function. Please find the code below.
[xs, ys, zs] = ellipsoid(0,0,0, 1.3103e+07, 1.2800e+07, 1.2473e+07, 30);
[xx, yy, zz] = ellipsoid(0,0,0, 7.8349e+06, 1.3571e+07, 7.8349e+06, 30);
figure
surf(xs, ys, zs)
alpha 0.3
axis equal
hold on
surf(xx, yy, zz)
alpha 0.3
xlabel('k_{x} (m^{-1})')
ylabel('k_{y} (m^{-1})')
zlabel('k_{z} (m^{-1})')

Answers (1)

David Goodmanson
David Goodmanson on 17 Aug 2020
Hi PA,
As has been alluded to, the intersection of these two surfaces results in lines. These are calculated below.
The abc coefficients are for the equation
a*x^2 + b*y^2 + c*z^2 = 1
Since there are two ellipsoids, abc is a 3x2 matrix, and that is solved to get a particular value of x^2,y^2,z^2. With two equations and three unknowns there is extra room to roam in the null space of abc (which is a single vector) and trace out a line, using [particular value] + lambda*null_vector. Here lambda is a parameter. We have
xyz2point =
0
152.8432
10.4421
n = % null vector
0.6856
0.1044
-0.7204
Since the values here are squares of x,y,z, they have to all be positive, which restricts allowed values of lambda. I used a scale factor initially so as to have values of lambda that are not that large.
The code below takes positive square roots of x^2,y^2,z^2 to find x,y,z, and traces out one quarter of a closed curve which you will see on the plot.. Taking all combinations of pos or neg square roots will complete two closed curves for the intersection.
I suppose one could do linear programming to find the solution as well.
scalefac = 1e6;
m = (1/scalefac)*[1.3103e+07, 1.2800e+07, 1.2473e+07;
7.8349e+06, 1.3571e+07, 7.8349e+06];
abc = 1./m.^2;
xyz2point = abc\[1;1] % xyz2 means x^2,y^2,z^2
n = -null(abc) % minus sign for convenience
lambdamax = xyz2point(3)/(-n(3))
lambda = 0:.1:lambdamax;
xyz2 = xyz2point+ n.*lambda;
xyz = scalefac*sqrt(xyz2); % all positive square roots, rescale back
[xs, ys, zs] = ellipsoid(0,0,0, 1.3103e+07, 1.2800e+07, 1.2473e+07, 30);
[xx, yy, zz] = ellipsoid(0,0,0, 7.8349e+06, 1.3571e+07, 7.8349e+06, 30);
figure(1)
surf(xs, ys, zs)
alpha 0.3
axis equal
hold on
surf(xx, yy, zz)
hold on
alpha 0.3
xlabel('k_{x} (m^{-1})')
ylabel('k_{y} (m^{-1})')
zlabel('k_{z} (m^{-1})')
plot3(xyz(1,:),xyz(2,:),xyz(3,:),'linewidth',2,'color','c')
hold off
  1 Comment
PhysicsAnalyst
PhysicsAnalyst on 17 Aug 2020
Edited: PhysicsAnalyst on 17 Aug 2020
Awesome! That is exactly what I've been looking for!
Thank you so much, David.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!