# Generating uniformly distributed points within rotated ellipsoid

I want to generate uniformly distributed points within an ellipsoid that has been rotated. I have the vectors for the major, minor and intermediate axes of the rotated ellipsoid. The ellipsoid is centered at the origin. There are two problems in my opinion:

1) Finding the rotation matrix that will move the major axis from its original position (say aligned along z) onto the new major axis vector

2) Rewriting the equation of the ellipsoid as

$\vec{x} ^{T} R^{T}AR \vec{x} = 1$

and finding x,y and z so that the LHS $\leq 1$

I have attempted finding the rotation matrix this way, https://math.stackexchange.com/a/476311/475188 , where I tried to align the major axis (along z) with some arbitrary vector (both normalized) and inserting the obtained matrix into equation above, but this did not seem to work. Am I missing something or is there a simpler way?

Edit: following G Cab's answer, I've tested the solution on a unit sphere in Python, this should be returning another unit sphere, but the shape looks off. Any ideas?

 import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

axes = np.array([[0.83977192, -0.39189744, -0.37576525],
[0.481634,    0.8571836,   0.18238689],
[0.25062285,  -0.3341447,  0.90858984]])

#Orthogonal
print np.dot(axes[0], axes[1])
print np.dot(axes[1], axes[2])
print np.dot(axes[0], axes[2])
print axes
print np.cross(axes[0],axes[1])

#Defining the matrix
r_m = np.zeros((3,3))

r_m[:,0] = axes[0,:]
r_m[:,1] = axes[1,:]
r_m[:,2] = axes[2,:]

#Taking transpose
r_m = r_m.T

x1 = np.random.uniform(-1,1,10000)
y1 = np.random.uniform(-1,1,10000)
z1 = np.random.uniform(-1,1,10000)

#Unit sphere
a = 1.
b = 1.
c = 1.

r = x1**2. + y1**2. + z1**2.
ch, = np.where(r <= 1)

x = x1[ch]*a
y = y1[ch]*b
z = z1[ch]*c

#Plotting original sphere
fig = plt.figure()
ax.set_xlim([-1,1])
ax.set_ylim([-1,1])
ax.set_zlim([-1,1])
ax.scatter(x1[ch],y1[ch],z1[ch], s= 0.5, color = 'red')
plt.show()

#Matrix multiplying
x = x*r_m[0,0] + y*r_m[0,1] + z*r_m[0,2]
y = x*r_m[1,0] + y*r_m[1,1] + z*r_m[1,2]
z = x*r_m[2,0] + y*r_m[2,1] + z*r_m[2,2]

#Shape of output sphere looks weird
fig = plt.figure()
ax.set_xlim([-1,1])
ax.set_ylim([-1,1])
ax.set_zlim([-1,1])
ax.scatter(x,y,z, s= 0.5, color = 'blue')
plt.show()

• Given a symmetric $A$, you need to find $R$ such that $R^\top A R$ is diagonal. Is this correct? If so lookup diagonalization with eigenvectors. – John Alexiou Aug 27 '17 at 14:58

We have two questions here:
1) the equation of the ellipsoid in the given rotated axes, and with the given semi-axes;
2) how to evenly randomly distribute points in its interior.

Question 1

Once you have the three (normalized) vectors $\bf v_1, \bf v_2, \bf v_3$, corresponding to the three ellipsoid's axes, and you construct a matrix with the vectors stacked vertically $${\bf C} = \left( {\matrix{ {v_{1,\,x} } & {v_{2,\,x} } & {v_{3,\,x} } \cr {v_{1,\,y} } & {v_{2,\,y} } & {v_{3,\,y} } \cr {v_{1,\,z} } & {v_{2,\,z} } & {v_{3,\,z} } \cr } } \right)$$ Then that is your conversion (rotation) matrix, in the sense that \eqalign{ & \left( {{\bf i'},\;{\bf j'},\;{\bf k'}} \right) = {\bf C}\left( {{\bf i},\;{\bf j},\;{\bf k}} \right) \cr & {\bf P} = \left( {\matrix{ {x'} \cr {y'} \cr {z'} \cr } } \right) = {\bf C^T}\left( {\matrix{ x \cr y \cr z \cr } } \right) \cr}

Now your ellipsoid has the equation $${{x'^{\,2} } \over {a^{\,2} }} + {{y'^{\,2} } \over {b^{\,2} }} + {{z'^{\,2} } \over {c^{\,2} }} \le 1$$ and once you distribute your points in it, you can convert from $(x',y',z')$ to $(x,y,z)$, by inverting the relation given before.

Concerning the second question, I suppose you know how to distribute the points in an ellipsoid, canonically oriented.

• The OP is using ${\bf A}$ for the ellipsoid coefficients and $\mathtt{R}$ for the rotation. I hope they don't get confused with your notation. – John Alexiou Aug 27 '17 at 15:19
• I have tried implementing the solution in numpy (please see edits), testing on a simple unit sphere. It appears the resultant shape is a bit off, have I made a mistake? – Anna Aug 27 '17 at 15:58
• @ja72: right remark, thanks, I changed notation – G Cab Aug 27 '17 at 22:14
• @Anna. I cannot follow you through the python code (I am not familiar with it). Put down the process you want to follow in mathematical terms, first of coding it. – G Cab Aug 27 '17 at 22:28

When you say "I have the vectors for the ... axes of the rotated ellipsoid", do they have unit length? And are they truly orthogonal? Yes? Good: Take the 3 vectors, make them have unit length if necessary, assemble them into matrix. That matrix, or its transpose, is your rotation matrix. I'd generate my points in the unit ball, then multiply the $x$, $y$, and $z$ coordinates by the lengths of the axes, to get uniform points in a non-rotated ellipsoid of the right shape. Finally, I'd multiply them by that matrix (or its transpose).

• When you stretch the points from the ball to the ellipse you loose uniform distribution points will be less dense along the semimajor axis the along the minor axis. – John F. Miller Oct 23 '17 at 4:08
• No. What you say is true of surface area of the ellipsoid but not volume. That is to say, the jacobian of the transformation is a constant. – kimchi lover Oct 23 '17 at 11:49