Explaining Quaternion Rotation Matrix

With a quaternion ${\bf q} = q_{0} + q_{1}i + q_{2}j + q_{3}k$, the rotation matrix is:

$$\begin{bmatrix}     q_{0}^{2}+q_{1}^{2}-q_{2}^{2}-q_{3}^{2} & 2q_{1}q_{2}-2q_{0}q_{3} & 2q_{1}q_{3}+2q_{0}q_{2} \\     2q_{1}q_{2}+2q_{0}q_{3} & q_{0}^{2}+q_{2}^{2}-q_{1}^{2}-q_{3}^{2} & 2q_{2}q_{3}-2q_{0}q_{1}\\     2q_{1}q_{3}-2q_{0}q_{2} & 2q_{2}q_{3}+2q_{0}q_{1} &  q_{0}^{2}+q_{3}^{2}-q_{1}^{2}-q_{2}^{2}  \\ \end{bmatrix}$$

A quaternion is determined by two factors: normalized rotation axis $(a_{x},a_{y},a_{z})$ and rotation angle $\theta$. ${\bf q} = cos(\theta/2) + (a_{x}i + a_{y}j + a_{z}k)sin(\theta/2)$ and ${\bf q^{-1}} = cos(\theta/2) - (a_{x}i + a_{y}j + a_{z}k)sin(\theta/2)$. This means that $||{\bf q}|| = 1$.

The rotation matrix in term of $(a_{x},a_{y},a_{z})$ and $\theta$ is:
$$\begin{bmatrix} 1+(a_{x}^{2}-1)(1-cos(\theta)) & -a_{z}sin(\theta)+a_{x}a_{y}(1-cos(\theta)) & a_{y}sin(\theta)+a_{x}a_{z}(1-cos(\theta))\\ a_{z}sin(\theta)+a_{x}a_{y}(1-cos(\theta)) & 1+(a_{y}^2-1)(1-cos(\theta)) & -a_{x}sin(\theta)+a_{y}a_{z}(1-cos(\theta))\\ -a_{y}sin(\theta)+a_{x}a_{z}(1-cos(\theta)) & a_{x}sin(\theta)+a_{y}a_{z}(1-cos(\theta)) & 1+(a_{z}^{2}-1)(1-cos(\theta)) \end{bmatrix}$$

From the expression of quaternion, we have
$$cos(\theta) = cos^{2}(\theta/2)-sin^{2}(\theta/2) =q_{0}^2-(q_{1}^2+q_{2}^2+q_{3}^2)$$
$$1-cos(\theta) = 1-cos^{2}(\theta/2)+sin^{2}(\theta/2) =2sin^{2}(\theta/2)$$
$$(a_{x},a_{y},a_{z})=(q_{1},q_{2},q_{3})/sin(\theta/2)$$

Based on these equations, we can derive the quaternion rotation matrix. For example,
$$1+(a_{x}^{2}-1)(1-cos(\theta)) = cos(\theta)+a_{x}^{2}(1-cos(\theta))\\ =q_{0}^2-(q_{1}^2+q_{2}^2+q_{3}^2)+2q_{1}^{2}\\ = q_{0}^2+q_{1}^2-q_{2}^2-q_{3}^2$$
$$-a_{z}sin(\theta)+a_{x}a_{y}(1-cos(\theta)) = \frac{-q_{3}}{sin(\theta/2)}sin(\theta)+\frac{q_{1}q_{2}}{sin^{2}(\theta/2)}2sin^{2}(\theta/2)\\ =2q_{1}q_{2}-2q_{0}q_{3}$$