22
$\begingroup$

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

|cite|improve this question|||||
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – jibounet Oct 20 '14 at 12:55
  • $\begingroup$ @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. $\endgroup$ – CiaPan Oct 20 '14 at 13:16
  • $\begingroup$ 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! $\endgroup$ – Georgy Oct 20 '14 at 13:18
  • $\begingroup$ @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 ? $\endgroup$ – jibounet Oct 20 '14 at 13:21
  • $\begingroup$ @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. $\endgroup$ – Georgy Oct 20 '14 at 13:30
13
+25
$\begingroup$

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.

|cite|improve this answer|||||
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Georgy Oct 21 '14 at 10:31
  • $\begingroup$ 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) $\endgroup$ – mercio Oct 21 '14 at 10:41
  • $\begingroup$ 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. $\endgroup$ – Georgy Oct 21 '14 at 12:32
  • $\begingroup$ 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 ? $\endgroup$ – Georgy Oct 21 '14 at 12:43
  • 1
    $\begingroup$ @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 $\endgroup$ – mercio Feb 1 '18 at 17:20
2
$\begingroup$

Not sure if this is the same answer as @jibounet, but what if you use the normal distribution method, as you propose, but use the respective radii as the variances and then scale using the 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?

|cite|improve this answer|||||
$\endgroup$
  • $\begingroup$ I agree that it does not look clustered $\endgroup$ – Makogan Mar 15 '18 at 22:36
0
$\begingroup$

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

|cite|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for?Browse other questions tagged or ask your own question.