cRotationMatrix就其特殊目的來說是稱職的，但也正因為如此，它的廣泛應用受到了限制。cMatrix4x3類是一個更加一般化的矩陣，它被用來處理更加復雜的變換。這個矩陣類保存了一個一般仿射變換矩陣。旋轉、縮放、鏡像、投影和平移變換它都支持，該矩陣還能求逆和組合。

因此，cMatrix4x3類的語義和cRotationMatrix類完全不同。cRotationMatrix僅應用于特殊的物體空間和慣性空間，而cMatrix4x3有更一般的應用，所以我們使用更一般化的術語"源"和"目標"坐標空間。和cRotationMatrix不一樣，它的變換方向是在矩陣創建時指定的，之后點只能向那個方向（源到目標）變換。如果要向相反的方向變換，須先計算逆矩陣。

這里使用線性代數的乘法記法，operator*()被同時用來變換點和組合矩陣。因為我們的約定是行向量不是列向量，變換的順序和讀句子一樣，從左向右。

Matrix4x3.h:

（1）identity()將矩陣設為單位矩陣。

（2）zero_translation()通過將最后一行設為[0, 0, 0]來取消矩陣的平移部分，線性變換部分(3x3部分)不受影響。

（3）set_translation() 將矩陣的平移部分設為指定值，不改變3x3部分。setup_translation()設置矩陣來執行平移，上面的3x3部分設為單位矩陣，平移部分設為指定向量。

（4）setup_local_to_parent()創建一個矩陣能將點從局部坐標空間變換到父坐標空間，需要給出局部坐標空間在父坐標空間中的位置和方向。最常用到該方法的可能是將點從物體坐標空間變換到世界坐標空間的時候。局部坐標空間的方位可以用歐拉角或旋轉矩陣定義，用旋轉矩陣更快一些，因為它沒有實數運算，只有矩陣元素的復制。setup_parent_to_local()設置用來執行相反變換的矩陣。

（5）setup_rotate()的兩個重載方法都創建一個繞軸旋轉的矩陣。如果軸是坐標軸，將使用一個數字來代表坐標軸。繞任意軸旋轉時，用第二個版本的setup_rotate()，它用一個單位向量代表旋轉軸。

（6）from_quat()將一個四元數轉換到矩陣形式，矩陣的平移部分為0。

（7）setup_scale()創建一個矩陣執行沿坐標軸的均勻或非均勻縮放。輸入向量包含沿x、y、z軸的縮放因子。對均勻縮放，使用一個每個軸的值都相同的向量。

（8）setup_scale_along_axis()創建一個矩陣執行沿任意方向縮放。這個縮放發生在一個穿過原點的平面上----該平面垂直于向量參數，向量是單位化的。

（9）setup_shear()創建一個切變矩陣。

（10）setup_project()創建一個投影矩陣，該矩陣向穿過原點且垂直于給定向量的平面投影。

（11）setup_reflect()創建一個沿平面鏡像的矩陣。第一個版本中，坐標軸用一個數字指定，平面不必穿過原點。第二個版本中指定任意的法向量，且平面必須穿過原點。（對于沿任意不穿過原點的平面鏡像，必須將該矩陣和適當的變換矩陣連接。）

determinant()函數計算矩陣的行列式。實際只使用了3x3部分，如果假設最后一列總為[0, 0, 0, 1]^T，那么最后一行（平移部分）會被最后一列的0消去。

inverse()計算并返回矩陣的逆，理論上不可能對4x3矩陣求逆，只有方陣才能求逆。所以再聲明一次，假設的最后一列[0, 0, 0, 1]^T保證了合法性。

get_translation()是一個輔助函數，幫助以向量形式提取矩陣的平移部分。

get_position_from_parent_to_local_matrix()和get_position_from_local_to_parent_matrix()是從父坐標空間中提取局部坐標空間位置的函數，需要傳入變換矩陣。在某種程度上，這兩個方法是setup_local_to_parent()和setup_parent_to_local()關于位置部分的逆向操作。當然，你可以對任意矩陣使用這兩個方法（假設它是一個剛體變換），而不僅僅限于setup_local_to_parent()和setup_parent_to_local()產生的矩陣。以歐拉角形式從變換矩陣中提取方位，需要使用cEulerAngles類的一個方法。

Matrix4x3.cpp：

    #include <assert.h>
    #include <math.h>
    #include "Vector3.h"
    #include "EulerAngles.h"
    #include "Quaternion.h"
    #include "RotationMatrix.h"
    #include "Matrix4x3.h"
    #include "MathUtil.h"
   
    /////////////////////////////////////////////////////////////////////////////
    //
    // MATRIX ORGANIZATION
    //
    // The purpose of this class is so that a user might perform transformations
    // without fiddling with plus or minus signs or transposing the matrix
    // until the output "looks right."  But of course, the specifics of the
    // internal representation is important.  Not only for the implementation
    // in this file to be correct, but occasionally direct access to the
    // matrix variables is necessary, or beneficial for optimization.  Thus,
    // we document our matrix conventions here.
    //
    // We use row vectors, so multiplying by our matrix looks like this:
    //
    //               | m11 m12 m13 |
    //     [ x y z ] | m21 m22 m23 | = [ x' y' z' ]
    //               | m31 m32 m33 |
    //               | tx  ty  tz  |
    //
    // Strict adherance to linear algebra rules dictates that this
    // multiplication is actually undefined.  To circumvent this, we can
    // consider the input and output vectors as having an assumed fourth
    // coordinate of 1.  Also, since we cannot technically invert a 4x3 matrix
    // according to linear algebra rules, we will also assume a rightmost
    // column of [ 0 0 0 1 ].  This is shown below:
    //
    //                 | m11 m12 m13 0 |
    //     [ x y z 1 ] | m21 m22 m23 0 | = [ x' y' z' 1 ]
    //                 | m31 m32 m33 0 |
    //                 | tx  ty  tz  1 |
    //
    /////////////////////////////////////////////////////////////////////////////
   
    //---------------------------------------------------------------------------
    // Set the matrix to identity
    //---------------------------------------------------------------------------
    void cMatrix4x3::identity()
    {
        m11 = 1.0f;        m12 = 0.0f;        m13 = 0.0f;
        m21 = 0.0f;        m22 = 1.0f;        m23 = 0.0f;
        m31 = 0.0f;        m32 = 0.0f;        m33 = 1.0f;
        tx  = 0.0f;        ty  = 0.0f;        tz  = 1.0f;
    }
   
    //---------------------------------------------------------------------------
    // Zero the 4th row of the matrix, which contains the translation portion.
    //---------------------------------------------------------------------------
    void cMatrix4x3::zero_translation()
    {
        tx = ty = tz = 0.0f;
    }
   
    //---------------------------------------------------------------------------
    // Sets the translation portion of the matrix in vector form.
    //---------------------------------------------------------------------------
    void cMatrix4x3::set_translation(const cVector3& d)
    {
        tx = d.x;    ty = d.y;    tz = d.z;
    }
   
    //---------------------------------------------------------------------------
    // Sets the translation portion of the matrix in vector form.
    //---------------------------------------------------------------------------
    void cMatrix4x3::setup_translation(const cVector3& d)
    {
        // Set the linear transformation portion to identity
            m11 = 1.0f;        m12 = 0.0f;        m13 = 0.0f;
        m21 = 0.0f;        m22 = 1.0f;        m23 = 0.0f;
        m31 = 0.0f;        m32 = 0.0f;        m33 = 1.0f;
   
        // Set the translation portion
        tx = d.x; ty = d.y; tz = d.z;
    }
   
    //---------------------------------------------------------------------------
    // Setup the matrix to perform a local -> parent transformation, given
    // the position and orientation of the local reference frame within the
    // parent reference frame.
    //
    // A very common use of this will be to construct a object -> world matrix.
    // As an example, the transformation in this case is straightforward.  We
    // first rotate from object space into inertial space, then we translate
    // into world space.
    //
    // We allow the orientation to be specified using either euler angles,
    // or a cRotationMatrix.
    //---------------------------------------------------------------------------
    void cMatrix4x3::setup_local_to_parent(const cVector3& pos, const cEulerAngles& orient)
    {
        // create a rotation matrix
        cRotationMatrix orient_matrix;
        orient_matrix.setup(orient);
   
        // Setup the 4x3 matrix.  Note: if we were really concerned with speed, we could 
        // create the matrix directly into these variables, without using the temporary 
        // cRotationMatrix object.  This would save us a function call and a few copy operations.
            setup_local_to_parent(pos, orient_matrix);
    }
   
    void cMatrix4x3::setup_local_to_parent(const cVector3& pos, const cRotationMatrix& orient)
    {
        // Copy the rotation portion of the matrix.  According to the comments in RotationMatrix.cpp, 
        // the rotation matrix is "normally" an inertial->object matrix, which is parent->local.  
        // We want a local->parent rotation, so we must transpose while copying.
   
        m11 = orient.m11;    m12 = orient.m21;    m13 = orient.m31;
        m21 = orient.m12;    m22 = orient.m22;    m23 = orient.m32;
        m31 = orient.m13;    m32 = orient.m23;    m33 = orient.m33;
   
        // Now set the translation portion.  Translation happens "after" the 3x3 portion, 
        // so we can simply copy the position field directly.
            tx = pos.x;        ty = pos.y;        tz = pos.z;
    }
   
    //---------------------------------------------------------------------------
    // Setup the matrix to perform a parent -> local transformation, given
    // the position and orientation of the local reference frame within the
    // parent reference frame.
    //
    // A very common use of this will be to construct a world -> object matrix.
    // To perform this transformation, we would normally FIRST transform
    // from world to inertial space, and then rotate from inertial space into
    // object space.  However, our 4x3 matrix always translates last.  So
    // we think about creating two matrices T and R, and then concatenating M = TR.
    //
    // We allow the orientation to be specified using either euler angles,
    // or a RotationMatrix.
    //---------------------------------------------------------------------------
    void cMatrix4x3::setup_parent_to_local(const cVector3& pos, const cEulerAngles& orient)
    {
        // create a rotation matrix
        cRotationMatrix orient_matrix;
        orient_matrix.setup(orient);
   
        // setup the 4x3 matrix
            setup_parent_to_local(pos, orient_matrix);
    }
   
    void cMatrix4x3::setup_parent_to_local(const cVector3& pos, const cRotationMatrix& orient)
    {
        // Copy the rotation portion of the matrix.  We can copy the elements directly 
        // (without transposing) according to the layout as commented in RotationMatrix.cpp.
   
        m11 = orient.m11;    m12 = orient.m12;    m13 = orient.m13;
        m21 = orient.m21;    m22 = orient.m22;    m23 = orient.m23;
        m31 = orient.m31;    m32 = orient.m32;    m33 = orient.m33;
   
        // Now set the translation portion.  Normally, we would translate by the negative of 
        // the position to translate from world to inertial space.  However, we must correct
        // for the fact that the rotation occurs "first."  So we must rotate the translation portion. 
        // This is the same as create a translation matrix T to translate by -pos,
        // and a rotation matrix R, and then creating the matrix as the concatenation of TR.
   
        tx = -(pos.x * m11 + pos.y * m21 + pos.z * m31);
        ty = -(pos.x * m12 + pos.y * m22 + pos.z * m32);
        tz = -(pos.x * m13 + pos.y * m23 + pos.z * m33);
    }
   
    //---------------------------------------------------------------------------
    // Setup the matrix to perform a rotation about a cardinal axis.
    //
    // theta is the amount of rotation, in radians.  The left-hand rule is used to 
    // define "positive" rotation.
    //
    // The translation portion is reset.
    //
    // Rotate around x axis:
    //                | 1            0        0      |
    //        R(a) =  | 0            cosa    sina  |
    //                | 0            -sina    cosa  |
    //
    // Rotate around y axis:
    //                | cosa        0        -sina |
    //        R(a) =  | 0            1        0      |
    //                | sina        0        cosa  |
    //
    // Rotate around z axis:
    //                | cosa        sina    0      |
    //        R(a) =  | -sina        cosa    0      |
    //                | 0            0        1      |
    //---------------------------------------------------------------------------
    void cMatrix4x3::setup_rotate(int axis, float theta)
    {
        // get sin and cosine of rotation angle
        float s, c;
        sin_cos(&s, &c, theta);
   
        // check which axis they are rotating about
        switch(axis)
        {
        case ROTATE_AROUND_X:
            m11 = 1.0f; m12 = 0.0f; m13 = 0.0f;
            m21 = 0.0f; m22 = c;    m23 = s;
            m31 = 0.0f; m32 = -s;   m33 = c;
            break;
   
        case ROTATE_AROUND_Y:
            m11 = c;    m12 = 0.0f; m13 = -s;
            m21 = 0.0f; m22 = 1.0f; m23 = 0.0f;
            m31 = s;    m32 = 0.0f; m33 = c;
            break;
   
        case ROTATE_AROUND_Z: 
            m11 = c;    m12 = s;    m13 = 0.0f;
            m21 = -s;   m22 = c;    m23 = 0.0f;
            m31 = 0.0f; m32 = 0.0f; m33 = 1.0f;
            break;
   
        default:    // bogus axis index
            assert(false);
        }
   
        // Reset the translation portion
        tx = ty = tz = 0.0f;
    }
   
    //-------------------------------------------------------------------------------------
    // Setup the matrix to perform a rotation about an arbitrary axis.
    // The axis of rotation must pass through the origin.
    //
    // axis defines the axis of rotation, and must be a unit vector.
    //
    // theta is the amount of rotation, in radians.  The left-hand rule is
    // used to define "positive" rotation.
    //
    // The translation portion is reset.
    //
    //           | x^2(1-cosa) + cosa        xy(1-cosa) + zsina        xz(1-cosa) - ysina |
    // R(n, a) = | xy(1-cosa) - zsina        y^2(1-cosa) + cosa        yz(1-cosa) + xsina |
    //             | xz(1-cosa) + ysina        yz(1-cosa) - xsina        z^2(1-cosa) + cosa |
    //-------------------------------------------------------------------------------------
    void cMatrix4x3::setup_rotate(const cVector3& axis, float theta)
    {
        // Quick sanity check to make sure they passed in a unit vector to specify the axis
            assert(fabs(axis * axis - 1.0f) < 0.01f);
   
        // Get sin and cosine of rotation angle
        float s, c;
        sin_cos(&s, &c, theta);
   
        // Compute 1 - cos(theta) and some common subexpressions
   
        float    a = 1.0f - c;
        float    ax = a * axis.x;
        float    ay = a * axis.y;
        float    az = a * axis.z;
   
        // Set the matrix elements.  There is still a little more opportunity for optimization 
        // due to the many common subexpressions.  We'll let the compiler handle that
   
        m11 = ax * axis.x + c;
        m12 = ax * axis.y + axis.z * s;
        m13 = ax * axis.z - axis.y * s;
   
        m21 = ay * axis.x - axis.z * s;
        m22 = ay * axis.y + c;
        m23 = ay * axis.z + axis.x * s;
   
        m31 = az * axis.x + axis.y * s;
        m32 = az * axis.y - axis.x * s;
        m33 = az * axis.z + c;
   
        // Reset the translation portion
        tx = ty = tz = 0.0f;
    }
   
    //---------------------------------------------------------------------------
    // Setup the matrix to perform a rotation, given the angular displacement
    // in quaternion form.
    //
    // The translation portion is reset.
    //
    //     | 1 - 2(y^2 + z^2)        2(xy + wz)            2(xz - wy)         |
    // M = | 2(xy - wz)                1 - 2(x^2 + z^2)    2(yz + wx)         |
    //       | 2(xz + wy)                2(yz - wx)            1 - 2(x^2 + y^2) |
    //---------------------------------------------------------------------------
    void cMatrix4x3::from_quat(const cQuaternion& q)
    {
        // Compute a few values to optimize common subexpressions
        float double_w = 2.0f * q.w;
        float double_x = 2.0f * q.x;
        float double_y = 2.0f * q.y;
        float double_z = 2.0f * q.z;
   
        // Set the matrix elements.  There is still a little more opportunity for optimization 
        // due to the many common subexpressions.  We'll let the compiler handle that
   
        m11 = 1.0f - double_y * q.y - double_z * q.z;
        m12 = double_x * q.y + double_w * q.z;
        m13 = double_x * q.z - double_w * q.x;
   
        m21 = double_x * q.y - double_w * q.z;
        m22 = 1.0f - double_x * q.x - double_z * q.z;
        m23 = double_y * q.z + double_w * q.x; 
   
        m31 = double_x * q.z + double_w * q.y;
        m32 = double_y * q.z - double_w * q.x;
        m33 = 1.0f - double_x * q.x - double_y * q.y;
   
        // Reset the translation portion
        tx = ty = tz = 0.0f;
    }
   
    //---------------------------------------------------------------------------
    // Setup the matrix to perform scale on each axis.  For uniform scale by k,
    // use a vector of the form Vector3(k,k,k)
    //
    // The translation portion is reset.
    //---------------------------------------------------------------------------
    void cMatrix4x3::setup_scale(const cVector3& s)
    {
        // Set the matrix elements.  Pretty straightforward
            m11 = s.x;        m12 = 0.0f;        m13 = 0.0f;
        m21 = 0.0f;        m22 = s.y;        m23 = 0.0f;
        m31 = 0.0f;        m32 = 0.0f;        m33 = s.z;
   
        // Reset the translation portion
        tx = ty = tz = 0.0f;
    }
   
    //---------------------------------------------------------------------------
    // Setup the matrix to perform scale along an arbitrary axis.
    //
    // The axis is specified using a unit vector, the translation portion is reset.
    //
    //           | 1 + (k-1)x^2        (k-1)xy            (k-1)xz         |
    // S(n, k) = | (k-1)xy            1 + (k-1)y^2    (k-1)yz         |
    //             | (k-1)xz            (k-1)yz            1 + (k-1)z^2 |
    //---------------------------------------------------------------------------
    void cMatrix4x3::setup_scale_along_axis(const cVector3& axis, float k)
    {
        // Quick sanity check to make sure they passed in a unit vector to specify the axis.
            assert(fabs(axis * axis - 1.0f) < 0.01f);
   
        // Compute k-1 and some common subexpressions
        float    a = k - 1.0f;
        float    ax = a * axis.x;
        float    ay = a * axis.y;
        float    az = a * axis.z;
   
        // Fill in the matrix elements.  We'll do the common subexpression optimization 
        // ourselves here, since diagonally opposite matrix elements are equal.
   
        m11 = ax * axis.x + 1.0f;
        m22 = ay * axis.y + 1.0f;
        m32 = az * axis.z + 1.0f;
   
        m12 = m21 = ax * axis.y;
        m13 = m31 = ax * axis.z;
        m23 = m32 = ay * axis.z;
   
        // Reset the translation portion
        tx = ty = tz = 0.0f;
    }
   
    //---------------------------------------------------------------------------
    // Setup the matrix to perform a reflection about a plane parallel to a 
    // cardinal plane.
    //
    // axis is a 1-based index which specifies the plane to project about:
    //
    //    REFLECT_ABOUT_X => reflect about the plane x=k
    //    REFLECT_ABOUT_Y => reflect about the plane y=k
    //    REFLECT_ABOUT_Z => reflect about the plane z=k
    //
    // The translation is set appropriately, since translation must occur if k != 0
    //
    // Reflect about x = 0:
    //
    //       | -1.0f        0.0f    0.0f |
    // M = | 0.0f        1.0f    0.0f |
    //       | 0.0f        0.0f    1.0f |
    //
    // Reflect about y = 0:
    //
    //       | 1.0f        0.0f    0.0f |
    // M = | 0.0f       -1.0f    0.0f |
    //       | 0.0f        0.0f    1.0f |
    //
    // Reflect about z = 0:
    //
    //       | 1.0f        0.0f    0.0f |
    // M = | 0.0f        1.0f    0.0f |
    //       | 0.0f        0.0f   -1.0f |
    //---------------------------------------------------------------------------
    void cMatrix4x3::setup_shear(int axis, float s, float t)
    {
        // Check which type of shear they want
        switch (axis) 
        {
        case SHERE_AROUND_X: // Shear y and z using x
                m11 = 1.0f; m12 = s;    m13 = t;
            m21 = 0.0f; m22 = 1.0f; m23 = 0.0f;
            m31 = 0.0f; m32 = 0.0f; m33 = 1.0f;
            break;
   
        case SHERE_AROUND_Y: // Shear x and z using y
                m11 = 1.0f; m12 = 0.0f; m13 = 0.0f;
            m21 = s;    m22 = 1.0f; m23 = t;
            m31 = 0.0f; m32 = 0.0f; m33 = 1.0f;
            break;
   
        case SHERE_AROUND_Z: // Shear x and y using z
                m11 = 1.0f; m12 = 0.0f; m13 = 0.0f;
            m21 = 0.0f; m22 = 1.0f; m23 = 0.0f;
            m31 = s;    m32 = t;    m33 = 1.0f;
            break;
   
        default:
            // bogus axis index
            assert(false);
        }
   
        // Reset the translation portion
        tx = ty = tz = 0.0f;
    }
   
    //---------------------------------------------------------------------------
    // Setup the matrix to perform a projection onto a plane passing
    // through the origin.  The plane is perpendicular to the unit vector n.
    //
    //          | 1-x^2    -xy        -xz   |
    // P(n) = |-xy        1-y^2    -yz   |
    //          |-xz        -zy        1-z^2 |
    //---------------------------------------------------------------------------
    void cMatrix4x3::setup_project(const cVector3& n) 
    {
        // Quick sanity check to make sure they passed in a unit vector to specify the axis
            assert(fabs(n * n - 1.0f) < 0.01f);
   
        // Fill in the matrix elements.  We'll do the common subexpression optimization 
        // ourselves here, since diagonally opposite matrix elements are equal.
   
        m11 = 1.0f - n.x * n.x;
        m22 = 1.0f - n.y * n.y;
        m33 = 1.0f - n.z * n.z;
   
        m12 = m21 = -n.x * n.y;
        m13 = m31 = -n.x * n.z;
        m23 = m32 = -n.y * n.z;
   
        // Reset the translation portion
        tx = ty = tz = 0.0f;
    }
   
    //---------------------------------------------------------------------------
    // Setup the matrix to perform a reflection about a plane parallel to a 
    // cardinal plane.
    //
    // axis is a 1-based index which specifies the plane to project about:
    //
    //    REFLECT_ABOUT_X => reflect about the plane x=k
    //    REFLECT_ABOUT_Y => reflect about the plane y=k
    //    REFLECT_ABOUT_Z => reflect about the plane z=k
    //
    // The translation is set appropriately, since translation must occur if k != 0
    //
    // Reflect about x = 0:
    //
    //       | -1.0f        0.0f    0.0f |
    // M = | 0.0f        1.0f    0.0f |
    //       | 0.0f        0.0f    1.0f |
    //
    // Reflect about y = 0:
    //
    //       | 1.0f        0.0f    0.0f |
    // M = | 0.0f       -1.0f    0.0f |
    //       | 0.0f        0.0f    1.0f |
    //
    // Reflect about z = 0:
    //
    //       | 1.0f        0.0f    0.0f |
    // M = | 0.0f        1.0f    0.0f |
    //       | 0.0f        0.0f   -1.0f |
    //---------------------------------------------------------------------------
    void cMatrix4x3::setup_reflect(int axis, float k) 
    {
        // Check which plane they want to reflect about
        switch (axis) 
        {
        case REFLECT_ABOUT_X: // Reflect about the plane x=k
                m11 = -1.0f; m12 =  0.0f; m13 =  0.0f;
            m21 =  0.0f; m22 =  1.0f; m23 =  0.0f;
            m31 =  0.0f; m32 =  0.0f; m33 =  1.0f;
   
            tx = 2.0f * k;
            ty = 0.0f;
            tz = 0.0f;
   
            break;
   
        case REFLECT_ABOUT_Y: // Reflect about the plane y=k
                m11 =  1.0f; m12 =  0.0f; m13 =  0.0f;
            m21 =  0.0f; m22 = -1.0f; m23 =  0.0f;
            m31 =  0.0f; m32 =  0.0f; m33 =  1.0f;
   
            tx = 0.0f;
            ty = 2.0f * k;
            tz = 0.0f;
   
            break;
   
        case REFLECT_ABOUT_Z: // Reflect about the plane z=k
                m11 =  1.0f; m12 =  0.0f; m13 =  0.0f;
            m21 =  0.0f; m22 =  1.0f; m23 =  0.0f;
            m31 =  0.0f; m32 =  0.0f; m33 = -1.0f;
   
            tx = 0.0f;
            ty = 0.0f;
            tz = 2.0f * k;
   
            break;
   
        default:
            // bogus axis index
            assert(false);
        }
    }
   
    //---------------------------------------------------------------------------
    // Setup the matrix to perform a reflection about an arbitrary plane
    // through the origin.  The unit vector n is perpendicular to the plane.
    //
    // The translation portion is reset.
    //
    //                     | 1 - 2x^2        -2xy        -2xz     |
    // P(n) = S(n, -1) = | -2xy            1 - 2y^2    -2yz     |
    //                     | -2xz            -2zy        1 - 2z^2 |
    //---------------------------------------------------------------------------
    void cMatrix4x3::setup_reflect(const cVector3& n)
    {
        // Quick sanity check to make sure they passed in a unit vector to specify the axis
            assert(fabs(n * n - 1.0f) < 0.01f);
   
        // Compute common subexpressions
        float    ax = -2.0f * n.x;
        float    ay = -2.0f * n.y;
        float    az = -2.0f * n.z;
   
        // Fill in the matrix elements.  We'll do the common subexpression optimization ourselves 
        // here, since diagonally opposite matrix elements are equal.
   
        m11 = 1.0f + ax * n.x;
        m22 = 1.0f + ay * n.y;
        m33 = 1.0f + az * n.z;
   
        m12 = m21 = ax * n.y;
        m13 = m31 = ax * n.z;
        m23 = m32 = ay * n.z;
   
        // Reset the translation portion
        tx = ty = tz = 0.0f;
    }
   
    //---------------------------------------------------------------------------
    // Transform the point.  This makes using the vector class look like it
    // does with linear algebra notation on paper.
    //
    // We also provide a *= operator, as per C convention.
    //---------------------------------------------------------------------------
    cVector3 operator *(const cVector3& p, const cMatrix4x3& m)
    {
        // Grind through the linear algebra.
        return cVector3(p.x * m.m11 + p.y * m.m21 + p.z * m.m31 + m.tx,
                        p.x * m.m12 + p.y * m.m22 + p.z * m.m32 + m.ty,
                        p.x * m.m13 + p.y * m.m23 + p.z * m.m33 + m.tz);
    }
   
    cVector3& operator *=(cVector3& p, const cMatrix4x3& m)
    {
        p = p * m;
        return p;
    }
   
    //---------------------------------------------------------------------------
    // Matrix concatenation.  This makes using the vector class look like it
    // does with linear algebra notation on paper.
    //
    // We also provide a *= operator, as per C convention.
    //---------------------------------------------------------------------------
    cMatrix4x3 operator *(const cMatrix4x3& a, const cMatrix4x3& b)
    {
        cMatrix4x3 r;
   
        // Compute the upper 3x3 (linear transformation) portion
   
        r.m11 = a.m11 * b.m11 + a.m12 * b.m21 + a.m13 * b.m31;
        r.m12 = a.m11 * b.m12 + a.m12 * b.m22 + a.m13 * b.m32;
        r.m13 = a.m11 * b.m13 + a.m12 * b.m23 + a.m13 * b.m33;
   
        r.m21 = a.m21 * b.m11 + a.m22 * b.m21 + a.m23 * b.m31;
        r.m22 = a.m21 * b.m12 + a.m22 * b.m22 + a.m23 * b.m32;
        r.m23 = a.m21 * b.m13 + a.m22 * b.m23 + a.m23 * b.m33;
     
        r.m31 = a.m31 * b.m11 + a.m32 * b.m21 + a.m33 * b.m31;
        r.m32 = a.m31 * b.m12 + a.m32 * b.m22 + a.m33 * b.m32;
        r.m33 = a.m31 * b.m13 + a.m32 * b.m23 + a.m33 * b.m33;
   
        // Compute the translation portion
   
        r.tx = a.tx * b.m11 + a.ty * b.m21 + a.tz * b.m31 + b.tx;
        r.ty = a.tx * b.m12 + a.ty * b.m22 + a.tz * b.m32 + b.ty;
        r.tz = a.tx * b.m13 + a.ty * b.m23 + a.tz * b.m33 + b.tz;
   
        // Return it.  Ouch - involves a copy constructor call.  If speed
        // is critical, we may need a seperate function which places the
        // result where we want it
        return r;
    }
   
    cMatrix4x3& operator *=(cMatrix4x3& a, const cMatrix4x3& b)
    {
        a = a * b;
   
        return a;
    }
   
    //---------------------------------------------------------------------------
    // Compute the determinant of the 3x3 portion of the matrix.
    //---------------------------------------------------------------------------
    float determinant(const cMatrix4x3& m) 
    {
        float result = m.m11 * (m.m22 * m.m33 - m.m23 * m.m32) + 
                       m.m12 * (m.m23 * m.m31 - m.m21 * m.m33) + 
                       m.m13 * (m.m21 * m.m32 - m.m22 * m.m31);
   
        return result;
    }
   
        //---------------------------------------------------------------------------
    // Compute the inverse of a matrix.  We use the classical adjoint divided
    // by the determinant method.
    //---------------------------------------------------------------------------
    cMatrix4x3 inverse(const cMatrix4x3& m) 
    {    
        float det = determinant(m);
   
        // If we're singular, then the determinant is zero and there's no inverse.
        assert(fabs(det) > 0.000001f);
   
        // Compute one over the determinant, so we divide once and can multiply per element.
        float one_over_det = 1.0f / det;
   
        // Compute the 3x3 portion of the inverse, by dividing the adjoint by the determinant.
   
        cMatrix4x3    r;
   
        r.m11 = (m.m22 * m.m33 - m.m23 * m.m32) * one_over_det;
        r.m12 = (m.m13 * m.m32 - m.m12 * m.m33) * one_over_det;
        r.m13 = (m.m12 * m.m23 - m.m13 * m.m22) * one_over_det;
   
        r.m21 = (m.m23 * m.m31 - m.m21 * m.m33) * one_over_det;
        r.m22 = (m.m11 * m.m33 - m.m13 * m.m31) * one_over_det;
        r.m23 = (m.m13 * m.m21 - m.m11 * m.m23) * one_over_det;
   
        r.m31 = (m.m21 * m.m32 - m.m22 * m.m31) * one_over_det;
        r.m32 = (m.m12 * m.m31 - m.m11 * m.m32) * one_over_det;
        r.m33 = (m.m11 * m.m22 - m.m12 * m.m21) * one_over_det;
   
        // Compute the translation portion of the inverse
            r.tx = -(m.tx * r.m11 + m.ty * r.m21 + m.tz * r.m31);
        r.ty = -(m.tx * r.m12 + m.ty * r.m22 + m.tz * r.m32);
        r.tz = -(m.tx * r.m13 + m.ty * r.m23 + m.tz * r.m33);
   
        // Return it.  Ouch - involves a copy constructor call.  If speed is critical, 
        // we may need a seperate function which places the result where we want it
        return r;
    }
   
    //---------------------------------------------------------------------------
    // Return the translation row of the matrix in vector form
    //---------------------------------------------------------------------------
    cVector3 get_translation(const cMatrix4x3& m) 
    {
        return cVector3(m.tx, m.ty, m.tz);
    }
   
    //---------------------------------------------------------------------------
    // Extract the position of an object given a parent -> local transformation
    // matrix (such as a world -> object matrix).
    //
    // We assume that the matrix represents a rigid transformation.  (No scale,
    // skew, or mirroring).
    //---------------------------------------------------------------------------
        cVector3 get_position_from_parent_to_local_matrix(const cMatrix4x3& m)
    {
        // Multiply negative translation value by the transpose of the 3x3 portion.  
        // By using the transpose, we assume that the matrix is orthogonal.  
        // (This function doesn't really make sense for non-rigid transformations)
   
        return cVector3(-(m.tx * m.m11 + m.ty * m.m12 + m.tz * m.m13),
                        -(m.tx * m.m21 + m.ty * m.m22 + m.tz * m.m23),
                        -(m.tx * m.m31 + m.ty * m.m32 + m.tz * m.m33));
    }
   
    //---------------------------------------------------------------------------
    // Extract the position of an object given a local -> parent transformation
    // matrix (such as an object -> world matrix).
    //---------------------------------------------------------------------------
        cVector3 get_position_from_local_to_parent_matrix(const cMatrix4x3& m)
    {
        // Position is simply the translation portion
        return cVector3(m.tx, m.ty, m.tz);
    }