Gimbal lock, despite using Matrices and Quaternions

[Davespice]

Hi guys;
I'm having some trouble animating a 3D model. I am just trying to do a simple rotation on x, y and z through 360 degrees. So my goal is to be able to pitch, roll and yaw the model through any rotation angle. The code below uses a combination of matrices and quaternions to accomplish this. But despite that the resulting animation appears to have gimbal lock. It starts to spin correctly but then suddenly slows down and appears to go backwards for a moment before carrying on.

You'll see that I am trying to use pitch, roll and yaw angles to do this. Maybe it is that I need to completely get away from using these kinds of angles and come at this from another perspective. Anyway, I thought I would post my code here and see if people can point out to me why it's not working and what I really ought to be doing.

Here is the sudo code in brief;

The model data that I load comes out upside-down initially so first thing I do is flip 180 degrees on z.

I then do the translation to push the model away from the viewport a little.

I then do three stages of rotation, I extract the look, up and right vectors from the matrix create quaternions based on those vectors and my stored values for pitch roll and yaw. I then convert the quaternion back into a matrix and multiply the current matrix by it. But I do it sequentially, I apply the multiplication for each quaternion-matrix before getting the vector values for the next one from the resulting matrix. I also have to invert the quaternions to make the pitch, roll and yaw angles move in the direct direction since I originally did a flip on z.

So despite all this I still have gimbal lock, any advice or help is welcome - even if it's extended use of the "don't bat"

Code: Select all

static GLfloat scale = 1.0f;
static GLfloat pitch = 0.0f;
static GLfloat yaw = 0.0f;
static GLfloat roll = 0.0f;

static void drawView(STATE_T *state, EGL_NSH_MODEL *model)
{
	glLoadIdentity();

	pitch = (pitch == 360) ? 0 : pitch + 1;
	yaw = (yaw == 360) ? 0 : yaw + 1;
	roll = (roll == 360) ? 0 : roll + 1;
		
	Matrix3D identityMatrix;
	Matrix3DSetIdentity(identityMatrix);
	
	Matrix3D flipMatrix;
    Matrix3DSetRotationByDegrees(flipMatrix, 180.0f, 0.0f, 0.0f, 1.0f); //Flip model on z
	
	glMultMatrixf(flipMatrix);
	
    Matrix3D transScaleMatrix;
    Matrix3DSetTranslationScaling(transScaleMatrix, scale, scale, scale, 0.0f, 0.0f, -25.0);
	
	glMultMatrixf(transScaleMatrix);
	
	Matrix3D transformedMatrix;
	glGetFloatv(GL_MODELVIEW_MATRIX, transformedMatrix);
			
	//-------------
	
	Vector3D rollVector;
	rollVector.x =	transformedMatrix[2];
	rollVector.y =	transformedMatrix[6];
	rollVector.z =	transformedMatrix[10];
	
	Quaternion3D rollQuat;
	rollQuat = Quaternion3DMakeWithAxisAndAngle(rollVector, DegreesToRadians(roll));	
	Quaternion3DInvert(&rollQuat);
	
	Matrix3D rollMatrix;
	Matrix3DSetUsingQuaternion3D(rollMatrix, rollQuat);	
	glMultMatrixf(rollMatrix);
	
	Matrix3D rollResultMatrix;	
	glGetFloatv(GL_MODELVIEW_MATRIX, rollResultMatrix);
	
	//-------------
	
	Vector3D yawVector;
	yawVector.x = rollResultMatrix[1];
	yawVector.y = rollResultMatrix[5];
	yawVector.z = rollResultMatrix[9];
	
	Quaternion3D yawQuat;
	yawQuat = Quaternion3DMakeWithAxisAndAngle(yawVector, DegreesToRadians(yaw));
	Quaternion3DInvert(&yawQuat);
	
	Matrix3D yawMatrix;
	Matrix3DSetUsingQuaternion3D(yawMatrix, yawQuat);	
	glMultMatrixf(yawMatrix);
	
	Matrix3D rollYawResultMatrix;	
	glGetFloatv(GL_MODELVIEW_MATRIX, rollYawResultMatrix);
	
	//-------------
	
	Vector3D pitchVector;
	pitchVector.x =	rollYawResultMatrix[0];
	pitchVector.y =	rollYawResultMatrix[4];
	pitchVector.z =	rollYawResultMatrix[8];

	Quaternion3D pitchQuat;
	pitchQuat = Quaternion3DMakeWithAxisAndAngle(pitchVector, DegreesToRadians(pitch));
	Quaternion3DInvert(&pitchQuat);
	
	Matrix3D pitchMatrix;
	Matrix3DSetUsingQuaternion3D(pitchMatrix, pitchQuat);	
	glMultMatrixf(pitchMatrix);
	
	Matrix3D rollYawPitchResultMatrix;	
	glGetFloatv(GL_MODELVIEW_MATRIX, rollYawPitchResultMatrix);

	//-------------
	
	glClearColor(0.0, 0.0, 0.0, 1.0);
	glClearDepthf(1.0f);
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	
	glEnableClientState(GL_VERTEX_ARRAY);
	glEnableClientState(GL_COLOR_ARRAY);

	glVertexPointer(3, GL_FLOAT, 0, model->vertices);
	glColorPointer(4, GL_FLOAT, 0, model->colors);
		
	glDrawElements(GL_TRIANGLES, model->polysCount, GL_UNSIGNED_SHORT, model->polys);
	
	glLineWidth(2.0f);
	glDrawElements(GL_LINES, model->linesCount, GL_UNSIGNED_SHORT, model->lines);

	glDisableClientState(GL_VERTEX_ARRAY);
	glDisableClientState(GL_COLOR_ARRAY);

	eglSwapBuffers(state->display, state->surface);
}

Thanks!

[jmacey]

Without seeing the maths in your Matrix and Quat code it is difficult to see what is going on, however a couple of things to check

1. In this code

Code: Select all

   
rollVector.x =   transformedMatrix[2];
rollVector.y =   transformedMatrix[6];
rollVector.z =   transformedMatrix[10];

It appears that you are getting the current transform matrix, then adding a new rotation to it, If you are doing this it will make the rotations etc increase each time by a larger amount each time you call the code (for example you would get something like 2, 4 ,8 16 etc which would give the impression of it speeding up / slowing down as you are setting uneven rotations per frame).

2. I'm slightly suspicious of mixing the quaternions and the matrices in the same code, usually you would use all quaternions then extract the matrix once to set the GL values.

3. check the matrix order of your library vs that of OpenGL they could be row / col flipped so you would need to transpose the matrix.

If you can post a link to the rest of the library / code you are using I will have a look.

Jon

[Davespice]

Sure the code I am using for the Matrix and Quat functions come from here and here. Thanks so much for looking at it. Would it help if I posted a video of how the animation looks now?

[jmacey]

Just about to head to work but will have a look later, a video would be good.

The code below shows how I do transformations using just matrices, the m_scale, m_rotation and m_translation are just 3tuple values holding the x,y,z rotation (degrees) scale and translations for each. It's more about the order of setting things to get the final matrix that is important, I then use m_matrix and load it directly to the shader.

Code: Select all

ngl::Matrix mScale;
ngl::Matrix mRotationX;
ngl::Matrix mRotationY;
ngl::Matrix mRotationZ;
ngl::Matrix mTranslate;

// rotation/scale matrix
ngl::Matrix mRotationScale;
// effectivly sets the diagonals of the matrix 
mScale.scale(m_scale.m_x, m_scale.m_y, m_scale.m_z);
// these methods create a matrix based on the rotation value such as (for x)
//Real beta=radians(_deg);
//  Real sr = sin( beta );
//	Real cr = cos( beta );
//	m_11 =  cr;
//	m_21 = -sr;
//	m_12 =  sr;
//	m_22 =  cr;
mRotationX.rotateX(m_rotation.m_x);
mRotationY.rotateY(m_rotation.m_y);
mRotationZ.rotateZ(m_rotation.m_z);
// now combine the matrices in a know way to avoid gimbal lock etc
mRotationScale = mScale * mRotationX * mRotationY * mRotationZ;

// transform matrix
m_matrix = mRotationScale;
// now set translation
m_matrix.m_m[3][0] = m_position.m_x;
m_matrix.m_m[3][1] = m_position.m_y;
m_matrix.m_m[3][2] = m_position.m_z;
m_matrix.m_m[3][3] = 1;

[blu]

Davespice, from the code you posted I think you've missed a small detail about the GL1x matrix stacks, namely:

Each new MultMatrix(m) does current = current * m. Combined with the fact that the final vertex transformation uses operator-leftism, i.e the matrix operator is on the left of the column-vector comprised of the vertex spatial attributes (i.e. in formal notation M * (x, y, z, w)T ), you need to supply the individual transformations in an order opposite to their logical application. In other words, if you transformations' logical order in time is T1, T2, ... Tn, you need to pass them to the GL stack as

Code: Select all

glLoadIdentity();
glMultMatrix(Tn);
...
glMultMatrix(T2);
glMultMatrix(T1);

The above would yield a final transform of Tn * ... * T2 * T1 * (x, y, z, w)T

[Davespice]

Hi blu, thanks for replying.
I'm sorry but what you describe here is well above my head and it doesn't actually mean anything to me in terms of how my code should be. Could you elaborate a bit and use extracts from my code to explain?

Thanks again for trying to help though.

[blu]

Well, let's say the logical order of your transformations is:

1. flip object from its original orientation to the desired 'base' orientation - let's call this transform FLIP
2. roll/yaw/pitch transform - let's call it CARDAN
3. offset of the object away from the view plane - OFFSET

then the proper order of feeding those to the MODELVIEW stack would be:

Code: Select all

glLoadIdentity();
glMultMatrix(OFFSET);
glMultMatrix(CARDAN);
glMutlMatrix(FLIP):

[Janq]

The reason you still have gimbal lock is that you are basically still using euler angles!

Any time you are storing an absolute pitch,yaw,roll (or any other combination of 3 axis rotations) you will have gimbal lock.

The way to avoid it is to store your current rotation as a matrix or quaternion and apply delta rotations to it every frame. Small rotations are order independent and so can be combined in any order without affecting the result (for the most part anyway).

So, you should do something like this (pseudo code):

static Quaternion s_CurrentOrientation;

...
float xr = .. x delta rotation for this frame ..;
float yr = .. y delta rotation for this frame ..;
float zr = .. z delta rotation for this frame ..;
s_CurrentOrientation = s_CurrentOrientation * MakeXRotation(xr) * MakeYRotation(yr) * MakeZRotation(zr);
s_CurrentOrientation = Normalize(s_CurrentOrientation);
...
.. make matrix from s_CurrentOrientation and use for rendering ...

You can probably get away without renormalizing the quaternion but over very long time periods it may drift away from unity so its worth doing.

[Davespice]

Hi all, thanks for the responses. I appreciate your patience, this is by far the most intense programming I have ever done. It’s taking me some time to grasp everything so again, thanks for bearing with me.

Well, let's say the logical order of your transformations is:

1. flip object from its original orientation to the desired 'base' orientation - let's call this transform FLIP
2. roll/yaw/pitch transform - let's call it CARDAN
3. offset of the object away from the view plane - OFFSET

then the proper order of feeding those to the MODELVIEW stack would be:

Code: Select all

glLoadIdentity();
glMultMatrix(OFFSET);
glMultMatrix(CARDAN);
glMutlMatrix(FLIP):

Hi Blu, thanks for this – I understand what you mean now. I’ll remember this, I also like the concept of the matrix stack – that helps my understanding even further. I now am thinking I might just set up all my matrices and then apply them sequentially like you have shown here (instead of after I create each one)… that might make the code easier to follow and read.

The reason you still have gimbal lock is that you are basically still using euler angles!

Any time you are storing an absolute pitch,yaw,roll (or any other combination of 3 axis rotations) you will have gimbal lock.

This is exactly what I have been wondering. Thanks for this! Now I know for definite I need to find another approach to the task.

So bear with me please. Okay so I just want to show this image now, because this is what I have used to work out how to extract the local angles/axis from a matrix.

from Diney Bomfim’s blog.

If you look at this you’ll see how I decided to make my roll, pitch and yaw vectors, which I make the Quats with, in the code from my first post. Anyway... moving on from that.

I take it what you’re saying is that the look vector is what would be the orientation quaternion? As a Quat needs a vector part… so how would I construct that vector? Same as the pic using look? Or where do I get it from? And would the angle part of the Quat then give me the roll? I have a feeling I am losing the plot at this point

Also, your MakeXRotation etc functions – this is basically quaternion multiplication isn’t it?
Sorry for so many questions but I really am starting from scratch and there may be a lot of knowledge that I am missing which stops things from being obvious. I appreciate everyone’s efforts to help me, I really do.

I’m also slightly worried about normalizing the values too often as I know the square root function is one of the slowest things the CPU can perform… so I probably don’t want to be doing it every game frame for instance.

[jmacey]

If you get a chance have a look at my lecture notes here http://nccastaff.bournemouth.ac.uk/jmac ... Camera.pdf This introduces the concept of a Virtual Camera as well as a transformation stack. This is done for Desktop Core Profile OpenGL using a shader to do the final vertex transformation but the principles are more or less the same for OpenGL ES2.0, I presume you are actually using ES1.2 so you just need to load the final matrix using the glMatrixf commands.

Most of the notes are based on the F.S. Hill book which I would recommend as a good introduction to both CG and OpenGL

http://books.google.co.uk/books/about/C ... edir_esc=y

If I get a chance later I will try and upload an example of doing basic transformations using a simple matrix class, I usually only use quaternions when I need to do interpolation of joints etc. I may have some old notes on Quaternion rotation lying around as well but need to find them, basically to do quat rotation you do something like

Code: Select all

void Quaternion::fromAxisAngle(
                               const Vector& _axis,
                               float _angle
															)
{
	Vector axis = _axis;
	axis.normalize();
	_angle=radians(_angle);
	Real sinAngle = static_cast<Real>(sin( _angle / 2.0f ));
	Real cosAngle = static_cast<Real>(cos( _angle / 2.0f ));
	m_s = cosAngle;
	m_x = axis.m_x * sinAngle;
	m_y = axis.m_y * sinAngle;
	m_z = axis.m_z * sinAngle;
}

This set the quaternion from an axis and an angle (in degrees), so for example I would use

Code: Select all

ngl::Quaternion q;
q.fromAxisAngle(ngl::Vector(1,0,0),m_xRot);

To rotate around the x Axis by m_xRot degrees.

Next I rotate each of the points in my mesh by using something like this

Code: Select all

for(int i=0; i<size; ++i)
  {
    p=m_points[i];
    q.rotatePoint(q,p);
    rotPoints.push_back(p.toVec3());
    // we could also use the following code
    //rotPoints.push_back(q*m_points[i]);
  }

Where the rotatePoint method looks like this

Code: Select all

void Quaternion::rotatePoint(
						const Quaternion& _r,
						ngl::Vector & io_p
					       )
{
Quaternion temp = -_r;
Quaternion point(0.0,io_p.m_x, io_p.m_y, io_p.m_z);
point = temp * point * _r;
io_p.set(point.m_x, point.m_y, point.m_z,1.0);
}

The main thing here is that to rotate a point with a quat we need to do

P'=-quat * P * quat; where -quat is the inverse ((aka conjugate) where we invert the vector part but leave the scalar alone)

Hope this helps a bit, quats are quite difficult to get your head around initially, If you really want to get into them I would recommend this book by my colleague John Vince http://www.amazon.com/Quaternions-Compu ... 0857297597 and also this one http://www.amazon.com/Mathematics-Compu ... 809&sr=1-1

Jon

[blu]

Davespice, I'd advise you to leave quaternions aside for now (you can always revisit them later) and focus on the Euler-based transforms. What Janq says about small, relative rotations avoiding the Gimbal Lock is true for any kind of transformation which implements rotations around the axes of a coordinate system.

So let's assume you have an arbitrary rotation (i.e. around an arbitrary axis and of arbitrary magnitude), which we'll call R. Let this rotation keep our object's current orientation in a given basis (i.e. a coordinate system). Using that R as a stepping stone, you can add small Euler-ian rotations on top of your R, and avoid any Gimbal Locks, like this:

R' = R * E3 * E2 * E1
v' = R' * v

where En are the small Euler-ian rotations (i.e. one around each axis), and v is your vertex subject to the final transformation.

Now, what the order of those En rotaions should be is entirely up to you and depends on the local coordinate system of the object, but for an object that 'looks' along the local Z-axis, and whose 'up' is along the local Y-axis, I normally find it useful to apply them in E(y) * E(x) * E(z), which gives a standard azimuth-declination-roll behavior.

[Davespice]

I know you all mean well, but this is really hard to understand. A lot of mathematics terms that I don’t understand are being used here. Just to give you some idea of where my head is at. I started doing this only a month or two ago (in my spare time in the evenings and some weekends) in that time I have managed to get a very simplistic model of the starship Enterprise, a Klingon cruiser and a Tie Fighter spinning around on the screen. I am coming from no previous experience so I think I’ve done okay so far... but I think you’re all assuming I know a lot of the background stuff or have studied this at university. I haven’t and a lot of the help here does baffle me a bit. Please understand I am not ungrateful, I think I just need you to explain everything really explicitly or I probably won’t understand.

Jon, thanks for those uni notes I will go through them.

[jmacey]

I've just uploaded an non graphical demo of a simple matrix class here http://nccastaff.bournemouth.ac.uk/jmac ... tation.tgz

If you extract the file using

tar vfxz MatrixRotation.tgz
cd MatrixRotation
make

It will build a simple demo, this uses a cut down version of my Mat4 matrix class to perform some rotations and scales, as well as setting the translation. The final matrix is printed out, but you would then just use this to load it to OpenGL instead.

the demo program looks like this

Code: Select all

#include "Mat4.h"

int main()
{
	// create some matrices =1 assigns it to be the ident matrix
	Mat4 xRot=1;
	Mat4 yRot=1;
	Mat4 zRot=1;
	Mat4 scale=1;

	Mat4 final;
	// or you can use the identity method
	final.identity();

	xRot.rotateX(45);
	yRot.rotateY(-90);
	zRot.rotateZ(25);
	scale.scale(1,2,1);
	final=scale*xRot*yRot*zRot;
	// now set translate -2 in z
	final.m_m[3][0]=0.0f;
	final.m_m[3][1]=0.0f;
	final.m_m[3][2]=-2.0f;
	final.m_m[3][3]=1.0f;
	std::cout<<final;

}

Hopefully you could use this in your own program to get your final matrix for the transformations, you will still need to setup your perspective projection matrix to make it properly 3D etc but this should make it easier for you.

[Davespice]

Ahh Jon... that is exactly what I need. I can tell you’re a university lecturer

[Davespice]

Hi Jon, I got your code to compile this morning. Works great, thanks for that.
I just needed to modify one include line for the GL library.

[jmacey]

great hope it helps