Slerp should take the shortest way between two quaternions

[mjoubert]

Hi all,

I am working on attitude interpolation for real-time scriptedOrbit and I had an issue with SLERP : it sometimes doesn't take the shortest path between two attitudes. I fixed that for my script.

I checked what is done in Celestia quaternion library for SLERP, and I realize that this issue is not solved yet.

For testing, use this SampledOrientation :

Code: Select all

2400010 -0.51 -0.51 -0.51 -0.51
2400020 0.49 0.51 0.51 0.51

You should see that your spacecraft is making a impressive looping! (where it should rotate for few degrees only)

In file quaternion.h, in slerp function :

Code: Select all

line 538:    T c = dot(q0, q1);

This statement gives c, which is cos(half Theta), where Theta is the angle between two interpolated attitudes.
We'd like Theta to be confined between -180 and 180.
So we have to add this protection:

Code: Select all

   if( c < 0 )
   {
      q0 = -q0 ;
      c  = -c  ;
   }

It's working well for me.
So I ask a Celestia developper to check this modification, and commit it if it's correct.

Mathieu (spacebel)

[rthorvald]

You should see that your spacecraft is making a impressive looping! (where it should rotate for few degrees only)

FYI, this behaviour is evident with the old Orientation declaration, too; if you multiply an object and give it different (progressive) orientation values, it will jump weirdly around itself in certain instances instead of rotating in small, even steps.

- rthorvald

[Rassilon]

Its probably due to a precision issue using floats in the mathematics... My guess would be that your calculations are correct but the core of celestia is based on float instead of double thus giving that jump in your steps...

[mjoubert]

Rassilon,

Thank you for your answer. My explaination was confusious. Let's do it again.

In slerp, the first step is the computation of the angle between attitudes. The dot product (a.w*b.w + a.x*b.x + ...) gives cosinus half this angle. Exemple (think in 2 dimensions) :

We want to go from Pi/3 to -Pi/3. Quaternion 1 is 0.85/0.5/0/0 (0.85 is sqrt(3)/2). Quaternion 2 is 0.85/-0.5/0/0. So, dot product (c in source code) is 0.5. So theta is 2*acos(0.5) = 2Pi/3 which is what we wanted.

Same thing between 2Pi/3 and -2Pi/3. The shortest path is direct rotation (through Pi). Q1 = 0.5/0.85/0/0 and Q2 = 0.5/-0.85/0/0. So c is now negative, -0.5. If you take 2*acos(-0.5), it's 4Pi/3. This is the looping!

Since we want to take the sortest path, we need c to be positive. Replacing Q1 by -Q1 (they are the same attitude), solve this problem. Now, c is positive, and so gives the shortest path.

That's why I suggess to add

Code: Select all

if( c < 0 )
   {
      q0 = -q0 ;
      c  = -c  ;
   }

to force Theta to be always between -180 and 180.

You can try the same test with this exemple values :

Code: Select all

2400010 0.5  0.85 0 0
2400020 0.5 -0.85 0 0


Without source code negative c protection : looping
With protection : no looping.

Hope it helps to understand.


Mathieu (spacebel)

[chris]

Mathieu,

Thanks for spotting this bug. I thought that the code to use the additive inverse was already there, but apparently not. I'll check in the fix to CVS.

--Chris

[mjoubert]

Thanks Chris.

I'll send you the interpolated version of rtorientation when it will be finished (some convertions from CNES defintions left).

I'll probably do a topic about that.

Mathieu (spacebel)