Hey guys,

I'm in the process of working out my math library, and so I was testing out a function that returns the angle between two quaternions..

{ return 2.0f * acos( Abs( Dot4( qa, qb ) ) ); }

.. but for some reason, I'm either getting a lot of floating point error in the result, or I'm not checking for a situation that I should be. While testing a quaternion (which was generated by a random axis+angle rotation and appears to be very close to normalized)..

{ x=0.0172970667 y=-0.0245058369 z=0.0205858145, w=-0.999337912 }

.. with itself, I'm getting a result angle of 0.00138106791 (or almost 0.1 degrees)..

I'm just wondering if this is acceptable error when working with float variables? And is there anything I can do to improve this issue other than switching to double type or something else as drastic?

edit note: After testing some more, the highest "error angle" I've been able to generate (through random axis-angles) is 0.001953125. And that was getting the angle (from itself) of a quaternion generated by the axis @ angle: { -0.833756,0.551120,-0.033417 @ 2.960138559341 } (quaternion result: { -0.830327,0.548853,-0.033279,0.090603 } )

Thank you

well you could check long double for it

just make something lie this

rewrite dot4 so it will get quaterions (like old did but i will operate on long doubles and will return long double)

ofc function returns long double

then return 2.0 * acosl( fabsl( DotDouble4( qa, qb ) ) ); }

this could use the best precision for floats

anyway try to rotate something by in ex 90 degrees and try to find that angle.

After further testing, it seems I was wrong when I said my quaternions were very close to normalized.

Apparently, my quaternion axis-angle rotation was not generating perfectly normalized results. Even though the lengths were around 0.9999+, it was enough to generate the angle differences above.

In the same test loops I used above, I normalized the quaternion after creating it from an axis+angle, and the angle between itself was exactly 0.

So what I've learned from this is that I have to normalize every quaternion I create from an axis+angle, even if the input axis was normalized..?

Sorry, I should have tested further before asking.

Thanks again.

yes many have such problems, they forget to normalize them :)

Numerical stability is a pain. Sometimes you have to really think through your algorithms and rewrite them to be somewhat less efficient but more stable. There's a _lot_ of work required to really understand numeric stability with floats. I don't myself pretend to really understand it. This is one of the (many) reasons why you should just use an existing well-tested math library written by math and CS experts rather than trying to create your own.

I just wanted to point a couple things out. I can't be 100% sure, but I believe the original quaternions were already normalized. After a lot of testing and examining, it looks as though the normalization simply traded error on one side for error on the other side, which resulted in a SEEMINGLY better angle result. I noticed this because I can literally renormalize a quaternion over and over again, and get different results each time.

I wanted to point this out because I am now no longer positive that a quaternion constructed from a normalized axis + angle will need normalized after being constructed. From what I've seen (from a lot of extreme random axis+angle generation), it is pretty much already normalized within a decent error tolerance.

edit: Here are the results:

q1: error-angle[0.00119] construct-length[0.999999880791], normalize-length[1.000000119209] quaternion[ -0.388950,-0.416175,-0.421772,0.705426 ]

q2: error-angle[0.00138] construct-length[0.999999880791], normalize-length[1.000000119209] quaternion[ -0.865744,0.019881,0.467512,0.177554 ]

q3: error-angle[0.00154] construct-length[0.999999880791], normalize-length[1.000000119209] quaternion[ -0.777527,-0.501608,-0.192569,0.326739 ]

(notice the error is always exactly 0.000000119209 distance from 1)

Numerical stability is a pain. Sometimes you have to really think through your algorithms and rewrite them to be somewhat less efficient but more stable. There's a _lot_ of work required to really understand numeric stability with floats. I don't myself pretend to really understand it. This is one of the (many) reasons why you should just use an existing well-tested math library written by math and CS experts rather than trying to create your own.

I agree. But from my limited experience, it seems next to impossible to create a 3D sim (regardless of a math library) without having to deal with these errors in some form or another anyway. Even something as simple as 1+1=2 can turn your game world upside-down if its not interpreted as 1+1=almost definitely 2.

Essentially, my engine relies on a very select subset of vital functions, so I really don't have a lot of work cut out for me. But getting this stumped on one function isn't really getting that point across, is it?

Thanks

I just wanted to point a couple things out. I can't be 100% sure, but I believe the original quaternions were already normalized. After a lot of testing and examining, it looks as though the normalization simply traded error on one side for error on the other side, which resulted in a SEEMINGLY better angle result. I noticed this because I can literally renormalize a quaternion over and over again, and get different results each time.

I think the quaternion you posted originally was not properly normalized. Do you care to post a particular quaternion for which the renormalization isn't stable?

Even something as simple as 1+1=2 can turn your game world upside-down if its not interpreted as 1+1=almost definitely 2.

That seems to indicate you don't understand floating-point numbers very much: 1 + 1 is most definitely exactly 2.

I just wanted to point a couple things out. I can't be 100% sure, but I believe the original quaternions were already normalized. After a lot of testing and examining, it looks as though the normalization simply traded error on one side for error on the other side, which resulted in a SEEMINGLY better angle result. I noticed this because I can literally renormalize a quaternion over and over again, and get different results each time.

I think the quaternion you posted originally was not properly normalized. Do you care to post a particular quaternion for which the renormalization isn't stable?

Well, it was as normalized as my normalize function will seem to allow it to be. And yeah, sure. Let me type something up..

Each line is of the same quaternion, renormalized again:

0.768565058708, 0.393124073744, 0.112918674946, 0.491945981979 (length: 1.000000119209)
0.768564939499, 0.393124014139, 0.112918660045, 0.491945922375 (length: 0.999999940395)
0.768565058708, 0.393124073744, 0.112918674946, 0.491945981979 (length: 0.999999940395)

another example:

-0.384467869997, -0.321283340454, 0.054897289723, 0.863682806492 (length: 1.000000119209)
-0.384467810392, -0.321283310652, 0.054897282273, 0.863682687283 (length: 0.999999940395)
-0.384467869997, -0.321283340454, 0.054897289723, 0.863682806492 (length: 0.999999940395)

one more:

0.335628807545, 0.205629050732, 0.067485921085, 0.916796505451 (length: 1.000000119209)
0.335628777742, 0.205629020929, 0.067485913634, 0.916796386242 (length: 0.999999940395)
0.335628807545, 0.205629050732, 0.067485921085, 0.916796505451 (length: 0.999999940395)

The error in length is always exactly 0.000000119209 away from 1. On the third attempt, I get the quaternion from the first with the length from the second. I assume this is because of digits beyond what I've printed out. I didn't try any further than 3 normalizes, but I assume it will jump back and forth.

Even something as simple as 1+1=2 can turn your game world upside-down if its not interpreted as 1+1=almost definitely 2.

That seems to indicate you don't understand floating-point numbers very much: 1 + 1 is most definitely exactly 2.

Sorry, I meant for that to be an exaggeration. You are correct, nonetheless.

I took your examples and tried to normalize them myself:

0.76856505870819091797, 0.39312407374382019043, 0.11291867494583129883, 0.49194598197937011719 (length: 1.00000011920928955078)
0.76856493949890136719, 0.39312401413917541504, 0.11291866004467010498, 0.49194592237472534180 (length: 1.00000000000000000000)
0.76856493949890136719, 0.39312401413917541504, 0.11291866004467010498, 0.49194592237472534180 (length: 1.00000000000000000000)

-0.38446786999702453613, -0.32128334045410156250, 0.05489728972315788269, 0.86368280649185180664 (length: 1.00000011920928955078)
-0.38446781039237976074, -0.32128331065177917480, 0.05489728227257728577, 0.86368268728256225586 (length: 1.00000000000000000000)
-0.38446781039237976074, -0.32128331065177917480, 0.05489728227257728577, 0.86368268728256225586 (length: 1.00000000000000000000)

0.33562880754470825195, 0.20562905073165893555, 0.06748592108488082886, 0.91679650545120239258 (length: 1.00000000000000000000)
0.33562880754470825195, 0.20562905073165893555, 0.06748592108488082886, 0.91679650545120239258 (length: 1.00000000000000000000)
0.33562880754470825195, 0.20562905073165893555, 0.06748592108488082886, 0.91679650545120239258 (length: 1.00000000000000000000)

Here's my code:

#include <cmath>
#include <cstdio>
#include <boost/math/quaternion.hpp>

typedef float Real;
typedef boost::math::quaternion<Real> Q;

void show_two_normalizations(Q q) {
  for (int i = 0; i < 3; ++i) {
    std::printf("%.20f, %.20f, %.20f, %.20f (length: %.20f)\n",
                q.R_component_1(),
                q.R_component_2(),
                q.R_component_3(),
                q.R_component_4(),
                abs(q));
    Real inv_length = 1 / abs(q);
    q *= inv_length;
  }
  std::puts("");
}

int main() {
  show_two_normalizations(Q(0.768565058708, 0.393124073744, 0.112918674946, 0.491945981979));
  show_two_normalizations(Q(-0.384467869997, -0.321283340454, 0.054897289723, 0.863682806492));
  show_two_normalizations(Q(0.335628807545, 0.205629050732, 0.067485921085, 0.916796505451));
}

Sorry, I'm not sure I understand what you're trying to show me?

To be clear, the floating point errors didn't concern me as much as wanting to make sure the angle function is doing its job. I mean if there's something obvious I'm doing to produce the errors, I will gladly work to correct it, but AFAIK, its inherent to basic math.

Thanks again :)