[Math] How to generate points uniformly distributed on the surface of an ellipsoid

analytic geometryellipsoidsprobabilityrandom variablesstochastic-processes

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

Best Answer

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.

Related Solutions

Uniform Distribution of Points on a Surface in $\mathbb{R}^n$

I did find the following ovservation. Let us denote $\rho = \phi \circ \psi$. We want to find $\psi$ such that for all open subsets $W \subseteq U$, the submanifold $\rho(W)$ satisfies that $\rho$ preserves the probabiliy:

$$\frac{\int_W \sqrt{\det (D\rho ^T D\rho)}\text{d}u}{\text{vol}_k(M)}=\frac{\int_W \text{d}u}{\text{vol}_k(U)}$$

Since this is for all $W$ we must have $$\frac{\sqrt{\det (D\rho ^T D\rho)}}{\text{vol}_k(M)}=\frac{1}{\text{vol}_k(U)}$$ Note that by the chain rule $$\sqrt{\det (D\rho ^T D\rho)}=\sqrt{\det ((D\phi \circ \psi D\psi)^T(D\phi \circ \psi D\psi)}=\sqrt{\det (D\psi^T (D\phi \circ \psi)^T (D\phi \circ \psi) D\psi}$$ $$=|\det D\psi|\sqrt{\det ((D\phi \circ \psi)^T (D\phi \circ \psi)}=\frac{\text{vol}_k(M)}{\text{vol}_k(U)}$$ From here I think it's pretty obvious that generally you can't solve explicitly for $\psi$, At least that's my guess.

[Math] Generate random points on perimeter of ellipse

The following may be of some use in this question. (Note: Some of these points have also been sketched out in the above comments-included here for completeness). In particular, the below code computes the transformation based on the following derivation:

The points on the ellipse are assumed to have the coordinates defined by $$ x=a\cos{\theta} \\ y=b\sin{\theta} \\ $$

The arclength differential $\mathrm{d}s$ along the perimeter of the ellipse is obtained from

$$ {\mathrm{d}s}^{2}={\mathrm{d}x}^{2}+{\mathrm{d}y}^{2} $$

$$ {\mathrm{d}s}^{2}=a^{2}\sin^{2}{\theta}{\mathrm{d}\theta}^{2}+b^{2}\cos^{2}{\theta}{\mathrm{d}\theta}^{2} $$

$$ {\mathrm{d}s}^{2}=\left(a^{2}\sin^{2}{\theta}+b^{2}\cos^{2}{\theta}\right){\mathrm{d}\theta}^{2} $$

$$ {\mathrm{d}s}=\sqrt{a^{2}\sin^{2}{\theta}+b^{2}\cos^{2}{\theta}}{\mathrm{d}\theta} $$

$$ \frac{{\mathrm{d}s}}{\mathrm{d}\theta}=\sqrt{a^{2}\sin^{2}{\theta}+b^{2}\cos^{2}{\theta}} $$

Now, the probability function is taken to be

$$ p\left(\theta\right)=\frac{{\mathrm{d}s}}{\mathrm{d}\theta} $$

with the interpretation that when the rate-of-change of arclength increases, we want a higher probability of sample points in that interval to keep the density of points uniform.

We can then set up the following expression:

$$ p\left(\theta\right){\mathrm{d}\theta}=p\left(x\right){\mathrm{d}x} $$

and assuming a uniform distribution for $x$:

$$ \int p\left(\theta\right){\mathrm{d}\theta}=x+K $$.

Some plots of the uncorrected and corrected ellipses are shown in the below figure, using the above derivation and code implementation below. I hope this helps.

Python code below:

import math
import matplotlib.pyplot as plt

# ellipse major (a) and minor (b) axis parameters
a=8
b=2

# num points for transformation lookup function
npoints = 1000
delta_theta=2.0*math.pi/npoints

theta=[0.0]
delta_s=[0.0]
integ_delta_s=[0.0]

# integrated probability density
integ_delta_s_val=0.0

for iTheta in range(1,npoints+1):
    # ds/d(theta):
    delta_s_val=math.sqrt(a**2*math.sin(iTheta*delta_theta)**2+ \
                          b**2*math.cos(iTheta*delta_theta)**2)

    theta.append(iTheta*delta_theta)
    delta_s.append(delta_s_val)
    # do integral
    integ_delta_s_val = integ_delta_s_val+delta_s_val*delta_theta
    integ_delta_s.append(integ_delta_s_val)
    
# normalize integrated ds/d(theta) to make into a scaled CDF (scaled to 2*pi)
integ_delta_s_norm = []
for iEntry in integ_delta_s:
    integ_delta_s_norm.append(iEntry/integ_delta_s[-1]*2.0*math.pi)    
    
#print('theta= ', theta)
#print('delta_theta = ', delta_theta)
#print('delta_s= ', delta_s)
#print('integ_delta_s= ', integ_delta_s)
#print('integ_delta_s_norm= ', integ_delta_s_norm)

# Plot tranformation function
x_axis_range=1.5*math.pi
y_axis_range=1.5*math.pi
plt.xlim(-0.2, x_axis_range)
plt.ylim(-0.2, y_axis_range)
plt.plot(theta,integ_delta_s_norm,'+')
# overplot reference line which are the theta values.
plt.plot(theta,theta,'.')
plt.show()

# Reference ellipse without correction.
ellip_x=[]
ellip_y=[]

# Create corrected ellipse using lookup function
ellip_x_prime=[]
ellip_y_prime=[]

npoints_new=40
delta_theta_new=2*math.pi/npoints_new

for theta_index in range(npoints_new):
    theta_val = theta_index*delta_theta_new
#    print('theta_val = ', theta_val)
    
# Do lookup:
    for lookup_index in range(len(integ_delta_s_norm)):
#        print('doing lookup: ', lookup_index)
#        print('integ_delta_s_norm[lookup_index]= ', integ_delta_s_norm[lookup_index])
        if theta_val >= integ_delta_s_norm[lookup_index] and theta_val < integ_delta_s_norm[lookup_index+1]:
#            print('value found in lookup table')
            theta_prime=theta[lookup_index]
#            print('theta_prime = ', theta_prime)
#            print('---')
            break
    
    # ellipse without transformation applied for reference
    ellip_x.append(a*math.cos(theta_val))
    ellip_y.append(b*math.sin(theta_val))
    
    # ellipse with transformation applied
    ellip_x_prime.append(a*math.cos(theta_prime))
    ellip_y_prime.append(b*math.sin(theta_prime))

# Plot reference and transformed ellipses  
x_axis_range=10
y_axis_range=10
plt.xlim(-x_axis_range, x_axis_range)
plt.ylim(-y_axis_range, y_axis_range)
plt.gca().set_aspect('equal', adjustable='box')
plt.plot(ellip_x, ellip_y, '.')
plt.title('Uncorrected')
plt.show()

plt.xlim(-x_axis_range, x_axis_range)
plt.ylim(-y_axis_range, y_axis_range)
plt.gca().set_aspect('equal', adjustable='box')
plt.plot(ellip_x_prime, ellip_y_prime, '.')
plt.title('Corrected')
plt.show()
```

Best Answer

Related Solutions

Uniform Distribution of Points on a Surface in $\mathbb{R}^n$

[Math] Generate random points on perimeter of ellipse

Related Question