Create a mesh of points on an ellipsoid

Question

I have given 4 vectors $m, a, b, c \in \mathbb{R}^3$ (the center of an ellipsoid and its 3 main axes). I am looking for an clever way to compute a mesh of points on the ellipse, which is not to computational heavy (I have thousands to plot in a single figure).

Wikipedia says I should just use the parameterform $$E(\phi, \theta) = m + a~\cos (\theta) \cos (\phi) + b~\cos (\theta) \sin (\phi) +c ~ \sin (\theta),~~~ \phi \in (0,2\pi),~ \theta\in (-\pi/2, \pi /2)$$ but it allways plots the wrong ellisoid.

I also read that there is an easier way by using the matrix in the alternative definition of an ellipsoid:

$$ (x - m)^T A (x-m) = 1 $$

where with $v_1, v_2, v_3$ parallel to $a,b,c$ and $\lambda_1, \lambda_2, \lambda_3$ equal to $1/||a||^2, 1/||b||^2, 1/||c||^2$.

I know there is the python package nestle which can draw ellipsoids using this matrix, but I won't find any documentaion on how it dows that (I would like to do it without the package if possible).

Thank you for any help :)

Answer 1

I found my answer by looking through the original code of the package.

You'll find the ellipsoid function and the code:

class Ellipsoid(object):
"""An N-ellipsoid.
Defined by::
    (x - v)^T A (x - v) = 1
where the vector ``v`` is the center of the ellipse and ``A`` is an N x N
matrix. Assumes that ``A`` is symmetric positive definite.
Parameters
----------
(...)
"""

def __init__(self, ctr, a):
    self.n = len(ctr)
    self.ctr = ctr    # center coordinates
    self.a = a        # ~ inverse of covariance of points contained
    self.vol = vol_prefactor(self.n) / np.sqrt(np.linalg.det(a))

    # eigenvalues (l) are a^-2, b^-2, ... (lengths of principle axes)
    # eigenvectors (v) are normalized principle axes
    l, v = np.linalg.eigh(a)
    self.axlens = 1. / np.sqrt(l)

    # Scaled eigenvectors are the axes: axes[:,i] is the i-th
    # axis.  Multiplying this matrix by a vector will transform a
    # point in the unit n-sphere into a point in the ellipsoid.
    self.axes = np.dot(v, np.diag(self.axlens))

Here they also explain how they calculate ellipsoid.axes

While I do not fully understand the math behind this, I now can use a simpler version of this to make my code more efficient.

So how do I plot the ellipsoid? My solution is something like this:

# ctr ... (3x1) numpy array
# ellipsoid ... (13x1) numpy arrays of the lenght a,b,c of each axis and its directions 
# v1, v2, v3 in the ellipsoid. The directions must be a orthonormal 
# system. The order is (a, v1, b, v2, c, v3),
def plot_ellipsoid(ctr, ellipsoid):
  # the axes of the ellipsoid (normalizes and orthogonal)
  e_vec = np.array([ellipsoid[1:4], ellipsoid[5:8], ellipsoid[9:12]])
  #the length of each axis of the ellipsoid
  (a,b,c) = [ellipsoid[0], ellipsoid[4], ellipsoid[8]]

  ## calculate a few point on the ellipsoid to create a simple mesh  
  # convert to matrix for the equation (x-ctr)^T * A * (x-ctr) = 1
  e_val = np.diag([1/a**2,1/b**2,1/c**2])
  A = np.linalg.inv(e_vec) @ e_val @ e_vec
  # stuff from nestle
  A = A.T @ A
  lam, vec = np.linalg.eigh(A)
  axes = np.dot(vec, np.diag(1. / np.sqrt(lam)))

  # points on unit sphere
  phi   = np.linspace(0.0, 2.0 * np.pi, 8)  # TODO should I mark this as something from the nestle python package site?
  theta = np.linspace(0.0, np.pi, 8)
  x = np.outer(np.ones_like(phi), np.cos(theta))
  y = np.outer(np.sin(phi),       np.sin(theta))
  z = np.outer(np.cos(phi),       np.sin(theta))
  # transform points to ellipsoid
  for i in range(len(x)):
      for j in range(len(x)):
          x[i,j], y[i,j], z[i,j] = ctr + axes @ [x[i,j],y[i,j],z[i,j]] 
return x, y, z

(might not work 100% since I had to change some things)

Have fun :)