search

czaloj / bullet

Automatically exported from code.google.com/p/bullet
0 stars 0 forks source link

btQuaternion angle and slerp do not deal with redundancy in quaternion representation correctly #140

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago 
There are two problems:
1. I identified in issue 98 comment 5
http://code.google.com/p/bullet/issues/detail?id=98#c5 
2. The same issue is a problem in the calculation of slerp.  

The general problem is that for any given quaternion q,  q == -q
The angle calculation needs to make sure that q and the other q are in the
same representation for otherwise the angle will be much too high.  

Likewise when interpolating between two quaternions with similar values but
the other representation the calculation will compute the weighted average
between a component and almost it's negative which will return close to
zero every time.  

If you apply either one of these patches it will get worse results.  (close
to zero if angle patch applied w/o slerp patch, or very large values if
only the slerp patch is applied)  

What steps will reproduce the problem?
1.
  btQuaternion q1 (M_PI/4,0,0);
  btQuaternion q2(-q1.x(), -q1.y(), -q1.z()-.0001, -q1.w()+.0001);
  q2.normalize();
  btQuaternion q3 = q2.slerp(q1, .5);

  printf("%f %f %f %f,%f %f %f %f\n", q2.x(), q2.y(), q2.z(), q2.w(),
q3.x(), q3.y(), q3.z(), q3.w());

What is the expected output? What do you see instead?
The second 4 values should be a valid quaternion very close(3-4 significant
figures) to the first four values. (Not transposed and partially negated)

(as in r1497)
-0.000000 -0.000000 -0.382804 -0.923830,0.000000 0.000000 -0.923855 0.382744

What version of the product are you using? On what operating system?
svn r1497 linux ubuntu fiesty.

Original issue reported on code.google.com by Tully.Foote on 11 Nov 2008 at 6:07

Attachments:

GoogleCodeExporter commented 9 years ago 

I suppose that is why we use the 'farthest' before using 'getAngle' in 
btTransformUtil:: calculateDiffAxisAngleQuaternion:

    btQuaternion orn1 = orn0.farthest(orn1a);
        btQuaternion dorn = orn1 * orn0.inverse();
        ///floating point inaccuracy can lead to w component > 1..., which 
breaks 
        dorn.normalize();
        angle = dorn.getAngle();

The two patches are conflicting. Can you provide a single patch that solves 
your 
issue, without breaking existing functionality in 
calculateDiffAxisAngleQuaternion?

Thanks,
Erwin

Original comment by erwin.coumans on 11 Nov 2008 at 4:13

GoogleCodeExporter commented 9 years ago 
Sorry I wasn't clear.  The two patches are alternatives.  The first patch 
applies
both the catch for the long direction and clamps the nan values from the patch 
for
ticket 98.  This is my preferred patch which fixes all my outstanding issues.  

The second patch does not apply the clamping from ticket 98. But does apply the 
logic
for the redundant representations to both angle and slerp.  I put this patch in 
if
the overhead of clamping is considered too much for common uses in bullet.  

These patches are inncompatible with previous patches in issue 98 because I
duplicated the line of code which I patched there.  

I will make sure that calculateDiffAxisAngleQuaternion continues to operate as
expected and get back to you.

Original comment by Tully.Foote on 11 Nov 2008 at 6:27

GoogleCodeExporter commented 9 years ago 
Attached is my full patch.  
It fixes the clamping from issue 98  as well as a few other instances of 
inverse trig
functions.  My old patch tested in the report for issue 101, which is reported 
to
have failed should also be taken care of since it is now a full clamp not a 
single
sided one.  
It cleans up calculateDiffAxisAngleQuaternion since clamping protection was 
added in
getAngle, and the farthest calculation was wrong.  
And it resolves the redundancy of slerp notation for this issue.

So I've checked that calculateDiffAxisAngleQuaternion works with the following 
code.
 It's not exhaustive.  Before removing the farthest it was failing this test.  After
the cleanup it passes.

TEST(Bullet, calculateDiffAxisAngleQuaternion)
{
  btVector3 vec;
  btScalar ang;
  for (unsigned int i = 0 ; i < 1000 ; i++)
  {
    btQuaternion q1(M_PI*2 *(double) i / 1000, 0, 0);
    btQuaternion q2(M_PI/2*0, 0,0);
    btTransformUtil::calculateDiffAxisAngleQuaternion(q1, q2, vec, ang);
    //    printf("%f %f %f, ang %f\n", vec.x(), vec.y(), vec.z(), ang);             

    EXPECT_NEAR(M_PI*2 *(double) i / 1000, ang, 0.001);
    btTransformUtil::calculateDiffAxisAngleQuaternion(q2, q1, vec, ang);
    //    printf("%f %f %f, ang %f\n", vec.x(), vec.y(), vec.z(), ang);             

    EXPECT_NEAR(M_PI*2 *(double) i / 1000, ang, 0.001);
  }
  for (unsigned int i = 0 ; i < 1000 ; i++)
  {
    btQuaternion q1(0, M_PI*2 *(double) i / 1000,1);
    btQuaternion q2(0, 0, 1);
    btTransformUtil::calculateDiffAxisAngleQuaternion(q1, q2, vec, ang);
    //    printf("%f %f %f, ang %f\n", vec.x(), vec.y(), vec.z(), ang);             

    EXPECT_NEAR(M_PI*2 *(double) i / 1000, ang, 0.001);
    btTransformUtil::calculateDiffAxisAngleQuaternion(q2, q1, vec, ang);
    //    printf("%f %f %f, ang %f\n", vec.x(), vec.y(), vec.z(), ang);             

    EXPECT_NEAR(M_PI*2 *(double) i / 1000, ang, 0.001);
  }
  for (unsigned int i = 0 ; i < 1000 ; i++)
  {
    btQuaternion q1(0, 0, M_PI*2 *(double) i / 1000);
    btQuaternion q2(0, 0,0);
    btTransformUtil::calculateDiffAxisAngleQuaternion(q1, q2, vec, ang);
    //    printf("%f %f %f, ang %f\n", vec.x(), vec.y(), vec.z(), ang);             

    EXPECT_NEAR(M_PI*2 *(double) i / 1000, ang, 0.001);
    btTransformUtil::calculateDiffAxisAngleQuaternion(q2, q1, vec, ang);
    //    printf("%f %f %f, ang %f\n", vec.x(), vec.y(), vec.z(), ang);             

    EXPECT_NEAR(M_PI*2 *(double) i / 1000, ang, 0.001);
  }
}

Original comment by Tully.Foote on 2 Dec 2008 at 7:45

Attachments:

GoogleCodeExporter commented 9 years ago 
Sorry I just noticed that I accidentally included an additional method I was
developing btQuaternion.getAxis in the patch.  

It is a function to match getAngle.  I'd prefer it to be added. Attached is full
patch with the cleaner implementation of getAngle. 

There is also attached a patch with that function removed since it is 
unrelated.

Original comment by Tully.Foote on 5 Dec 2008 at 3:37

Attachments:

GoogleCodeExporter commented 9 years ago 
I'm working toward releasing my code with bullet built in, and I'd like not to
diverge from bullet mainline too much.  Is there any change of getting this 
change
merged in.  I've attached an updated patch against head.

Original comment by Tully.Foote on 4 Sep 2009 at 11:22

Attachments:

GoogleCodeExporter commented 9 years ago 
We applied the slerp_redundency.patch in latest trunk

Thanks for the fix!

Original comment by erwin.coumans on 3 Nov 2009 at 6:22

GoogleCodeExporter commented 9 years ago 
Hello, I updated the bullet sources from repo today, but the problem is actual.
There are assert(x <= btScalar(1.f)) fails on call btQuaternion::angle(const 
btQuaternion&) const;

At the start of slerp method :
280     btQuaternion slerp(const btQuaternion& q, const btScalar& t) const
281     {
282         btScalar theta = angle(q);

2 quaternions are generated by bullet and stay correct.

Now I added the temp first-aid code in btScalar::btACos :
if (x > 1.f) { x = 1.f; } this works,
but I hope the next update will fix the problem.

Original comment by elwood...@gmail.com on 9 Nov 2009 at 2:54

GoogleCodeExporter commented 9 years ago 

I'm confused about your report: btScalar::btACos already has this fix in latest 
trunk.

Can you explain more in detail, and provide a new patch against latest trunk?

Thanks a lot,
Erwin

Original comment by erwin.coumans on 9 Nov 2009 at 10:28

GoogleCodeExporter commented 9 years ago 
Hello, I breaked on this line in debugger and got dump of 2 transforms' 
openGLMatrix :

2009-11-10 16:13:03.163 bowlmoney[1585:20b] value1_transform OpenGLMatrix :
2009-11-10 16:13:03.164 bowlmoney[1585:20b] m[0] = -0.642792522907
2009-11-10 16:13:03.165 bowlmoney[1585:20b] m[1] = 0.766040325165
2009-11-10 16:13:03.166 bowlmoney[1585:20b] m[2] = -0.000225860480
2009-11-10 16:13:03.166 bowlmoney[1585:20b] m[3] = 0.000000000000
2009-11-10 16:13:03.167 bowlmoney[1585:20b] m[4] = -0.766040325165
2009-11-10 16:13:03.167 bowlmoney[1585:20b] m[5] = -0.642792463303
2009-11-10 16:13:03.167 bowlmoney[1585:20b] m[6] = 0.000271065888
2009-11-10 16:13:03.168 bowlmoney[1585:20b] m[7] = 0.000000000000
2009-11-10 16:13:03.168 bowlmoney[1585:20b] m[8] = 0.000062465995
2009-11-10 16:13:03.168 bowlmoney[1585:20b] m[9] = 0.000347257330
2009-11-10 16:13:03.168 bowlmoney[1585:20b] m[10] = 0.999999940395
2009-11-10 16:13:03.169 bowlmoney[1585:20b] m[11] = 0.000000000000
2009-11-10 16:13:03.169 bowlmoney[1585:20b] m[12] = 20.570852279663
2009-11-10 16:13:03.170 bowlmoney[1585:20b] m[13] = 739.996643066406
2009-11-10 16:13:03.170 bowlmoney[1585:20b] m[14] = 5.903697013855
2009-11-10 16:13:03.170 bowlmoney[1585:20b] m[15] = 1.000000000000
2009-11-10 16:13:03.171 bowlmoney[1585:20b] value2_transform OpenGLMatrix :
2009-11-10 16:13:03.171 bowlmoney[1585:20b] m[0] = -0.642792999744
2009-11-10 16:13:03.171 bowlmoney[1585:20b] m[1] = 0.766040086746
2009-11-10 16:13:03.172 bowlmoney[1585:20b] m[2] = -0.000226838863
2009-11-10 16:13:03.172 bowlmoney[1585:20b] m[3] = 0.000000000000
2009-11-10 16:13:03.172 bowlmoney[1585:20b] m[4] = -0.766040086746
2009-11-10 16:13:03.173 bowlmoney[1585:20b] m[5] = -0.642792940140
2009-11-10 16:13:03.173 bowlmoney[1585:20b] m[6] = 0.000282453373
2009-11-10 16:13:03.173 bowlmoney[1585:20b] m[7] = 0.000000000000
2009-11-10 16:13:03.174 bowlmoney[1585:20b] m[8] = 0.000070560200
2009-11-10 16:13:03.174 bowlmoney[1585:20b] m[9] = 0.000355326687
2009-11-10 16:13:03.174 bowlmoney[1585:20b] m[10] = 0.999999940395
2009-11-10 16:13:03.175 bowlmoney[1585:20b] m[11] = 0.000000000000
2009-11-10 16:13:03.175 bowlmoney[1585:20b] m[12] = 20.570899963379
2009-11-10 16:13:03.176 bowlmoney[1585:20b] m[13] = 739.996704101562
2009-11-10 16:13:03.177 bowlmoney[1585:20b] m[14] = 5.903925895691
2009-11-10 16:13:03.177 bowlmoney[1585:20b] m[15] = 1.000000000000

ratioBetweenValue1AndValue2 = 0.999999225

Calling code :

btTransform value1_transform;
btTransform value2_transform;
value1.getTransform(value1_transform);
value2.getTransform(value2_transform);

DZLog(@"value1_transform OpenGLMatrix :")
btScalar m[16];
value1_transform.getOpenGLMatrix(m);
for (int i = 0; i < 16; i++) {
    DZLog(@"m[%u] = %.12f", i, m[i])
}
DZLog(@"value2_transform OpenGLMatrix :")
value2_transform.getOpenGLMatrix(m);
for (int i = 0; i < 16; i++) {
    DZLog(@"m[%u] = %.12f", i, m[i])
}

value1.setRotation(value1_transform.getRotation().slerp(value2_transform.getRota
tion(),
btScalar(ratioBetweenValue1AndValue2)));

Oh, I see. My code calling the btAcos from 251 row in btScalar.h
251 SIMD_FORCE_INLINE btScalar btAcos(btScalar x) { 
    btAssert(x <= btScalar(1.));  // <- cause an error
    if (x<btScalar(-1)) 
        x=btScalar(-1); 
    if (x>btScalar(1))  
        x=btScalar(1);
    return acosf(x); 
}

because condition (defined(BT_USE_DOUBLE_PRECISION) || 
defined(BT_FORCE_DOUBLE_FUNCTIONS)
) doesn't take true.

Original comment by elwood...@gmail.com on 10 Nov 2009 at 1:24

GoogleCodeExporter commented 9 years ago 
ok, with clamping in place we removed the assert.
http://code.google.com/p/bullet/source/detail?r=1863

Thanks for the report.
Erwin

Original comment by erwin.coumans on 21 Dec 2009 at 11:44