Angle between two Euler Angles or Quaternions, Gimbal Lock

Question

I have inertial measurement unit sensors that can output data in either quaternions or Euler angles. As a biomechanist, Euler angles make more sense to me, but I sorta understand quaternions as well, but I've never really studied them. I do have a math background, so I'm not completely lost and I understand the Gimbal lock effect in Euler Angles.

I'm looking to calculate angles between two vectors whether they're quaternions or Euler angles, specifically on a human subject. I basically want to find the axis of rotations and calculate the difference in angles in the three basic components (x,y,z), and it seems incredibly unlikely that a person would be able to contort their body and reach Gimbal-lock.

I've read this paper and it seems like the way you choose to approach a rotation (x->y->z gets you to the same point as x->z->y but are different paths in terms of angles taken) is where Gimbal-lock comes into play, but the proposed XZ'Y'' sequence seems to avoid Gimbal-lock altogether.

I've read that quaternions are simply easier for computers to calculate which is where I would like to keep using quaternions since I am using a Pi, but I just don't quite fully understand how to go from quaternions to your basic x,y,z components. So I guess my questions are:

Are quaternions necessary for human movement?
Would maintaining numbers in quaternions until the very final angle calculations and converting to Euler angle in the last step avoid Gimbal lock?

BitTickler BitTickler · Accepted Answer · 2017-01-08T08:07:18

Basically, you have 2 main options to work e.g. with a skeleton.

Use 4x4 matrix (which allows rotation and translation)
Use (unit) quaternions for rotation and offsets for translation.

If you look at a typical implementation of a function taking 2 vectors and returning a quaternion, giving the rotation between them, you will see that it is not just a simple formula. Edge cases are being identified and taken care of.

let rotFromVectors (v1 : vec3) (v2 : vec3) : quat =
    let PI = System.Math.PI
    let PI_BY_TWO = PI / 2.0
    let TWO_PI = 2.0 * PI
    let ZERO_ROTATION = quat(0.0f,0.0f,0.0f,1.0f)
    let aabb = sqrt (float (vec3.dot(v1, v1)) * float (vec3.dot(v2,v2)))
    if aabb <> 0.0
    then
        let ab = float (vec3.dot(v1,v2)) / aabb
        let c = 
            vec3
                ( float32 ((float v1.y * float v2.z - float v1.z * float v2.y) / aabb)
                , float32 ((float v1.z * float v2.x - float v1.x * float v2.z) / aabb)
                , float32 ((float v1.x * float v2.y - float v1.y * float v2.x) / aabb)
                )
        let cc = float (vec3.dot(c, c))
        if cc <> 0.0
        then
            let s =
                match ab > -sin (PI_BY_TWO) with //0.707107f
                | true -> 1.0 + ab
                | false -> cc / (1.0 + sqrt (1.0-cc))
            let m = sqrt (cc + s * s)
            quat(float32 (float c.x / m), float32 (float c.y / m), float32 (float c.z / m), float32(s / m))
        else
            if ab > 0.0
            then 
                ZERO_ROTATION
            else
                let m = sqrt (v1.x * v1.x + v1.y * v1.y)
                if(m <> 0.0f)
                then
                    quat(v1.y / m, (-v1.x) / m, 0.0f, 0.0f)
                else
                    quat(1.0f,0.0f,0.0f,0.0f)
    else
        ZERO_ROTATION

Where quat is the type for a quaternion and vec3 the type of a 3D vector in the code above.

The code to rotate a vector by a quaternion is just about as straightforward as the math suggests:

let rotateVector (alpha : quat) (v:vec3) : vec3 =
    let s = vec3.length v
    quat.inverse alpha * (vecToPureQuat v) * alpha |> pureQuatToVec |> fun v' -> v' * s

And last not least the conversion functions between (some... Euler angles - there are actually 24 different versions of Euler angles, 12 with fixed angle rotations and 12 with consecutive rotations) uses the half angle approach.

let eulerToRot (v:vec3) : quat =
    let d = 0.5F
    let t0 = cos (v.z * d)
    let t1 = sin (v.z * d)
    let t2 = cos (v.y * d)
    let t3 = sin (v.y * d)
    let t4 = cos (v.x * d)
    let t5 = sin (v.x * d)
    quat
        (   t0 * t3 * t4 - t1 * t2 * t5
        ,   t0 * t2 * t5 + t1 * t3 * t4
        ,   t1 * t2 * t4 - t0 * t3 * t5
        ,   t0 * t2 * t4 + t1 * t3 * t5
        )
    |> quat.normalize

let rotToEuler (q:quat) : vec3 =
    let ysqr = q.y * q.y
    // roll (x-axis rotation)
    let t0 = +2.0f * (q.w * q.x + q.y * q.z)
    let t1 = +1.0f - 2.0f * (q.x * q.x + ysqr)
    let roll = atan2 t0 t1

    // pitch (y-axis rotation)
    let t2 = 
        let t2' = +2.0f * (q.w * q.y - q.z * q.x)
        match t2' with
        | _ when t2' > 1.0f -> 1.0f
        | _ when t2' < -1.0f -> -1.0f
        | _ -> t2'
    let pitch = asin t2

    // yaw (z-axis rotation)
    let t3 = +2.0f * (q.w * q.z + q.x *q.y)
    let t4 = +1.0f - 2.0f * (ysqr + q.z * q.z)
    let yaw = atan2 t3 t4
    vec3(roll,pitch,yaw)

The final trick to know is, that to turn a vector into a (pure) quaternion comes in handy for the rotateVector function.

let vecToPureQuat (v:vec3) : quat =
    quat(v.x,v.y,v.z,0.0f)

let pureQuatToVec (q:quat) : vec3 =
    vec3(q.x,q.y,q.z)

So, to answer your main question: Are quaternions necessary? No. You can as well use 4x4 matrices.

And you can go from one to the other if it deems you useful:

    let offsetAndRotToMat (offset:vec3) (q:quat) : mat4 =
        let ux = v3 1 0 0
        let uy = v3 0 1 0
        let uz = v3 0 0 1
        let rx = rotateVector q ux 
        let ry = rotateVector q uy 
        let rz = rotateVector q uz 
        mat4
            (
                rx.x, rx.y, rx.z, 0.0f,
                ry.x, ry.y, ry.z, 0.0f,
                rz.x, rz.y, rz.z, 0.0f,
                offset.x,offset.y,offset.z,1.0f
            )

Angle between two Euler Angles or Quaternions, Gimbal Lock

1 Answers