yade-dev team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

Hi,

you shouldn't be running wild through the code adding
quaternion.normalize() every second line in the code.

I've went through this stage already, and let me explain what is the
problem and how to proceed.


Quaternions:

Quaternions, are just like complex numbers, except that instead of
single imaginary "i", they have three imaginary numbers: "i", "j", "k".

We are using them to perform rotations in 3D space. Namely, we are
using them as a representation of an SO(3) rotation group:
http://en.wikipedia.org/wiki/SO(3). Fortunately multiplying the
quaternions is just like composition of two rotations.
(BTW: orthogonal rotation matrices 3x3 are also a representation of SO(3) ). 

Quaternions, are quite like complex numbers, and can be of ANY
magnitude (length) that you could imagine. You can calculate sin(q),
ln(q), e^q, solve polynomial equations with them, and other stuff.
Which we are not interested in.


The problem:

The problem is that any given quaternion belongs to the SE(3)
rotation group only, and only, if its length equals ONE.
If any quaternion suddenly doesn't have length=1, it does not belong
any more to SE(3).

Analytically multiplying quaternions with length=1 will always
produce a quaternion with length=1.

Numerically multiplying two quaternions requires 12 additions and 16
multiplications. There is a numerical error of order 1e-17 which can
build up, making the resultant quaternion<double> to be not unitary
any more.


The way to solve this:

Analytically + theoretically if a quaternion's length is not 1, then
this is for us a critical error since it is NOT a rotation any more.

Numerically + practically there is a numerical error introduced by
the computer. It is of magnitude std::numeric_limit<TYPE>::epsilon()

http://www.cplusplus.com/reference/std/limits/numeric_limits/
http://msdn.microsoft.com/en-us/library/6x7575x3%28VS.71%29.aspx

(where used TYPE is usually float, double or long double, for double
IIRC it was 1e-17).

We have to cope with this error, and due to this error it is
impossible to have a quaternion's length always equal to 1.


How this epsilon() propagates through calculations (12 additions, 16
multiplications) depends exactly on the way in which it was
implemented. The worst case would be: a single epsilon() becomes
12*16*epsilon() larger after single multiplication. I am sure that
this is never the case. But you get the picture. We can call that
12*16 a MAX_ERROR and use 200 instead of 12*16=192.


The problem again:

So, we want to use UNITARY quaternions, but the library (wm3 and
eigen) provide quaternions ONLY. Not the unitary quaternions.

Our problem is inherent to the mathematics of SO(3) and computer's
numerical errors. The end user must never worry about it.


The solution:

We need a library that provides UNITARY quaternions.

If we don't have such library, then we need to normalize quaternions
every second line in the code. And suddenly YADE's code base is
bigger by 5000 lines.

The only correct way to cope with this, is to extend existing
quaternion library with a simple function that checks if quaternion
is UNITARY.

if quaternion's length differs from 1, then:

- if it differs less that MAX_ERROR*epsilon(). Then the error is not
  critical. We have to normalize that quaternion (divide it by its
  length) and silently move on.

- if the difference from 1 is bigger than MAX_ERROR*epsilon()
  then YADE has to report a critical error, throw exception and instantly
  terminate all calculations. Since whatever calculations there were,
  they ARE wrong.

This simple function has to be called at the end of each non-const
member function in quaternion class.



Let me repeat the solution:

1) Do NOT run wild through the code normalizing quaternions everywhere.

2) Do check if quaternion is unitary at the end of each non-const
member function of quaternion class.

3) quaternions are in eigen lib and we can not modify it? Then we
have to either:
   - send them a patch, witch enables the library to provide unitary
     quaternions, or
   - derive from their class and extend it.


I hope that it is clear now.

This error is NOT of type: OMG all our calculations were wrong!


The error build-up exists, usually is negligible. It shows itself only if:

- same quaternion is used a lot in a single loop. lading to lots
  of multiplications in single iteration. This is not the case in
  DEM. Such multiplication in DEM with spheres is done just once per
  iteration per quaternion.

- a quaternion is mistakenly initialized to be non-unitary. Which CAN
  be the case quite often, if you don't read what you are writing.


I stumbled on this error with lattice model, because rotating a
single rod required about 10 to 20 quaternion multiplications per
iteration per quaternion. And I thought that it should be fixed in miniwm3, but:
- I was too lazy to modify miniwm3
- I was scared of computation overhead that this is going to introduce
- I saw that in fact this problem manifests itself rarely.


In fact that function instead of doing unitarity check on every call,
could count the number of times that it was called, and perform
normalization only every 10th call. That solution doesn't seem to be
"perfect", but I tend to like it.

Maybe it is better to do all those checks ALWAYS (each time) in DEBUG mode.

And do them every 10th time in non debug mode. You cannot skip them
totally in non-debug mode, because in non-debug mode the numerical
buildup error occurs too!


well, I hope that now you understand what's going on here, and aren't
scared anymore
-- 
Janek Kozicki                               http://janek.kozicki.pl/