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

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)}$$

• I might have the wrong idea... you could generate a set of points uniformly on the sphere $\mathbb{S}^{2} \subset \mathbb{R}^{3}$. Let $\left\{ c_{1},\ldots,c_{N} \right\}$ be this set of points. Then the mapping $\displaystyle \begin{bmatrix} x \\ y \\ z \end{bmatrix} \, \longmapsto \, \begin{bmatrix} ax \\ by \\ cz \end{bmatrix}$ maps a point $c_{i}$ on the sphere to a point on the ellipsoid. Commented Oct 20, 2014 at 12:55
• @jibounet Your solution would transform a uniform distribution over the ball volume into a uniform distribution over the elipsoid volume, however on surfaces it will fail. Consider very, very long cigar-like elipsoid ($a\gg b = c = 1$) - the density at the cigar's tip ($x\approx \pm a$) will be close to that of the unit sphere, but the density at $x\approx 0$ will decrease as $1/a$ relative to that on the sphere. Commented Oct 20, 2014 at 13:16
• Hm... the sphere can be parametrized as $\mathbf{F}(u,v),$ polar coordinates or whatever and the surface element can be calculated using the first fundamental form as $dS^2=H(u,v)\ du\ dv$. How would that transform under your diagonal transformation? If it's something simpler we are getting close! Commented Oct 20, 2014 at 13:18
• @CiaPan : You're right. I guess we could use the same method as here : mathworld.wolfram.com/SpherePointPicking.html with the parametrization of the ellipsoid instead on the parametrization of the sphere, right ? Commented Oct 20, 2014 at 13:21
• @jibounet: I don't think this calc can go through because in the case of the sphere the surface element is proportional to the solid angle element which is not the case for the ellipsoid. Commented Oct 20, 2014 at 13:30

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.

• If e.g. $b>>a$ and $c>>a$, we only keep a small fraction of points generated. I was originally hoping to come up with a solution without rejecting points, but now it seems an impossible task. And from all methods presented so far your's is the most simple and direct. Commented Oct 21, 2014 at 10:31
• A direct method would need to compute the function $x \mapsto$area to the left of the $X=x$ plane, and then invert that function. I haven't checked at all if either step is easy or not (I would guess they are not) Commented Oct 21, 2014 at 10:41
• I have looked into that and it comes down to inverting a function that involves elliptic functions, so it can only be done numerically by some root finding procedure. Commented Oct 21, 2014 at 12:32
• @mercio I am trying to generalize your approach to sample points on an ellipsoid in arbitrary dimension $n$, but I have never seen your notation using the $\wedge$. Does this easily generalize to arbitrary dimension and if so, what would be the probability to keep a point in dimension $n$? Commented Feb 1, 2018 at 16:26
• @CookieMaster It is called a wedge product, or exterior product. The important computation rule is $a \wedge b = - b \wedge a$. And then you still use the euclidean norm, except that now you have $n(n-1)/2$ coordinates to sum. I would expect that the probability in the end still is (smallest axis)/(largest axis) but I could be wrong Commented Feb 1, 2018 at 17:20

Here's the code, This works in ANY dimension :

import numpy as np

dim = 2
r=1
a = 10
b = 4
A = np.array([[1/a**2, 0],
[0, 1/b**2]])
L = np.linalg.cholesky(A).T

x = np.random.normal(0,1,(200,dim))

z = np.linalg.norm(x, axis=1)
z = z.reshape(-1,1).repeat(x.shape[1], axis=1)
y = x/z * r
y_new = np.linalg.inv(L) @ y.T

plt.figure(figsize=(15,15))
plt.plot(y[:,0],y[:,1], linestyle="", marker='o', markersize=2)
plt.plot(y_new.T[:,0],y_new.T[:,1], linestyle="", marker='o', markersize=5)
plt.gca().set_aspect(1)
plt.grid()


And here's the theory : https://freakonometrics.hypotheses.org/files/2015/11/distribution_workshop.pdf

• But it seems that this will not give us a uniform distribution on the surface of a ellipsoid. Commented Apr 1 at 10:56

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.

\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:

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