Dejemos que $\mathcal{A} \subset \mathcal{M}_{3\times 3}(\mathbb{R})$ sea el espacio de $3 \times 3$ matrices simétricas reales. Es un tramo de tres matrices
$$L_x = \begin{bmatrix}0 & 0 & 0\\0 & 0 & -1\\0 & 1 & 0\end{bmatrix},\quad L_y = \begin{bmatrix}0 & 0 & 1\\0 & 0 & 0\\-1& 0 & 0\end{bmatrix},\quad L_z = \begin{bmatrix}0 & -1 & 0\\1 & 0 & 0\\0 & 0 & 0\end{bmatrix} $$ Para cada $A \in \mathcal{A}$ podemos ampliar $A$ como $x L_x + y L_y + z L_z$ y asociado un vector $\vec{A} \in \mathbb{R}^3$ y un cuaternión $\tilde{A} \in \mathbb{H}$ a ella.
$$A = x L_x + y L_y + z L_z \quad\leftrightarrow\quad \vec{A} = (x,y,z) \quad\leftrightarrow\quad \tilde{A} = x{\bf i} + y{\bf j} + z{\bf k}$$
Dejemos que $\mathcal{A}_u = \{\; A \in \mathcal{A} : |\vec{A}| = 1\; \}$ . Dada cualquier matriz de rotación $M \in SO(3)$ podemos encontrar un $\theta \in [0,\pi]$ y $L \in \mathcal{A}_u$ tal que
$$M = e^{\theta L} = I_3 + \sin\theta L + (1-\cos\theta) L^2$$
El $\theta$ es el ángulo de rotación asociado a $M$ y $\vec{L}$ será un vector unitario en la dirección del eje de rotación.
Cuando $\theta \ne 0$ ni $\pi$ podemos extraer $L$ tomando la parte antisimétrica de $M$ y luego "normalizar" el vector correspondiente porque
$$ \sin\theta L = \frac12 \left(M - M^T\right)$$
Para fijar el valor de $\theta$ podemos utilizar la relación $\text{Tr}(M) = 1 + 2\cos\theta$ .
Una vez $\theta$ y $L$ es conocida, el cuaternión correspondiente a la matriz de rotación $M$ viene dada por
$$e^{\frac{\theta}{2} \tilde{L}} = \cos\frac{\theta}{2} + \sin\frac{\theta}{2} \tilde{L} = \frac{\sqrt{1+\text{Tr}(M)}}{2}\left[ 1 + \frac{\widetilde{M - M^{T}}}{1 + \text{Tr}(M)}\right] $$