rotation_quat fails in some cases

[lucieb]

With rotation matrix {0,-1,0,0,-1,0,0,0,-1}, the function rotation_quat returns {0.7071068,0.7071068,0,0} instead of {-0.7071068,0.7071068,0,0}.

[melax]

check the matrix first, since the above is not just non-orthonormal, its singular

[lucieb]

Typo, I meant {{0,-1,0},{-1,0,0},{0,0,-1}} (row-ordered)

[melax]

off topic: note that if the input has a negative determinant. its is a reflection. thus cant be expressed as just a rotation.

[melax]

testing now I get the right answer:

    float3x3 m({ 0,-1,0 }, { -1,0,0 }, { 0,0,-1 });
    float4 q = quatfrommat(m);
    std::cout << q << std::endl; 

gives

-0.707107 0.707107 0 0

I'll test your copy of linalg.h to make sure its the same as mine.

[melax]

OOps, just realized my test was using my quatfrommat() code instead. , looks like linalg's rotation_quat might be broken. in the meantime, here's some code that works.

inline float4 quatfrommat(const float3x3 &m)
{
    float magw = m[0][0] + m[1][1] + m[2][2];
    float magxy;
    float magzw;
    float3 pre;
    float3 prexy;
    float3 prezw;

    bool wvsz = magw  > m[2][2];
    magzw = wvsz ? magw : m[2][2];
    prezw = wvsz ? float3(1, 1, 1) : float3(-1, -1, 1);
    auto postzw = wvsz ? float4(0, 0, 0, 1) : float4(0, 0, 1, 0);

    bool xvsy = m[0][0] > m[1][1];
    magxy = xvsy ? m[0][0] : m[1][1];
    prexy = xvsy ? float3(1, -1, -1) : float3(-1, 1, -1);
    auto postxy = xvsy ? float4(1, 0, 0, 0) : float4(0, 1, 0, 0);

    bool zwvsxy = magzw > magxy;
    pre = zwvsxy ? prezw : prexy;
    auto post = zwvsxy ? postzw : postxy;

    float t = pre.x * m[0][0] + pre.y * m[1][1] + pre.z * m[2][2] + 1;
    float s = 1 / sqrt(t) / 2;
    float4 qp
    {
        (pre.y * m[1][2] - pre.z * m[2][1]) * s,
        (pre.z * m[2][0] - pre.x * m[0][2]) * s,
        (pre.x * m[0][1] - pre.y * m[1][0]) * s,
        t * s
    };
    return qmul(qp, post);
}

[sgorsten]

Good catch. This is a numerical stability issue in my implementation related to exact 90 degree rotations.

I've added a failing test case for now while I investigate my options. Stan's function works perfectly, to the best of my knowledge. I'd like to see if I can come up with something a little terser, but definitely make use of his function if you have an immediate need.

[sgorsten]

And it looks like we've got it! Do let me know if you encounter any other corner cases, and I apologize for the inconvenience.

[melax]

wow that implementation is definitely clean and compact. very nice. thx