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

Question

I am trying to find a way to generate random points uniformly distributed on the surface of an ellipsoid.

If it was a sphere there is a neat way of doing it: Generate three $N(0,1)$ variables $\{x_1,x_2,x_3\}$, calculate the distance from the origin

$$d=\sqrt{x_1^2+x_2^2+x_3^2}$$

and calculate the point

$$\mathbf{y}=(x_1,x_2,x_3)/d.$$

It can then be shown that the points $\mathbf{y}$ lie on the surface of the sphere and are uniformly distributed on the sphere surface, and the argument that proves it is just one word, "isotropy". No prefered direction.

Suppose now we have an ellipsoid

$$\frac{x_1^2}{a^2}+\frac{x_2^2}{b^2}+\frac{x_3^2}{c^2}=1$$

How about generating three $N(0,1)$ variables as above, calculate

$$d=\sqrt{\frac{x_1^2}{a^2}+\frac{x_2^2}{b^2}+\frac{x_3^2}{c^2}}$$

and then using $\mathbf{y}=(x_1,x_2,x_3)/d$ as above. That way we get points guaranteed on the surface of the ellipsoid but will they be uniformly distributed? How can we check that?

Any help greatly appreciated, thanks.

PS I am looking for a solution without accepting/rejecting points, which is kind of trivial.

EDIT:

Switching to polar coordinates, the surface element is $dS=F(\theta,\phi)\ d\theta\ d\phi$ where $F$ is expressed as $$\frac{1}{4} \sqrt{r^2 \left(16 \sin ^2(\theta ) \left(a^2 \sin ^2(\phi )+b^2 \cos ^2(\phi )+c^2\right)+16 \cos ^2(\theta ) \left(a^2 \cos ^2(\phi )+b^2 \sin ^2(\phi )\right)-r^2 \left(a^2-b^2\right)^2 \sin ^2(2 \theta ) \sin ^2(2 \phi )\right)}$$

Answer 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.

Answer 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)

Answer 3

What if you use the normal distribution method, as you propose, but use the respective radii as the variances? Then scale with d, as proposed. Empirically this does not look more bunched (at least to me) at "cigar poles" than a brute force selection from uniform distributed points in $\mathbb{R}^3$.

So, the method is generate random variables:

$$x_1 \sim N(0,a^2)$$ $$x_2 \sim N(0,b^2)$$ $$x_3 \sim N(0,c^2)$$

Then the rest is as suggested:

$$d=\sqrt{\frac{x_1^2}{a^2}+\frac{x_2^2}{b^2}+\frac{x_3^2}{c^2}}$$

and use points $\mathbf{y}=(x_1,x_2,x_3)/d$

The proof of the original spherical selection is simply that normal distributions are radially symmetric? Does that not still apply with unequal variances?

NOTE: This is different from the answer in the comment from @jibounet that proposes scaling the points that were generated on the surface of a sphere, because the terms in $d$ are different, i.e. for

$$x \sim N(0,n^2)$$ $$ \frac{x^2}{n^2} \neq N(0,1) $$

Answer 4

Idea for an approximate solution: divide the ellipsoid in small enough, flat enough quasi-rectangles $P_i$, choose an almost-area-preserving parametrization of each piece: $$f_i:R_i\longrightarrow P_i$$ with $R_i\subset\Bbb R^2$ a rectangle. Choose randomly an index $j$ with probability $\text{area}(R_j)/\sum_i\text{area}(R_i)$, choose a point in $R_j$ and finally apply $f_j$.

Answer 5

• Simulate some points on the sphere (with arbitrary dimension) $x_1$, $x_2$, $...$

• Take the upper Cholesky factor $U$ of the shape matrix $A$ of the ellipsoid ($x'Ax = r^2$)

• Then the solution $Y$ of $UY = (x_1 \, x_2 \, ...)$ is a matrix of uniformly random points on the surface of the ellipsoid.