Updating a quaternion orientation by a vector of euler angles

I'm trying to understand why this formula works to update an orientation with an angular velocity represented as a vector of rotations in $radians/{second}$.

I understand that two quaternions multiplied together will form a corresponding quaternion with each rotation applied, but what I don't understand is

• What exactly a pure imaginary quaternion multiplied with an orientation represents
• How adding this quaternion to the original orientation results in updating the original quaternion by the rotations stored in the angular velocity

// angular velocity represented by a vector of radians
$angular_{velocity}= vector(x_{axis\ rotation}, y_{axis\ rotation}, z_{axis\ rotation})$

// quaternion(w, x, y, z) or quaternion(r, i, j, k)
$w_{quat} = quaternion(0, angular.x, angular.y, angular.z)$

// t = time in seconds
$new\_orientation = orientation + {\Delta t\over2} (w_{quat} * orientation)$

$new\_orientation = normalize(new\_orientation)$

The space of unit quaternions is isomorphic to $SO(3)$ (the space of possible rotations in 3D space), plus a double cover (which is not important in this context). It's easier to picture a sphere for visualization purposes, though. Now imagine you want to change the quaternion vector (which is a vector from the origin to the surface of the (hyper-)sphere). Since you want to keep the length of the vector the same, you have to move it orthogonally to its direction (for a sphere, this means up/down and left/right, but not in/out of the sphere). Thus, you need a vector whose dot product with the $orientation$ quaternion is zero.

But pure imaginary quaternions are orthogonal to pure real ones, and the real axis represents the identity rotation (no rotation at all). Additionally, if $q_1$ and $q_2$ are quaternions, $u$ is a unit quaternion, and $q_1\cdot q_2=0$, then $(q_1u)\cdot(q_2u)=0$ as well (this is because under the interpretation of quaternions as rotations, they are isometries, so they preserve distance and angles), so $orientation\cdot (q_i\ orientation)=0$, if $q_i$ is any pure-imaginary quaternion.

It's a bit more work that the specific choice $q_i=\frac{\Delta t}2(\omega_x{\bf i}+\omega_y{\bf j}+\omega_z{\bf k})$ gives you the amount and axis of rotation you want in your mapping from quaternions to actual rotations, but accepting that, the sum $orientation+(q_i\ orientation)$ gives you a new vector, which is the same length as the old one (to first order), but pointing in a different direction, which is what you wanted. Since there is a second order term that will tend to increase the distance to the the center (which should be intuitively obvious; imagine moving tangent to a sphere from a point on the surface—you stay near the surface for a while, but as you continue straight, you will drift away from the surface), you presumably want a step after those shown to renormalize the vector back onto the sphere.

The answer I just gave is more a programmer's answer than a mathematician's, since your question seemed quite concrete and grounded in some numerical computation. That said, the math angle here has to do with the isomorphism from the space of unit quaternions to the group $\operatorname{Spin}(3)$, which is the simply connected double cover of $SO(3)$ (as I mentioned. Moreover, the Lie algebra $\mathfrak{so}(3)$ generated by (either of) these groups is exactly the space of skew-symmetric $3\times3$ matrices $[\omega]_\times$ for all possible 3-vectors $\omega$, which is exactly what you are doing with your own angular velocity vector. (Note that $[\omega]_\times$ is the 3D Hodge star, defined here.)

Let $u$ be some unit, pure imaginary quaternion.

A quaternion $q$ can be expressed as the exponential of some pure imaginary quaternion. That is, $q=e^{u\theta/2}$, so that a vector (which is also a pure imaginary quaternion) $r$ rotates to

$$r' = qrq^{-1} = e^{u\theta/2} r e^{-u\theta/2}$$

corresponds to rotating the vector $r$ by $\theta$ about the axis defined by $u$. This correspondence is key; otherwise, we could build quaternions to have almost any form we want, but to get the right correspondence, there is this need to use half-angles in the quaternions to get the full-angle rotation.

Now, for our convenience, we'll call $u\theta = \Theta$, the (vector) angular orientation. Then, the time derivative of $q$ is given by

$$\frac{dq}{dt} = \frac{d}{dt} e^{\Theta/2} = \frac{1}{2} \frac{d\Theta}{dt} e^{\Theta/2} = \frac{\omega}{2} q$$

where $\omega = d\Theta/dt$.

From here, you're just using $q(t + \Delta t) \approx q(t) + q'(t)\Delta t$. The key is, why 1/2 in our exponential form of $q$? That's a question a bit larger in scope, but the need to use half angles in building up quaternions is well-established.