Rotate model around x,y,z axes, without gimbal lock, with input data always as x,y,z axes angle rotations

Question

I have an input device that gives me 3 angles -- rotation around x,y,z axes.

Now I need to use these angles to rotate the 3D space, without gimbal lock. I thought I could convert to Quaternions, but apparently since I'm getting the data as 3 angles this won't help?

If that's the case, just how can I correctly rotate the space, keeping in mind that my input data simply is x,y,z axes rotation angles, so I can't just "avoid" that. Similarly, moving around the order of axes rotations won't help -- all axes will be used anyway, so shuffling the order around won't accomplish anything. But surely there must be a way to do this?

If it helps, the problem can pretty much be reduced to implementing this function:

void generateVectorsFromAngles(double &lastXRotation,
                                                      double &lastYRotation,
                                                      double &lastZRotation,
                                                      JD::Vector &up,
                                                      JD::Vector &viewing) {
    JD::Vector yaxis = JD::Vector(0,0,1);
    JD::Vector zaxis = JD::Vector(0,1,0);
    JD::Vector xaxis = JD::Vector(1,0,0);
    up.rotate(xaxis, lastXRotation);
    up.rotate(yaxis, lastYRotation);
    up.rotate(zaxis, lastZRotation);
    viewing.rotate(xaxis, lastXRotation);
    viewing.rotate(yaxis, lastYRotation);
    viewing.rotate(zaxis, lastZRotation);
}

in a way that avoids gimbal lock.

Answer 1

If your device is giving you absolute X/Y/Z angles (which implies something like actual gimbals), it will have some specific sequence to describe what order the rotations occur in.

Since you say that "the order doesn't matter", this suggests your device is something like (almost certainly?) a 3-axis rate gyro, and you're getting differential angles. In this case, you want to combine your 3 differential angles into a rotation vector, and use this to update an orientation quaternion, as follows:

given differential angles (in radians):
  dXrot, dYrot, dZrot

and current orientation quaternion Q such that:
  {r=0, ijk=rot(v)} = Q {r=0, ijk=v} Q*

construct an update quaternion:
  dQ = {r=1, i=dXrot/2, j=dYrot/2, k=dZrot/2}

and update your orientation:
  Q' = normalize( quaternion_multiply(dQ, Q) )

Note that dQ is only a crude approximation of a unit quaternion (which makes the normalize() operation more important than usual). However, if your differential angles are not large, it is actually quite a good approximation. Even if your differential angles are large, this simple approximation makes less nonsense than many other things you could do. If you have problems with large differential angles, you might try adding a quadratic correction to improve your accuracy (as described in the third section).

However, a more likely problem is that any kind of repeated update like this tends to drift, simply from accumulated arithmetic error if nothing else. Also, your physical sensors will have bias -- e.g., your rate gyros will have offsets which, if not corrected for, will cause your orientation estimate Q to precess slowly. If this kind of drift matters to your application, you will need some way to detect/correct it if you want to maintain a stable system.

---

If you do have a problem with large differential angles, there is a trigonometric formula for computing an exact update quaternion dQ. The assumption is that the total rotation angle should be linearly proportional to the magnitude of the input vector; given this, you can compute an exact update quaternion as follows:

given differential half-angle vector (in radians):
  dV = (dXrot, dYrot, dZrot)/2

compute the magnitude of the half-angle vector:
  theta = |dV| = 0.5 * sqrt(dXrot^2 + dYrot^2 + dZrot^2)

then the update quaternion, as used above, is:
  dQ = {r=cos(theta), ijk=dV*sin(theta)/theta}
     = {r=cos(theta), ijk=normalize(dV)*sin(theta)}

Note that directly computing either sin(theta)/theta ornormalize(dV) is is singular near zero, but the limit value of vector ijk near zero is simply ijk = dV = (dXrot,dYrot,dZrot), as in the approximation from the first section. If you do compute your update quaternion this way, the straightforward method is to check for this, and use the approximation for small theta (for which it is an extremely good approximation!).

---

Finally, another approach is to use a Taylor expansion for cos(theta) and sin(theta)/theta. This is an intermediate approach -- an improved approximation that increases the range of accuracy:

cos(x)    ~  1 - x^2/2 + x^4/24 - x^6/720 ...
sin(x)/x  ~  1 - x^2/6 + x^4/120 - x^6/5040 ...

So, the "quadratic correction" mentioned in the first section is:

dQ = {r=1-theta*theta*(1.0/2), ijk=dV*(1-theta*theta*(1.0/6))}
Q' = normalize( quaternion_multiply(dQ, Q) )

Additional terms will extend the accurate range of the approximation, but if you need more than +/-90 degrees per update, you should probably use the exact trig functions described in the second section. You could also use a Taylor expansion in combination with the exact trigonometric solution -- it may be helpful by allowing you to switch seamlessly between the approximation and the exact formula.

Answer 2

I think that the 'gimbal lock' is not a problem of computations/mathematics but rather a problem of some physical devices.

Given that you can represent any orientation with XYZ rotations, then even at the 'gimbal lock point' there is a XYZ representation for any imaginable orientation change. Your physical gimbal may be not able to rotate this way, but the mathematics still works :).

The only problem here is your input device - if it's gimbal then it can lock, but you didn't give any details on that.

---

EDIT: OK, so after you added a function I think I see what you need. The function is perfectly correct. But sadly, you just can't get a nice and easy, continuous way of orientation edition using XYZ axis rotations. I haven't seen such solution even in professional 3D packages.

The only thing that comes to my mind is to treat your input like a steering in aeroplane - you just have some initial orientation and you can rotate it around X, Y or Z axis by some amount. Then you store the new orientation and clear your inputs. Rotations in 3DMax/Maya/Blender are done the same way.

If you give us more info about real-world usage you want to achieve we may get some better ideas.