20

Does anyone have sample code for plotting ellipsoids? There is one for sphere on matplotlib site, but nothing for ellipsoids. I am trying to plot

x**2 + 2*y**2 + 2*z**2 = c

where c is a constant (like 10) that defines an ellipsoid. I tried the meshgrid(x,y) route, reworked the equation so z is on one side, but the sqrt is a problem. The matplotlib sphere example works with angles, u,v, but I am not sure how to work that for ellipsoid.

Hooked
  • 77,871
  • 38
  • 181
  • 253
BBSysDyn
  • 4,129
  • 8
  • 43
  • 60

3 Answers3

33

Here is how you can do it via spherical coordinates:

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure(figsize=plt.figaspect(1))  # Square figure
ax = fig.add_subplot(111, projection='3d')

coefs = (1, 2, 2)  # Coefficients in a0/c x**2 + a1/c y**2 + a2/c z**2 = 1 
# Radii corresponding to the coefficients:
rx, ry, rz = 1/np.sqrt(coefs)

# Set of all spherical angles:
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)

# Cartesian coordinates that correspond to the spherical angles:
# (this is the equation of an ellipsoid):
x = rx * np.outer(np.cos(u), np.sin(v))
y = ry * np.outer(np.sin(u), np.sin(v))
z = rz * np.outer(np.ones_like(u), np.cos(v))

# Plot:
ax.plot_surface(x, y, z,  rstride=4, cstride=4, color='b')

# Adjustment of the axes, so that they all have the same span:
max_radius = max(rx, ry, rz)
for axis in 'xyz':
    getattr(ax, 'set_{}lim'.format(axis))((-max_radius, max_radius))

plt.show()

The resulting plot is similar to

enter image description here

The program above actually produces a nicer looking "square" graphics.

This solution is strongly inspired from the example in Matplotlib's gallery.

Eric O Lebigot
  • 86,421
  • 44
  • 210
  • 254
  • thanks for this answer. The axis adjustment lines gave me an error however, "AttributeError: 'Axes3DSubplot' object has no attribute 'set_zlim'". Is it because I am on older version of matplotlib (I am on 1.0.1) – BBSysDyn Oct 19 '11 at 12:28
  • 3
    yep that was the problem. I upgraded to 1.1.0 and it all works now. – BBSysDyn Oct 19 '11 at 12:42
  • 1
    rx, ry, rz = 1/np.sqrt(coef) Should be rx, ry, rz = 1/np.sqrt(coefs) – derchambers Jun 10 '15 at 15:38
  • If I have semi-major axis as a and c, semi minor axes as b, will the coeff be (1/a**2,1/b**2,1/c**2) ? – Alex Punnen Oct 28 '15 at 12:40
  • Yes. Of course in that case you can directly do `rx, ry, rz = a, b, c` (no need for coefficients `coefs`). – Eric O Lebigot Oct 28 '15 at 20:35
  • For those, like me, who didn't initially understand how `np.outer` works here, use this instead: `u = np.linspace(0, 2 * np.pi, 100); v = np.linspace(0, np.pi, 100); u,v = np.meshgrid(u,v,indexing='ij'); x = rx * np.cos(u) * np.sin(v); y = ry * np.sin(u) * np.sin(v); z = rz * np.cos(v)` – mhdadk Jan 22 '22 at 18:34
14

Building on EOL's answer. Sometimes you have an ellipsoid in matrix format:

A and c Where A is the ellipsoid matrix and c is a vector representing the centre of the ellipsoid.

import numpy as np
import numpy.linalg as linalg
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# your ellispsoid and center in matrix form
A = np.array([[1,0,0],[0,2,0],[0,0,2]])
center = [0,0,0]

# find the rotation matrix and radii of the axes
U, s, rotation = linalg.svd(A)
radii = 1.0/np.sqrt(s)

# now carry on with EOL's answer
u = np.linspace(0.0, 2.0 * np.pi, 100)
v = np.linspace(0.0, np.pi, 100)
x = radii[0] * np.outer(np.cos(u), np.sin(v))
y = radii[1] * np.outer(np.sin(u), np.sin(v))
z = radii[2] * np.outer(np.ones_like(u), np.cos(v))
for i in range(len(x)):
    for j in range(len(x)):
        [x[i,j],y[i,j],z[i,j]] = np.dot([x[i,j],y[i,j],z[i,j]], rotation) + center

# plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(x, y, z,  rstride=4, cstride=4, color='b', alpha=0.2)
plt.show()
plt.close(fig)
del fig

So, not too much new here, but helpful if you've got an ellipsoid in matrix form which is rotated and perhaps not centered at 0,0,0 and want to plot it.

minillinim
  • 630
  • 6
  • 10
  • 2
    Hi, why you need to compute the inverse square root of the singular values, instead of just the square root, to get the radii? I mean why radii = 1.0/np.sqrt(s) instead of radii = np.sqrt(s). I think the later gives the correct radii. – isti_spl Jan 08 '14 at 11:56
  • You don't need any of this complicated math at all, least of all the SVD. The whole point of the covariance matrix is that if you multiply a sphere by it, you get the ellipse. – Mad Physicist Apr 23 '21 at 10:54
  • I've posted an answer demonstrating what I mean. – Mad Physicist Apr 23 '21 at 11:04
3

If you have an ellipsoid specified by an arbitrary covariance matrix cov and offset bias, you do not need to figure out the intuitive parameters of the ellipsoid to get the shape. Specifically, you don't need the individual axes or rotations. The whole point of the matrix is that it transforms a unit sphere (represented by the identity matrix) into your ellipse.

Starting with

u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)

Make a unit sphere

x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones_like(u), np.cos(v))

Now transform the sphere:

ellipsoid = (cov @ np.stack((x, y, z), 0).reshape(3, -1) + bias).reshape(3, *x.shape)

You can plot the result pretty much as before:

ax.plot_surface(*ellipsoid, rstride=4, cstride=4, color='b', alpha=0.75)
Mad Physicist
  • 95,415
  • 23
  • 151
  • 231
  • When I do this, I get a sample covariance ellipsoid which is quite a bit larger than the sampled points. Is there a way to scale it to a confidence interval? – Translunar Aug 11 '21 at 20:13
  • @Translunar. Any scalar multiplier is valid on the matrix – Mad Physicist Aug 11 '21 at 23:29
  • This answer is a little bit of a hack Take zero mean random vector X with covariance matrix is C = < X X’>. C can be decomposed as RLR’ where R is the rotation matrix of eigenvectors and L is the diagonal matrix of eigenvalues. With S = \sqrt{L} we have C = RSS’R’$. Suppose N is a normal random vector so that < NN’> = I. We can transform N into X by X = RSN. This is confirmed because = = RS < N N’>S’R’ = RSIS’ R’ = C. So the “proper” transformation to be done is X = RSN. – Jagerber48 Aug 17 '21 at 02:47
  • However, in this case you get away with it because our “dataset” is the points on a sphere. When you act R’ on that you get the same set of points back, so your transformation C = RSS'R' is effectively the transformation RSS’. This is almost RS like it should be, you just scale the data twice by S instead of once. This is why Translunar complains that the covariance ellipsoid is not scaled appropriately to the data. So this method does give you an ellipsoid with almost the correct shape, but like I say it's a little bit of a hack. – Jagerber48 Aug 17 '21 at 02:48