Possible bug in Matrix3::ToEulerAnglesZYX

[samapico]

Here's a sample of a GTEST unit test I made that easily reproduces the problem:

Code: Select all

    Ogre::Matrix3 m;

    Ogre::Radian y,p,r;

    m.FromEulerAnglesZYX(Ogre::Degree(0), Ogre::Degree(-90), Ogre::Degree(0));

    m.ToEulerAnglesZYX(y,p,r);

    EXPECT_NEAR(y.valueDegrees(),   0, 0.00001);
    EXPECT_NEAR(p.valueDegrees(), -90, 0.00001);
    EXPECT_NEAR(r.valueDegrees(),   0, 0.00001);

The result of this is:
y == 180
p == -90
r == 0

Now, I know about the gimbal lock thing with yaw-pitch-roll values, but no matter how you look at this, the result is at a 180 degree rotation from the expected result. I'm also aware of the not-so-intuitive notion that the 'y' parameter in this case is actually the Z rotation (roll), since the order of the parameters match the order of the rotations.
Note that the rotation matrix obtained with 'FromEulerAnglesZYX' is correct, so the problem is definitely in 'ToEulerAnglesZYX'. I don't know if similar issues exist with different rotation orders as I have not done any tests with them.

Using Ogre 1.8.1, but there are no changes to that part of the code since that version, so a newer version would certainly not fix the issue.
The issue happens with p == -90, no matter the values of y and r. It does not happen with p == +90.

This looks as a bug to me, unless I'm missing something. Any ideas?

For reference, here's Matrix3::ToEulerAnglesZYX with some annotations:

Code: Select all

bool Matrix3::ToEulerAnglesZYX (Radian& rfYAngle, Radian& rfPAngle,
        Radian& rfRAngle) const
    {
        // rot =  cy*cz           cz*sx*sy-cx*sz  cx*cz*sy+sx*sz
        //        cy*sz           cx*cz+sx*sy*sz -cz*sx+cx*sy*sz
        //       -sy              cy*sx           cx*cy

        rfPAngle = Math::ASin(-m[2][0]);   ///<<< This gives -1.57..., or -PI/2
        if ( rfPAngle < Radian(Math::HALF_PI) )
        {
            if ( rfPAngle > Radian(-Math::HALF_PI) )
            {
                rfYAngle = Math::ATan2(m[1][0],m[0][0]);
                rfRAngle = Math::ATan2(m[2][1],m[2][2]);
                return true;
            }
            else
            {
                // WARNING.  Not a unique solution.                     ///< We enter this case
                Radian fRmY = Math::ATan2(-m[0][1],m[0][2]);     ///< This gives -PI; it would need to be 0 somehow for the result to make sense
                rfRAngle = Radian(0.0);  // any angle works
                rfYAngle = rfRAngle - fRmY;
                return false;
            }
        }
        else
        {
            // WARNING.  Not a unique solution.
            Radian fRpY = Math::ATan2(-m[0][1],m[0][2]);
            rfRAngle = Radian(0.0);  // any angle works
            rfYAngle = fRpY - rfRAngle;
            return false;
        }
    }

Edit:
I think that:

Code: Select all

Radian fRmY = Math::ATan2(-m[0][1],m[0][2]);

should be

Code: Select all

Radian fRmY = Math::ATan2(-m[1][2],m[1][1]);

With P = 90, both are equivalent, but not with P = -90.

I'm having a hard time finding sources that show how to properly deal with the gimbal lock. I got that from some old LabView code I have, I'm not sure if it's right...

[samapico]

Some quick similar tests seem to indicate that the XYZ and YXZ orders might have similar issues.

Pseudo code of the tests I ran:

Code: Select all

for yaw = -90 to 90    (can use steps of 5 or 10 degrees, exceptions show mostly at -90 and 90 anyway)
    for pitch = -90 to 90
        for roll = -90 to 90
            Matrix3 m1.FromEulerAngles(yaw, pitch, roll)

            m1.ToEulerAngles(yaw2,pitch2,roll2)

            Matrix3 m2.FromEulerAngles(yaw2,pitch2,roll2)

            for i = 0 to 2
                for j = 0 to 2
                     test abs(m1[i][j] - m2[i][j]) < 0.0001

//Repeat with all 6 possible rotation orders

[Kojack]

Ogre's math (at lest some of it) was originally taken from Wild Magic. Here's their version of the matrix 3 code: http://www.geometrictools.com/LibMathem ... atrix3.inl
Their ExtractEulerZYX is the equivalent to our ToEulerAnglesZYX.
The first part of the code is roughly identical in functionality (we do the asin before the two ifs, they do it after, but otherwise the same).
But the two edge cases (pitch 90 and pitch -90) are different.

Ours for -90 pitch:

Code: Select all

// rfPAngle is already Math::ASin(-m[2][0]);
Radian fRmY = Math::ATan2(-m[0][1],m[0][2]);     ///< This gives -PI; it would need to be 0 somehow for the result to make sense
rfRAngle = Radian(0.0);  // any angle works
rfYAngle = rfRAngle - fRmY;

Theirs:

Code: Select all

yAngle = Math<Real>::HALF_PI;
zAngle = -Math<Real>::ATan2(mEntry[1], mEntry[2]);
xAngle = (Real)0;

(m[0][1] is the same element as mEntry[1] and m[0][2] is the same element as mEntry[2])
Ours keeps the pitch angle from the beginning, theirs replaces it.

I haven't actually run the numbers yet to see if that's enough to fix it, I'm a little distracted at the moment waiting for a courier delivery that I missed yesterday.