matrix4.inl

Go to the documentation of this file.

/*
    Sheep - A Rigid Body Dynamics Engine
    Copyright (C) 2001-2004 Francois Beaune
    Contact: http://toxicengine.sourceforge.net/

    This file is part of Sheep.

    Sheep is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    Sheep is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sheep; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

inline
Matrix4::Matrix4() {}

inline
Matrix4::Matrix4(Real r) {
    m_v[0] = m_v[1] = m_v[2] = m_v[3] =
    m_v[4] = m_v[5] = m_v[6] = m_v[7] =
    m_v[8] = m_v[9] = m_v[10] = m_v[11] =
    m_v[12] = m_v[13] = m_v[14] = m_v[15] = r;
}

inline
Matrix4::Matrix4(const Matrix4 &m) {
    m_v[0] = m.m_v[0];
    m_v[1] = m.m_v[1];
    m_v[2] = m.m_v[2];
    m_v[3] = m.m_v[3];
    m_v[4] = m.m_v[4];
    m_v[5] = m.m_v[5];
    m_v[6] = m.m_v[6];
    m_v[7] = m.m_v[7];
    m_v[8] = m.m_v[8];
    m_v[9] = m.m_v[9];
    m_v[10] = m.m_v[10];
    m_v[11] = m.m_v[11];
    m_v[12] = m.m_v[12];
    m_v[13] = m.m_v[13];
    m_v[14] = m.m_v[14];
    m_v[15] = m.m_v[15];
}

inline
Matrix4 Matrix4::Identity() {
    Matrix4 res;

    res.m_v[0] = res.m_v[5] = res.m_v[10] = res.m_v[15] = 1.0;
    res.m_v[1] = res.m_v[2] = res.m_v[3] = res.m_v[4] =
    res.m_v[6] = res.m_v[7] = res.m_v[8] = res.m_v[9] =
    res.m_v[11] = res.m_v[12] = res.m_v[13] = res.m_v[14] = 0.0;

    return res;
}

inline
Matrix4 Matrix4::Rotation(const Quaternion &q) {    // q must be normalized
    Matrix4 res;

    res.m_v[0] = 1.0 - 2.0 * (q.m_v.m_y * q.m_v.m_y + q.m_v.m_z * q.m_v.m_z);
    res.m_v[1] = 2.0 * (q.m_v.m_x * q.m_v.m_y - q.m_s * q.m_v.m_z);
    res.m_v[2] = 2.0 * (q.m_v.m_x * q.m_v.m_z + q.m_s * q.m_v.m_y);
    res.m_v[4] = 2.0 * (q.m_v.m_x * q.m_v.m_y + q.m_s * q.m_v.m_z);
    res.m_v[5] = 1.0 - 2.0 * (q.m_v.m_x * q.m_v.m_x + q.m_v.m_z * q.m_v.m_z);
    res.m_v[6] = 2.0 * (q.m_v.m_y * q.m_v.m_z - q.m_s * q.m_v.m_x);
    res.m_v[8] = 2.0 * (q.m_v.m_x * q.m_v.m_z - q.m_s * q.m_v.m_y);
    res.m_v[9] = 2.0 * (q.m_v.m_y * q.m_v.m_z + q.m_s * q.m_v.m_x);
    res.m_v[10] = 1.0 - 2.0 * (q.m_v.m_x * q.m_v.m_x + q.m_v.m_y * q.m_v.m_y);

    res.m_v[3] = res.m_v[7] = res.m_v[11] =
    res.m_v[12] = res.m_v[13] = res.m_v[14] = 0.0;
    res.m_v[15] = 1.0;

    return res;
}

inline
Matrix4 Matrix4::Rotation(Real angle, const Vector3 &axis) {
    //!\todo implement.
    assert(!"Not implemented yet.");
    return Matrix4();
}

inline
Matrix4 Matrix4::Rotation(Real rotx, Real roty, Real rotz) {
    const Real a = cos(rotx);
    const Real b = sin(rotx);
    const Real c = cos(roty);
    const Real d = sin(roty);
    const Real e = cos(rotz);
    const Real f = sin(rotz);

    const Real ad = a * d;
    const Real bd = b * d;

    Matrix4 res;

    res.m_v[0] = c * e;
    res.m_v[1] = -a * f + bd * e;
    res.m_v[2] = b * f + ad * e;
    res.m_v[3] = 0.0;
    res.m_v[4] = c * f;
    res.m_v[5] = a * e + bd * f;
    res.m_v[6] = -b * e + ad * f;
    res.m_v[7] = 0.0;
    res.m_v[8] = -d;
    res.m_v[9] = b * c;
    res.m_v[10] = a * c;
    res.m_v[11] = 0.0;
    res.m_v[12] = 0.0;
    res.m_v[13] = 0.0;
    res.m_v[14] = 0.0;
    res.m_v[15] = 1.0;

    return res;
}

inline
Matrix4 Matrix4::Scale(const Vector3 &s) {
    Matrix4 res;

    res.m_v[0] = s.m_x;
    res.m_v[5] = s.m_y;
    res.m_v[10] = s.m_z;
    res.m_v[15] = 1.0;

    res.m_v[1] = res.m_v[2] = res.m_v[3] = res.m_v[4] =
    res.m_v[6] = res.m_v[7] = res.m_v[8] = res.m_v[9] =
    res.m_v[11] = res.m_v[12] = res.m_v[13] = res.m_v[14] = 0.0;

    return res;
}

inline
Matrix4 Matrix4::Translation(const Vector3 &t) {
    Matrix4 res;

    res.m_v[0] = res.m_v[5] = res.m_v[10] = res.m_v[15] = 1.0;
    res.m_v[1] = res.m_v[2] = res.m_v[4] =
    res.m_v[6] = res.m_v[8] = res.m_v[9] =
    res.m_v[12] = res.m_v[13] = res.m_v[14] = 0.0;

    res.m_v[3] = t.m_x;
    res.m_v[7] = t.m_y;
    res.m_v[11] = t.m_z;

    return res;
}

inline
Matrix4 &Matrix4::operator=(Real r) {
    m_v[0] = m_v[1] = m_v[2] = m_v[3] =
    m_v[4] = m_v[5] = m_v[6] = m_v[7] =
    m_v[8] = m_v[9] = m_v[10] = m_v[11] =
    m_v[12] = m_v[13] = m_v[14] = m_v[15] = r;

    return *this;
}

inline
Matrix4 &Matrix4::operator=(const Matrix4 &m) {
    m_v[0] = m.m_v[0];
    m_v[1] = m.m_v[1];
    m_v[2] = m.m_v[2];
    m_v[3] = m.m_v[3];
    m_v[4] = m.m_v[4];
    m_v[5] = m.m_v[5];
    m_v[6] = m.m_v[6];
    m_v[7] = m.m_v[7];
    m_v[8] = m.m_v[8];
    m_v[9] = m.m_v[9];
    m_v[10] = m.m_v[10];
    m_v[11] = m.m_v[11];
    m_v[12] = m.m_v[12];
    m_v[13] = m.m_v[13];
    m_v[14] = m.m_v[14];
    m_v[15] = m.m_v[15];

    return *this;
}

inline
Matrix4 Matrix4::operator+(const Matrix4 &m) const {
    Matrix4 res;

    res.m_v[0] = m_v[0] + m.m_v[0];
    res.m_v[1] = m_v[1] + m.m_v[1];
    res.m_v[2] = m_v[2] + m.m_v[2];
    res.m_v[3] = m_v[3] + m.m_v[3];
    res.m_v[4] = m_v[4] + m.m_v[4];
    res.m_v[5] = m_v[5] + m.m_v[5];
    res.m_v[6] = m_v[6] + m.m_v[6];
    res.m_v[7] = m_v[7] + m.m_v[7];
    res.m_v[8] = m_v[8] + m.m_v[8];
    res.m_v[9] = m_v[9] + m.m_v[9];
    res.m_v[10] = m_v[10] + m.m_v[10];
    res.m_v[11] = m_v[11] + m.m_v[11];
    res.m_v[12] = m_v[12] + m.m_v[12];
    res.m_v[13] = m_v[13] + m.m_v[13];
    res.m_v[14] = m_v[14] + m.m_v[14];
    res.m_v[15] = m_v[15] + m.m_v[15];

    return res;
}

inline
Matrix4 Matrix4::operator-(const Matrix4 &m) const {
    Matrix4 res;

    res.m_v[0] = m_v[0] - m.m_v[0];
    res.m_v[1] = m_v[1] - m.m_v[1];
    res.m_v[2] = m_v[2] - m.m_v[2];
    res.m_v[3] = m_v[3] - m.m_v[3];
    res.m_v[4] = m_v[4] - m.m_v[4];
    res.m_v[5] = m_v[5] - m.m_v[5];
    res.m_v[6] = m_v[6] - m.m_v[6];
    res.m_v[7] = m_v[7] - m.m_v[7];
    res.m_v[8] = m_v[8] - m.m_v[8];
    res.m_v[9] = m_v[9] - m.m_v[9];
    res.m_v[10] = m_v[10] - m.m_v[10];
    res.m_v[11] = m_v[11] - m.m_v[11];
    res.m_v[12] = m_v[12] - m.m_v[12];
    res.m_v[13] = m_v[13] - m.m_v[13];
    res.m_v[14] = m_v[14] - m.m_v[14];
    res.m_v[15] = m_v[15] - m.m_v[15];

    return res;
}

inline
Matrix4 Matrix4::operator-() const {
    Matrix4 res;

    res.m_v[0] = - m_v[0];
    res.m_v[1] = - m_v[1];
    res.m_v[2] = - m_v[2];
    res.m_v[3] = - m_v[3];
    res.m_v[4] = - m_v[4];
    res.m_v[5] = - m_v[5];
    res.m_v[6] = - m_v[6];
    res.m_v[7] = - m_v[7];
    res.m_v[8] = - m_v[8];
    res.m_v[9] = - m_v[9];
    res.m_v[10] = - m_v[10];
    res.m_v[11] = - m_v[11];
    res.m_v[12] = - m_v[12];
    res.m_v[13] = - m_v[13];
    res.m_v[14] = - m_v[14];
    res.m_v[15] = - m_v[15];

    return res;
}

inline
Matrix4 Matrix4::operator*(Real r) const {
    Matrix4 res;

    res.m_v[0] = m_v[0] * r;
    res.m_v[1] = m_v[1] * r;
    res.m_v[2] = m_v[2] * r;
    res.m_v[3] = m_v[3] * r;
    res.m_v[4] = m_v[4] * r;
    res.m_v[5] = m_v[5] * r;
    res.m_v[6] = m_v[6] * r;
    res.m_v[7] = m_v[7] * r;
    res.m_v[8] = m_v[8] * r;
    res.m_v[9] = m_v[9] * r;
    res.m_v[10] = m_v[10] * r;
    res.m_v[11] = m_v[11] * r;
    res.m_v[12] = m_v[12] * r;
    res.m_v[13] = m_v[13] * r;
    res.m_v[14] = m_v[14] * r;
    res.m_v[15] = m_v[15] * r;

    return res;
}

inline
Matrix4 Matrix4::operator/(Real r) const {
    Real inv_r = 1.0 / r;

    Matrix4 res;

    res.m_v[0] = m_v[0] * inv_r;
    res.m_v[1] = m_v[1] * inv_r;
    res.m_v[2] = m_v[2] * inv_r;
    res.m_v[3] = m_v[3] * inv_r;
    res.m_v[4] = m_v[4] * inv_r;
    res.m_v[5] = m_v[5] * inv_r;
    res.m_v[6] = m_v[6] * inv_r;
    res.m_v[7] = m_v[7] * inv_r;
    res.m_v[8] = m_v[8] * inv_r;
    res.m_v[9] = m_v[9] * inv_r;
    res.m_v[10] = m_v[10] * inv_r;
    res.m_v[11] = m_v[11] * inv_r;
    res.m_v[12] = m_v[12] * inv_r;
    res.m_v[13] = m_v[13] * inv_r;
    res.m_v[14] = m_v[14] * inv_r;
    res.m_v[15] = m_v[15] * inv_r;

    return res;
}

inline
Matrix4 &Matrix4::operator+=(const Matrix4 &m) {
    m_v[0] += m.m_v[0];
    m_v[1] += m.m_v[1];
    m_v[2] += m.m_v[2];
    m_v[3] += m.m_v[3];
    m_v[4] += m.m_v[4];
    m_v[5] += m.m_v[5];
    m_v[6] += m.m_v[6];
    m_v[7] += m.m_v[7];
    m_v[8] += m.m_v[8];
    m_v[9] += m.m_v[9];
    m_v[10] += m.m_v[10];
    m_v[11] += m.m_v[11];
    m_v[12] += m.m_v[12];
    m_v[13] += m.m_v[13];
    m_v[14] += m.m_v[14];
    m_v[15] += m.m_v[15];

    return *this;
}

inline
Matrix4 &Matrix4::operator-=(const Matrix4 &m) {
    m_v[0] -= m.m_v[0];
    m_v[1] -= m.m_v[1];
    m_v[2] -= m.m_v[2];
    m_v[3] -= m.m_v[3];
    m_v[4] -= m.m_v[4];
    m_v[5] -= m.m_v[5];
    m_v[6] -= m.m_v[6];
    m_v[7] -= m.m_v[7];
    m_v[8] -= m.m_v[8];
    m_v[9] -= m.m_v[9];
    m_v[10] -= m.m_v[10];
    m_v[11] -= m.m_v[11];
    m_v[12] -= m.m_v[12];
    m_v[13] -= m.m_v[13];
    m_v[14] -= m.m_v[14];
    m_v[15] -= m.m_v[15];

    return *this;
}

inline
Matrix4 &Matrix4::operator*=(Real r) {
    m_v[0] *= r;
    m_v[1] *= r;
    m_v[2] *= r;
    m_v[3] *= r;
    m_v[4] *= r;
    m_v[5] *= r;
    m_v[6] *= r;
    m_v[7] *= r;
    m_v[8] *= r;
    m_v[9] *= r;
    m_v[10] *= r;
    m_v[11] *= r;
    m_v[12] *= r;
    m_v[13] *= r;
    m_v[14] *= r;
    m_v[15] *= r;

    return *this;
}

inline
Matrix4 &Matrix4::operator/=(Real r) {
    Real inv_r = 1.0 / r;

    m_v[0] *= inv_r;
    m_v[1] *= inv_r;
    m_v[2] *= inv_r;
    m_v[3] *= inv_r;
    m_v[4] *= inv_r;
    m_v[5] *= inv_r;
    m_v[6] *= inv_r;
    m_v[7] *= inv_r;
    m_v[8] *= inv_r;
    m_v[9] *= inv_r;
    m_v[10] *= inv_r;
    m_v[11] *= inv_r;
    m_v[12] *= inv_r;
    m_v[13] *= inv_r;
    m_v[14] *= inv_r;
    m_v[15] *= inv_r;

    return *this;
}

inline
Matrix4 Matrix4::operator*(const Matrix4 &m) const {
    Matrix4 res;

    res.m_v[0] = m_v[0] * m.m_v[0] + m_v[1] * m.m_v[4] + m_v[2] * m.m_v[8] + m_v[3] * m.m_v[12];
    res.m_v[1] = m_v[0] * m.m_v[1] + m_v[1] * m.m_v[5] + m_v[2] * m.m_v[9] + m_v[3] * m.m_v[13];
    res.m_v[2] = m_v[0] * m.m_v[2] + m_v[1] * m.m_v[6] + m_v[2] * m.m_v[10] + m_v[3] * m.m_v[14];
    res.m_v[3] = m_v[0] * m.m_v[3] + m_v[1] * m.m_v[7] + m_v[2] * m.m_v[11] + m_v[3] * m.m_v[15];

    res.m_v[4] = m_v[4] * m.m_v[0] + m_v[5] * m.m_v[4] + m_v[6] * m.m_v[8] + m_v[7] * m.m_v[12];
    res.m_v[5] = m_v[4] * m.m_v[1] + m_v[5] * m.m_v[5] + m_v[6] * m.m_v[9] + m_v[7] * m.m_v[13];
    res.m_v[6] = m_v[4] * m.m_v[2] + m_v[5] * m.m_v[6] + m_v[6] * m.m_v[10] + m_v[7] * m.m_v[14];
    res.m_v[7] = m_v[4] * m.m_v[3] + m_v[5] * m.m_v[7] + m_v[6] * m.m_v[11] + m_v[7] * m.m_v[15];

    res.m_v[8] = m_v[8] * m.m_v[0] + m_v[9] * m.m_v[4] + m_v[10] * m.m_v[8] + m_v[11] * m.m_v[12];
    res.m_v[9] = m_v[8] * m.m_v[1] + m_v[9] * m.m_v[5] + m_v[10] * m.m_v[9] + m_v[11] * m.m_v[13];
    res.m_v[10] = m_v[8] * m.m_v[2] + m_v[9] * m.m_v[6] + m_v[10] * m.m_v[10] + m_v[11] * m.m_v[14];
    res.m_v[11] = m_v[8] * m.m_v[3] + m_v[9] * m.m_v[7] + m_v[10] * m.m_v[11] + m_v[11] * m.m_v[15];

    res.m_v[12] = m_v[12] * m.m_v[0] + m_v[13] * m.m_v[4] + m_v[14] * m.m_v[8] + m_v[15] * m.m_v[12];
    res.m_v[13] = m_v[12] * m.m_v[1] + m_v[13] * m.m_v[5] + m_v[14] * m.m_v[9] + m_v[15] * m.m_v[13];
    res.m_v[14] = m_v[12] * m.m_v[2] + m_v[13] * m.m_v[6] + m_v[14] * m.m_v[10] + m_v[15] * m.m_v[14];
    res.m_v[15] = m_v[12] * m.m_v[3] + m_v[13] * m.m_v[7] + m_v[14] * m.m_v[11] + m_v[15] * m.m_v[15];

    return res;
}

inline
Matrix4 &Matrix4::operator*=(const Matrix4 &m) {
    Real u, v, w;

    u = m_v[0] * m.m_v[0] + m_v[1] * m.m_v[4] + m_v[2] * m.m_v[8] + m_v[3] * m.m_v[12];
    v = m_v[0] * m.m_v[1] + m_v[1] * m.m_v[5] + m_v[2] * m.m_v[9] + m_v[3] * m.m_v[13];
    w = m_v[0] * m.m_v[2] + m_v[1] * m.m_v[6] + m_v[2] * m.m_v[10] + m_v[3] * m.m_v[14];

    m_v[3] = m_v[0] * m.m_v[3] + m_v[1] * m.m_v[7] + m_v[2] * m.m_v[11] + m_v[3] * m.m_v[15];
    m_v[0] = u;
    m_v[1] = v;
    m_v[2] = w;

    u = m_v[4] * m.m_v[0] + m_v[5] * m.m_v[4] + m_v[6] * m.m_v[8] + m_v[7] * m.m_v[12];
    v = m_v[4] * m.m_v[1] + m_v[5] * m.m_v[5] + m_v[6] * m.m_v[9] + m_v[7] * m.m_v[13];
    w = m_v[4] * m.m_v[2] + m_v[5] * m.m_v[6] + m_v[6] * m.m_v[10] + m_v[7] * m.m_v[14];
    m_v[7] = m_v[4] * m.m_v[3] + m_v[5] * m.m_v[7] + m_v[6] * m.m_v[11] + m_v[7] * m.m_v[15];

    m_v[4] = u;
    m_v[5] = v;
    m_v[6] = w;

    u = m_v[8] * m.m_v[0] + m_v[9] * m.m_v[4] + m_v[10] * m.m_v[8] + m_v[11] * m.m_v[12];
    v = m_v[8] * m.m_v[1] + m_v[9] * m.m_v[5] + m_v[10] * m.m_v[9] + m_v[11] * m.m_v[13];
    w = m_v[8] * m.m_v[2] + m_v[9] * m.m_v[6] + m_v[10] * m.m_v[10] + m_v[11] * m.m_v[14];
    m_v[11] = m_v[8] * m.m_v[3] + m_v[9] * m.m_v[7] + m_v[10] * m.m_v[11] + m_v[11] * m.m_v[15];

    m_v[8] = u;
    m_v[9] = v;
    m_v[10] = w;

    u = m_v[12] * m.m_v[0] + m_v[13] * m.m_v[4] + m_v[14] * m.m_v[8] + m_v[15] * m.m_v[12];
    v = m_v[12] * m.m_v[1] + m_v[13] * m.m_v[5] + m_v[14] * m.m_v[9] + m_v[15] * m.m_v[13];
    w = m_v[12] * m.m_v[2] + m_v[13] * m.m_v[6] + m_v[14] * m.m_v[10] + m_v[15] * m.m_v[14];
    m_v[15] = m_v[12] * m.m_v[3] + m_v[13] * m.m_v[7] + m_v[14] * m.m_v[11] + m_v[15] * m.m_v[15];

    m_v[12] = u;
    m_v[13] = v;
    m_v[14] = w;

    return *this;
}

/*inline
Vector3 Matrix4::operator*(const Vector3 &v) const {
    // We treat v as the 4-dimensional vector (v.x, v.y, v.z, 1.0), and by
    // multiplying this vector by the 4x4 matrix, we obtain a 4-dimensional
    // vector (res.x, res.y, res.z, w). Finally, we divide the resulting
    // vector res by w to obtain a 3-dimensional vector.

    Real w = m_v[12] * v.m_x + m_v[13] * v.m_y + m_v[14] * v.m_z + m_v[15];

    assert(w != 0.0);

    Vector3 res;

    res.m_x = m_v[0] * v.m_x + m_v[1] * v.m_y + m_v[2] * v.m_z + m_v[3];
    res.m_y = m_v[4] * v.m_x + m_v[5] * v.m_y + m_v[6] * v.m_z + m_v[7];
    res.m_z = m_v[8] * v.m_x + m_v[9] * v.m_y + m_v[10] * v.m_z + m_v[11];

    return res / w;
}*/

inline
Vector4 Matrix4::operator*(const Vector4 &v) const {
    return Vector4(
        m_v[0] * v.m_x + m_v[1] * v.m_y + m_v[2] * v.m_z + m_v[3] * v.m_w,
        m_v[4] * v.m_x + m_v[5] * v.m_y + m_v[6] * v.m_z + m_v[7] * v.m_w,
        m_v[8] * v.m_x + m_v[9] * v.m_y + m_v[10] * v.m_z + m_v[11] * v.m_w,
        m_v[12] * v.m_x + m_v[13] * v.m_y + m_v[14] * v.m_z + m_v[15] * v.m_w);
}

inline
Vector3 Matrix4::operator*(const Vector3 &v) const {
    return Vector3(
        m_v[0] * v.m_x + m_v[1] * v.m_y + m_v[2] * v.m_z,
        m_v[4] * v.m_x + m_v[5] * v.m_y + m_v[6] * v.m_z,
        m_v[8] * v.m_x + m_v[9] * v.m_y + m_v[10] * v.m_z);
}

inline
Point3 Matrix4::operator*(const Point3 &v) const {
    return Point3(
        m_v[0] * v.m_x + m_v[1] * v.m_y + m_v[2] * v.m_z + m_v[3],
        m_v[4] * v.m_x + m_v[5] * v.m_y + m_v[6] * v.m_z + m_v[7],
        m_v[8] * v.m_x + m_v[9] * v.m_y + m_v[10] * v.m_z + m_v[11]);
}

inline
Real &Matrix4::operator[](int i) {
    assert(i >= 0 && i < 16);

    return m_v[i];
}

inline
const Real &Matrix4::operator[](int i) const {
    assert(i >= 0 && i < 16);

    return m_v[i];
}

inline
Real &Matrix4::operator()(int row, int column) {
    assert(row >= 0 && row < 4);
    assert(column >= 0 && column < 4);

    return m_v[row + row + row + column];
}

inline
const Real &Matrix4::operator()(int row, int column) const {
    assert(row >= 0 && row < 4);
    assert(column >= 0 && column < 4);

    return m_v[row + row + row + column];
}

inline
Matrix4 Matrix4::Transpose() const {
    Matrix4 res;

    res.m_v[0] = m_v[0];
    res.m_v[1] = m_v[4];
    res.m_v[2] = m_v[8];
    res.m_v[3] = m_v[12];
    res.m_v[4] = m_v[1];
    res.m_v[5] = m_v[5];
    res.m_v[6] = m_v[9];
    res.m_v[7] = m_v[13];
    res.m_v[8] = m_v[2];
    res.m_v[9] = m_v[6];
    res.m_v[10] = m_v[10];
    res.m_v[11] = m_v[14];
    res.m_v[12] = m_v[3];
    res.m_v[13] = m_v[7];
    res.m_v[14] = m_v[11];
    res.m_v[15] = m_v[15];

    return res;
}

#define DET(c0, c1, c2, c3, c4, c5, c6, c7, c8) \
    (c0 * (c4 * c8 - c5 * c7) + \
     c1 * (c5 * c6 - c3 * c8) + \
     c2 * (c3 * c7 - c4 * c6))

inline
Matrix4 Matrix4::Inverse() const {
    // Compute the matrix inverse using Cramer's rule.

    const Real c0 = m_v[0];
    const Real c1 = m_v[1];
    const Real c2 = m_v[2];
    const Real c3 = m_v[3];
    const Real c4 = m_v[4];
    const Real c5 = m_v[5]; 
    const Real c6 = m_v[6];
    const Real c7 = m_v[7];
    const Real c8 = m_v[8];
    const Real c9 = m_v[9];
    const Real c10 = m_v[10];
    const Real c11 = m_v[11];
    const Real c12 = m_v[12];
    const Real c13 = m_v[13]; 
    const Real c14 = m_v[14];
    const Real c15 = m_v[15];

    const Real u = DET(c5, c9, c13, c6, c10, c14, c7, c11, c15);
    const Real v = - DET(c4, c8, c12, c6, c10, c14, c7, c11, c15);
    const Real w = DET(c4, c8, c12, c5, c9, c13, c7, c11, c15);
    const Real x = - DET(c4, c8, c12, c5, c9, c13, c6, c10, c14);

    Matrix4 res;

    // Compute cofactors.
    res.m_v[0] = u;
    res.m_v[1] = - DET(c1, c9, c13, c2, c10, c14, c3, c11, c15);
    res.m_v[2] = DET(c1, c5, c13, c2, c6, c14, c3, c7, c15);
    res.m_v[3] = - DET(c1, c5, c9, c2, c6, c10, c3, c7, c11);
    res.m_v[4] = v;
    res.m_v[5] = DET(c0, c8, c12, c2, c10, c14, c3, c11, c15);
    res.m_v[6] = - DET(c0, c4, c12, c2, c6, c14, c3, c7, c15);
    res.m_v[7] = DET(c0, c4, c8, c2, c6, c10, c3, c7, c11);
    res.m_v[8] = w;
    res.m_v[9] = - DET(c0, c8, c12, c1, c9, c13, c3, c11, c15);
    res.m_v[10] = DET(c0, c4, c12, c1, c5, c13, c3, c7, c15);
    res.m_v[11] = - DET(c0, c4, c8, c1, c5, c9, c3, c7, c11);
    res.m_v[12] = x;        
    res.m_v[13] = DET(c0, c8, c12, c1, c9, c13, c2, c10, c14);
    res.m_v[14] = - DET(c0, c4, c12, c1, c5, c13, c2, c6, c14);
    res.m_v[15] = DET(c0, c4, c8, c1, c5, c9, c2, c6, c10);

    // Compute matrix determinant.
    const Real det = m_v[0] * u + m_v[1] * v + m_v[2] * w + m_v[3] * x;

    //!\todo Throw an exception if the matrix to be inverted
    //!\todo is singular (i.e. if fabs(det) <= EPSILON).
    assert(det != 0.0);

    // Multiply the obtained matrix by the reciprocal of the
    // input matrix determinant.
    const Real inv_det = 1.0 / det;
    res.m_v[0] *= inv_det;
    res.m_v[1] *= inv_det;
    res.m_v[2] *= inv_det;
    res.m_v[3] *= inv_det;
    res.m_v[4] *= inv_det;
    res.m_v[5] *= inv_det;
    res.m_v[6] *= inv_det;
    res.m_v[7] *= inv_det;
    res.m_v[8] *= inv_det;
    res.m_v[9] *= inv_det;
    res.m_v[10] *= inv_det;
    res.m_v[11] *= inv_det;
    res.m_v[12] *= inv_det;
    res.m_v[13] *= inv_det;
    res.m_v[14] *= inv_det;
    res.m_v[15] *= inv_det;

    return res;
}

inline
Real Matrix4::Determinant() const {
    const Real c4 = m_v[4];
    const Real c5 = m_v[5]; 
    const Real c6 = m_v[6];
    const Real c7 = m_v[7];
    const Real c8 = m_v[8];
    const Real c9 = m_v[9]; 
    const Real c10 = m_v[10];
    const Real c11 = m_v[11];
    const Real c12 = m_v[12];
    const Real c13 = m_v[13]; 
    const Real c14 = m_v[14];
    const Real c15 = m_v[15];

    return
          m_v[0] * DET(c5, c9, c13, c6, c10, c14, c7, c11, c15)
        - m_v[1] * DET(c4, c8, c12, c6, c10, c14, c7, c11, c15)
        + m_v[2] * DET(c4, c8, c12, c5, c9, c13, c7, c11, c15)
        - m_v[3] * DET(c4, c8, c12, c5, c9, c13, c6, c10, c14);
}

#undef DET

inline
Real Matrix4::Trace() const {
    return m_v[0] + m_v[5] + m_v[10] + m_v[15];
}

inline
Matrix4 operator*(Real r, const Matrix4 &m) {
    return m * r;
}

---