Plotting Ellipsoid in Python

Question

I have a 3x3 positive definite matrix (a diffusion tensor to be specific) and am hoping to plot an ellipsoid who's axis directions correspond to the eigenvalues of the tensor with scale proportional to the associated eigenvalues. I am working in Python/jupyter notebooks. I have successfully used the code in the answer to this question to plot an ellipsoid with scaled axes, but I can not figure out how to alter the axis orientations to be along the directions of the eigenvectors.

Answer 1

If S is the central matrix of the ellipsoid, then this ellipsoid is the set of points x such that (x-c)' S (x-c) <= 1 where c is the center. Therefore its boundary is the isosurface (x-c)' S (x-c) = 1. Here is the way to create the isosurface with PyVista and to plot it with plotly. Maybe you need to invert S in fact (try both and check).

I think it's possible to directly plot an isosurface with plotly but I don't know how.

import numpy as np
import pyvista as pv
import plotly.graph_objects as go

# central matrix of the ellispoid
S = np.array([
  [14, 12, 8],
  [12, 12, 8],
  [8, 8, 6]    
])

def f(x, y, z):
    v = np.array([x, y, z])
    return (
        np.dot(v, np.matmul(S, v))
    )
vecf = np.vectorize(f)
# generate data grid for computing the values
X, Y, Z = np.mgrid[(-1):1:250j, (-1.5):1.5:250j, (-1.5):1.5:250j]
# create a structured grid
grid = pv.StructuredGrid(X, Y, Z)
# compute and assign the values
values = vecf(X, Y, Z)
grid.point_data["values"] = values.ravel(order = "F")
# compute the isosurface f(x, y, z) = 0
isosurf = grid.contour(isosurfaces = [1])
mesh = isosurf.extract_geometry()

# check the mesh is triangle
mesh.is_all_triangles # True
# otherwise, do: mesh.triangulate()

# extract the vertices and the triangles
points = mesh.points
triangles = mesh.faces.reshape(-1, 4)

# Now it’s a child game to plot the isosurface with Plotly:

fig = go.Figure(data=[
    go.Mesh3d(
        x=points[:, 0],
        y=points[:, 1],
        z=points[:, 2],
        colorscale = [[0, 'gold'],
                     [0.5, 'mediumturquoise'],
                     [1, 'magenta']],
        # Intensity of each vertex, which will be interpolated and color-coded
        intensity = np.linspace(0, 1, len(mesh.points)),
        # i, j and k give the vertices of the triangles
        i = triangles[:, 1],
        j = triangles[:, 2],
        k = triangles[:, 3],
        showscale = False
    )
])

fig.show()
fig.write_html("ellipsoid.html")