Fixed point maths, quaternions, rounding errors and normalisation.
 

I’ve just updated the 3D game thing I’m writing to use quaternions to hold the orientation / attitude of the camera and the objects I’m rendering.

It all seems to work and was not difficult to code. I’m reasonably happy with it.

But, I know there are rounding errors accumulating each time I pre-multiply the current quaternion with one that represents the change in orientation. All the literature says to expect this in the form of the quaternion drifting away from unit length, and to rectify it, once it has reached a certain point, by dividing the quaternion by its length, thus returning it to a unit length.

I’m using fixed point maths for everything. Currently 16 bits, 2 before the point, 14 after.

My rounding errors so far are small, but they do exist. They don’t become a problem in the small demo / play about code that I have at the moment, but I do suspect that they’ll cause issues over the time it takes to play the game for real.

I have questions about the use of quaternions in this way.

• Is 2:14 fixed point format enough precision for this? Should I keep my precision when creating the rotation matrix that results from them? If so, how high? 2:30? And should I embiggen my sine table to match?
  
  
• How far should I let the quaternions drift before I normalise them? Inversely, how do I calculate when it’s possible for the normalisation to be beneficial given the limitations of the fixed point maths, including the calculation of the length which involves a square root.

Without knowing your target precision, nobody can answer. But why do you use quaternions for this and not ordinary rotation matrices? It seems to overly complicate the thing (plus, having to take care about both real and imaginary part). For computing the square root in fixed point, I provided an algorithm. It is not hard to do that.

Thomas, I'm not sure how you're trying to help. Would you like me to explain the issues that "ordinary rotation matrices" bring, and how the use of quaternions overcome them, while leading to more efficient code? Or are you confused about the supposedly "imaginary" nature of quaternions? And with regard to square root algorithms, I can point you to one from a trusted, published, source.

I just wonder if ultimately this (quaternions) would be slower than standard 3x3 matrix math on Amiga?

Btw, I work in Unity on a daily basis, and most of the maths is "under the hood" so to speak...

Standard approach is to convert to a matrix before you do all your vertices, so, no. Instead of creating and maintaining a matrix, you just convert, and the number of calculations is the same for the bit that matters.

And "under the hood" Unity use quaternions :)

I can't help you with this because your 3D math is sure much advanced than mine (I can only rotate some stupid cubes, some polygons or some stars in space :D)
And in any case I never used quaternions.
But for the various matrix calculations I used 2:14 and the results always seemed good to me.
If you want to continue using a bare 68k, 2:30 it's not an option, the internal ALU is 16-bit, so even a simple sum would require many more cycles.
I am very curious to see your results with quaternions.

I'm definitely sticking with bare 68k.

Maybe 2:30 can be an option if it's restricted to a tiny section of code where the performance isn't critical though? Maybe it's necessary for the normalisation step, which may turn out to be very infrequent. These are things I just don't know.

I will try and upload an executable and log files later.

In a Doom-type game you could just reset the vectors each time the player view more or less aligns with one of the base vectors. Who would notice, when turning around, that one turn had a slightly smaller or larger angle to make the camera vector hit a right angle?

By asking the right questions and providing code. How can one answer the question whether a given precision is good enough without knowing what good enough means? If the peak error (l^infty) of the rotated vector is to be estimated, a hard error bound is given by the l^1 row norm of the errors of the rotation matrix.



It is the obvious approach, isn't it? I do not see any issues, I just wonder why this rather obvious approach is not used.



No, I just wonder why one picks quaternions when rotation matrices do the job. Even more so as quaternions (aka SU(2)) are not fully equivalent to rotation matrices (SO(3)), but only up to a sign the former can pick up, see https://en.wikipedia.org/wiki/3D_rot...SO(3)_and_SU(2)



Is an international standard good enough for you? The described algorithm is an integer implementation of it - actually, it is pretty much known and not novel. It is also the binary version of the manual "pen and paper" digit-by-digit square root.



If this makes you feel any better, mathieeedoubbas uses the same algorithm, for the IEEE mantissa of double precision numbers.

Look, once again, without knowing the precision you want to reach in the application (as in, "how precise do the vectors you attempt to rotate have to be over how many iterations"), nobody would be able to answer.

I don't think I have any opportunities like that, unfortunately. Possibly when the player changes camera angle or performs an action like shooting which would disguise any adjustments that the normalisation would cause, but I think I need to be more regular than that.

But, they don't though, do they? It sounds like you're in a better position than me to explain why.

It might surprise you, but I didn't actually come here to discuss square root algorithms.

I came here looking for advice from people with experience.

Perhaps someone with experience could say, "well, game engine x renormalises its quaternions when they drift more than y% away from unit length". That would be a worthwhile contribution to the conversation.

Or they could explain how to calculate the limitations of the particular fixed point size used in order to work out when it's not even worth trying to normalise. That would be helpful too.

We could even discuss the cost of normalisation to try and work out how often we can afford to normalise.

But to refuse to engage until your question is answered is unhelpful, particularly as it's a question that, a) can only be answered through experience, b) you have provided no help or guidance is how to answer, and c) may not even be relevant once the limitations of the fixed point maths are taken into account.

I've attached the latest executable and the log output it generates.

My main loop currently does this:

The output is something like this:

The last column is the square of the length of the quaternion without the right shift that would normally be done after multiplying fixed point numbers.

You can see it drifts away from 1 (<<28).

Edit: Very slowly, 0x0FFF7544 = 0.999868, and if I run it for longer, I get 0x0FE88AC1 = 0.994273.

Is the following too slow?

I'll use this 'rounding' when more precision is needed in my code (using partial results if possible a little better..).

Thanks, but I think you're answering the wrong question. There's nothing wrong with the 4.28 numbers, it's the rest of them that need to be corrected.

The square root of the 4.28 number is the length of the quaternion (think vector) described by the 2.14 numbers. If there were no rounding errors then the 4.28 number would always be 1<<28, but because we're always moving the quaternion (by multiplying it with another one representing where the control stick is pointed), the rounding errors build up, leading to a quaternion that is slightly off unit length. I need to divide the four 2:14 numbers by the square root of the 4.28 number (or multiply by the inverse square root, but I don't think that's possible in fixed point).

I need to work out, first, when is the division actually possible with the fixed point maths, and second, when do I want to do it? I see these two questions as interlinked, as the answer to the second may be forced by the first.

I doubt there is an answer like this. The answer is more likely: If your screen has a resolution of W, and you multiply a vector with a matrix to obtain a position on the screen, what is the resolution to ensure that the placement of such a vector is pixel-precise after N multiplications. That is an answer I can give, but for that, the screen dimension needs to be known, and the number of multiplications needs to be known.


If you have a f-fractional bit matrix, then clearly, the resulting relative error of a multiplication with a 4x4 matrix can not be higher than 4*2^-f. If you have a screen resolution of W pixel, then you get a maxmal absolute error of W*4*2^-f. Approximately, after N iterations, you get an error of N*W*4*2^-f pixels.


This is not really "high mathematics", just some basic estimates.



Concerning quaternions: After thinking about it, I believe you must be using homogeneous coordinates. In such a case, the sign collected from the SU(2) rotations (instead of SO(3) = rotation matrices) does not matter. But without seeing your formulae, it is hard to tell. I'm also not clear why this should be better than the typical 4x4 matrices you use for coordinate computation in projective space (aka "homogeneous coordinates).

That's swell, but it would be an answer to a question that wasn't asked.

Except we haven't moved from quaternion to matrix yet. And given that the quaternion is the only point where we can attempt to rectify the rounding errors, and the 14 bits of precision and effective [-2, +2) range we have will certainly constrain us, can we calculate the earliest we can normalise, knowing that we need to divide by the norm (or multiply by it's inverse) and that involves a square root?

Nope. We are concerned only with rotation here. Translation and the perspective transform occur separately (and must do, as they rely on a different level of precision in their fixed point maths), so it's all normal, not homogeneous coordinates.

Just for understanding, when I use fixed point I do with 8:8. Which is for movement of BOBs or Sprites on a 2d screen.

How can you just use 2:14? I mean, two bits for the coordinates of the whole number seems a bit small, isn't it?
How does this work?