Bug in dQuaternion::CalcAverageOmega

Report any bugs here and we'll post fixes

Moderators: Sascha Willems, Thomas

Bug in dQuaternion::CalcAverageOmega

Postby JoeJ » Thu Jan 03, 2013 4:41 am

I continued helping PJani via PM, and we found the reason why is powering did not work, while mine did:

dFloat omegaMag = dFloat(2.0f) * dAtan2 (dirMag, dq.m_q0) / dt;

should be:

dFloat omegaMag = dFloat(2.0f) * acos (dq.m_q0) / dt;

I may be wrong but i don't understand the use of atan2 here.
It may be the reason this is commented out in dgQuat, seems it's only used by Kinematic Joint.

EDIT: fixed the suggested fix
User avatar
JoeJ
 
Posts: 1453
Joined: Tue Dec 21, 2010 6:18 pm

Re: Bug in dQuaternion::CalcAverageOmega

Postby JoeJ » Thu Jan 03, 2013 6:36 am

Hmm, after some testing there's no difference, so sorry - no bug :oops:
I expected different signs in some cases but it seems not so for the kinematic joint usecase.
Maybe PJani had different signs in the input quats or something like that.
User avatar
JoeJ
 
Posts: 1453
Joined: Tue Dec 21, 2010 6:18 pm

Re: Bug in dQuaternion::CalcAverageOmega

Postby Julio Jerez » Thu Jan 03, 2013 8:51 am

the atan2 comes form this

quat = [cos (angle/2), dir.Scale (sin (angle/2)]

basically a quat is a pair of a [cos sin] values
it jus happen that teh sin is a vector compenemet

if you simple use acos to get teh angle you will always get the magniti of the angle but not the sign.
Julio Jerez
Moderator
Moderator
 
Posts: 12249
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Bug in dQuaternion::CalcAverageOmega

Postby PJani » Thu Jan 03, 2013 10:23 am

I think i know what was happening...

The only problem i see is how the shortest path is being found...
Code: Select all
//this one is from kinematic joint
      Ogre::Quaternion q0_tmp (q0);
      if (q0_tmp.Dot(q1) < 0.0f) {
         q0_tmp.w *= -1;
         q0_tmp.x *= -1;
         q0_tmp.y *= -1;
         q0_tmp.z *= -1;
      }
      Ogre::Quaternion dq (q0_tmp.Inverse() * q1);

Code: Select all
//this is JoeJ version
   static inline Ogre::Quaternion quaternionAtoB(const Ogre::Quaternion& qA, const Ogre::Quaternion& qB)
   {
      
      Ogre::Quaternion q;
      if (qA.Dot(qB) < 0.0f)
      {
         q.x = qA.x; q.y = qA.y; q.z = qA.z;
         q.w = -qA.w;
      }
      else
      {
         q.x = -qA.x; q.y = -qA.y; q.z = -qA.z;
         q.w = qA.w;
      }

      return qB * q;

   }
| i7-5930k@4.2Ghz, EVGA 980Ti FTW, 32GB RAM@3000 |
| Dell XPS 13 9370, i7-8550U, 16GB RAM |
| Ogre 1.7.4 | VC++ 9 | custom OgreNewt, Newton 300 |
| C/C++, C# |
User avatar
PJani
 
Posts: 448
Joined: Mon Feb 02, 2009 7:18 pm
Location: Slovenia

Re: Bug in dQuaternion::CalcAverageOmega

Postby Julio Jerez » Thu Jan 03, 2013 11:53 am

the firs one is tehg correct way.

the secund is using the conjugate of a matrix, wich mean it is mutiplying by teh inverse of the target matrix each time angle between teh two quat is laraler that 180 degrees. and that makes no sence.


what you want t o avoid is the effect that because

matrix (Q) = matrix (-Q)
each time is time you get a matrix for quat Q, you will so get the same matrix from -Q

edit:
actuallly after looking at Joe version I thonk his is correct too.
he try to get teh sortes patt and teh conjugate in teh same step. I am no sure is that works every time but it may.
I prefer to do simnpler first check that the shoprtest path is the alone the positive projection of the target quat over the origin quat
the mutiply the current quat by the conjugatye of the target quat, is cleaner and less prone to error.
Julio Jerez
Moderator
Moderator
 
Posts: 12249
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Bug in dQuaternion::CalcAverageOmega

Postby JoeJ » Thu Jan 03, 2013 12:12 pm

the atan2 comes form this

quat = [cos (angle/2), dir.Scale (sin (angle/2)]

basically a quat is a pair of a [cos sin] values
it jus happen that teh sin is a vector compenemet

if you simple use acos to get teh angle you will always get the magniti of the angle but not the sign.



You're right, but after you do this:

dQuaternion q0 (*this);
if (q0.DotProduct (q1) < 0.0f) {
q0.Scale(-1.0f);
}
dQuaternion dq (q0.Inverse() * q1);

dq.m_q0 MUST be positive, or am i wrong?
So atan2 is not necessary but it can't be a bug. I'm still confused how it can be a difference for PJani - his code looks right.
Maybe ogre gives him corrupt quats.

P.S.: Me fusion of two things is just an optimisation - i can't see a difference between our stuff until someone finds how dq.m_q0 cen become negative.
User avatar
JoeJ
 
Posts: 1453
Joined: Tue Dec 21, 2010 6:18 pm

Re: Bug in dQuaternion::CalcAverageOmega

Postby JoeJ » Thu Jan 03, 2013 1:30 pm

I wrore some testcode to proof.
The breakpoint is never reached, but it is if the quats are unnormalized.
PJani, try to normalize them or grab them from newton instead of ogre, as you planned to do already.
Never trust graphics engine when it's about physics :)

Code: Select all
dQuaternion n[8] =
      {
         dQuaternion( 1, 1, 1, 1),
         dQuaternion( 1,-1, 1, 1),
         dQuaternion( 1,-1,-1, 1),
         dQuaternion( 1,-1,-1,-1),
         dQuaternion(-1, 1, 1, 1),
         dQuaternion(-1,-1, 1, 1),
         dQuaternion(-1,-1,-1, 1),
         dQuaternion(-1,-1,-1,-1)
      };
      for (int i=0; i<8; i++) for (int j=0; j<8; j++)
      {
         dQuaternion q0 = n[i]; q0.Normalize();
         dQuaternion q1 = n[j]; q1.Normalize();

         if (q0.DotProduct (q1) < 0.0f) {
            q0.Scale(-1.0f);
         }
         dQuaternion dq (q0.Inverse() * q1);

         if (dq.m_q0 < 0.0)
         {
            int breakpoint = 1;
         }
      }
User avatar
JoeJ
 
Posts: 1453
Joined: Tue Dec 21, 2010 6:18 pm

Re: Bug in dQuaternion::CalcAverageOmega

Postby Julio Jerez » Thu Jan 03, 2013 2:58 pm

JoeJ wrote:You're right, but after you do this:

dQuaternion q0 (*this);
if (q0.DotProduct (q1) < 0.0f) {
q0.Scale(-1.0f);
}
dQuaternion dq (q0.Inverse() * q1);

dq.m_q0 MUST be positive, or am i wrong?
So atan2 is not necessary but it can't be a bug.

where are you getting that q0 must be positive? this is not so.
a quetrnion is a unit circle in four dimension. It is represented in a cartesian space as
Code: Select all
q0^2  + q1^2  + q2^2  + q3^2   = 1.0

any of the q's have any especial meaning, there is not reason to think that q0 have to have any special value, or must be a positive value.
the expression
q0 = [cos(W/s)
[q1, q2, q3] = unit vector.Scale (sin (w/2) )

is just one convenient way to initiliaze a quat, but you can also use any permuation, for example I could use if I want
q2 = cos(w/2)
[q0, q2, q3] = unit vector.Scale (sin (w/2) )

and it will be just as valid.
It is difficut in four dimensions to see why the better way multiply two quat is by testing the dot product, but you can make a projection to a lower dimension and apply the idea.
for example if you take a unit sphere
x^2 + y ^ 2 + z ^2 = 1

and you make teh z value zero, what you are doint is a prjection oven teh x-y plane, you get
x^2 + y ^ 2 = 1

that is a sphere too, in the sence that every point on the circle of that expression satisfy all of the property of every point on the surafce of the more general sphere. You can do the same with a quat, and force two of the cordinate to be zero
q0 = cos(W/s)
[q1, q2, q3] = unit vector.Scale (sin (w/2) )

let unit vector be [1, 0, 0] and you get

q0 = cos(W/s)
q1 = sin(W/s)
q2 = 0
q3 = 0

now the equation will be

q0^2 + q1 ^2 = 1

which is a circle, we can apply the principle that I applied to the quaternion delta rotation to this special quat by making an image
quat.Rotationpng.png
quat.Rotationpng.png (10.36 KiB) Viewed 9902 times


in the image say the red line is the source quat r0 which we kknow have q2 and q3 = 0
we want to rotate to quaternion r1 represented by the cyan line.
we want to calculate the delta rotation matrix that will rotate r0 into r1.

there are two equally valid possible solutions to the problem, because we already say that the rotation matrix that maps to quat Q, also maps to -Q
in this case r0 is given, therfore we have to chose between quat r1 and quat -r1
you can resolve this by any other geomtrical or analical method and you will always get r1 and -r1 as two equal valid solution, because this is a secund order problem equation.
To resolve the ambiguity a we need a third peice of informations. in this case we always want the solution tha yield the shortest arch distance from r0 to r1,
that information is determined by projecting r1 onto r0 and -r1 onto r0, and as you can see in the image r1 is the one yielding the shortest distance
so you alwys do

dot =dotproduct (r0, r1)

if (dot < 0) then is must be that r1 is the reflection of what I wanted, so you reflect r1
and then you do teh rest of the operation.
when all the other dimension are not zero, the problem is identical because the two quat determine a hiperplane that devide the hypersphere into two hyper circles.
I let you do the sustitution of yor method and see if it yeild the same result.
Bascially you place r1, and teh you place the conjugate of r1, see i fteh product ot q1 and it scinjugate yeled teh result you want in all cases, I am alomost 100% certain it will not. But I did not did the test, It may be but I am not going to do it myselft.





to resolve the ambigutu
Julio Jerez
Moderator
Moderator
 
Posts: 12249
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Bug in dQuaternion::CalcAverageOmega

Postby JoeJ » Thu Jan 03, 2013 3:33 pm

hmm, i must look closer to that spending more time, i feel it'll help me to understand how quats work...

On the thing where PJani postet two versions of our code - don't think about it, it's not worth the time.
Just look at the instructions and ignore math background. There's an if where you might flip all signs, followed by an inversion where you (again) flip 3 signs,
you can optimise that to either flip 3 or one sign(s), which is what i do but it's the same.

Julio Jerez wrote:wher are you getiing that q0 must be positive? this is no so.


I mean that only the 'w component', or 'dq.m_q0', or 'angular component' - no matter how we call it, is positive in this case which is special because:
1. q0 and *this are assumed to be rotation quaternions
2. There dot product is => 0
3. dq = q0.Inverse() * q1

I can't explain why this is, but my test code seems to proof that this is right.
I'm pretty sure you can get a mathematical proof, if you're interested.
It's just that guys like my observe such stuff over time and become a feeling for it (but never be sure), because math beckground is weak.
User avatar
JoeJ
 
Posts: 1453
Joined: Tue Dec 21, 2010 6:18 pm

Re: Bug in dQuaternion::CalcAverageOmega

Postby Julio Jerez » Thu Jan 03, 2013 4:16 pm

q0 does not have to be positive. and q0 is not any "angular component" of a quat
that is the equivalent to say "in a circle to say that x is the angular component", just because is polar cordenate the x componet is given by
x = cos (angle)

in a quaternion all value of q from 0 to 3 have equal weight, neither is more or less important that the other and neither is special
Julio Jerez
Moderator
Moderator
 
Posts: 12249
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Bug in dQuaternion::CalcAverageOmega

Postby JoeJ » Thu Jan 03, 2013 5:07 pm

Julio Jerez wrote:q0 is not any "angular component" of a quat

This is very hard for me to get (note that i'm not ready to leave the 3D-rotation use case for now).

Is this also true for rotation use case, where we build the quat from angle and axis?
And is the quaternion multiplication special for the rotation use-case, where q0 is build only from subtractions, while q1-q3 have additions too?
I always assumed the 'standart' multiplaction we use is constructed only for the rotation usage (i don't know other practical quat usecases).

I actually understand some of the typical quat usage, like how and why the dot product works, axis/angle conversion and why Q maps to -Q,
but i absolutely do not understand the muliplication.
I'll try to get a better understanding by analysing how vector rotation works, that might be the key for me.
I actually failed to create a 'triple-value-quat' for 2D rotations :roll:
User avatar
JoeJ
 
Posts: 1453
Joined: Tue Dec 21, 2010 6:18 pm

Re: Bug in dQuaternion::CalcAverageOmega

Postby Julio Jerez » Thu Jan 03, 2013 6:46 pm

JoeJ wrote:Is this also true for rotation use case, where we build the quat from angle and axis?

yes,
when you buidl a quat form a( cos(w/2), unitVector.Scale (sin(w/2)]

q0 = cos(w/2)

you get q0 as a positive valube because we never pass angle that are larger than pi, this is because we knwo that teh quad can no hold a roation larger that 180 degrees, what happens if you pass a larger 180 degree anagle is that you get and alias rotataion, meaning the rotation wrapt around

in this scase q0 = be negative, but if you multiply all q's by -1 then you will ge a q0 tha is positive but teh axis is flipped, but tso are all valid rotation.
is is that by defintion we linmit teh angle to be bewteen 0 and 180.

Regard to quatenion mutiplication, I have to say I do no undertand it either, but shis is not because you an I are stupid, it is because quation are element of four dimensions, the is not way that a human mind can visualize geometrically how object works in higher dimension.

The only way to underntand it is my the algebra of quartion aritmetic.

The trditiona theory of Quaternion rotation was writen and explened bu a man named Ken Shoemake,
hwoever I find confusing and no intutive al all, he try to make intuitive and in my opbnion he fail at it.
I have a book on Clasical mechanic that was part of my reading material in college,
that have a much more elegane explanation of the algebra of quaterion rotation, I do not remember the name now, I will post it later.

But bascially he formulate the quaternion as a complex number with only two componet in the frequency domain.
(silimilarto the euler theorem, and the furier thansforms) them he apply the laws of derivation, mutiplication, series explansion and so on,
and from there you can see how everything fall into place with an austanding elegancy.
Basically the quaternion algebra is correct not because it is intuitive, but because the mathematic says so.
This is the case for every almost everything, once you deal with objects in more than three dimensions.
Julio Jerez
Moderator
Moderator
 
Posts: 12249
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Bug in dQuaternion::CalcAverageOmega

Postby JoeJ » Mon Jan 07, 2013 5:56 pm

I'm through with it and understand now. That took me some time, but now it's easy and everybody with knowledge of 3D vectors and trigonometry can understand quaternions.
I'll show how vector rotation and quaternion multiplication works. More than that, i'll show that it's possible to reinvent quaternions without using complex numbers or 4D space.
I don't wanna fight against anything with that, it's just that i think this might be useful for anyone interested in quaternions (except those who already know, because there's surely nothing new).


1. Quaternions - definition, usecase and visualisation

quatvis.png
quatvis.png (12.62 KiB) Viewed 9845 times

I'll discuss only rotation quaternions and define them as:
A 3D vector representing the axis of rotation times the sine of the rotation halfangle,
and a scalar, which is the cosine of the rotation halfangle.

We can treat it as a 'optimized axis / angle rotation' - optimized because encoding sin and cos has advantages, as we'll see...
As visualisation i choose a cone that shows the angle of rotation, and a vector to show the rotation axis and it's length.
User avatar
JoeJ
 
Posts: 1453
Joined: Tue Dec 21, 2010 6:18 pm

Re: Bug in dQuaternion::CalcAverageOmega

Postby JoeJ » Mon Jan 07, 2013 6:20 pm

2. How to rotate a vector using a quaternion

xform_vec.png
xform_vec.png (14.37 KiB) Viewed 9844 times


To transform a vector by a quaternion, we build two orthonormal directions aligned to the rotation plane
and calculate the difference vector to the rotated vector.

Code: Select all
vec &x = (sVec3&)q[0]; // map x to the quaternion axis
vec y = x.Cross(v); // v is the vector we want to rotate


Let's look at the length properties of the y vector we've got:
|y| = |v| * sin(angle(x,v)) * sin(halfangle)
note: |v| * sin(angle(x,v)) = Distance from v to rotation axis, we call that term 'axisDist' in the following equations.

Code: Select all
y *= 2.0;
vec z = x.Cross(y);


Z is the second vector, and it has a length that we can already use:
|z| = axisDist * 2 * sin(halfangle)^2
|z| = axisDist * 2 * (cos(angle+PI)/2 + 0.5)
|z| = axisDist * ((cos(angle+PI)) + 1)
note: We go from a point on a circle (not it's center), thus the offset of 1 in the cosine term.

Code: Select all
y *= q[3]; // multiply with the scalar part from the quaternion


|y| = axisDist * 2 * sin(halfangle) * cos(halfangle)
|y| = axisDist * sin(angle)
note: Because y is perpendicular to circle-point - circle-center, no offset in sine term.

Code: Select all
vec result = (v + y + z);


Add the vectors and we're done.
It's interesting to see that so much trigonometry is done using so few instructions.
This shows how useful it is to encode sine and cos in the quaternion.
There's not much to explain, just follow how the two orthonormal directions and their trigonometric length properties build a circle.
User avatar
JoeJ
 
Posts: 1453
Joined: Tue Dec 21, 2010 6:18 pm

Re: Bug in dQuaternion::CalcAverageOmega

Postby JoeJ » Mon Jan 07, 2013 6:54 pm

3. How to rotate a orientation with a quaternion ( = the hard to understand quaternion multiplication)

quatmult.png
quatmult.png (19.1 KiB) Viewed 9844 times


We do some mappings to access vector and scalar parts from our quats individually:
Code: Select all
sQuat q, r, qr; // we want qr = q * r

vec &qrV = (sVec3&)qr;
float &qrS = qr[3];
vec &qV = (sVec3&)q;
float &qS = q[3];
vec &rV = (sVec3&)r;
float &rS = r[3];


First (ignoring the image) we'll look at an easy special case: Both rotations have their axis at the same line, so we can treat this as an 2D-problem.
We detect that easy case by taking axis cross product and check it's length:

Code: Select all
sVec3 qcv = qV.Cross(rV);
if (qcv.Length()<0.001)
{


It's clear that we just need to sum (or subtract) both angles,
because rotating first 20, then 5 degrees over the same axis is the same as rotating degrees 25 at once.
(If the second axis points in the opposite direction, we should end with 15 degrees.)

If we want to sum two angles, where sines and cosines are known,
we can calculate the cosine of the summed angle by:

cosSummedAngle = (cos1 * cos2) - (sin1 * sin2)

Thus the scalar part (cosine of final halfangle) is:

Code: Select all
   qrS = qS*rS - qV.Dot(rV);


Note: The dot product also includes the factor cos(angle(q.axis, r.axis)),
which in this case is either 1 or -1 and handles axis pointing in opposite directions nicely.

Further, we can calculate the sine of the summed angle by:

sinSummedAngle = (sin1 * cos2) + (cos1 * sin2)

We build the vector part (sine times axis) accordingly:

Code: Select all
   qrV = sVec3(qV*rS + rV*qS);
}
else
{


Now the general case, which we treat as a 3D Problem (now look at the image).

One dimension more means one more angle to look at:
The angle between the two rotations axis that we wanna combine.
We already have their sine and cosine in our dot and cross products,
So we'll see how to use them...

We'll use following example:
Put a cube in your hand and remember it's initial orientation.
Rotate it 90 degrees over one of its local axis, followed by another 90 degree rotation around another local axis.
To rotate it back to the initial state in one go, you need pick it at two opposite corner points and rotate 120 degrees.

How can we calculate this angle?

It turns out, the cosine factor from the dot product has another nice property:
In the example situation it zeroes the sine term,
which gives us exactly the angle we want (cos = 0,5 -> 60 degrees halfangle).
Thus the 3D case of our scalar part becomes:

cosSummedAngle = (cos1 * cos2) - (sin1 * sin2) * cos(angle(axis1, axis))

So no need to change anything about the scalar part.

Code: Select all
   qrS = qS*rS - qV.Dot(rV);


In the image you see three quaternions (red * green = blue).
I added a cube so you can see the blue axis goes through opponent cube corners to verify the example.
The orange line shows the axis result from our formula we already have.

Code: Select all
   qrV = sVec3(qV*rS + rV*qS);


We see it is at the right place on the plane formed by red and green axis,
and this is how we get the missing distance - we just add the cross produkt (yellow line):

Code: Select all
   qrV += qcv;
}


This time the sine factor in the cross product comes to our advance,
because it gives the cross product the right length and we are done.

At this point, it may seem still a bit magical,
but it is not more magical than saying: "cos(a) + sin(s) = 1"


To be true, the above code has some instructions more than the standart quaternion multiplication,
Maybe i'm in the mood to figure out what terms cancel each other out, but i reached my goal:
Understand, how quaternions work.

So please let me know what you think guys, but please don't say 'What ???' :)
User avatar
JoeJ
 
Posts: 1453
Joined: Tue Dec 21, 2010 6:18 pm


Return to Bugs and Fixes

Who is online

Users browsing this forum: No registered users and 4 guests