[Math] Quaternion Differentiation

ordinary differential equationsquaternions

I have an application that tracks an image and estimates its position and orientation. The orientation is given by a quaternion, and it is modified by an angular velocity every frame.

To predict the orientation I calculate the differential quaternion basing on the angular rate $\vec \omega$ and the previous quaternion $\vec q$. I found these equations.

$$q_x=\frac{1}{2}(w_x q_w+w_y q_z-w_z q_y) $$
$$q_y=\frac{1}{2}(w_y q_w+w_z q_x-w_x q_z) $$
$$q_z=\frac{1}{2}(w_z q_w+w_x q_y-w_y1 q_x)$$
$$q_w=-\frac{1}{2}(w_x q_x+w_y q_y-w_z q_z)$$

Is this approach correct? Should I use $\vec \omega$ or do I need to take into account the time interval between frames $\vec \omega\Delta t $?

After this, the predicted quaternion $\hat q$ would be the sum of the previous one and the differentiation, wouldn't it?

Best Answer

The derivative of a function $q: \mathbb{R} \rightarrow \mathbb{H}$ representing a rotation with angular velocity $\omega: \mathbb{R} \rightarrow \mathbb{R}^3$ around some axis by unit quaternions is

$$ \frac{\mathrm{d}}{\mathrm{d}t}q(t) = \frac{1}{2} * (\omega_x(t) \mathrm{i} + \omega_y(t) \mathrm{j} + \omega_z(t) \mathrm{k}) * q(t),$$

where $q(t) = q_w(t) + q_x(t) \mathrm{i} + q_y(t) \mathrm{j} + q_z(t) \mathrm{k}$, $\omega(t) = \left(\begin{smallmatrix} \omega_x(t) \\ \omega_y(t) \\ \omega_z(t) \end{smallmatrix}\right)$ and $*$ represents the quaternion multiplication operator. A derivation can be found at Euclidean Space. By evaluating the quaternion multiplications the derivatives of the real and imaginary components of $q(t)$ can be identified to be

$$ \left(\begin{matrix}\frac{\mathrm{d}}{\mathrm{d}t} q_w(t) \\ \frac{\mathrm{d}}{\mathrm{d}t} q_x(t) \\ \frac{\mathrm{d}}{\mathrm{d}t} q_y(t) \\ \frac{\mathrm{d}}{\mathrm{d}t} q_z(t) \end{matrix}\right) = \frac{1}{2} \left[\begin{matrix} -q_x(t) & -q_y(t) & -q_z(t) \\ q_w(t) & q_z(t) & -q_y(t) \\ -q_z(t) & q_w(t) & q_x(t) \\ q_y(t) & -q_x(t) & q_w(t) \end{matrix} \right] \left(\begin{matrix} \omega_x(t) \\ \omega_y(t) \\ \omega_z(t) \end{matrix} \right).$$

Note that the left-hand sides of your equations miss the derivatives and the derivative of the real part $q_w(t)$ has a sign error. To numerically integrate these equations in time you need to discretize them. You can probably use a simple explicit Euler scheme for your purpose, where e.g.

$$ \frac{\mathrm{d}}{\mathrm{d}t} q_w(t) \approx \frac{q_w(t + \delta t) - q_w(t)}{\delta t}.$$

Thus to "predict" the orientation you can use

$$ \left(\begin{matrix}q_w(t+\delta t) \\ q_x(t+\delta t) \\ q_y(t+\delta t) \\ q_z(t+\delta t) \end{matrix}\right) = \left(\begin{matrix}q_w(t) \\ q_x(t) \\ q_y(t) \\ q_z(t) \end{matrix}\right) + \frac{1}{2} \delta t \left[\begin{matrix} -q_x(t) & -q_y(t) & -q_z(t) \\ q_w(t) & q_z(t) & -q_y(t) \\ -q_z(t) & q_w(t) & q_x(t) \\ q_y(t) & -q_x(t) & q_w(t) \end{matrix} \right] \left(\begin{matrix} \omega_x(t) \\ \omega_y(t) \\ \omega_z(t) \end{matrix} \right).$$

So if your equations are read as update formulas, then no they are not quite correct (see note above), yes you should take the time interval into account and yes you should add the previous quaternion. Also you need to make sure that you normalize your quaternion after updating since it can loose its unit length when numerically integrating it.