Certain unit quaternions produce invalid rotation matrices

[lierdakil]

Case in point, take this quaternion:

let q = UnitQuaternion::new_normalize(Quaternion::new(
    0.1045170328727042,
    -0.6993398242910853,
    0.10451703287270421,
    0.6993398242910854,
));

then convert it to rotation matrix:

dbg!(q.to_rotation_matrix());

You'll likely get:

[
   [
        -1.6653345369377348e-16,
        0.0,
        -1.0000000000000002,
    ],
    [
        -0.2923716936184902,
        -0.9563047593579446,
        2.7755575615628914e-17,
    ],
    [
        -0.9563047593579445,
        0.29237169361849025,
        1.6653345369377348e-16,
    ],
]

Notice that the last component of the first row is less than -1. This will create issues (of the NAN variety) if, say, this matrix is then converted to Euler angles.

[mikebauer]

I reproduced the same behavior, but I think that matrix is a valid rotation (it's orthonormal and its determinant is 1). Here's the code I'm running to quickly verify:

use nalgebra::{UnitQuaternion, Quaternion};

fn main() {
    let q = UnitQuaternion::new_normalize(Quaternion::new(
        0.1045170328727042,
        -0.6993398242910853,
        0.10451703287270421,
        0.6993398242910854,
    ));

    let R = q.to_rotation_matrix();

    dbg!(R);

    dbg!(R.matrix().determinant());

    dbg!(R.transpose() * R);
}

Which yields:

[tmp/src/main.rs:15:5] R = [
    [
        -1.6653345369377348e-16,
        0.0,
        -1.0000000000000002,
    ],
    [
        -0.2923716936184902,
        -0.9563047593579446,
        2.7755575615628914e-17,
    ],
    [
        -0.9563047593579445,
        0.29237169361849025,
        1.6653345369377348e-16,
    ],
]
[tmp/src/main.rs:17:5] R.matrix().determinant() = 1.0000000000000007
[tmp/src/main.rs:19:5] R.transpose() * R = [
    [
        1.0000000000000004,
        2.0934092284956047e-17,
        -7.276719334102089e-18,
    ],
    [
        2.0934092284956047e-17,
        1.0000000000000007,
        -5.551115123125782e-17,
    ],
    [
        -7.276719334102089e-18,
        -5.551115123125782e-17,
        1.0000000000000004,
    ],
]

Because the set of orthonormal matrices with determinant 1 is a unique representation of SO(3), there should be a unique set of Euler angles for it. Furthermore, I don't think there's any other such matrix which represents the same orientation. If you let me know which code you're running that results in a NaN or other numerical issues, I'm happy to have a look. Feel free to let me know if I've missed anything.

[lierdakil]

Taking asin/acos of a value larger than 1.0 produces NAN on most architectures. So...

A pure rotation matrix, by construction, simply can't have the absolute value of its elements be greater than 1. -1.000...02 is a rounding error -- it should've been -1.0 (the quaternion in question is pointing at nadir IIRC)

Here's a specific example where this breaks:

use nalgebra::*;

fn main() {
    let q = UnitQuaternion::new_normalize(Quaternion::new(
        0.1045170328727042,
        -0.6993398242910853,
        0.10451703287270421,
        0.6993398242910854,
    ));
    dbg!(q.to_rotation_matrix().euler_angles_ordered(
        [
            UnitVector3::new_normalize(Vector3::z()),
            UnitVector3::new_normalize(Vector3::y()),
            UnitVector3::new_normalize(Vector3::x())
        ],
        false
    ));
}

outputs

(
    [
        2.677945044588987,
        NaN,
        0.0,
    ],
    false,
)

[Ralith]

euler_angles_ordered should probably internally clamp into the valid range. In general, any rotation matrix (or other normalized structure) you might handle is going to be prone to that type of rounding error.

[lierdakil]

There's an argument that UnitQuaternion::to_rotation_matrix should clamp to the valid range, as Rotation3 is, in concept, a rotation matrix, so, in principle, it shouldn't have elements' modulus > 1. Way I see it, that's one of the main points of distinguishing Rotation3 from just Matrix3. Doing this any other way is just asking for more issues like this to crop up, because of the implicit assumption that Rotation3 is indeed a valid rotation matrix... which doesn't necessarily hold, mind, because _unchecked exists, so it's not a perfect fix. But still, I really don't expect to encounter NANs (or having to check the validity of each -- possibly hidden -- intermediate result manually) if I don't use _unchecked in my code and there are no bogus inputs.

[lierdakil]

That being said, I guess there may be performance reasons not to do it this way (as it would have to be done any time the matrix is touched, not just when it's constructed), so I'm fine with just adding a clamp to euler_angles_ordered and moving on. If anyone's going to have a look at euler_angles_ordered though, consider also checking #1482, I think nadir handling there is a bit wonky in general (whether it's a bug in the implementation or in the original algorithm I have no idea).

[Ralith]

I really don't expect to encounter NANs (or having to check the validity of each -- possibly hidden -- intermediate result manually) if I don't use _unchecked in my code

Then you're going to have a bad time. Rounding errors will happen as a result of any operation, e.g. composing rotations together. It's generally your responsibility to renormalize and clamp as needed.

[lierdakil]

Frankly, that's a very strange stance to take. If this "A valid Rotation3 isn't guaranteed to remain valid via valid operations, and will eventually throw NANs at you if you don't police it" is the official stance, then I suggest it should be written in large letters all over its docs, because I would guess most people don't expect it to behave that way.

[Ralith]

Doc PRs would be very welcome! This is the behavior throughout nalgebra. Rotation3::renormalize and Unit::renormalize may be useful references.

[lierdakil]

FTR, Rotation3::renormalize doesn't, in fact, clamp to -1..1.

[Ralith]

As I said in #1481 (comment), that's the responsibility of the logic that has that constraint, i.e. euler_angles_ordered.