How to obtain axis-angle from rotation matrix?

Question

I need to obtain some data from an openGL rotation matrix. I need to obtain the equivalent euler angles (already did it), the equivalent quaternion (did it, but just copying it from the Internet) and the equivalent axis-angle.

I dont know if a rotation matrix can be expresed as a single rotation of a certain angle around an certain vector. Are these equivalent? If they are, how can I obtain one from the other?

Also, i would like to understand better the meaning of a quaternion, and the insides of a rotation matrix. Where should i go to learn about this?

Answer 1

Yes any rotation matrix/unit quaternion is equivalent to a rotation around a single axis. If we call this axis n and the angle theta then the quaternion for this rotation is:

[n * sin(theta / 2) cos(theta / 2)]

To reconstruct this use acos on the w element of the quaternion to get theta / 2. After you have theta you can divide x,y and z component with sin(theta / 2) to reconstruct the axis.

Answer 2

Here's a function which converts a 3x3 matrix into an axis, angle (using a quatention, so perhaps theres a more efficient way which bypasses that step).

void axis_angle_from_mat3(float r_axis[3], float *r_angle, float mat[3][3])
{
    float q[4];

    /* -------------------------------------------------------------------- */
    /* matrix to quaternion */
    double tr, s;
    float tmat[3][3];

    /* work on a copy */
    memcpy(tmat, mat, sizeof(tmat));

    /* normalize the matrix */
    int i;
    for (i = 0; i < 3; i++) {
        float d = (tmat[i][0] * tmat[i][0] + tmat[i][1] * tmat[i][1] + tmat[i][2] * tmat[i][2]);

        if (d > 1.0e-35f) {
            d = sqrtf(d);
            tmat[i][0] /= d;
            tmat[i][1] /= d;
            tmat[i][2] /= d;
        }
        else {
            tmat[i][0] = 0.0f;
            tmat[i][1] = 0.0f;
            tmat[i][2] = 0.0f;
            d = 0.0f;
        }
    }


    tr = 0.25 * (double)(1.0f + tmat[0][0] + tmat[1][1] + tmat[2][2]);

    if (tr > (double)1e-4f) {
        s = sqrt(tr);
        q[0] = (float)s;
        s = 1.0 / (4.0 * s);
        q[1] = (float)((double)(tmat[1][2] - tmat[2][1]) * s);
        q[2] = (float)((double)(tmat[2][0] - tmat[0][2]) * s);
        q[3] = (float)((double)(tmat[0][1] - tmat[1][0]) * s);
    }
    else {
        if (tmat[0][0] > tmat[1][1] && tmat[0][0] > tmat[2][2]) {
            s = 2.0f * sqrtf(1.0f + tmat[0][0] - tmat[1][1] - tmat[2][2]);
            q[1] = (float)(0.25 * s);

            s = 1.0 / s;
            q[0] = (float)((double)(tmat[1][2] - tmat[2][1]) * s);
            q[2] = (float)((double)(tmat[1][0] + tmat[0][1]) * s);
            q[3] = (float)((double)(tmat[2][0] + tmat[0][2]) * s);
        }
        else if (tmat[1][1] > tmat[2][2]) {
            s = 2.0f * sqrtf(1.0f + tmat[1][1] - tmat[0][0] - tmat[2][2]);
            q[2] = (float)(0.25 * s);

            s = 1.0 / s;
            q[0] = (float)((double)(tmat[2][0] - tmat[0][2]) * s);
            q[1] = (float)((double)(tmat[1][0] + tmat[0][1]) * s);
            q[3] = (float)((double)(tmat[2][1] + tmat[1][2]) * s);
        }
        else {
            s = 2.0f * sqrtf(1.0f + tmat[2][2] - tmat[0][0] - tmat[1][1]);
            q[3] = (float)(0.25 * s);

            s = 1.0 / s;
            q[0] = (float)((double)(tmat[0][1] - tmat[1][0]) * s);
            q[1] = (float)((double)(tmat[2][0] + tmat[0][2]) * s);
            q[2] = (float)((double)(tmat[2][1] + tmat[1][2]) * s);
        }
    }


    /* normalize the quat */
    float len;
    len = sqrtf(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]);
    if (len != 0.0f) {
        q[0] /= len;
        q[1] /= len;
        q[2] /= len;
        q[3] /= len;
    }
    else {
        q[1] = 1.0f;
        q[0] = q[2] = q[3] = 0.0f;
    }


    /* -------------------------------------------------------------------- */
    /* quaternion to axis angle */

    float ha, si;

    ha = acosf(q[0]);
    si = sinf(ha);

    *r_angle = ha * 2;

    if (fabsf(si) < FLT_EPSILON)
        si = 1.0f;

    r_axis[0] = q[1] / si;
    r_axis[1] = q[2] / si;
    r_axis[2] = q[3] / si;
}