Let $\mathcal{A} \subset \mathcal{M}_{3\times 3}(\mathbb{R})$ be the space of $3 \times 3$ real skew symmetric matrices. It is span by three
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}
$$
For each $A \in \mathcal{A}$, we can expand $A$ as $x L_x + y L_y + z L_z$ and
associated a vector $\vec{A} \in \mathbb{R}^3$ and a quaternion $\tilde{A} \in \mathbb{H}$ to it.

$$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}$$

Let $\mathcal{A}_u = \{\; A \in \mathcal{A} : |\vec{A}| = 1\; \}$.
Given any rotation matrix $M \in SO(3)$, we can find a $\theta \in [0,\pi]$ and
$L \in \mathcal{A}_u$ such that

$$M = e^{\theta L} = I_3 + \sin\theta L + (1-\cos\theta) L^2$$

The $\theta$ is the angle of rotation associated with $M$ and $\vec{L}$ will be a unit vector in the direction of the rotational axis.

When $\theta \ne 0$ nor $\pi$, we can extract $L$ out by taking the antisymmetric part
of $M$ and then "normalize" the corresponding vector because

$$ \sin\theta L = \frac12 \left(M - M^T\right)$$

To fix the value of $\theta$, we can use the relation $\text{Tr}(M) = 1 + 2\cos\theta$.

Once $\theta$ and $L$ is known, the quaternion corresponding to the rotation matrix $M$ is then given by

$$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]
$$