Generating uniformly distributed points within rotated ellipsoid

Question

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 = fig.add_subplot(111, projection='3d')
 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 = fig.add_subplot(111, projection='3d')
 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()

Answer 1

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.

Answer 2

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