Quaternion instant numeric drift

Question

I'm propagating a quaternion with little changes created from user mouse input where I create a quaternion for the added rotation in the X-axis and one for the Y-axis and combine these rotations into one quaternion. I keep adding (appending) these changes to a camera quaternion and so created a glimbal-lock free camera (which is still rolling like crazy :S)

The moment I start rotating my camera the rotation matrix starts to scale and skew out of control. Everything works fine when I normalize the camera quaternion every frame, but normalizing every frame should not be needed right? If this IS numeric drift, it couldn't be visible so soon could it?

In this question: Quaternion and normalization

David Hammen posted this comment:

Normalization is a kludge to drag the quaternion back to the unit 3-sphere. The reason quaternions drift off the manifold is because that Euler step is in a sense invalid mathematically: The unit quaternions are a group, not an algebra. Not using unit quaternions is also invalid; now you don't have a chart on so(3). Do it right (e.g., Lie group integration techniques) and the quaternion will pretty much stay on the manifold. You might still need to normalize occasionally

Could it be that because I'm creating two quaternions from euler rotations that this is the cause of my drastic rotation errors?

Should quaternion generated from euler rotations always be normalized before use?

Answer 1

First, find if your quaternion is become NONunit? and problem source in that.

I recommend use robust conversion to matrix of nonunit quaternions, in pseudocode

inline void to_matrix( matrix33& m  )const  {
        float wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
        float s  = 2.0f/norm_2_sqr(); 
        x2 = x * s;    y2 = y * s;    z2 = z * s;
        xx = x * x2;   xy = x * y2;   xz = x * z2;
        yy = y * y2;   yz = y * z2;   zz = z * z2;
        wx = w * x2;   wy = w * y2;   wz = w * z2;

        m.m[0][0] = 1.0f - (yy + zz);
        m.m[0][1] = xy - wz;
        m.m[0][2] = xz + wy;

        m.m[1][0] = xy + wz;
        m.m[1][1] = 1.0f - (xx + zz);
        m.m[1][2] = yz - wx;

        m.m[2][0] = xz - wy;
        m.m[2][1] = yz + wx;
        m.m[2][2] = 1.0f - (xx + yy);

        //1 division    16 multiplications    15 addidtions    12 variables
    }