I have coordinates corresponding to a set of 2D contours, each corresponding to different heights. These contours do not draw out a perfect ellipsoid in 3D, and instead what I would like to do is to find the best fitting ellipsoid. I do not have any knowledge on the origin of this ellipsoid.

My first thought was to incorporate some type of least squares algorithm, where I find the ellipsoid parameters that minimize the distance between points. I imagine this would be quite expensive and not too far from a brute force approach. I am convinced there is a more elegant and efficient way of doing this. If there is an existing library that handles this (preferably in Python) that would be even better.

I have already seen a related question (

The best fit ellipse is described by the general equation

A x^2 B y^2 2C x y 2D x 2E y = 1

and once the coefficients (A,B,C,D,E) are found, the ellipse of fully described. See below in how to find the the curve coordinates (x,y) from the coefficients and a parameter t=0 .. 1.

To find the coefficients of the ellipse, form 5 vectors, each a column of a n×5 matrix Q

for i = 1 to n
  Q(i,1) = x(i)^2
  Q(i,2) = y(i)^2
  Q(i,3) = 2*x(i)*y(i)
  Q(i,4) = 2*x(i)
  Q(i,5) = 2*y(i)
next i

and a vector K filled with 1 for the right-hand side

for i = 1 to n
  K(i) = 1.0
next i

Find the coefficients using a least-squares fit with some linear algebra

[A,B,C,D,E] = inv(tr(Q)*Q)*tr(Q)*K

where tr(Q) is the transpose of Q and * is matrix/vector product

Now we need to extract the geometric properties of the ellipse from the coefficient. I want to have a the semi-major axis, b the semi-minor axis, φ the rotation angle, xc the x-axis center, yc the y-axis center.

xc = -(B*D-C*E)/(A*B-(C^2))
yc = -(A*E-C*D)/(A*B-(C^2))
φ = atan( 2*C/(A-B) )/2
a = SQRT(2*(A*(B E^2) B*D^2-C*(C 2*D*E))/((A*B-C^2)*(A B-SQRT((A-B)^2 4*C^2))))
b = SQRT(2*(A*(B E^2) B*D^2-C*(C 2*D*E))/((A*B-C^2)*(A B SQRT((A-B)^2 4*C^2))))

Finally to plot the ellipse you need to generate a set of points (x,y) from the curve parameter t=0..1 using the above 5 coefficients.

1. Generate the centered aligned coordinates (u,v) with

   u = a*cos(2*π*t)
   v = b*sin(2*π*t)
   

2. Generate the centered rotated coordinates (x',y') with

   x' = u*cos(φ) - v*sin(φ)
   y' = u*sin(φ)   v*cos(φ)
   

3. Generate the ellipse coordinates (x,y) with

   x = x'   xc
   y = y'   yc
   

The result is observed above in the first picture.

Now for the total solution, each 2D slice would have its own ellipse. But all the slices would not generate an ellipsoid this way.

Extending the above into 3D coordinates (x,y,z) is doable, but the math is quite involved and I feel [SO] is not a good place to develop such an algorithm. You can hack it together, by finding the average center for each slice (weighted by the ellipse area π*a*b). Additionally, the rotation angle should be the same for all contours, and so another averaging is needed. Finally, the major and minor axis values would fall on an elliptical curve along the z-axis and it would require another least-fit solution. This one is driven by the equation

(x/a)^2 (y/b)^2 (z/c)^2 = 1

but rather in the aligned coordinates (u,v,w)

(u/a)^2 (v/b)^2 (w/c)^2 = 1