numerical drift in quaternion rotations

hi all,

I'm using quaternions to make a lot of small rotations, but after quite a few (10's) rotations the returned vector starts taking silly values.

The quaternion itself is generated fresh each time, but the resulting rotated vector is preserved and repeatedly rotated.

I'm considering the possibility that the rotated vector is drifting of from its unit length and eventually runs away to crazy values - does anyone have any experience with how many times you can rotate a vector before its numerical drift becomes significant?
I'm using double precision (64bit) reals as variables (its actually a physics simulation rather than a game, but game programmers seem to know the most about quaternions!)

thanks in advance :)

Maybe calculate the total rotation from time=0 to now, rather than rely on a sum of deltas.

That way, you only have one calculation error rather than accumulating n calculation errors.

I'm accruing scattering events through the simulation, which are generated randomly from a scatter distribution. So the final oritentation is history dependant and that history is rather too long to store (>570 rotations for possibly over 1M particles).

I think it probably is a numerical drift error - i've watched the code exectute in the debugger and noticed if I normalise the direction after each rotation sometimes the very last decimal place or two will change.

Obviously, normalising the direction prevents run away by itself even if the rotation code is completely wrong, but since the change is only ever in the last few decimal places (and the simulation is behaving reasonably now) i'm going to assume it is numerical drift.

I would have to see some of your code. I assume you are rotating the vectors using a homomorphism, eg:

q' = q * v * q-1

where q is the axis/angle representation of the quaternion, with components:

<cos(angle/2), axis * sin(angle/2)>

'v' is the quaternion made from the vector you are trying to rotate, made of components:

<0, vx,vy,vz>

A better approach is to store the total angle about each axis, and rotate it fully from the initial vector

If you're using angular velocities, this is what I do:

Make sure to Normalize the quaternion every frame.

I use quaternions to integrate the orientation with respect to time. Note that a typical mathbook will tell you to integrate quaternions with respect to the formula:

q' = q*w*.5*dt

where w is the quaternion with components <0,wx,wy,wz>

and w is the angular velocity vector

I do something a little differently, rather just concatenating the rotations into a single quaternion, instead of using that formula (I found it to be much more accurate)

Code:

---

        //PREDICTED ORIENTATION (integrate quaternion)
        Vector3        Axis = CurrentState.mAngularVelocity;
       
        double        ang_vel_mag = Axis.GetLength();
        double        ang_accel_mag = (mJInverseGlobal * Axis).GetLength();
        Axis.Normalize(ang_vel_mag);
       
#if        USE_ACCEL_TERM
        double        Angle = ang_vel_mag * DT + (.5f * ang_accel_mag * DT * DT);
#else
        double        Angle = ang_vel_mag * DT;
#endif
       
        Quaternion        dQ(Axis,Angle);

        CurrentState.mQuat_Orientation.Normalize();

        PredictedStateVector.mQuat_Orientation        =        CurrentState.mQuat_Orientation * dQ;
        PredictedStateVector.mQuat_Orientation.Normalize();

---

You must re-normalize normals and/or vectors at times during rotations and/or vectors used in light calculations because the FPU is very subject to rounding errors thus causing large numerical drift. I'm not sure how you are approaching this but you should be only rotating from the original amount to the current amount. Cumulative rotations on numbers that have already been processed through matrices will result in very strange values.

You cannot do cumulative rotations in 3D- that is, on values that have already been through the pipeline previously. You must start from the origin of the vertex, change the rotation value, and then rotate according to that amount. Wash, rinse, repeat.
The same holds true for translation. Everything starts at the origin and then is transformed into different spaces.

In other words you cannot rotate an object by say .65 radians and then use those values to rotate another .65 radians. Rather you start at 0 radians, rotate at .65 first and render. Then you go back to 0 radians and now rotate 1.70 radians and then render.

This is a common misconception in 3D and in computations. The world always starts as-is and then you transform that world to represent the orientation of the camera/user.

Object/Model space->World space->Camera space->Clip space->Screen space

Some pipelines add other spaces for effects but this is the general approach.

Note every frame goes through this process. You cannot transform once to World space and then use those transformed coordinates in the next frame. You must always start in Object/Model space or you will end up with a gigantic mess.

Thanks for the suggestions. As I said, this is a physics simulation not a game so I don't have to transform through world->view->projection each time, I just work in world space with physically real units defined.


I've had to rewrite my code to use the Quaternions to time integrate the scattering process (thanks BobMcGee123) - applying consecutive quaternion rotations resulted in QNAN errors in a small number of cases that I couldn't trace. It seems to work nicely now.


I don't really disagree with what you say Bubba about applying consecutive rotations, although I think it should be possible in principle since that is essentially what you do when you concatenate matrices provided you recognise the non-commutivity of 3D rotations.
Further more, in principle, the vector (0,1,0) should be indentical to the vector (1,0,0) rotate 90 degrees about Z. 3D math is not history dependant (I suppose you could argue non-commutivity means it is, but lets put that aside) so at a later time rotating (0,1,0) by any rotation should be no different to rotating the vector that started life as (1,0,0) and was rotated.

Having said that, consecutive rotations do regularly seem to blow up on me!

Since at the heart of most modern games are physics simulations and since physics simulations also follow the same principles as games in that they have a frame rate, you should be transforming in exactly the same way per frame.

Consecutive rotations on transformed objects or objects that have been through the simulation code once, will result in a mess.

Floating point is extremely inaccurate by nature. Double data types do a bit better but eventually the rounding errors present in all calculations will bring your rotation code to it's knees, and possible into the realm of NaN's.

I suggest you google this to see what the rest of the internet has to say about it. As well I would recommend buying some books on physics, not for the equations per se, but for the transformation process.

I've got an excess of physics books, especially since I just finished my physics masters degree today :D

(Accumulating in the quaternion worked, accumulating in the vector did not. guess thats the bottom line for here)

Congrats on the degree.

But we judge people here based on what they know and demonstrate, not by how many letters are behind their name.

Thats fine. All I was saying is I already have more physics books than I care to own :)

Essentially Bob and I are saying the same thing (again) from different perspectives.

Quote:

---

A better approach is to store the total angle about each axis, and rotate it fully from the initial vector

---

I agree with Bob's suggestion, and the method Bubba explained is the only way to insure that you keep away from rounding since you recalculate the points each time.

Based on your first post though, maybe it would help to make sure you set the vector back to its original magnitude after the transform. Like normalizing, but not keeping the length to 1.

Code:

---

double mag = vec.GetMagnitude();
quaternion.Transform( vec );
vec.SetMagnitude( mag );

---

Quote:

---

Originally Posted by reanimated

I'm using quaternions to make a lot of small rotations, but after quite a few (10's) rotations the returned vector starts taking silly values.

---

It's strange that you get silly values after so few rotations.

It did seem odd to me too, which was why i didn't want to put it down to numrical drift immediately.
The reason could be that every other rotation was only a few milli-radians in magnitude - its a small rotation so shouldn't matter so much to the end result, but conversely the smaller numbers need a higher degree of accuracy to get the same end accuracy compared to a rotation of a radian.

Data type usage would be my guess.

Floats are highly inaccurate with small values. Double's are less so but still prone.

You may need to download a thirdy party data type library which has greater accuracy than can be achieved with the standard data types.

You wouldn't want to compute the flight path of the shuttle from here to the moon using floats or even using doubles. You would need something with more precision.