matrix4x4.h

Go to the documentation of this file.

// Copyright © 2008-2021 Pioneer Developers. See AUTHORS.txt for details
// Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
 
#ifndef _MATRIX4X4_H
#define _MATRIX4X4_H
 
#include "matrix3x3.h"
#include "vector3.h"
#include <math.h>
#include <stdio.h>
#include <type_traits>
 
template <typename T>
class matrix4x4 {
private:
    T cell[16];
    using other_float_t = typename std::conditional<std::is_same<T, float>::value, double, float>::type;
 
public:
    matrix4x4() {}
    matrix4x4(T val)
    {
        cell[0] = cell[1] = cell[2] = cell[3] = cell[4] = cell[5] = cell[6] =
            cell[7] = cell[8] = cell[9] = cell[10] = cell[11] = cell[12] = cell[13] =
                cell[14] = cell[15] = val;
    }
    matrix4x4(const T *vals)
    {
        memcpy(cell, vals, sizeof(T) * 16);
    }
    matrix4x4(const matrix3x3<T> &m)
    {
        LoadFrom3x3Matrix(m.Data());
    }
    matrix4x4(const matrix3x3<T> &m, const vector3<T> &v)
    {
        LoadFrom3x3Matrix(m.Data());
        SetTranslate(v);
    }
    explicit matrix4x4(const matrix4x4<other_float_t> &m)
    {
        for (int i = 0; i < 16; i++)
            cell[i] = T(m[i]);
    }
 
    void SetTranslate(const vector3<T> &v)
    {
        cell[12] = v.x;
        cell[13] = v.y;
        cell[14] = v.z;
    }
    vector3<T> GetTranslate() const { return vector3<T>(cell[12], cell[13], cell[14]); }
    void SetRotationOnly(const matrix4x4 &m)
    {
        for (int i = 0; i < 12; i++)
            cell[i] = m.cell[i];
    }
    matrix3x3<T> GetOrient() const
    {
        matrix3x3<T> m;
        m[0] = cell[0];
        m[1] = cell[4];
        m[2] = cell[8];
        m[3] = cell[1];
        m[4] = cell[5];
        m[5] = cell[9];
        m[6] = cell[2];
        m[7] = cell[6];
        m[8] = cell[10];
        return m;
    }
    // row-major 3x3 matrix
    void LoadFrom3x3Matrix(const T *r)
    {
        cell[0] = r[0];
        cell[4] = r[1];
        cell[8] = r[2];
        cell[12] = 0;
        cell[1] = r[3];
        cell[5] = r[4];
        cell[9] = r[5];
        cell[13] = 0;
        cell[2] = r[6];
        cell[6] = r[7];
        cell[10] = r[8];
        cell[14] = 0;
        cell[3] = 0;
        cell[7] = 0;
        cell[11] = 0;
        cell[15] = 1;
    }
    // row-major
    void SaveTo3x3Matrix(T *r) const
    {
        r[0] = cell[0];
        r[1] = cell[4];
        r[2] = cell[8];
        r[3] = cell[1];
        r[4] = cell[5];
        r[5] = cell[9];
        r[6] = cell[2];
        r[7] = cell[6];
        r[8] = cell[10];
    }
    static matrix4x4 Identity()
    {
        matrix4x4 m = matrix4x4(0.0);
        m.cell[0] = m.cell[5] = m.cell[10] = m.cell[15] = 1.0f;
        return m;
    }
    //glscale equivalent
    void Scale(T x, T y, T z)
    {
        *this = (*this) * ScaleMatrix(x, y, z);
    }
    void Scale(T s)
    {
        *this = (*this) * ScaleMatrix(s, s, s);
    }
    static matrix4x4 ScaleMatrix(T x, T y, T z)
    {
        matrix4x4 m;
        m[0] = x;
        m[1] = m[2] = m[3] = 0;
        m[5] = y;
        m[4] = m[6] = m[7] = 0;
        m[10] = z;
        m[8] = m[9] = m[11] = 0;
        m[12] = m[13] = m[14] = 0;
        m[15] = 1;
        return m;
    }
    static matrix4x4 ScaleMatrix(T scale)
    {
        matrix4x4 m;
        m[0] = scale;
        m[1] = m[2] = m[3] = 0;
        m[5] = scale;
        m[4] = m[6] = m[7] = 0;
        m[10] = scale;
        m[8] = m[9] = m[11] = 0;
        m[12] = m[13] = m[14] = 0;
        m[15] = 1;
        return m;
    }
    static matrix4x4 MakeRotMatrix(const vector3<T> &rx, const vector3<T> &ry, const vector3<T> &rz)
    {
        matrix4x4 m;
        m[0] = rx.x;
        m[4] = rx.y;
        m[8] = rx.z;
        m[12] = 0;
        m[1] = ry.x;
        m[5] = ry.y;
        m[9] = ry.z;
        m[13] = 0;
        m[2] = rz.x;
        m[6] = rz.y;
        m[10] = rz.z;
        m[14] = 0;
        m[3] = 0;
        m[7] = 0;
        m[11] = 0;
        m[15] = 1;
        return m;
    }
    static matrix4x4 MakeInvRotMatrix(const vector3<T> &rx, const vector3<T> &ry, const vector3<T> &rz)
    {
        matrix4x4 m;
        m[0] = rx.x;
        m[4] = ry.x;
        m[8] = rz.x;
        m[12] = 0;
        m[1] = rx.y;
        m[5] = ry.y;
        m[9] = rz.y;
        m[13] = 0;
        m[2] = rx.z;
        m[6] = ry.z;
        m[10] = rz.z;
        m[14] = 0;
        m[3] = 0;
        m[7] = 0;
        m[11] = 0;
        m[15] = 1;
        return m;
    }
 
    // Matrix Construction Functions
    // NOTE: all matrix functions here are optimized for reverse-Z depth buffers.
    // Compared to "standard" DirectX or OpenGL matricies they invert the Z value
    // so it ranges from 1.0 at the near plane to 0.0 at the far plane.
 
    // Construct a perspective projection matrix based on arbitrary left/right/top/bottom
    // plane positions.
    // This method is slower than the others, but supports view frustrums that are not
    // aligned with the Z-axis. Unless you know what you're doing, you shouldn't use this.
    //
    // @param left - the minimum x-value of the view volume at the near plane
    // @param right - the maxiumum x-value of the view volume at the near plane
    // @param bottom - the maxiumum y-value of the view volume at the near plane
    // @param top - the maxiumum y-value of the view volume at the near plane
    // @param znear - the near clipping plane
    // @param zfar - the far clipping plane
    static matrix4x4 FrustumMatrix(T left, T right, T bottom, T top, T znear, T zfar)
    {
        assert((znear > T(0)) && (zfar > T(0)));
        // these expressions come from the documentation for glFrustum
        const T sx = (T(2) * znear) / (right - left);
        const T sy = (T(2) * znear) / (top - bottom);
        const T A = (right + left) / (right - left);
        const T B = (top + bottom) / (top - bottom);
        const T C = (zfar) / (zfar - znear) - 1;
        const T D = (zfar * znear) / (zfar - znear);
        matrix4x4 m;
 
        // http://glprogramming.com/red/appendixf.html
        // OpenGL 'Red Book' on Perspective Projection
        // Presented here in row-major notation (because that's what matrix4x4f uses internally)
        T perspective[16] = {
            sx, 0, 0, 0,
            0, sy, 0, 0,
            A, B, C, -1,
            0, 0, D, 0
        };
        return matrix4x4(&perspective[0]);
    }
 
    // Construct a perspective projection matrix based field of view and aspect ratio.
    // This method is the optimized case when you know your screen aspect ratio and
    // field of view and aren't interested in fancy math. Use this function or
    // InfinitePerspectiveMatrix if at all possible.
    //
    // @param fovR - the camera FOV in radians
    // @param aspect - the aspect ratio (width / height) of the viewport
    // @param znear - the near clipping plane
    // @param zfar - the far clipping plane
    // @param fovX - whether the field of view is horizontal or vertical (default)
    static matrix4x4 PerspectiveMatrix(T fovR, T aspect, T znear, T zfar, bool fovX = false)
    {
        assert((znear > T(0)) && (zfar > znear));
 
        const T e = 1 / tan(fovR / T(2));
        const T x = fovX ? e : e / aspect;
        const T y = fovX ? e / aspect : e;
        const T z = (znear) / (zfar - znear);
        const T w = (zfar * znear) / (zfar - znear);
 
        // Based on: http://www.terathon.com/gdc07_lengyel.pdf
        // Unlike gluProject / FrustumMatrix, this projection matrix can only be
        //  symmetric about the Z axis.
        // This is what you want in 99% of cases, and simplifies the math a good deal.
        T perspective[16] = {
            x, 0, 0, 0,
            0, y, 0, 0,
            0, 0, z, -1,
            0, 0, w, 0
        };
 
        return matrix4x4(&perspective[0]);
    }
 
    // Construct an infinite far-plane perspective projection matrix.
    // Unless you specifically want to clip objects beyond a specific distance,
    // this projection will work for any object at any distance.
    //
    // @param fovR - the camera FOV in radians
    // @param aspect - the aspect ratio (width / height) of the viewport
    // @param znear - the near clipping plane
    // @param fovX - whether the field of view is horizontal or vertical (default)
    static matrix4x4 InfinitePerspectiveMatrix(T fovR, T aspect, T znear, bool fovX = false)
    {
        assert(znear > T(0));
 
        const T e = 1 / tan(fovR / T(2));
        const T x = fovX ? e : e / aspect;
        const T y = fovX ? e / aspect : e;
        const T w = znear;
 
        // Based on: http://dev.theomader.com/depth-precision/
        // An 'infinite far-plane' projection matrix. There is no concept of a zFar value,
        // and it can handle everything up to and including homogeneous coordinates with w=0.
        T perspective[16] = {
            x, 0, 0, 0,
            0, y, 0, 0,
            0, 0, 0, -1,
            0, 0, w, 0
        };
 
        return matrix4x4(&perspective[0]);
    }
 
    // set a orthographic frustum with 6 params similar to glOrtho()
    // (left, right, bottom, top, near, far)
    //
    // Derived from:
    // [1] https://docs.microsoft.com/en-us/windows/win32/direct3d9/d3dxmatrixorthorh
    // [2] https://thxforthefish.com/posts/reverse_z/
    //
    // Specifically, the `tz` and `c` terms are based on a right-handed DirectX-style
    // (Z=0..1) matrix [1], multiplied by the given reversing matrix [2].
    // If looking at the sources, keep in mind that [1] is presented in row-major order
    // and [2] is presented in column-major order, and thus multiplication occurs in
    // column-column fashion.
    static matrix4x4 OrthoFrustum(T left, T right, T bottom, T top, T znear, T zfar)
    {
        assert((znear >= T(-1)) && (zfar > T(0)));
        T a = T(2) / (right - left);
        T b = T(2) / (top - bottom);
        T c = T(1) / (zfar - znear);
 
        T tx = (right + left) / (left - right);
        T ty = (top + bottom) / (bottom - top);
        T tz = (zfar) / (zfar - znear);
 
        T ortho[16] = {
            a, 0, 0, 0,
            0, b, 0, 0,
            0, 0, c, 0,
            tx, ty, tz, 1
        };
        matrix4x4 m(&ortho[0]);
        return m;
    }
 
    // Optimized form that takes width/height and near/far planes
    static matrix4x4 OrthoMatrix(T width, T height, T znear, T zfar)
    {
        assert((znear >= T(-1)) && (zfar > T(0)));
        T a = T(2) / width;
        T b = T(2) / height;
        T c = T(1) / (zfar - znear);
 
        T tz = (zfar) / (zfar - znear);
 
        T ortho[16] = {
            a, 0, 0, 0,
            0, b, 0, 0,
            0, 0, c, 0,
            0, 0, tz, 1
        };
        matrix4x4 m(&ortho[0]);
        return m;
    }
 
    //glRotate equivalent (except radians instead of degrees)
    void Rotate(T ang, T x, T y, T z)
    {
        *this = (*this) * RotateMatrix(ang, x, y, z);
    }
    // (x,y,z) must be normalized
    static matrix4x4 RotateMatrix(T ang, T x, T y, T z)
    {
        matrix4x4 m;
        T c = cos(ang);
        T s = sin(ang);
        m[0] = x * x * (1 - c) + c;
        m[1] = y * x * (1 - c) + z * s;
        m[2] = x * z * (1 - c) - y * s;
        m[3] = 0;
        m[4] = x * y * (1 - c) - z * s;
        m[5] = y * y * (1 - c) + c;
        m[6] = y * z * (1 - c) + x * s;
        m[7] = 0;
        m[8] = x * z * (1 - c) + y * s;
        m[9] = y * z * (1 - c) - x * s;
        m[10] = z * z * (1 - c) + c;
        m[11] = 0;
        m[12] = 0;
        m[13] = 0;
        m[14] = 0;
        m[15] = 1;
        return m;
    }
    void RotateZ(T radians) { *this = (*this) * RotateZMatrix(radians); }
    void RotateY(T radians) { *this = (*this) * RotateYMatrix(radians); }
    void RotateX(T radians) { *this = (*this) * RotateXMatrix(radians); }
    static matrix4x4 RotateXMatrix(T radians)
    {
        matrix4x4 m;
        T cos_r = cosf(float(radians));
        T sin_r = sinf(float(radians));
        m[0] = 1.0f;
        m[1] = 0;
        m[2] = 0;
        m[3] = 0;
 
        m[4] = 0;
        m[5] = cos_r;
        m[6] = -sin_r;
        m[7] = 0;
 
        m[8] = 0;
        m[9] = sin_r;
        m[10] = cos_r;
        m[11] = 0;
 
        m[12] = 0;
        m[13] = 0;
        m[14] = 0;
        m[15] = 1.0f;
        return m;
    }
    static matrix4x4 RotateYMatrix(T radians)
    {
        matrix4x4 m;
        T cos_r = cosf(float(radians));
        T sin_r = sinf(float(radians));
        m[0] = cos_r;
        m[1] = 0;
        m[2] = sin_r;
        m[3] = 0;
 
        m[4] = 0;
        m[5] = 1;
        m[6] = 0;
        m[7] = 0;
 
        m[8] = -sin_r;
        m[9] = 0;
        m[10] = cos_r;
        m[11] = 0;
 
        m[12] = 0;
        m[13] = 0;
        m[14] = 0;
        m[15] = 1.0f;
        return m;
    }
    static matrix4x4 RotateZMatrix(T radians)
    {
        matrix4x4 m;
        T cos_r = cosf(float(radians));
        T sin_r = sinf(float(radians));
        m[0] = cos_r;
        m[1] = -sin_r;
        m[2] = 0;
        m[3] = 0;
 
        m[4] = sin_r;
        m[5] = cos_r;
        m[6] = 0;
        m[7] = 0;
 
        m[8] = 0;
        m[9] = 0;
        m[10] = 1.0f;
        m[11] = 0;
 
        m[12] = 0;
        m[13] = 0;
        m[14] = 0;
        m[15] = 1.0f;
        return m;
    }
    void Renormalize()
    {
        vector3<T> x(cell[0], cell[4], cell[8]);
        vector3<T> y(cell[1], cell[5], cell[9]);
        vector3<T> z(cell[2], cell[6], cell[10]);
        x = x.Normalized();
        z = x.Cross(y).Normalized();
        y = z.Cross(x).Normalized();
        cell[0] = x.x;
        cell[4] = x.y;
        cell[8] = x.z;
        cell[1] = y.x;
        cell[5] = y.y;
        cell[9] = y.z;
        cell[2] = z.x;
        cell[6] = z.y;
        cell[10] = z.z;
    }
    void ClearToRotOnly()
    {
        cell[12] = 0;
        cell[13] = 0;
        cell[14] = 0;
    }
    T &operator[](const size_t i) { return cell[i]; }
    const T &operator[](const size_t i) const { return cell[i]; }
    const T *Data() const { return cell; }
    T *Data() { return cell; }
    friend matrix4x4 operator+(const matrix4x4 &a, const matrix4x4 &b)
    {
        matrix4x4 m;
        for (int i = 0; i < 16; i++)
            m.cell[i] = a.cell[i] + b.cell[i];
        return m;
    }
    friend matrix4x4 operator-(const matrix4x4 &a, const matrix4x4 &b)
    {
        matrix4x4 m;
        for (int i = 0; i < 16; i++)
            m.cell[i] = a.cell[i] - b.cell[i];
        return m;
    }
    friend matrix4x4 operator-(const matrix4x4 &a)
    {
        matrix4x4 m;
        for (int i = 0; i < 16; ++i) {
            m.cell[i] = -a.cell[i];
        }
        return m;
    }
    friend matrix4x4 operator*(const matrix4x4 &a, const matrix4x4 &b)
    {
        matrix4x4 m;
        m.cell[0] = a.cell[0] * b.cell[0] + a.cell[4] * b.cell[1] + a.cell[8] * b.cell[2] + a.cell[12] * b.cell[3];
        m.cell[1] = a.cell[1] * b.cell[0] + a.cell[5] * b.cell[1] + a.cell[9] * b.cell[2] + a.cell[13] * b.cell[3];
        m.cell[2] = a.cell[2] * b.cell[0] + a.cell[6] * b.cell[1] + a.cell[10] * b.cell[2] + a.cell[14] * b.cell[3];
        m.cell[3] = a.cell[3] * b.cell[0] + a.cell[7] * b.cell[1] + a.cell[11] * b.cell[2] + a.cell[15] * b.cell[3];
 
        m.cell[4] = a.cell[0] * b.cell[4] + a.cell[4] * b.cell[5] + a.cell[8] * b.cell[6] + a.cell[12] * b.cell[7];
        m.cell[5] = a.cell[1] * b.cell[4] + a.cell[5] * b.cell[5] + a.cell[9] * b.cell[6] + a.cell[13] * b.cell[7];
        m.cell[6] = a.cell[2] * b.cell[4] + a.cell[6] * b.cell[5] + a.cell[10] * b.cell[6] + a.cell[14] * b.cell[7];
        m.cell[7] = a.cell[3] * b.cell[4] + a.cell[7] * b.cell[5] + a.cell[11] * b.cell[6] + a.cell[15] * b.cell[7];
 
        m.cell[8] = a.cell[0] * b.cell[8] + a.cell[4] * b.cell[9] + a.cell[8] * b.cell[10] + a.cell[12] * b.cell[11];
        m.cell[9] = a.cell[1] * b.cell[8] + a.cell[5] * b.cell[9] + a.cell[9] * b.cell[10] + a.cell[13] * b.cell[11];
        m.cell[10] = a.cell[2] * b.cell[8] + a.cell[6] * b.cell[9] + a.cell[10] * b.cell[10] + a.cell[14] * b.cell[11];
        m.cell[11] = a.cell[3] * b.cell[8] + a.cell[7] * b.cell[9] + a.cell[11] * b.cell[10] + a.cell[15] * b.cell[11];
 
        m.cell[12] = a.cell[0] * b.cell[12] + a.cell[4] * b.cell[13] + a.cell[8] * b.cell[14] + a.cell[12] * b.cell[15];
        m.cell[13] = a.cell[1] * b.cell[12] + a.cell[5] * b.cell[13] + a.cell[9] * b.cell[14] + a.cell[13] * b.cell[15];
        m.cell[14] = a.cell[2] * b.cell[12] + a.cell[6] * b.cell[13] + a.cell[10] * b.cell[14] + a.cell[14] * b.cell[15];
        m.cell[15] = a.cell[3] * b.cell[12] + a.cell[7] * b.cell[13] + a.cell[11] * b.cell[14] + a.cell[15] * b.cell[15];
        return m;
    }
    friend vector3<T> operator*(const matrix4x4 &a, const vector3<T> &v)
    {
        vector3<T> out;
        out.x = a.cell[0] * v.x + a.cell[4] * v.y + a.cell[8] * v.z + a.cell[12];
        out.y = a.cell[1] * v.x + a.cell[5] * v.y + a.cell[9] * v.z + a.cell[13];
        out.z = a.cell[2] * v.x + a.cell[6] * v.y + a.cell[10] * v.z + a.cell[14];
        return out;
    }
    // scam for doing a transpose operation
    friend vector3<T> operator*(const vector3<T> &v, const matrix4x4 &a)
    {
        vector3<T> out;
        out.x = a.cell[0] * v.x + a.cell[1] * v.y + a.cell[2] * v.z;
        out.y = a.cell[4] * v.x + a.cell[5] * v.y + a.cell[6] * v.z;
        out.z = a.cell[8] * v.x + a.cell[9] * v.y + a.cell[10] * v.z;
        return out;
    }
 
    // Transform a vector by the affine inverse of a matrix4x4
    // internally this does a transpose operation, and thus only works on
    // Euclidean (translation + rotation) matricies.
    vector3<T> InvTransform(const vector3<T> &inVec)
    {
        // Formula derivation from songho (https://songho.ca/opengl):
        //
        // M = [ R | T ]
        //     [ --+-- ]    (R denotes 3x3 rotation/reflection matrix)
        //     [ 0 | 1 ]    (T denotes 1x3 translation matrix)
        //
        // y = M*x  ->  y = R*x + T  ->  x = R^-1*(y - T)  ->  x = R^T*y - R^T*T
        // (R is orthogonal,  R^-1 = R^T)
 
        // thanks to https://stackoverflow.com/a/2625420
        // "Depending on your situation, it may be faster to compute the result of
        // inv(A) * x instead of actually forming inv(A)..."
        //
        // inv(A) * [x] = [ inv(M) * (x - b) ]
        //          [1] = [        1         ]
        //
 
        vector3<T> v = inVec - GetTranslate();
        vector3<T> out;
        out.x = cell[0] * v.x + cell[1] * v.y + cell[2] * v.z;
        out.y = cell[4] * v.x + cell[5] * v.y + cell[6] * v.z;
        out.z = cell[8] * v.x + cell[9] * v.y + cell[10] * v.z;
        return out;
    }
 
    friend matrix4x4 operator*(const matrix4x4 &a, T v)
    {
        matrix4x4 m;
        for (int i = 0; i < 16; i++)
            m[i] = a.cell[i] * v;
        return m;
    }
    friend matrix4x4 operator*(T v, const matrix4x4 &a)
    {
        return (a * v);
    }
    vector3<T> ApplyRotationOnly(const vector3<T> &v) const
    {
        vector3<T> out;
        out.x = cell[0] * v.x + cell[4] * v.y + cell[8] * v.z;
        out.y = cell[1] * v.x + cell[5] * v.y + cell[9] * v.z;
        out.z = cell[2] * v.x + cell[6] * v.y + cell[10] * v.z;
        return out;
    }
    //gltranslate equivalent
    void Translate(const vector3<T> &t)
    {
        Translate(t.x, t.y, t.z);
    }
    void Translate(T x, T y, T z)
    {
        matrix4x4 m = Identity();
        m[12] = x;
        m[13] = y;
        m[14] = z;
        *this = (*this) * m;
    }
    static matrix4x4 Translation(const vector3<T> &v)
    {
        return Translation(v.x, v.y, v.z);
    }
    static matrix4x4 Translation(T x, T y, T z)
    {
        matrix4x4 m = Identity();
        m[12] = x;
        m[13] = y;
        m[14] = z;
        return m;
    }
    matrix4x4 Inverse() const
    {
        matrix4x4 m;
        // this only works for matrices containing only rotation and transform
        m[0] = cell[0];
        m[1] = cell[4];
        m[2] = cell[8];
        m[4] = cell[1];
        m[5] = cell[5];
        m[6] = cell[9];
        m[8] = cell[2];
        m[9] = cell[6];
        m[10] = cell[10];
        m[12] = -(cell[0] * cell[12] + cell[1] * cell[13] + cell[2] * cell[14]);
        m[13] = -(cell[4] * cell[12] + cell[5] * cell[13] + cell[6] * cell[14]);
        m[14] = -(cell[8] * cell[12] + cell[9] * cell[13] + cell[10] * cell[14]);
        m[3] = m[7] = m[11] = 0;
        m[15] = 1.0f;
 
        return m;
    }
    matrix4x4 Transpose() const
    {
        matrix4x4 m;
        m[0] = cell[0];
        m[1] = cell[4];
        m[2] = cell[8];
        m[3] = cell[12];
        m[4] = cell[1];
        m[5] = cell[5];
        m[6] = cell[9];
        m[7] = cell[13];
        m[8] = cell[2];
        m[9] = cell[6];
        m[10] = cell[10];
        m[11] = cell[14];
        m[12] = cell[3];
        m[13] = cell[7];
        m[14] = cell[11];
        m[15] = cell[15];
        return m;
    }
    void Print() const
    {
        for (int i = 0; i < 4; i++) {
            printf("%.12f %.12f %.12f %.12f\n", cell[i], cell[i + 4], cell[i + 8], cell[i + 12]);
        }
        printf("\n");
    }
 
    //convenience accessors for getting right/up/back vectors
    //from rotation matrices
    vector3<T> Right() const
    {
        return vector3<T>(cell[0], cell[4], cell[8]);
    }
 
    vector3<T> Up() const
    {
        return vector3<T>(cell[1], cell[5], cell[9]);
    }
 
    vector3<T> Back() const
    {
        return vector3<T>(cell[2], cell[6], cell[10]);
    }
};
 
typedef matrix4x4<float> matrix4x4f;
typedef matrix4x4<double> matrix4x4d;
 
static const matrix4x4f matrix4x4fIdentity(matrix4x4f::Identity());
static const matrix4x4d matrix4x4dIdentity(matrix4x4d::Identity());
 
#endif /* _MATRIX4X4_H */