# 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. 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. 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! 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 ? 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. Oct 20, 2014 at 13:30

## 5 Answers

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. 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) 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. Oct 21, 2014 at 12:32
• One other thing is your choice of $v_1$ and $v_2$. They are not orthogonal (do they have to be?) and not unit vectors. They are just two vectors in the tangent plane. So in the cross product you get the direction right but what about the magnitude? That depends also on the angle between them. If you had chosen another set of vectors spanning the tangent plane we 'd have a different result, no ? Oct 21, 2014 at 12:43
• @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 Feb 1, 2018 at 17:20

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:

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

• I agree that it does not look clustered Mar 15, 2018 at 22:36
• Take the limit $b \to 0$ and $c \to 0$. This will only generate values at the two endpoints of the (degenerated) ellipsoid.
– Marc
Aug 19, 2020 at 13:41
• @Marc, I'm not sure what you mean. In the degenerate case you describe $𝑥_1$ remains a normally distributed RV. And the limits for the $x_2$ and $x_3$ terms in the equation for d are not zero, but are simply scaled back to unit normal RVs. How would this only generate the end points? Aug 24, 2020 at 8:01
• @elhuhdron I'm sorry. I didn't read carefully enough and wrote before thinking. So, yes, now that I read your solution again, I agree that the points won't necessarily just be at the endpoints. But now I am wondering what the difference between your suggestion and the first comment by jibounet below the original question is.
– Marc
Aug 24, 2020 at 16:55
• Basically in @jibounet the scaling happens after the points are on the unit sphere. In my answer they are also scaled before the transformation because of the non-unit variance of the sampled variables, and since the transformation is non-linear, this is not the same. Jan 27 at 9:51

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

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