Hi,

Generating an angleAxis object form a quaternion and getting the angle of it sometimes returns NaN. This result is obtained after relatively simple operations on quaternions representing near-zero rotations, without normalization (example below).
Does it means we have to normalize the quaternions all the time, before converting them to angleAxis?

Best regards.

Bruno Chareyre

For instance, lets take those quaternions :
q1=0.745315 -0.0308916 -0.665997 -1.14489e-10 (orientation of body 1)
q2=0.745315 -0.0308916 -0.665997 0 (initial orientation of body 1)
q3=0.746671 -0.166278 -0.644077 0 (initial orientation of body 2)
q4=0.746671 -0.166278 -0.644077 1.47527e-18 (orientation of body 2)

1. compute delta=(q1*q2.conjugate())*(q3*q4.conjugate()) (the relative rotation between two bodies, after very small movements),
2. define AngleAxisr aa(delta),
3. get aa.angle() -> NaN

This result seems to appear only if the rotations are very small, else the angle is correct. Here we have the rotations :
(q1*q2.conjugate())=1 -1.56132e-11 7.47311e-11 -1.7066e-10
(q3*q4.conjugate())=1 9.50188e-19 -2.45305e-19 0

and the relative rotation :
delta = 1 -1.56132e-11 7.47311e-11 -1.7066e-10

does that still occur if you normalize your input quaternions (q1 to q4) ??

anyway, though our rule is to assume the use gives normalized quaternions, here a small numerical error can indeed lead to NaN (because of the acos), so I've enabled auto-normalization for this function.

For the record, it is line 179 in src/Geometry/AngleAxis.h:

m_angle = Scalar(2)*std::acos(q.w()/internal::sqrt(n2+q.w()*q.w()));

Thank you for the reply and code update.
Yes, the same thing occurs if q1,...,q4 are normalized. Only normalizing delta fixes it.
Actually, we didn't want to normalize delta all the time for speed reasons, although it was the obvious way to avoid nan's.
A workaround we used was to check nan's on the result, and set angle=zero if NaN was found.
Isn't this fix almost as expensive as normalizing delta each time?

The problem sits on the last few decimal places, and has roughly the same magnitude as numerical precision epsilon 'e', you know, when 1.0+e==1.0 occurs. 1.0 + 2*e != 1.0, for us it is 1.0, but for acos it is not, and gives NaN.

I would say that this line:

m_angle = Scalar(2)*std::acos(q.w()/internal::sqrt(n2+q.w()*q.w()));

isn't perfect, beucase you are already calculating magnitude of quaternion and perform division, then why not divide other components too? It's only three divisions more. If dividing is too expensive, then make it multiplications. Then you would be normalizing whole quaternion and not just one component.

And... in fact if we are already spending CPU cycles to normalize it, then perhaps even better would be to normalize the original variable, not just its local copy in the function? I mean - pass the argument by reference and normalize it?

This breaks convention of having a `const` argument, so I don't know.

Bruno, when you say that only normalizing delta solves the problem - it is exactly what we are trying to do above in AxisAngle, right? Or you meant normalizing in some other place?

There is no need to apply this normalization to other components since the others are already normalized by there respective norm to get the axis. In other words, no need to normalize twice.

There is no way we modify the content of the input Quaternion!

However, we are considering to revert this change, and simply clamp the w component to [-1..1], and let the responsibility of the user to construct and pass normalized quaternion as we do everywhere else.

In this precise example the problem was that the input quaternions q1 to q4 were not normalized (maybe this is only because you printed them with limited precision), and then the numerical errors in delta were too high to be detected with double precision accuracy.

So if you write/load quaternions to /from files make sure you use full precision, and/or renormalize the input yourself.

there lies the problem: in our code we want the library to make sure that our quaternions are unit quaternions. The library finds it difficult, because of speed overhead and passes it back onto the user.

It is difficult. Maybe clamping indeed is a solution.

Numerically ideal solution, and a very slow one, would be to normalize quaterions after each multiplication. Not before, because AxisAngle uses output of multiplication, not its input. And multiplying normalized quaterion will introduce numerical error anyway.

But this ideal solution is too slow.

There could be a flag for instance, which allows to enable such behaviour in the library, and in case of numerical problems it could be turned on for debugging.

just some wild ideas. I know this problem is diffult, can we solve this is a way that is comfortable for everyone?

for reference, there our discussion about it:

https://lists.launchpad.net/yade-dev/msg05979.html

and also, there is an attachment with C++ code that is testing this behavior:

https://lists.launchpad.net/yade-dev/msg05988.html

@ggael
>In this precise example the problem was that the input quaternions q1 to q4 were not normalized (maybe this is only because you printed them with limited precision)

I confirm that in the pasted example, q1,...,q4 were not normalized, but I tried normalizing them today, and I still found NaN angles (no file input/output involved either).
I tend to like the clamp idea.

@Janek
"Normalizing" -> I meant q.normalize() before angleDelta(). Yes, it is mathematically similar, but (very) slightly slower than internal normalization since in the latest x^2+y^2+z^2 is already defined.
But ok, with q.normalize() we have the advantage of keeping a normalized quaternion...

The worst situation would be normalizing once before angleDelta, and then normalize again internaly...

ok, so please do the clamping.

consider adding an option that enables normalizing at the end of each multiplication. This could be a compile time flag for speed reasons, or a runtime flag for flexibility reasons.

Or both! Compile time flag which adds the capability of a runtime flag