Find the rotation angles in $x, y$ and $z$ to convert between axis systems I have 2 axis systems. One is just the natural axis system:
$e_{x} = \{1,0,0\}$,
$e_{y} = \{0,1,0\}$,
$e_{z} = \{0,0,1\}$
The other axis system is defined by 3 orthogonal vectors: $a$, $b$ and $c$.
What is the simplest method to find the Euler angles for this conversion, such that I can convert between axis systems with using a rotation matrix?
 A: You don't need the Euler angles to obtain the rotation matrix (and its inverse) between the two systems.  (There is no reason to work with Euler angles; they will only confuse you.  Versors, or unit quaternions, are both better and easier to use to describe orientations.)
Let's say the other coordinate system basis — the three orthogonal unit vectors corresponding to the axes — are
$$\begin{aligned}
\hat{u} &= ( x_u , y_u , z_u ) \\
\hat{v} &= ( x_v , y_v , z_v ) \\
\hat{w} &= ( x_w , y_w , z_w ) \end{aligned}$$
Note that these are unit vectors, i.e. $\left\lVert \hat{u} \right\rVert = 1$, $\left\lVert\hat{v}\right\rVert = 1$, $\left\lVert\hat{w}\right\rVert = 1$. If your orthogonal vectors are of arbitrary (but nonzero) length, then simply normalize them. If you have some vector $\vec{a}$,
$$\vec{a} = \left [ \begin{matrix} x_a \\ y_a \\ z_a \end{matrix} \right ] = ( x_a , y_a , z_a )$$
divide its components by its length, to normalize them:
$$\hat{a} = \frac{\vec{a}}{\left\lVert\vec{a}\right\rVert} = \frac{\vec{a}}{\sqrt{\vec{a}\cdot\vec{a}}} = \left [ \begin{matrix} 
\frac{x_a}{\sqrt{x_a^2 + y_a^2 + z_a^2}} \\
\frac{y_a}{\sqrt{x_a^2 + y_a^2 + z_a^2}} \\
\frac{z_a}{\sqrt{x_a^2 + y_a^2 + z_a^2}} \end{matrix} \right ]$$
The rotation matrix $\mathbf{R}$ to rotate $\vec{p}$ in the natural coordinate system to $\vec{q}$ in the other coordinate system, via $\vec{q} = \mathbf{R} \vec{p}$, is
$$\mathbf{R} = \left [ \begin{matrix} \hat{u} & \hat{v} & \hat{w} \end{matrix} \right ] = \left [ \begin{matrix}
x_u & x_v & x_w \\
y_u & y_v & y_w \\
z_u & z_v & z_w \end{matrix} \right ]$$
As a pure rotation matrix, $\mathbf{R}$ is an orthonormal matrix, and its inverse is its transpose,
$$\mathbf{R}^{-1} = \mathbf{R}^T = \left [ \begin{matrix}
x_u & y_u & z_u \\
x_v & y_v & z_v \\
x_w & y_w & z_w \end{matrix} \right ]$$
This latter matrix rotates $\vec{q}$ back to $\vec{p}$ in the original coordinate system. In other words:
$$\vec{q} = \mathbf{R} \vec{p} \quad \iff \quad \vec{p} = \mathbf{R}^T \vec{q}$$
