Quaternion based point rotations using GLM

Question

I am trying to rotate a point using quaternions as implemented in GLM. The end goal is to use this code to create an orbital camera but this is a side note to help understand the motivation behind the code.

To understand quaternion based rotations better I have written a bit of code that contains two loops. The first loop will incrementally change the orientation of the quaternion by rotating it in steps around the X axis all the way to 90 degrees, and the second loop will continue to apply a rotation all the way to 90 degrees in steps around the Z axis. The loops do 4 steps each. So each loop incrementally rotates for 90 / 4 = 22.5 degrees around their respective axis. The change of orientation is applied using quaternion multiplication and tracked using Euler angles. The loops should end with a quaternion that will rotate a point at (0, 0, 3) to (3, 0, 0). Please note, I am not just trying to determine the quaternion that will do this rotation. The goal is to perform a series of incremental rotations.

If we look at the below picture the transition from C to I happens in the first loop and then the transition from I to R in the second (pardon the sparse point naming).

Rotation of a point is defined as (see here and here):

v' = q * v * q^-1

where v should be considered a pure quaternion (with a zero scalar term w) and q needs to be a unit quaternion (of length 1). And from what I understand the right hand multiplication with the inverse of the quaternion is needed to keep the resulting v' in 3D space and not end up with a 4D vector. So v' needs to be a pure quaternion as well.

Then there is the doubling effect of the rotation where the left hand multiplication with q contributes half of the desired rotation and the right hand side multiplication with the inverse adds another half of the desired rotation.

There is an excellent interactive visualisation and explanation of quaternions by Ben Eater and Grant Sanderson, that I used as a cross-reference. It can be found here.

So we first need to use a quaternion that rotates by 11.25 degrees around the X axis and GLM returns this quaternion for Euler angles (quaternion notation [w, [x, y, z]] is used):

Rotation of [ 11.25, 0.00,  0.00] deg => Q: [ 0.9952, [ 0.0980,  0.0000,  0.0000]]

According to this, and since we are rotating purely around the X axis, we could verify the rotation amount in the GLM calculated quaternion by performing an acos on the w component of the quaternion:

float angle = acosf(q.w)

then:

acos(0.9952) = 0.0980 rad / 5.6 degrees

Which is half the desired angle... And this is also confirmed in a cross check with the interactive animation (pardon the rounding):

So the quaternion returned by GLM for 11.25 degrees actually rotates for half the desired angle... If we look at the GLM code the calculation of the w component from Euler angles is a little more complex because rotation can happen around an arbitrary axis of rotation... But there is a distinct halving of the Euler angles:

template <typename T, precision P>
GLM_FUNC_QUALIFIER tquat<T, P>::tquat(tvec3<T, P> const & eulerAngle)
{
    tvec3<T, P> c = glm::cos(eulerAngle * T(0.5));
    tvec3<T, P> s = glm::sin(eulerAngle * T(0.5));
    
    this->w = c.x * c.y * c.z + s.x * s.y * s.z;
    this->x = s.x * c.y * c.z - c.x * s.y * s.z;
    this->y = c.x * s.y * c.z + s.x * c.y * s.z;
    this->z = c.x * c.y * s.z - s.x * s.y * c.z;
}

My first question is why is GLM halving the angle?

Despite the difference in the desired rotation angle I went ahead to check the rotation results with the two loops. And the results were... Unexpected.

If I used the "incorrect form" of rotation (suggested by some OpenGL online tutorials) and rotated the point by a left hand multiplication only (but for a full step of 22.5 degrees):

v' = q * v

I got the result I was hoping for. The point was following all the intermediate steps correctly and went from (0, 0, 3) to (3, 0, 0). Also the w component was 0 at all intermediate steps.

But if I used the "correct form" of rotation and rotated the point by a left hand multiplication with q and right hand multiplication with the inverse of q (for a half step of 11.25 degrees to account for the doubling of rotation):

v' = q * v * q^-1

I start getting wrong results as soon as the second loop starts to rotate the point around the Z axis. A small but distinct Z component starts to creep in and the rotation is just short of the full step of 22.5 degrees. This is visible in the green dots in the image below.

The w component of the rotated point remains 0 for both methods of rotation...

Can anyone explain why the GLM rotation works correctly with a single multiplication from the left?

Is this some kind of optimisation to reduce the number of operations to a minimum?

Can I use the v' = q * v rotation in GLM to get consistent and correct results for all rotations?

Code:

const int rotSteps = 4;
// Rotate around X axis in steps to 90deg
vec3 eulerState = vec3(0.0f);
// point we want to rotate (use vec4 to track the w component during rotations)
vec4 v = vec4(0.0f, 0.0f, 3.0f, 0.0f);

// Full Euler steps for q * v rotation
quat orientF   = quat(1.0f, 0.0f, 0.0f, 0.0f);
vec3 euler     = vec3(RAD(90.0f), RAD(0.0f), RAD(0.0f));
vec3 eulerStep = euler / (float)rotSteps;
quat qEulerF   = quat(eulerStep); // GetRotQuat(eulerStep);

vec4 qa          = ToAngularForm(qEulerF);
vec3 orientEuler = eulerAngles(qEulerF);
CLogD(TAG, "Rot Full Step    Q [W, X, Y, Z]: " FMT_Q(4)  " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerF), PAR_V3(degrees(orientEuler)), PAR_QA(qa));

// Half Euler steps for q * v * q^-1 rotation
quat orientH    = quat(1.0f, 0.0f, 0.0f, 0.0f);
vec3 eulerStepH = eulerStep / 2.0f;
quat qEulerH    = quat(eulerStepH); // GetRotQuat(eulerStepH);

qa          = ToAngularForm(qEulerH);
orientEuler = eulerAngles(qEulerH);
CLogD(TAG, "Rot Half Step    Q [W, X, Y, Z]: " FMT_Q(4) " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerH), PAR_V3(degrees(orientEuler)), PAR_QA(qa));

quat qEulerHI = inverse(qEulerH);
vec4 qai      = ToAngularForm(qEulerHI);
orientEuler   = eulerAngles(qEulerHI);
CLogD(TAG, "Rot Half Step Q^-1 [W, X, Y, Z]: " FMT_Q(4) " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerHI), PAR_V3(degrees(orientEuler)), PAR_QA(qai));


for (int rotStep = 1; rotStep <= rotSteps; ++rotStep)
{
    // Track the absolute Euler rotation
    eulerState += eulerStep;
    // Rotate by incremental rotation as defined by Euler angles
    orientH = qEulerH * orientH;
    orientEuler = eulerAngles(orientH);
    CLogI(TAG, "Rot Step %d. Curr Abs Q: " FMT_Q(4) "/" FMT_V3(2) "deg, Abs Euler: " FMT_V3(2) "deg",
          rotStep, PAR_Q(orientH), PAR_V3(degrees(orientEuler)), PAR_V3(degrees(eulerState)));

    // Transform the point using the correct q * v * q^-1 rotation and multiply from Left and Right
    quat orientHI = inverse(orientH);
    qa  = ToAngularForm(orientH);
    qai = ToAngularForm(orientHI);

    vec4 rotV = orientH * v * orientHI;
    CLogD(TAG, "Rot      QL: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientH), PAR_QA(qa));
    CLogD(TAG, "Rot      QR: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientHI), PAR_QA(qai));
    CLogD(TAG, "Rot LR   -> " FMT_V4(1), PAR_V4(rotV));

    // Transform the point using the incorrect q * v rotation and multiply from Left only
    orientF = qEulerF * orientF;
    qa      = ToAngularForm(orientF);

    rotV = orientF * v;
    CLogD(TAG, "Rot      QR: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientF), PAR_QA(qa));
    CLogD(TAG, "Rot L    -> " FMT_V4(1), PAR_V4(rotV));
}


// Rotate for 90 degrees around the Z axis
// Full Euler steps for q * v rotation
euler = vec3(RAD(0.0f), RAD(0.0f), RAD(90.0f));
eulerStep = euler / (float)rotSteps;
qEulerF = quat(eulerStep); // GetRotQuat(eulerStep);

qa = ToAngularForm(qEulerF);
orientEuler = eulerAngles(qEulerF);
CLogD(TAG, "Rot Full Step    Q [W, X, Y, Z]: " FMT_Q(4)  " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerF), PAR_V3(degrees(orientEuler)), PAR_QA(qa));

// Half Euler steps for q * v * q^-1 rotation
eulerStepH = eulerStep / 2.0f;
qEulerH = quat(eulerStepH); // GetRotQuat(eulerStepH);

qa = ToAngularForm(qEulerH);
orientEuler = eulerAngles(qEulerH);
CLogD(TAG, "Rot Half Step    Q [W, X, Y, Z]: " FMT_Q(4) " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerH), PAR_V3(degrees(orientEuler)), PAR_QA(qa));

qEulerHI = inverse(qEulerH);
qai = ToAngularForm(qEulerHI);
orientEuler = eulerAngles(qEulerHI);
CLogD(TAG, "Rot Half Step Q^-1 [W, X, Y, Z]: " FMT_Q(4) " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerHI), PAR_V3(degrees(orientEuler)), PAR_QA(qai));


for (int rotStep = 1; rotStep <= rotSteps; ++rotStep)
{
    // Track the absolute Euler rotation
    eulerState += eulerStep;
    // Rotate by incremental rotation as defined by Euler angles
    orientH = qEulerH * orientH;
    orientEuler = eulerAngles(orientH);
    CLogI(TAG, "Rot Step %d. Curr Abs Q: " FMT_Q(4) "/" FMT_V3(2) "deg, Abs Euler: " FMT_V3(2) "deg",
        rotStep, PAR_Q(orientH), PAR_V3(degrees(orientEuler)), PAR_V3(degrees(eulerState)));

    // Transform the point using the correct q * v * q^-1 rotation and multiply from Left and Right
    quat orientHI = inverse(orientH);
    qa = ToAngularForm(orientH);
    qai = ToAngularForm(orientHI);

    vec4 rotV = orientH * v * orientHI;
    CLogD(TAG, "Rot      QL: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientH), PAR_QA(qa));
    CLogD(TAG, "Rot      QR: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientHI), PAR_QA(qai));
    CLogD(TAG, "Rot LR   -> " FMT_V4(1), PAR_V4(rotV));

    // Transform the point using the incorrect q * v rotation and multiply from Left only
    orientF = qEulerF * orientF;
    qa = ToAngularForm(orientF);

    rotV = orientF * v;
    CLogD(TAG, "Rot      QR: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientF), PAR_QA(qa));
    CLogD(TAG, "Rot L    -> " FMT_V4(1), PAR_V4(rotV));
}

Output:

Rot Full Step    Q [W, X, Y, Z]: [ 0.9808, [ 0.1951,  0.0000,  0.0000]] / [ 22.50, -0.00,  0.00]deg / cos( 11.25) + sin( 11.25)( 1.00i +  0.00j +  0.00k)
Rot Half Step    Q [W, X, Y, Z]: [ 0.9952, [ 0.0980,  0.0000,  0.0000]] / [ 11.25, -0.00,  0.00]deg / cos( 5.63) + sin( 5.63)( 1.00i +  0.00j +  0.00k)
Rot Half Step Q^-1 [W, X, Y, Z]: [ 0.9952, [-0.0980, -0.0000, -0.0000]] / [-11.25, -0.00,  0.00]deg / cos( 5.63) + sin( 5.63)(-1.00i + -0.00j + -0.00k)
Rot Step 1. Curr Abs Q: [ 0.9952, [ 0.0980,  0.0000,  0.0000]]/[ 11.25, -0.00,  0.00]deg, Abs Euler: [ 22.50,  0.00,  0.00]deg
Rot      QL: [ 0.9952, [ 0.0980,  0.0000,  0.0000]] / cos( 5.6) + sin( 5.6)( 1.0i +  0.0j +  0.0k)
Rot      QR: [ 0.9952, [-0.0980, -0.0000, -0.0000]] / cos( 5.6) + sin( 5.6)(-1.0i + -0.0j + -0.0k)
Rot LR   -> [ 0.0, -1.1,  2.8,  0.0]
Rot      QR: [ 0.9808, [ 0.1951,  0.0000,  0.0000]] / cos( 11.3) + sin( 11.3)( 1.0i +  0.0j +  0.0k)
Rot L    -> [ 0.0, -1.1,  2.8,  0.0]
Rot Step 2. Curr Abs Q: [ 0.9808, [ 0.1951,  0.0000,  0.0000]]/[ 22.50, -0.00,  0.00]deg, Abs Euler: [ 45.00,  0.00,  0.00]deg
Rot      QL: [ 0.9808, [ 0.1951,  0.0000,  0.0000]] / cos( 11.3) + sin( 11.3)( 1.0i +  0.0j +  0.0k)
Rot      QR: [ 0.9808, [-0.1951, -0.0000, -0.0000]] / cos( 11.2) + sin( 11.2)(-1.0i + -0.0j + -0.0k)
Rot LR   -> [ 0.0, -2.1,  2.1,  0.0]
Rot      QR: [ 0.9239, [ 0.3827,  0.0000,  0.0000]] / cos( 22.5) + sin( 22.5)( 1.0i +  0.0j +  0.0k)
Rot L    -> [ 0.0, -2.1,  2.1,  0.0]
Rot Step 3. Curr Abs Q: [ 0.9569, [ 0.2903,  0.0000,  0.0000]]/[ 33.75, -0.00,  0.00]deg, Abs Euler: [ 67.50,  0.00,  0.00]deg
Rot      QL: [ 0.9569, [ 0.2903,  0.0000,  0.0000]] / cos( 16.9) + sin( 16.9)( 1.0i +  0.0j +  0.0k)
Rot      QR: [ 0.9569, [-0.2903, -0.0000, -0.0000]] / cos( 16.9) + sin( 16.9)(-1.0i + -0.0j + -0.0k)
Rot LR   -> [ 0.0, -2.8,  1.1,  0.0]
Rot      QR: [ 0.8315, [ 0.5556,  0.0000,  0.0000]] / cos( 33.8) + sin( 33.8)( 1.0i +  0.0j +  0.0k)
Rot L    -> [ 0.0, -2.8,  1.1,  0.0]
Rot Step 4. Curr Abs Q: [ 0.9239, [ 0.3827,  0.0000,  0.0000]]/[ 45.00, -0.00,  0.00]deg, Abs Euler: [ 90.00,  0.00,  0.00]deg
Rot      QL: [ 0.9239, [ 0.3827,  0.0000,  0.0000]] / cos( 22.5) + sin( 22.5)( 1.0i +  0.0j +  0.0k)
Rot      QR: [ 0.9239, [-0.3827, -0.0000, -0.0000]] / cos( 22.5) + sin( 22.5)(-1.0i + -0.0j + -0.0k)
Rot LR   -> [ 0.0, -3.0,  0.0,  0.0]
Rot      QR: [ 0.7071, [ 0.7071,  0.0000,  0.0000]] / cos( 45.0) + sin( 45.0)( 1.0i +  0.0j +  0.0k)
Rot L    -> [ 0.0, -3.0,  0.0,  0.0]

Rot Full Step    Q [W, X, Y, Z]: [ 0.9808, [ 0.0000,  0.0000,  0.1951]] / [ 0.00, -0.00,  22.50]deg / cos( 11.25) + sin( 11.25)( 0.00i +  0.00j +  1.00k)
Rot Half Step    Q [W, X, Y, Z]: [ 0.9952, [ 0.0000,  0.0000,  0.0980]] / [ 0.00, -0.00,  11.25]deg / cos( 5.63) + sin( 5.63)( 0.00i +  0.00j +  1.00k)
Rot Half Step Q^-1 [W, X, Y, Z]: [ 0.9952, [-0.0000, -0.0000, -0.0980]] / [ 0.00, -0.00, -11.25]deg / cos( 5.63) + sin( 5.63)(-0.00i + -0.00j + -1.00k)
Rot Step 1. Curr Abs Q: [ 0.9194, [ 0.3808,  0.0375,  0.0906]]/[ 45.00,  0.00,  11.25]deg, Abs Euler: [ 90.00,  0.00,  22.50]deg
Rot      QL: [ 0.9194, [ 0.3808,  0.0375,  0.0906]] / cos( 23.2) + sin( 23.2)( 1.0i +  0.1j +  0.2k)
Rot      QR: [ 0.9194, [-0.3808, -0.0375, -0.0906]] / cos( 23.2) + sin( 23.2)(-1.0i + -0.1j + -0.2k)
Rot LR   -> [ 1.0, -2.8,  0.0,  0.0]
Rot      QR: [ 0.6935, [ 0.6935,  0.1379,  0.1379]] / cos( 46.1) + sin( 46.1)( 1.0i +  0.2j +  0.2k)
Rot L    -> [ 1.1, -2.8,  0.0,  0.0]
Rot Step 2. Curr Abs Q: [ 0.9061, [ 0.3753,  0.0747,  0.1802]]/[ 45.00, -0.00,  22.50]deg, Abs Euler: [ 90.00,  0.00,  45.00]deg
Rot      QL: [ 0.9061, [ 0.3753,  0.0747,  0.1802]] / cos( 25.0) + sin( 25.0)( 0.9i +  0.2j +  0.4k)
Rot      QR: [ 0.9061, [-0.3753, -0.0747, -0.1802]] / cos( 25.0) + sin( 25.0)(-0.9i + -0.2j + -0.4k)
Rot LR   -> [ 1.9, -2.4,  0.1,  0.0]
Rot      QR: [ 0.6533, [ 0.6533,  0.2706,  0.2706]] / cos( 49.2) + sin( 49.2)( 0.9i +  0.4j +  0.4k)
Rot L    -> [ 2.1, -2.1,  0.0,  0.0]
Rot Step 3. Curr Abs Q: [ 0.8841, [ 0.3662,  0.1111,  0.2682]]/[ 45.00,  0.00,  33.75]deg, Abs Euler: [ 90.00,  0.00,  67.50]deg
Rot      QL: [ 0.8841, [ 0.3662,  0.1111,  0.2682]] / cos( 27.9) + sin( 27.9)( 0.8i +  0.2j +  0.6k)
Rot      QR: [ 0.8841, [-0.3662, -0.1111, -0.2682]] / cos( 27.9) + sin( 27.9)(-0.8i + -0.2j + -0.6k)
Rot LR   -> [ 2.5, -1.6,  0.3,  0.0]
Rot      QR: [ 0.5879, [ 0.5879,  0.3928,  0.3928]] / cos( 54.0) + sin( 54.0)( 0.7i +  0.5j +  0.5k)
Rot L    -> [ 2.8, -1.1,  0.0,  0.0]
Rot Step 4. Curr Abs Q: [ 0.8536, [ 0.3536,  0.1464,  0.3536]]/[ 45.00,  0.00,  45.00]deg, Abs Euler: [ 90.00,  0.00,  90.00]deg
Rot      QL: [ 0.8536, [ 0.3536,  0.1464,  0.3536]] / cos( 31.4) + sin( 31.4)( 0.7i +  0.3j +  0.7k)
Rot      QR: [ 0.8536, [-0.3536, -0.1464, -0.3536]] / cos( 31.4) + sin( 31.4)(-0.7i + -0.3j + -0.7k)
Rot LR   -> [ 2.9, -0.7,  0.4,  0.0]
Rot      QR: [ 0.5000, [ 0.5000,  0.5000,  0.5000]] / cos( 60.0) + sin( 60.0)( 0.6i +  0.6j +  0.6k)
Rot L    -> [ 3.0,  0.0,  0.0,  0.0]

omahena omahena · Accepted Answer · 2020-07-22T08:21:47

I have the answers to my question and a working orbital camera, but haven't had the time to double check if the sample code now works correctly - it should.

First question was why GLM is halving the angle during the quaternion conversion and it looks like according to the extended Euler's formula... It has to. This part could stand a little more investigation but due to lack of time I'll have to accept it.

The vector rotation in GLM was implemented using the multiplication operator. This means that when multiplying a vec3 with a quaternion will not do a conversion of the the vector into a quaternion and then performing a multiplication, it will do a vector rotation instead:

template <typename T, precision P>
GLM_FUNC_QUALIFIER tvec3<T, P> operator*(tquat<T, P> const & q, tvec3<T, P> const & v)
{
    tvec3<T, P> const QuatVector(q.x, q.y, q.z);
    tvec3<T, P> const uv(glm::cross(QuatVector, v));
    tvec3<T, P> const uuv(glm::cross(QuatVector, uv));

    return v + ((uv * q.w) + uuv) * static_cast<T>(2);
}

So, yes, the correct way to rotate a vector using a quaternion is to use the multiplication operator between a quaternion and a vector like this:

v' = q * v

or in C++:

vec3 posOrigin;
quat rotQ;
...
vec3 posRot = rotQ * posOrigin;

This code does not actually do a direct quaternion multiplication. It does a rotation. Personally I'd prefer GLM offered a function call like rotate(quat, vec)... But I am sure there is a reason for operator overloading obfuscation.

Please also note that the operand order matters, since the multiplication between a vector and a quaternion is defined like this:

template <typename T, precision P>
GLM_FUNC_QUALIFIER tvec3<T, P> operator*(tvec3<T, P> const & v, tquat<T, P> const & q)
{
    return glm::inverse(q) * v;
}

and will therefore rotate the vector in the inverse sense.

Note that GLM also implements multiplication between quaternions, but for this the multiplication operator between two quaternions needs to be used:

template <typename T, precision P>
template <typename U>
GLM_FUNC_QUALIFIER tquat<T, P> & tquat<T, P>::operator*=(tquat<U, P> const & r)
{
    tquat<T, P> const p(*this);
    tquat<T, P> const q(r);

    this->w = p.w * q.w - p.x * q.x - p.y * q.y - p.z * q.z;
    this->x = p.w * q.x + p.x * q.w + p.y * q.z - p.z * q.y;
    this->y = p.w * q.y + p.y * q.w + p.z * q.x - p.x * q.z;
    this->z = p.w * q.z + p.z * q.w + p.x * q.y - p.y * q.x;
    return *this;
}

Since GLM has precious little documentation I could find, such operator overloading leads to bad assumptions and a significant loss of time. So I suppose I should have read the GLM code instead of assuming what it does...

Quaternion based point rotations using GLM

1 Answers