How to find vector for the quaternion from X Y Z rotations

Question

I am creating a very simple project on OpenGL and I'm stuck with rotations. I am trying to rotate an object indepentdently in all 3 axes: X, Y, and Z. I've had sleepless nights due to the "gimbal lock" problem after rotating about one axis. I've then learned that quaternions would solve my problem. I've researched about quaternions and implementd it, but I havent't been able to convert my rotations to quaternions. For example, if I want to rotate around Z axis 90 degrees, I just create the {0,0,1} vector for my quaternion and rotate it around that axis 90 degrees using the code here:

http://iphonedevelopment.blogspot.com/2009/06/opengl-es-from-ground-up-part-7_04.html (the most complicated matrix towards the bottom)

That's ok for one vector, but, say, I first want to rotate 90 degrees around Z, then 90 degrees around X (just as an example). What vector do I need to pass in? How do I calculate that vector. I am not good with matrices and trigonometry (I know the basics and the general rules, but I'm just not a whiz) but I need to get this done. There are LOTS of tutorials about quaternions, but I seem to understand none (or they don't answer my question). I just need to learn to construct the vector for rotations around more than one axis combined.

UPDATE: I've found this nice page about quaternions and decided to implement them this way: http://www.cprogramming.com/tutorial/3d/quaternions.html

Here is my code for quaternion multiplication:

void cube::quatmul(float* q1, float* q2, float* resultRef){
 float w = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3];
 float x = q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2];
 float y = q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1];
 float z = q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0];

 resultRef[0] = w;
 resultRef[1] = x;
 resultRef[2] = y;
 resultRef[3] = z;
}

Here is my code for applying a quaternion to my modelview matrix (I have a tmodelview variable that is my target modelview matrix):

void cube::applyquat(){
 float& x = quaternion[1];
 float& y = quaternion[2];
 float& z = quaternion[3];
 float& w = quaternion[0];
 float magnitude = sqrtf(w * w + x * x + y * y + z * z);
 if(magnitude == 0){
     x = 1;
     w = y = z = 0;
 }else
 if(magnitude != 1){
     x /= magnitude;
     y /= magnitude;
     z /= magnitude;
     w /= magnitude;
  }

 tmodelview[0] = 1 - (2 * y * y) - (2 * z * z);
 tmodelview[1] = 2 * x * y + 2 * w * z;
 tmodelview[2] = 2 * x * z - 2 * w * y;
 tmodelview[3] = 0; 

 tmodelview[4] = 2 * x * y - 2 * w * z;
 tmodelview[5] = 1 - (2 * x * x) - (2 * z * z);
 tmodelview[6] = 2 * y * z - 2 * w * x;
 tmodelview[7] = 0;

 tmodelview[8] = 2 * x * z + 2 * w * y;
 tmodelview[9] = 2 * y * z + 2 * w * x;
 tmodelview[10] = 1 - (2 * x * x) - (2 * y * y);
 tmodelview[11] = 0;

 glMatrixMode(GL_MODELVIEW);
 glPushMatrix();
 glLoadMatrixf(tmodelview);
 glMultMatrixf(modelview);
 glGetFloatv(GL_MODELVIEW_MATRIX, tmodelview);
 glPopMatrix();
}

And my code for rotation (that I call externally), where quaternion is a class variable of the cube:

void cube::rotatex(int angle){
 float quat[4];
 float ang = angle * PI / 180.0;
 quat[0] = cosf(ang / 2);
 quat[1] = sinf(ang/2);
 quat[2] = 0;
 quat[3] = 0;
 quatmul(quat, quaternion, quaternion);
 applyquat();
}

void cube::rotatey(int angle){
 float quat[4];
 float ang = angle * PI / 180.0;
 quat[0] = cosf(ang / 2);
 quat[1] = 0;
 quat[2] = sinf(ang/2);
 quat[3] = 0;
 quatmul(quat, quaternion, quaternion);
 applyquat();
}

void cube::rotatez(int angle){
 float quat[4];
 float ang = angle * PI / 180.0;
 quat[0] = cosf(ang / 2);
 quat[1] = 0;
 quat[2] = 0;
 quat[3] = sinf(ang/2);
 quatmul(quat, quaternion, quaternion);
 applyquat();
}

I call, say rotatex, for 10-11 times for rotating only 1 degree, but my cube gets rotated almost 90 degrees after 10-11 times of 1 degree, which doesn't make sense. Also, after calling rotation functions in different axes, My cube gets skewed, gets 2 dimensional, and disappears (a column in modelview matrix becomes all zeros) irreversibly, which obviously shouldn't be happening with a correct implementation of the quaternions.

Answer 1

You're approaching this the wrong way. If you have three Euler angle rotations for the X, Y, and Z axes, turning them into a quaternion and then into a matrix won't help you. Gimbal lock arises because of the representation of your desired rotation. If you store the rotation that you want as X-Y-Z Euler angles, then you will get Gimbal lock.

You need to store your desired orientation as a quaternion to get the advantages. That is, it's possible to take a current orientation as a quaternion and then ask the question "how do I rotate that around the Z axis by 90 degrees and use that as my new orientation?", but it's not useful to ask "my current orientation is defined by these X-Y-Z Euler angles, how do I turn that into a quaternion?".

A full treatment of the relevant parts of quaternions would be pretty lengthy. It's worth noting that the site you linked to appears to really be talking about axis-angle rotations, not quaternions.

Edit: The code you posted is correct except that the signs for tmodelview[6] and tmodelview[9] are wrong.