[Math] Plotting an ellipsoid using eigenvectors and eigenvalues

conic sectionseigenvalues-eigenvectorslinear algebraquaternionsrotations

I have an ellipse represented by a positive definite symmetric matrix, $E$.
My goal is to plot this ellipsoid, $E$, by scaling and rotating a unit sphere.

I do not have access to the list of points constituting the mesh of the sphere.
Therefore, simply transforming the individual points of the unit sphere mesh by multiplying them with the square root matrix, $\sqrt{E}$, is not a solution.

As such, here is my alternative process:

1.

I started by computing the eigenvalues $e$ and right eigenvectors $x$ of $E$.
I am using NumPy, so I used numpy.linalg.eig for the job:

e, x = np.linalg.eig(E)

Given $e$, I can get the radii of the ellipsoid with:

r = np.sqrt(e)

The eigenvectors $x$ are unit vectors in the principal semi-axes directions and can be scaled by the radii $r$ to yield radius vectors that can be plotted:

  • x[:, 0] (red)
  • x[:, 1] (green)
  • x[:, 2] (blue)

01-semi-axes

Because the unit eigenvectors $x$ seem to flip at some point, I decided to elongate them in each corresponding negative direction in order to get a better visualisation:

02-axes

These are effectively the principal axes of the ellipsoid (I like to think of them as the "skeleton" of $E$).
And, for what's worth, they are correct. They display a continuous motion of what I expect to be the correct representation of what I am trying to plot.

Important note
OK, something that by now already seemed a little bit off was the order of the plotted axes. You see, my world frame axes look like this:

world frame axes

But let's continue…

2.

Given $r$, I can scale my unit sphere in $x$, $y$, $z$ with r[0], r[1], and r[2], respectively.

Afterwards, I just need to rotate it to match the "skeleton", i.e., the principal axes.
For that I need to use a quaternion (again, this is a constraint imposed by the tool I have available to plot the scaled sphere).

So, how do I get the quaternion?
I just built a $3$x$3$ rotation matrix using the unit eigenvectors $x$:

rot_matrix = np.array([x[:, 0],
                       x[:, 1],
                       x[:, 2]])

…and then I am converting this rotation matrix into a quaternion using a library called KDL.

3.

Finally, now that I have the orientation in quaternion format, I can set the orientation for the scaled ellipsoid, and voilla!

Right? Not really… Here is what I get as a result:

03-ellipsoid

There is clearly some rotation going on, but it does not match the one of the principal axes… The axes and the ellipsoid should overlay perfectly and move as one single object.

And basically, this is where I am stuck at…
Any help or suggestion would be more than welcome!

Best Answer

Turns out I was indeed doing two things wrong!

1. Axes order
Like I suspected — and mentioned in the OP as an important note — the frame defined by the axes of unit eigenvectors $x$ was sometimes switching direction, and thus not respecting the right-hand rule and world frame convention.

To overcome this, I just reassigned the 3rd eigenvector as the cross product between the 1st and 2nd eigenvectors:

x[:, 2] = np.cross(x[:, 0], x[:, 1])

Since the eigenvectors $x$ returned by numpy.linalg.eig are already normalized (i.e., unit eigenvectors), there is no need to normalise the result of the cross product.

2. Rotation matrix
The second thing I was not doing properly was during the construction of the $3$x$3$ rotation matrix. I was missing a transpose operation.

To solve this I just changed the indexing of the eigenvectors $x$ to:

rot_matrix = np.array([x[0, :],
                       x[1, :],
                       x[2, :]])

Using numpy.transpose would have worked too.

Below is the correct expected final result:

final-result

Special thanks to Vladimir Ivan (a.k.a. "Slovak")!

Related Question