Direction of flight

Hi all,
I’m writing a little OGL app that simulates flocking behaviour. I got them moving nicely, so I then wanted to get them pointing in the direction they’re flying. Enter “Computer Graphics: Principles and Practice” by Foley, van Dam, et al. I turn to page 221, and lo! there’s a nice looking way of doing exactly what I want. It tells me to create a rotation matrix thus:

| | | | | | 0
|y x DOF| |DOF x (y x DOF)| |DOF| 0
| | | | | | 0
0 0 0 1

where y is the “up” vector (in this case [0,1,0]), DOF is the vector in which the object has just moved, and x means the cross product.
Now, it says that if DOF is near or equal to y then the matrix degenerates, fair enough (what appears to happen is that the matrix becomes a shear matrix, rather than a rotation one). However, I seem to get shearing way before DOF appoaches y, and because normal vectors don’t like shear matrices, I get some odd lighting at steeper angles.
Is anyone familiar with this way of figuring a rotation matrix? If so can you shed any light on why I get a shear matrix and not a rotation one? I know that in a standard OpenGL matrix, elements 4 and 8 control shearing as well as being part of rotation, so I guess these numbers are getting too big (with the other rotational parts getting too small), but is there anything I can do about it?
On a related note, it also says in CG:PP that this method doesn’t handle banking. Whilst I’m not too bothered about this, it would be nice to have my birds banking as well as pitching and yawing. Any ideas?
I’ll post some sample code if people need it.

John

[This message has been edited by JohnD (edited 05-17-2000).]

Remember an article I read some time ago. Was a quite interesting one. It was about eulerangles and what kind of problems that may occur (and how to solve them too ) when using rotations about more than one axis. The author had an example about a planes direction/rotation. Think you should read it and see if you can come up with something.
http://www.3dgamedev.com/articles/eulers_are_evil.htm

Thanks for the response, Bob,
you’re right it is an interesting article, but I’m already doing just what he says! Well, sort of. He assumes that the object’s DOF comes from Euler angles, whereas mine is just a vector. It does appear, though, that in priciple I’m doing it right.
I am interested in the quaternion solution (sounds like a title for a cheesy thriller movie ). I’ve read about them before, and got a little confused about halfway through the explanation (my maths isn’t as good as it should be). Does anyone have any URLs for an OpenGL oriented idiot’s guide to quaternions?

Big thanks,
John

Ah! No need for quat URLs, I just found one on another thread.
Sorry for any inconvenience

He Who Now Looks Slightly Embarassed
-John

I came up with a similar problem in my app. I just figured out a quick solution a few days ago. Note, its not perfect, and I’m ironing out details still, but it seems to work.

Here’s a bit of code…

glRotatef(glAngleX, yAxisTilt, xAxisTilt, 0.0f);
glMultMatrixf(mvmatrix);
glGetFloatv(GL_MODELVIEW_MATRIX, mvmatrix);

First, load the identity matrix into modelview (otherwise you’re doing some really funky matrixmults).
glTranslatef(0.0f, 0.0f, -200.0f);
glMultMatrixf(mvmatrix);

Make your rotation (based on your vector) from that. Note, your vector must be broken down into components, which is fairly easy to do.

Next, multiply that rotation buy the previous rotation matrix (stored in mvmatrix here). This will “accumulate” the previous rotations with this one. Finally, store the current rotation matrix for the next rotation rendering pass.

You’ll note the second identity and translation, this simply is post rotational positioning, so the object is viewable (for my app… its just an example here).

The thing to note is that all multiplication must start from the identity matrix. Otherwise you’re multiplying ALL accumulations (rotations, translations, scalings) together, which does some pretty freaky stuff.

Its simple, and easy.

Hope this helps.

D.

Siwko, aren’t you effectively using Euler angles there?

Besides, that glGetFloatv(GL_MODELVIEW_MATRIX,mvmatrix) is probably going to hurt your performance quite a bit.

no, because that process rotates incrementally by small amounts every frame, and don’t exhibit any axial constrains.
in other words, it’s isotropic.

it can be translated to euler angle, but it’s a bit tricky.

i use basically the same method, but instead of querying opengl for the mv matrix, i keep the status into a class, wich is surely better.

hey siwko and the other users: remember to re-orthogonalize the mvmatrix sometime…

Dolo//\ightY

[This message has been edited by dmy (edited 05-17-2000).]

Thanks for the responses,
I have just one question about your method, Siwko, what are the three parameters in your glRotatef(glAngleX, yAxisTilt, xAxisTilt, 0.0f)?
You said to get the vector’s components, but I can’t understand (not enough caffiene, perhaps?) what glAngleX is, nor can I see why you have yAxisTilt before xAxisTilt (I assume these are just the x and y components of the vector).

He Who Wants A glRotateInTheDirectionIWant() function.

-John

[This message has been edited by JohnD (edited 05-18-2000).]

Originally posted by JohnD:
Thanks for the responses,
what are the three parameters in your glRotatef(glAngleX, yAxisTilt, xAxisTilt, 0.0f)?

Well, I can get a bit vague here. The variable names may be misleading because I was reusing them from a failed previous attempt at getting this working correctly.

glAngleX (which should just be glAngle now) is a constant angle value in which you want to rotate in. IE: set (@ dimensionally speaking for simplicity), you’re heading is 0 degrees, you want to turn 90 degrees to your right, it would be +90.0, turn 90 degrees left, -90.0, etc.

Next, I’m assuming you understand the glRotatef() function. The next three parameters specify the vector axis you wish to apply the rotation about. In otherwords, its specified by the vector (x,y,z) which when multiplied together gives you your vector --> this those are your vector components. So:

glRotatef(45.0f, 1.0f, 1.0f, 0.0f)

will rotate you 45.0 degrees clockwise looking down a vector that can be described by (1,1,0), or a vector extending out 45 degrees (no relevance to the 45.0f angle specified) up from the positive x and right from the positive y axis.

You’ll notice the 0.0f for the z component, but thats just because I don’t care to include any z axis rotations. For my app, the user interface is only 2D anyways.

I don’t know if that makes anything clearer.

And as far as it being euler angles to the above post (not yours), no. This has the cumulative effect in all 3 axises, because it is specified as a rotation about a single vector, not multiple rotations about individual vectors.

As for the performance of the glMultMatrix call, from what I understand, doing your own rotations this way should be faster. If you’re using euler angles, you have to call glRotate 3 times, which essentially translates in to 3 matrix mults. In my app, the performance hit, if any, is completely unnoticable. (I’m moving usually 100k triangles around at any given time, and I’m managing respectable frame rates).

Any more questions on how I got it working, let me know. I’m happy to help, especially since I’ve only been doing opengl for 2 months.

Originally posted by dmy:
[b]it can be translated to euler angle, but it’s a bit tricky.

i use basically the same method, but instead of querying opengl for the mv matrix, i keep the status into a class, wich is surely better.
[/b]

I believe that any type of rotational algorithms are capable of being turned into eulers, but once you start accumulating them, euler breaks down, correct?

Also, I’m interested in how you maintain the matrix status in the class. Can you let me in on that?

And finally,

hey siwko and the other users: remember to re-orthogonalize the mvmatrix sometime…

Enlighten me. What mean you by this?

Siwko

(Apologies in advance for long post.)
From what I can tell, Siwko, your code gets a vector and rotates an amount around that. However, I have my object lying flat, pointing in the -z direction, and I want to rotate it to point in the direction of my vector, which is calculated from scratch every frame. So at that point, I don’t know by how much to rotate it. Because the DOF vector is recalculated every time, I don’t want to store the previous rotation matrix, then multiply it with the current one.
BTW, my DOF vector is calculated by adding the vector from the boid to the centre of the flock and the vector from the boid to the “goal”, and others such as collision avoidance vectors.

void CBoid::Move() {
double px, py, pz;
double vx, vy, vz;

// Stay with the centre of the flock.
PerceiveCentre();
m_seekCentre -= m_position;
m_seekCentre.Scale(1.0 / 128.0);

// Try to keep the same speed as the rest of the flock.
AverageVelocity();
m_matchVelocity -= m_velocity;
m_matchVelocity.Scale(1.0 / 40.0);

// Move away from any other boids.
AvoidBoid();
m_avoidBoid.Scale(1.0 / 32.0);

m_seekGoal = m_theFlock->GetGoal();
m_seekGoal.Scale(1.0 / 128.0);

m_velocity += (m_seekCentre + m_matchVelocity + m_avoidBoid + m_seekGoal);
m_velocity.Limit(20.0);
m_position += m_velocity;
}

void CBoid::PerceiveCentre() {
m_seekCentre = m_theFlock->GetCentre() - m_position;
}

void CBoid::AverageVelocity() {
m_matchVelocity = m_theFlock->GetVelocity() - m_velocity;
}

void CBoid::AvoidBoid() {
std::vector::iterator iter;

for (iter = m_theFlock->m_boidList.begin(); iter != m_theFlock->m_boidList.end(); ++iter) {
if (m_position.RectangularDistance(iter->GetPosition()) <= 150.0) {
m_avoidBoid += (m_position - iter->GetPosition());
}
}
}

Here’s my interpretation of the algorithm I got from “CG:P&P”:

// For each boid…
glPushMatrix();
// Move the boid.
glTranslated(px, py, pz);

// Caclulate the DOF vector (see “Computer Graphics: Priciples and Practice” pp. 221-222).
c1 = m_up.Cross(vec);
c2 = vec.Cross(c1);

c1.Get(rot[0], rot[1], rot[2]);
c2.Get(rot[4], rot[5], rot[6]);
vec.Get(rot[8], rot[9], rot[10]);

// Multiply the current matrix by our DOF matrix.
glMultMatrixd(rot);

// Draw me!
glCallList(iter->m_geometryId);
glPopMatrix();

I’m still sure quats are the way to go, but I’m not sure how to go from a vector to a quat (it doesn’t have that in the quat FAQ). I’d also be interested in hearing about re-orthogonalising my matrix.

He Who Thought Re-orthogonalising Involved Surgery.
-John

[This message has been edited by JohnD (edited 05-19-2000).]

I think reorthogonalising is making sure your x, y and z axis are perpicular to each other.

So, its basically 2 crossproducts and some normalizing.

I.e take the crossproduct of x and y axis to get a new z axis. Then take the crossproduct of x and new z axis to get a new y axis… And then normalize new y and new z axis…

Well, thats the way I do it.

Hmm, okay… Ii understand your pain JohnD. I didn’t think your simulation could really involve an about-face though (180 degrees), which essentially is what your saying.

Hmm, you do have a valid point though, my code may not work for you in that case.

As far as quats go, someone posted a url for me on that subject a couple weeks ago, I’ll see if I can find it for you. It had the full code to do quat calculations.

Get back to you with that sometime today.

I.e take the crossproduct of x and y axis to get a new z axis. Then take the crossproduct of x and new z axis to get a new y axis… And then normalize new y and new z axis…

I was under the impression that taking the cross product of two normalised vectors also produced a normalised vector. Am I wrong? If so, then that could be the cause of my problem.
BTW, thanks Siwko (and everyone else) for your help.

He Who’s Impressions Have Got Him Into Trouble

-John

*** STOP PRESS *** STOP PRESS *** STOP PRESS *** STOP PRESS ***
Major apologies for my lack of prowess in the brain department
I tried that normalisation of each cross product vector, and lo! it works a treat!

Thanks to dmy for suggesting the ortho-normalisation thing, AndersO for telling us less intelligent people how to do it, and Siwko for answering in the first place.

He Who Now Feels REALLY Silly.

-John

[This message has been edited by JohnD (edited 05-19-2000).]

Now if only I could understand the cross-product thing myself. I don’t really understand how it applies and what its for.

Enlighten me?

Also, I don’t quite know if it’s too relevant to my situation, as I really am only concerned with two axises for my purposes anyways (x and y).

My goal is always rotate about the coordinate system’s axises, not the model’s.

Am I doing this right? Or is it going to break in the long run?

Well, I can tell you what the cross product is for.
I’m sure you know that one way of defining a plane in 3-space is by using 3 points. Now, imagine that one of those points is the origin, with a vector heading out to each of the two other points (so you have a sort of V shape). The cross product of those two vectors is another vector perpendicular them (this is the “normal” vector), so it points straight out from the plane. It is usually used for lighing calculations in OpenGL, but it can also be used for back-face culling. That is, if we know the normal is pointing away from the camera, we don’t need to draw the polygon.
In the code that I had, I assumed that a cross product of two unit-length vectors would result in a third unit-length vector. Of course, when I think about it, that’s a silly idea, but then who ever said that developers were clever? In my defense, I would like to point out that in the book I got it from, it didn’t say anything about renormalisation, so it’s the book’s fault, not mine.
Are you developing a 3D application, Siwko? If so, why are you only interested in the x & y axis? Without knowing more, I would say that to rotate about the world axis, you’d do your rotations first, then your translations, i.e.

glRotatef(anAngle, 1.0f, 0.0f, 0.0f); // Rotate about the x axis.
glRotatef(anotherAngle, 0.0f, 1.0f, 0.0f); // Rotate about the y axis.
glTranslatef(x, y, z);
DrawObject();
glFlush();

He Who Gets Dizzy Just Thinking About Transformations.
-John

if we were avatars into a perfect digital world the IEEE floating point format would be error-less.

but instead we are all in 2K and floating point operations exhibits little roundoff errors: if i have, say, a vector, and i rotate it incrementally:
1: take the vector components
2: multiply them with a rotation matrix
3: store the results as the new vector
4: goto

…by a small amount, the errors involved will sum up over and over, and after some (actually many) rotations i’ll end up with a vector with a modulus (lenght of the vector) wich is different than the original.

what will happen to my objects, then? scaling and shearing, like the starship enterprise when it entered the shunt space distortion made by mr barkley from the holodeck.

so, by re-orthonormalizing the matrix wich holds the values for my three unit vectors (basically the modelview matrix without the fourth row and column) will then be correct again, every vector will be of lenght 1 and they will be perpendicular each other.

the right way should be to evaluate the current error into the matrix, and if it is over a threshold, re-normalize…
i read this into Foley et al CG:PnP.

but also to do normalization on a regular basis of 30-60-100 frames will do good.
it depends on the application.

basically i keep into a class two matrices, one for translations (xlt), one for rotations (rot).
when i have to translate, i modify only the xlt matrix, so into xlt there will always the world location. this is useful to me because one of my projects is a space simulator.
when i have to model rotations, i rotate the rot matrix with custom functions.
then, it’s time to send it all to opengl, and i do it by:

glMatrixMode(GL_MODELVIEW);
glMultMatrix(xlt);
glMultMatrix(rot);
<render>

when you reach the singularities you get the gimbal lock, and things will start to shake… and you won’t get out safe of it without some more work!
yes, euler angles are good for things like camera-radar-cannon pointing, but not for things wich float free into space (like spaceships )
If i’m not wrong, i remember the first BBC Elite by David Braeben and Ian Bell: your wireframe cobra ship, as well the others, was all affected by euler angles problems.

whew, that was long!
over-n-out.

-He Who Scored Two Kills Into A Soft-Air Combat Yesterday, But He Was Also Shot Down By An Allied Because His Clan Is Made Of Dead-Meat-Buddies Without A Common-Band Radio
-DMY

Dolo//\ightY

PS: why not create a game about soft-air? (let’s not count Extreme Paintball…

Originally posted by JohnD:
Are you developing a 3D application, Siwko? If so, why are you only interested in the x & y axis? Without knowing more, I would say that to rotate about the world axis, you’d do your rotations first, then your translations

I do rotate then translate. I don’t know if I posted the code right then. Mind you, I look at the from the perspective that I’m not rotating the world, just the model. I know that this is the same as the inverse rotation of the world. But now its beginning to sound perverse! j/k

And yes, I am developing an app. Really, all the rotations serve are to rotation the model in response to user input. Just a single model, not the entire world. So that section is encapsulated in a push/pop for me.

As far as normalization goes, err, how long does it take before the effects of non-normalization start showing up? I’m curious, because I have not once seen them over the past two weeks.

Just curious.

Siwko

PS - I’ll probably pester you guys about normalization in the next couple weeks. Or at least how to.

On the subject of when rounding errors creep in, you need to know what the epsilon value of a float is. If you look under FLT_EPSILON in MSVC 5 docs, you’ll see it says: 1.192092896e-07
This is the value such that 1.0 + FLT_EPSILON != 1.0, in other words this is the smallest difference between two floats.
Now I could imagine in one frame, a float could be modified (added to, multiplied by, and so on) 5 or more times. Lets say you’re running at 30fps, this equates to 150 modifications per second. So in 10 seconds this equates to roughly 2 thousandths difference between what the number is and what it should be. Which is starting to become significant. This is a fairly conservative estimate, because in one matrix there are 16 such floats, all having a knock-on effect on each other.
Obviously, doubles are more accurate, with an epsilon of 2.2204460492503131e-016 which is a very small number, but still could be significant in, say, a flight sim which is running for 18 hours a day.
Added to the epsilon problem, is that many numbers cannot be stored accurately in a float or a double, so you end up with rounding errors before you even start!

He Who Blames His Increasing Waistline On Rounding Errors.
-JohnD