How to generate a random sample of points from a 3-D ellipsoid using Python?

Question

I am trying to sample around 1000 points from a 3-D ellipsoid, uniformly. Is there some way to code it such that we can get points starting from the equation of the ellipsoid?

I want points on the surface of the ellipsoid.

Answer 1

Depending on what "uniformly" refers to, different methods are applicable. In any case, we can use the parametric equations using spherical coordinates (from Wikipedia):

where s = 1 refers to the ellipsoid given by the semi-axes a > b > c. From these equations we can derive the relevant volume/area element and generate points such that their probability of being generated is proportional to that volume/area element. This will provide constant volume/area density across the surface of the ellipsoid.

1. Constant volume density

This method generates points on the surface of an ellipsoid such that their volume density across the surface of the ellipsoid is constant. A consequence of this is that the one-dimensional projections (i.e. the x, y, z coordinates) are uniformly distributed; for details see the plot below.

The volume element for a triaxial ellipsoid is given by (see here):

and is thus proportional to sin(theta) (for 0 <= theta <= pi). We can use this as the basis for a probability distribution that indicates "how many" points should be generated for a given value of theta: where the area density is low/high, the probability for generating a corresponding value of theta should be low/high, too.

Hence, we can use the function f(theta) = sin(theta)/2 as our probability distribution on the interval [0, pi]. The corresponding cumulative distribution function is F(theta) = (1 - cos(theta))/2. Now we can use Inverse transform sampling to generate values of theta according to f(theta) from a uniform random distribution. The values of phi can be obtained directly from a uniform distribution on [0, 2*pi].

Example code:

import matplotlib.pyplot as plt
import numpy as np
from numpy import sin, cos, pi


rng = np.random.default_rng(seed=0)

a, b, c = 10, 3, 1
N = 5000

phi = rng.uniform(0, 2*pi, size=N)
theta = np.arccos(1 - 2*rng.random(size=N))

x = a*sin(theta)*cos(phi)
y = b*sin(theta)*sin(phi)
z = c*cos(theta)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(x, y, z, s=2)
plt.show()

which produces the following plot:

The following plot shows the one-dimensional projections (i.e. density plots of x, y, z):

import seaborn as sns

sns.kdeplot(data=dict(x=x, y=y, z=z))
plt.show()

2. Constant area density

This method generates points on the surface of an ellipsoid such that their area density is constant across the surface of the ellipsoid.

Again, we start by calculating the corresponding area element. For simplicity we can use SymPy:

from sympy import cos, sin, symbols, Matrix

a, b, c, t, p = symbols('a b c t p')

x = a*sin(t)*cos(p)
y = b*sin(t)*sin(p)
z = c*cos(t)

J = Matrix([
    [x.diff(t), x.diff(p)],
    [y.diff(t), y.diff(p)],
    [z.diff(t), z.diff(p)],
])
print((J.T @ J).det().simplify())

This yields

-a**2*b**2*sin(t)**4 + a**2*b**2*sin(t)**2 + a**2*c**2*sin(p)**2*sin(t)**4 - b**2*c**2*sin(p)**2*sin(t)**4 + b**2*c**2*sin(t)**4

and further simplifies to (dividing by (a*b)**2 and taking the sqrt):

sin(t)*np.sqrt(1 + ((c/b)**2*sin(p)**2 + (c/a)**2*cos(p)**2 - 1)*sin(t)**2)

Since for this case the area element is more complex, we can use rejection sampling:

import matplotlib.pyplot as plt
import numpy as np
from numpy import cos, sin


def f_redo(t, p):
    return (
        sin(t)*np.sqrt(1 + ((c/b)**2*sin(p)**2 + (c/a)**2*cos(p)**2 - 1)*sin(t)**2)
        < rng.random(size=t.size)
    )


rng = np.random.default_rng(seed=0)
N = 5000
a, b, c = 10, 3, 1

t = rng.uniform(0, np.pi, size=N)
p = rng.uniform(0, 2*np.pi, size=N)

redo = f_redo(t, p)
while redo.any():
    t[redo] = rng.uniform(0, np.pi, size=redo.sum())
    p[redo] = rng.uniform(0, 2*np.pi, size=redo.sum())
    redo[redo] = f_redo(t[redo], p[redo])

x = a*np.sin(t)*np.cos(p)
y = b*np.sin(t)*np.sin(p)
z = c*np.cos(t)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(x, y, z, s=2)
plt.show()

which yields the following distribution:

The following plot shows the corresponding one-dimensional projections (x, y, z):

Answer 2

One way of doing this whch generalises for any shape or surface is to convert the surface to a voxel representation at arbitrarily high resolution (the higher the resolution the better but also the slower). Then you can easily select the voxels randomly however you want, and then you can select a point on the surface within the voxel using the parametric equation. The voxel selection should be completely unbiased, and the selection of the point within the voxel will suffer the same biases that come from using the parametric equation but if there are enough voxels then the size of these biases will be very small.

You need a high quality cube intersection code but with something like an elipsoid that can optimised quite easily. I'd suggest stepping through the bounding box subdivided into voxels. A quick distance check will eliminate most cubes and you can do a proper intersection check for the ones where an intersection is possible. For the point within the cube I'd be tempted to do something simple like a random XYZ distance from the centre and then cast a ray from the centre of the elipsoid and the selected point is where the ray intersects the surface. As I said above, it will be biased but with small voxels, the bias will probably be small enough.

There are libraries that do convex shape intersection very efficiently and cube/elipsoid will be one of the options. They will be highly optimised but I think the distance culling would probably be worth doing by hand whatever. And you will need a library that differentiates between a surface intersection and one object being totally inside the other.

And if you know your elipsoid is aligned to an axis then you can do the voxel/edge intersection very easily as a stack of 2D square intersection elipse problems with the set of squares to be tested defined as those that are adjacent to those in the layer above. That might be quicker.

One of the things that makes this approach more managable is that you do not need to write all the code for edge cases (it is a lot of work to get around issues with floating point inaccuracies that can lead to missing or doubled voxels at the intersection). That's because these will be very rare so they won't affect your sampling.

It might even be quicker to simply find all the voxels inside the elipse and then throw away all the voxels with 6 neighbours... Lots of options. It all depends how important performance is. This will be much slower than the opther suggestions but if you want ~1000 points then ~100,000 voxels feels about the minimum for the surface, so you probably need ~1,000,000 voxels in your bounding box. However even testing 1,000,000 intersections is pretty fast on modern computers.

Answer 3

Theory

Using this excellent answer to the MSE question How to generate points uniformly distributed on the surface of an ellipsoid? we can

generate a point uniformly on the sphere, apply the mapping f : (x,y,z) -> (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).

Assuming that the 3 axes of the ellipsoid are named such that

0 < a < b < c

We discard a generated point with

p(x,y,z) = 1 - mu(x,y,y)/mu_max

probability, ie we keep it with mu(x,y,y)/mu_max probability where

mu(x,y,z) = ((acy)^2 + (abz)^2 + (bcx)^2)^0.5

and

mu_max = bc

---

Implementation

import numpy as np
np.random.seed(42)

# Function to generate a random point on a uniform sphere
# (relying on https://stackoverflow.com/a/33977530/8565438)

def randompoint(ndim=3):
    vec = np.random.randn(ndim,1)
    vec /= np.linalg.norm(vec, axis=0)
    return vec

# Give the length of each axis (example values):

a, b, c = 1, 2, 4

# Function to scale up generated points using the function `f` mentioned above:

f = lambda x,y,z : np.multiply(np.array([a,b,c]),np.array([x,y,z]))

# Keep the point with probability `mu(x,y,z)/mu_max`, ie

def keep(x, y, z, a=a, b=b, c=c):
    mu_xyz = ((a * c * y) ** 2 + (a * b * z) ** 2 + (b * c * x) ** 2) ** 0.5
    return mu_xyz / (b * c) > np.random.uniform(low=0.0, high=1.0)

# Generate points until we have, let's say, 1000 points:

n = 1000
points = []
while len(points) < n:
    [x], [y], [z] = randompoint()
    if keep(x, y, z):
        points.append(f(x, y, z))

---

Checks

Check if all points generated satisfy the ellipsoid condition (ie that x^2/a^2 + y^2/b^2 + z^2/c^2 = 1):

for p in points:
    pscaled = np.multiply(p,np.array([1/a,1/b,1/c]))
    assert np.allclose(np.sum(np.dot(pscaled,pscaled)),1)

Runs without raising any errors. Visualize the points:

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
points = np.array(points)
ax.scatter(points[:, 0], points[:, 1], points[:, 2])
# set aspect ratio for the axes using https://stackoverflow.com/a/64453375/8565438
ax.set_box_aspect((np.ptp(points[:, 0]), np.ptp(points[:, 1]), np.ptp(points[:, 2])))
plt.show()

These points seem evenly distributed.

---

Problem with currently accepted answer

Generating a point on a sphere and then just reprojecting it without any further corrections to an ellipse will result in a distorted distribution. This is essentially the same as setting this posts's p(x,y,z) to 0. Imagine an ellipsoid where one axis is orders of magnitude bigger than another. This way, it is easy to see, that naive reprojection is not going to work.

Answer 4

J.F. Williamson, "Random selection of points distributed on curved surfaces", Physics in Medicine & Biology 32(10), 1987, describes a general method of choosing a uniformly random point on a parametric surface. It is an acceptance/rejection method that accepts or rejects each candidate point depending on its stretch factor (norm-of-gradient). To use this method for a parametric surface, several things have to be known about the surface, namely—

• x(u, v), y(u, v) and z(u, v), which are functions that generate 3-dimensional coordinates from two dimensional coordinates u and v,

• The ranges of u and v,

• g(point), the norm of the gradient ("stretch factor") at each point on the surface, and

• gmax, the maximum value of g for the entire surface.

The algorithm is then:

• Generate a point on the surface, xyz.
• If g(xyz) >= RNDU01()*gmax, where RNDU01() is a uniform random variate in [0, 1), accept the point. Otherwise, repeat this process.

Chen and Glotzer (2007) apply the method to the surface of a prolate spheroid (one form of ellipsoid) in "Simulation studies of a phenomenological model for elongated virus capsid formation", Physical Review E 75(5), 051504 (preprint).

Answer 5

Here is a generic function to pick a random point on a surface of a sphere, spheroid or any triaxial ellipsoid with a, b and c parameters. Note that generating angles directly will not provide uniform distribution and will cause excessive population of points along z direction. Instead, phi is obtained as an inverse of randomly generated cos(phi).

    import numpy as np
    def random_point_ellipsoid(a,b,c):
        u = np.random.rand()
        v = np.random.rand()
        theta = u * 2.0 * np.pi
        phi = np.arccos(2.0 * v - 1.0)
        sinTheta = np.sin(theta);
        cosTheta = np.cos(theta);
        sinPhi = np.sin(phi);
        cosPhi = np.cos(phi);
        rx = a * sinPhi * cosTheta;
        ry = b * sinPhi * sinTheta;
        rz = c * cosPhi;
        return rx, ry, rz

This function is adopted from this post: https://karthikkaranth.me/blog/generating-random-points-in-a-sphere/

Answer 6

Consider using Monte-Carlo simulation: generate a random 3D point; check if the point is inside the ellipsoid; if it is, keep it. Repeat until you get 1,000 points.

P.S. Since the OP changed their question, this answer is no longer valid.