How to generate points uniformly distributed on the surface of an ellipsoid?

Solution 1:

One way to proceed is to generate a point uniformly on the sphere, apply the mapping $f : (x,y,z) \mapsto (x'=ax,y'=by,z'=cz)$ and then correct the distortion created by the map by discarding the point randomly with some probability $p(x,y,z)$ (after discarding you restart the whole thing).

When we apply $f$, a small area $dS$ around some point $P(x,y,z)$ will become a small area $dS'$ around $P'(x',y',z')$, and we need to compute the multiplicative factor $\mu_P = dS'/dS$.

I need two tangent vectors around $P(x,y,z)$, so I will pick $v_1 = (dx = y, dy = -x, dz = 0)$ and $v_2 = (dx = z,dy = 0, dz=-x)$

We have $dx' = adx, dy'=bdy, dz'=cdz$ ; $Tf(v_1) = (dx' = adx = ay = ay'/b, dy' = bdy = -bx = -bx'/a,dz' = 0)$, and similarly $Tf(v_2) = (dx' = az'/c,dy' = 0,dz' = -cx'/a)$

(we can do a sanity check and compute $x'dx'/a^2+ y'dy'/b^2+z'dz'/c^2 = 0$ in both cases)

Now, $dS = v_1 \wedge v_2 = (y e_x - xe_y) \wedge (ze_x-xe_z) = x(y e_z \wedge e_x + ze_x \wedge e_y + x e_y \wedge e_z)$ so $|| dS || = |x|\sqrt{x^2+y^2+z^2} = |x|$

And $dS' = (Tf \wedge Tf)(dS) = ((ay'/b) e_x - (bx'/a) e_y) \wedge ((az'/c) e_x-(cx'/a) e_z) = (x'/a)((acy'/b) e_z \wedge e_x + (abz'/c) e_x \wedge e_y + (bcx'/a) e_y \wedge e_z)$

And finally $\mu_{(x,y,z)} = ||dS'||/||dS|| = \sqrt{(acy)^2 + (abz)^2 + (bcx)^2}$.

It's quick to check that when $(x,y,z)$ is on the sphere the extrema of this expression can only happen at one of the six "poles" ($(0,0,\pm 1), \ldots$). If we suppose $0 < a < b < c$, its minimum is at $(0,0,\pm 1)$ (where the area is multiplied by $ab$) and the maximum is at $(\pm 1,0,0)$ (where the area is multiplied by $\mu_{\max} = bc$)

The smaller the multiplication factor is, the more we have to remove points, so after choosing a point $(x,y,z)$ uniformly on the sphere and applying $f$, we have to keep the point $(x',y',z')$ with probability $\mu_{(x,y,z)}/\mu_{\max}$.

Doing so should give you points uniformly distributed on the ellipsoid.

Solution 2:

I have been working on this problem this week and I have found the following approach effective.

The ellipsoid is parameterized by

$$\vec f(t, u) = \begin{bmatrix}a \cos u \cos t \\ b \cos u \sin t \\ c \sin u\end{bmatrix}, \quad \begin{align}t &\in [0, 2 \pi] \\ u &\in [0, \pi]\end{align}$$

Naively, we might pick $T \sim U(0, 2 \pi)$ and $U \sim U(0, \pi)$ and choose $\vec f(T, U)$ as our random point, but this is equivalent to stretching the sphere, which creates distortion as others have pointed out.

Instead, define the regions

$$\begin{align} D_{t'} :\quad & t \in [0, t'] , && u \in [0, \pi] \\ D_{u'} :\quad & t \in [0, 2 \pi] , && u \in [0, u'] \end{align}$$

And the functions

$$\begin{align} S_t(t') \equiv & \iint_{D_{t'}} dS \\ S_u(u') \equiv & \iint_{D_{u'}} dS \end{align}$$

That is, $S_t$ yields the cumulative surface area on $t\in [0,t']$ and for all $u$ , while $S_u$ does vice-versa.

Next, pick two numbers $X_1, X_2 \sim U(0,A)$ where $A$ is the total surface area $$A = S_t(2\pi) = S_u(\pi)$$

Since $S_t$ and $S_u$ are monotonically increasing over the parameter space, they can be evaluated numerically and inverted through interpolation.

Our random parameters are then

$$\begin{align} T = S_t^{-1}(X_1) \\ U = S_u^{-1}(X_2) \\ \end{align}$$

And a random point on the surface is $$\vec f(T, U)$$

The advantage of this method is that because we don't have to discard any points, we don't need to compute a multiplicative factor particular to this ellipse (as another answerer does). In fact, this method picks points uniformly distributed on any parametric surface whose parameter region is square.

It is also computationally efficient, because we can use the same data (namely, a grid of parameters and the values of the function at those parameters) to determine the overall surface area and, by interpolation, the cumulative surface area functions and their inverses.

I've developed a small Python module that implements this approach as r_surface(). Example 3 in the module uses an ellipsoid:

This implementation is indebted to the responses to these SE questions:

• Generate random points on perimeter of ellipse (Math SE)
• Generate random points on perimeter of ellipse (Code Review SE)