Fast quaternion math in shaders

Introduction

I was recently looking around on the internet for fast ways to compute $q_1 q_2$ and $q v q^{-1}$ in shader languages like GLSL or HLSL. To my surprise, I could find no such expressions for $q_1 q_2$, and I only found the final expression of $q v q^{-1}$, but no derivation. Hence this page.

Conventions and notation

We use the standard convention that a regular 3d vector is the quaternion $v = v_x \mathbf{i} + v_y \mathbf{j} + v_z \mathbf{k}$. A general quaternion $q$ may also have a scalar part: $q = q_x \mathbf{i} + q_y \mathbf{j} + q_z \mathbf{k} + q_w$.

Quaternion product

In the fully general case, we have two quaternions $p$ and $q$: $$ p q = (p_x \mathbf{i} + p_y \mathbf{j} + p_z \mathbf{k} + p_w)(q_x \mathbf{i} + q_y \mathbf{j} + q_z \mathbf{k} + q_w) $$ $$ \begin{matrix} p q = & p_x & \mathbf{i} &(q_x \mathbf{i} + q_y \mathbf{j} + q_z \mathbf{k} + q_w) \\ & p_y & \mathbf{j} &(q_x \mathbf{i} + q_y \mathbf{j} + q_z \mathbf{k} + q_w) \\ & p_z & \mathbf{k} &(q_x \mathbf{i} + q_y \mathbf{j} + q_z \mathbf{k} + q_w) \\ & p_w & &(q_x \mathbf{i} + q_y \mathbf{j} + q_z \mathbf{k} + q_w) \\ \end{matrix} $$ $$ \begin{matrix} p q & = & p_x ( & - & q_x & & + & q_y & \mathbf{k} & - & q_z & \mathbf{j} & + & q_w & \mathbf{i} &) \\ & & p_y ( & - & q_x & \mathbf{k} & - & q_y & & + & q_z & \mathbf{i} & + & q_w & \mathbf{j} &) \\ & & p_z ( & + & q_x & \mathbf{j} & - & q_y & \mathbf{i} & - & q_z & & + & q_w & \mathbf{k} &) \\ & & p_w ( & + & q_x & \mathbf{i} & + & q_y & \mathbf{j} & + & q_z & \mathbf{k} & + & q_w & &) \\ \end{matrix} $$ $$ \begin{matrix} p q & = & - & p_xq_x & & + & p_xq_y & \mathbf{k} & - & p_xq_z & \mathbf{j} & + & p_xq_w & \mathbf{i} & \\ & & - & p_yq_x & \mathbf{k} & - & p_yq_y & & + & p_yq_z & \mathbf{i} & + & p_yq_w & \mathbf{j} & \\ & & + & p_zq_x & \mathbf{j} & - & p_zq_y & \mathbf{i} & - & p_zq_z & & + & p_zq_w & \mathbf{k} & \\ & & + & p_wq_x & \mathbf{i} & + & p_wq_y & \mathbf{j} & + & p_wq_z & \mathbf{k} & + & p_wq_w & & \\ \end{matrix} $$ $$ \begin{matrix} p q & = & ( & + p_x q_w + p_y q_z - p_z q_y + p_w q_x & ) & \mathbf{i} \\ & + & ( & - p_x q_z + p_y q_w + p_z q_x + p_w q_y & ) & \mathbf{j} \\ & + & ( & + p_x q_y - p_y q_x + p_z q_w + p_w q_z & ) & \mathbf{k} \\ & + & ( & - p_x q_x - p_y q_y - p_z q_z + p_w q_w & ) & \\ \end{matrix} $$ $$ \begin{matrix} p q & = & p_v \times q_v & + & ( & + p_x q_w + p_w q_x & ) & \mathbf{i} \\ & & & + & ( & + p_y q_w + p_w q_y & ) & \mathbf{j} \\ & & & + & ( & + p_z q_w + p_w q_z & ) & \mathbf{k} \\ & & & + & ( & - p_v \cdot q_v + p_w q_w & ) & \\ \end{matrix} $$

We finally get the formula $$ \boxed{ p q = p_v \times q_v + q_w p_v + p_w q_v + p_w q_w - p_v \cdot q_v } $$

In HLSL notation, this becomes

                
float4
quaternion_product(float4 p, float4 q)
{
  return
    float4(
      cross(p.xyz, q.xyz) + q.w*p.xyz + p.w*q.xyz,
      p.w*q.w - dot(p.xyz, q.xyz)
    );
}
              

Quaternion-vector product

We can specialize the above to the case where one of the quternions is a vector: $$ \boxed{ q v = q_v \times v + q_w v - q_v \cdot v } $$

                
float4
quaternion_vector_product(float4 q, float3 v)
{
  return
    float4(
      cross(q.xyz, v) + q.w*v,
      -dot(q.xyz, v)
    );
}
              

Quaternion conjugation of a vector

In this section, we will assume that $\left| q \right| = 1$ so that $q^{-1} = \frac{q^*}{\left| q \right|} = q^* = -q_x \mathbf{i} - q_y \mathbf{j} - q_z \mathbf{k} + q_w = q_w - q_v$.

Of special interest is the conjugation of $p$ by $q$, namely the expression $q p q^{-1}$. We are interested in this expression because when $q$ is a unit quaternion and $p=v$ is a vector, then $q v q^{-1}$ is the vector which is $v$ rotated by the rotation represented by $q$.

$$ qvq^{-1} = (q_w + q_v)v(q_w - q_v) = q_w^2 v + q_w(q_v v - v q_v) - q_v p q_v $$ Let's simplify the expression $q_v v - v q_v$. Using the formula we derived for the quaternion product, we can see that most of the terms cancel, leaving us with just $2q_v \times v$.

Now let's focus on the very last term, namely $q_v v q_v$. We use our product formula again: $$ q_v v q_v = (q_v \times v - q_v \cdot v) q_v = (q_v \times v) \times q_v - (q_v \times v) \cdot q_v - (q_v \cdot v) q_v = -q_v \times (q_v \times v) - (q_v \cdot v) q_v $$

Substituting these expressions back, we get $$ qvq^{-1} = q_w^2 v + 2 q_w (q_v \times v) + q_v \times (q_v \times v) + (q_v \cdot v) q_v $$

We can now spot the vector triple product and simplify $q_v(q_v \cdot v) = q_v \times (q_v \times v) + v q_v^2$. Using this as well as $ 1 = \left| q \right|^2 = q_v^2 + q_w^2$, we obtain: $$ \boxed{ qvq^{-1} = v + 2 q_w (q_v \times v) + 2 q_v \times (q_v \times v) } $$

In HLSL notation, this can be written as

                
float4
quaternion_vector_conjugate(float4 q, float3 v)
{
  float3 t = 2 * cross(q.xyz, v);
  return v + q.w*t + cross(q.xyz, t);
}
              

This function can be found here and there on the internet, but I think it's nice to also see a derivation!