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

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})')

[KSSV]

Not four points, there will be many intersection points. Refer this: https://www.mathworks.com/matlabcentral/answers/93623-how-do-i-plot-the-line-of-intersection-between-two-surfaces

[PhysicsAnalyst]

Many thanks, SSV!

I will check out this link.

Answer 1

[David Goodmanson]

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

[PhysicsAnalyst]

Awesome! That is exactly what I've been looking for!

Thank you so much, David.