/*
 * convert a rotation matrix to quaternion
 * 
 * Quat class member variables:
 * float w; //angle of displacement: cos(theta/2)
 * float x; // x-axis of rotation: sin(theta/2)nx
 * float y; // y-axis of rotation: sin(theta/2)ny
 * float z; // z-axis of rotation: sin(theta/2)nz
 *
 * Mat4 class member variables:
 * float elem[16];
 *
 */
Quat Mat2Quat(const Mat4& m)
{
    // matrix to quaternion conversion as defined in
    //Dunn/Parberry's "3D Math Primer for Graphics and Game Development"
    //p. 284 - 287
    
    Quat m2q;
    float absVal[4];
    float m11 = m.elem[0];
    float m12 = m.elem[1];
    float m13 = m.elem[2];
    float m21 = m.elem[4];
    float m22 = m.elem[5];
    float m23 = m.elem[6];
    float m31 = m.elem[8];
    float m32 = m.elem[9];
    float m33 = m.elem[10];
    
    //determine which component has the largest absolute value
    absVal[0] = m11 + m22 + m33; //W component
    absVal[1] = m11 - m22 - m33;
    absVal[2] = m22 - m11 - m33;
    absVal[3] = m33 - m11 - m22;
    
    int idx = 0; //index of component with the largest value
    float bigVal = absVal[0]; //holds the largest value
    for (int i = 1; i < 4; i++)
    {
        if (absVal[i] > bigVal)
        {
            bigVal = absVal[i];
            idx = i;
        }
    }
    
    //perform square root and division to extract the
    //quaternion value from the matrix
    bigVal = sqrtf(bigVal + 1.0f) * 0.5f;
    float div = 0.25f / bigVal;
    
    //apply the table to calculate the quaternion
    switch (idx)
    {
    case 0:    //W component is largest
        m2q.w = bigVal;
        m2q.x = (m23 - m32) * div;
        m2q.y = (m31 - m13) * div;
        m2q.z = (m12 - m21) * div;
        break;
    case 1:    //X component is largest
        m2q.w = (m23 - m32) * div;
        m2q.x = bigVal;
        m2q.y = (m12 + m21) * div;
        m2q.z = (m31 + m13) * div;
        break;
    case 2:    //Y component is largest
        m2q.w = (m31 - m13) * div;
        m2q.x = (m12 + m21) * div;
        m2q.y = bigVal;
        m2q.z = (m23 + m32) * div;
        break;
    case 3:    //Z component is largest
        m2q.w = (m12 - m21) * div;
        m2q.x = (m31 + m13) * div;
        m2q.y = (m23 + m32) * div;
        m2q.z = bigVal;
        break;
    }
    return m2q;
}