Applying angular velocity to quaternion

Question

I am reading Game Physics Engine Development by Ian Millington currently. In his book, he represent an orientation of rigid body by quaternion. I don't understand the formula where he apply angular velocity to the quaternion. The formula is

$$ q_{new} = q_{0} + \frac{t}{2} * w * q_{0} $$

where w is the angular velocity in quaternion representation.

$$ w = \begin{bmatrix} 0\\ w_{x}\\ w_{y}\\ w_{z} \end{bmatrix} $$

He doesn't explain where the formula come from. Currently I apply this method in my physics engine where I compute the angular rotation then convert it to quaternion and multiply it. It is slower than the above operation because of the sine, cosine and square root operation when converting the rotation to quaternion. Can someone explain to me where did the formula above come from? I understand basic quaternion operation like representing axis angle by quaternion and rotating a vector by quaternion.

Answer 1

In the following I use a vector hat to represent a quaternion and a unit vector hat to represent a unit quaternion.

Let q(t) represent the quaternion rotation at any point in time. Let r be a quaternion representing the rotation traveled through in the time period delta t. Then we have:

$$\hat q(0)=\hat q_0$$ and $$\hat q(\Delta t)=\hat r\hat q_0$$

If we realize that the change from time zero to time delta t is linear we can get:

$$\hat q(t)=\hat r^{t/\Delta t}\hat q_0$$

Using Euler's Formula we obtain:

$$\hat r=e^{\frac{\theta}{2}\hat u}$$

For some angle theta and axis of rotation u (Note that u has no real part and is a "pure" quaternion). This leads us to:

$$\hat r^{t/\Delta t}=exp(\frac{\theta}{2}\hat u\frac t {\Delta t})$$

Plugging this in we obtain:

$$\hat q(t)=exp(\frac{\theta}{2}\hat u\frac t {\Delta t})\hat q_0$$

Taking the derivative of both sides we obtain:

$$\hat q'=\frac {\theta \hat u} {2 \Delta t} exp(\frac{\theta}{2}\hat u\frac t {\Delta t})\hat q_0$$ $$\hat q'=\frac {\theta \hat u} {2 \Delta t}\hat q$$

If we note that: $$\frac \theta {\Delta t} = w$$

For some scalar w. This represents the magnitude of the angular velocity. The axis of rotation is the vector given by u (u has no real part).

Or we could put it more bluntly if we say that:

$$\vec w=w\hat u=w_x i+w_y j+ w_z k$$

This gives us:

$$\hat q'=\frac 1 2 \vec w \hat q$$

Now the formula that you are citing is actually not exact but only approximate......You can say for "small" delta t's that:

$$\frac {\Delta \hat {q}} {\Delta t} \approx \frac 1 2 \vec w \hat q$$

And since delta t is small we can choose any value within that time period for q, so we choose the one we already have:

$$\frac {\Delta \hat {q}} {\Delta t} \approx \frac 1 2 \vec w \hat q_{old}$$ $${\Delta \hat {q}} \approx \frac 1 2 \vec w \hat q_{old} {\Delta t}$$ $$\hat q_{new} \approx \frac 1 2 \vec w \hat q_{old} {\Delta t} + \hat q_{old}$$ $$\hat q_{new} \approx (\frac 1 2 \vec w {\Delta t} + 1)\hat q_{old}$$

The formula you link to comes from solving the differential equation and is an exact formula

Answer 2

I finally found a very clear derivation, here: https://fgiesen.wordpress.com/2012/08/24/quaternion-differentiation/

It's by Fabian Giesen (ryg) of demo-scene fame.

Answer 3

The formula from the book looks wrong.

Try putting t = 1, W = [ 0, 0, 2*Pi ] into the formula. That input represents 1 complete revolution around Z axis, and the output should be identity quaternion, yet the formula returns something very far from identity.

Here’s my attempt. It’s in C# but I hope it’s simple enough to adopt to any other language.

/// <summary>Apply angular velocity to the quaternion</summary>
public Quaternion rotate( Vector3 angularVelocity, float deltaTime )
{
    Vector3 vec = angularVelocity * deltaTime;
    float length = vec.Length();
    if( length < 1E-6F )
        return this;    // Otherwise we'll have division by zero when trying to normalize it later on

    // Convert the rotation vector to quaternion. The following 4 lines are very similar to CreateFromAxisAngle method.
    float half = length * 0.5f;
    float sin = MathF.Sin( half );
    float cos = MathF.Cos( half );
    // Instead of normalizing the axis, we multiply W component by the length of it. This method normalizes result in the end.
    Quaternion q = new Quaternion( vec.X * sin, vec.Y * sin, vec.Z * sin, length * cos );

    q = Multiply( q, this );
    q.Normalize();
    // The following line is not required, only useful for people. Computers are fine with 2 different quaternion representations of each possible rotation.
    if( q.W < 0 ) q = Negate( q );
    return q;
}

Answer 4

Sort of in response to Soonts answer, but also for other people that may wonder, like I did, I found that the formula does work after an implementation in Unity, C#. I am sorry if "reviving" this is considered bad input.

Inputting W = [0, 0, Pi/2, 0] ([x, y, z, w] format) and letting t be the delta time between screen refreshes (Time.deltaTime), does yield a rotation of 90° in 1 second around the Z axis in the correct, trigonometric, direction.

An input of W = [0, 0, Pi, 0] does yield the same rotation, but with an angle of 180°, for instance.

Finally, inputting W = [Pi, Pi, Pi, 0] yields rotations in the correct directions on all axes, and made the test game object return to its original starting orientation after 1 second, as it should.

My attempt is very naive, and doesn't account for any specific scenarios where the angular velocity is too small or the such (but as far as I tested, down to 1e-06 at least, the rotation still executes normally. I did not test lower values, as it would take forever):

void UpdateTransform()
{   // transform.rotation is the game object's orientation quaternion in the game's worldspace. 'angular_velocity' is defined as a Quaternion of (x, y, z, w) parameters.

    transform.rotation = SummedQuaternion(ScaledQuaternion(angular_velocity, Time.deltaTime / 2) * transform.rotation, transform.rotation);
}

Quaternion ScaledQuaternion(Quaternion quaternion, float scale)
{   
    Quaternion scaled_quaternion = new Quaternion(quaternion.x * scale, quaternion.y * scale, quaternion.z * scale, quaternion.w * scale);

    return scaled_quaternion;
}

Quaternion SummedQuaternion(Quaternion quaternion_alpha, Quaternion quaternion_beta)
{
    Quaternion summed_quaternion = new Quaternion(quaternion_alpha.x + quaternion_beta.x, quaternion_alpha.y + quaternion_beta.y, quaternion_alpha.z + quaternion_beta.z, quaternion_alpha.w + quaternion_beta.w);

    return summed_quaternion;
}