天行健 君子當(dāng)自強而不息

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

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

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

Matrix4x3.h:

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

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

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

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

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

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

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

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

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

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

（11）setup_reflect()創(chuàng)建一個沿平面鏡像的矩陣。第一個版本中，坐標(biāo)軸用一個數(shù)字指定，平面不必穿過原點。第二個版本中指定任意的法向量，且平面必須穿過原點。（對于沿任意不穿過原點的平面鏡像，必須將該矩陣和適當(dāng)?shù)淖儞Q矩陣連接。）

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

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

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

get_position_from_parent_to_local_matrix()和get_position_from_local_to_parent_matrix()是從父坐標(biāo)空間中提取局部坐標(biāo)空間位置的函數(shù)，需要傳入變換矩陣。在某種程度上，這兩個方法是setup_local_to_parent()和setup_parent_to_local()關(guān)于位置部分的逆向操作。當(dāng)然，你可以對任意矩陣使用這兩個方法（假設(shè)它是一個剛體變換），而不僅僅限于setup_local_to_parent()和setup_parent_to_local()產(chǎn)生的矩陣。以歐拉角形式從變換矩陣中提取方位，需要使用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);
    }