Any rotation in 3D space can be viewed as three successive rotations about predefined axes, giving the Euler or Tait-Bryan angles). It can also be viewed as a pair of rotations, one about a predefined axis (the roll axis) and another that aligns the new roll axis to the original. In the aeronautical context, the angles of these rotations are termed the total angle of attack and the bank angle. Every rotation can be described by a unit quaternion, so the question of how to find these two rotations (as quaternions) from a given rotation (also as a quaternion) is one of factorizing a quaternion into two quaternions with some specific relationship. This factorization is the swing-twist decomposition. So named because the object “swings” to align the new and old roll axes, and then “twists” about the aligned axes to restore the original orientation. Personally I don’t like this name, because “twist” doesn’t make clear the fact that one particular axis is being chosen. “Swing-roll” would be better, and if you’re willing to think of the swing motion as akin to that of a rocking chair then you could call it the “rock-and-roll” decomposition. But I digress.
The defining feature of the “swing” motion is that it is a rotation about an axis that is perpendicular to the roll axis. So given our original quaternion $q$ and a unit vector $\hat{r}$ for the roll axis, we want to find unit quaternions $s = s_0 + \vec{s}$ and $r = r_0 + \vec{r}$ such that:
$$\displaystyle \begin{aligned} q &= sr \\ \vec{s}\cdot\hat{r} &= 0\end{aligned} $$
Since the roll axis direction is fixed, we can write $\vec{r} = \lambda\hat{r}$ and solve for the scalar $\lambda$ that satisfies the orthogonality condition. Since $q$ is given, we can write $s$ in terms of $q$ and the mostly-unknown $r$:
$$\displaystyle s= s_0 + \vec{s} = qr^{-1} = (q_0 + \vec{q})(r_0 – \lambda\hat{r}) $$
So we can find the vector part $\vec{s}$ of $s$ by multiplying out the quaternions on the right-hand side:
$$\displaystyle \vec{s} = -q_0\lambda\hat{r} + r_0\vec{q} – \lambda\vec{q}\times\hat{r} $$
We require $\vec{s}$ to be perpendicular to $\hat{r}$ so:
$$\displaystyle 0 = \vec{s}\cdot\hat{r} = -q_0\lambda + r_0\vec{q}\cdot\hat{r} $$
The cross product is orthogonal to both of its input vectors, one of which is $\hat{r}$ already. We can now put $\lambda$ in terms of $r_0$ and the known quantities, but since $r$ is a unit quaternion we know that:
$$\displaystyle r_0^2 + \lambda^2 = 1 $$
so we can write:
$$\displaystyle q_0\lambda = r_0\vec{q}\cdot\hat{r} = \pm\sqrt{1-\lambda^2}(\vec{q}\cdot\hat{r}) $$
It is convenient to square both sides:
$$\displaystyle q_0^2 \lambda^2 = (1-\lambda^2)(\vec{q}\cdot\hat{r})^2 $$
from which we get:
$$\displaystyle \begin{aligned} \lambda^2 &= \frac{(\vec{q}\cdot\hat{r})^2}{q_0^2 + (\vec{q}\cdot\hat{r})^2} \\ \lambda &= \pm \frac{\vec{q}\cdot\hat{r}}{\sqrt{q_0^2 + (\vec{q}\cdot\hat{r})^2}} \\ r_0^2 &= 1-\lambda^2 = \frac{q_0^2}{q_0^2 + (\vec{q}\cdot\hat{r})^2} \\ r_0 &= \pm\frac{q_0}{\sqrt{q_0^2 + (\vec{q}\cdot\hat{r})^2}} \end{aligned}$$
The fact from earlier that $q_0\lambda = r_0\vec{q}\cdot\hat{r}$ means that the $\pm$ for $\lambda$ and for $r_0$ must be the same so the quaternion $r$ is:
$$\displaystyle r = r_0 + \vec{r} = r_0 + \lambda\hat{r} = \pm\frac{q_0 + (\vec{q}\cdot\hat{r})\hat{r}}{\sqrt{q_0^2 + (\vec{q}\cdot\hat{r})^2}} $$
The denominator is just the norm of the numerator, so this is indeed a unit quaternion. The $\pm$ remains because we could have either $q = sr$ or $q = (-s)(-r)$ and which one you prefer depends on the sign convention you want for your angles. Recall that if $r$ represents a rotation through angle $\phi$ about the $\hat{r}$ axis, then $-r$ is a rotation through $\phi + \pi$ (radians) about the same $\hat{r}$ axis.If you need your $\phi$ to be within $[0, \pi]$, for example, you might need to choose $-r$ instead of $+r$ and this choice will affect the $s$ factor also.
To sum up, the swing-twist decomposition of a quaternion $q$ is a factorization $q = sr$ where the axis of $r$ is specified and orthogonal to that of $s$, which is not specified. It can be found by forming $r$ as
$$\displaystyle r = \frac{v}{\| v\| }$$
where $v$ is a quaternion formed from $q$ by:
$$\displaystyle v = q_0 + (\vec{q}\cdot\hat{r})\hat{r} $$
and $\|v\|$ is its norm. Then the $s$ factor is found by $s = qr^{-1}$.
The $r$ factor is the component of $q$ in the $\hat{r}$ direction, loosely speaking, so it is not surprising that $(\vec{q}\cdot\hat{r})\hat{r}$, the component of $q$ along the axis vector is involved.