1. 程式人生 > >3D遊戲之投影矩陣演算法技術實現

3D遊戲之投影矩陣演算法技術實現

筆者介紹:姜雪偉,IT公司技術合夥人,IT高階講師,CSDN社群專家,特邀編輯,暢銷書作者,國家專利發明人,已出版書籍:《手把手教你架構3D遊戲引擎》電子工業出版社 和《Unity3D實戰核心技術詳解》電子工業出版社等書籍。

在3D遊戲開發中,固定流水線是不可迴避的問題,大部分人都是像被公式一樣的記住了固定流水線的步驟,但是對其實現原理並不是很清楚,這樣對其深入學習3D技術是非常不利的。筆者在此給讀者對相機的投影矩陣做了一下推導,幫助讀者理解矩陣的實現過程。這樣以後自己封裝底層介面時也能很容易的做到。下面主要實現了透視投影矩陣和正交投影矩陣的推導過程,如果大家對於矩陣的計算不是很清楚,可以翻閱一下大學課本《線性代數》。這個矩陣的推導也是為了我正在CSDN講的視訊課程《

演算法與遊戲實戰技術》做的補充內容。

在OpenGL或者DirectX中,在3D攝像機空間中的點需要對映到視錐體近裁剪面上(也就是投影面)才能在螢幕上看到。假設在相機空間中的點(xe, ye, ze)對映在近裁剪面後的座標是(xp,yp,zp),從頂檢視的角度顯示的效果如下圖所示:


在x軸方向看的示意圖如下:


從視錐體的頂檢視看,在相機空間的x軸,xe對映後的值為xp,可以通過相似三角形的比值計算得到,公式如下所示:


同理,通過側檢視,yp同樣可以求得:


通過這兩個算式,很容易發現,xp和yp都依賴ze的值,換句話說,它們都是除以-ze ,這樣非常容易構建投影矩陣,相機裁剪座標可以通過相機座標與投影矩陣相乘得到,裁剪座標也是一個齊次座標,它最終會通過與w相除被轉化成標準裝置座標(NDC),這也就是我們通常說的視口(0,0,1,1)轉換公式如下:



在這裡我們可以設定裁剪座標的w分量為-ze ,從而投影矩陣的第四行數值為(0,0,-1,0),將上述投影矩陣公式轉化如下:


接下來,我們設定投影xp 和yp到NDC的xn和yn的關係為:[l, r] ⇒[-1, 1] 和[b, t] ⇒[-1, 1]。先看一下從xp對映到xn的圖示如下:


根據上圖進行計算公式如下:


對映yp到yn的效果圖如下所示:


從而得到計算公式如下所示:


我們將已經計算好的xp和yp算式替換掉上面公式將得到新的公式,如下所示:



通過以上公式計算,我們可以得到投影矩陣的第一行和第二行,繼續投影矩陣公式的填充如下:


截止到目前,我們只剩了投影矩陣的第三行還沒解決,我們發現zn

有一點不同於其他的計算,因為ze在相機空間一直投影在-n的近裁剪面上。但是我們需要獨一無二的z值對於裁剪和深度測試。我們知道z值不依賴於x或者y值, 這樣我們就能指定投影矩陣的第三行如下所示:


在相機空間,we等於1,因而方程可以表示如下:



為了找到係數A和B,我們使用(ze,zn),分別表示為(-n,-1)和(-f,1),並且把它們帶入到上述公式中:


為了解方程將上述公式(1)用B表示為如下公式:


將公式(1’)帶入到(2)中後得出A的表示式:


再將A帶入(1)求得B公式如下:


求得A和B後,我們可以得到ze和zn的關係如下所示:


將其帶入公式可以得到投影矩陣的所有元素,矩陣如下所示:


這個投影矩陣終於完成了,當然我們也可以考慮到特殊情況,如果這個視口是對稱的,這樣矩陣還能簡化,對稱滿足的條件是r = -l 和t = -b,簡化的公式如下:


下面再計算一下正交攝像機的投影矩陣,相對透視攝像機矩陣,這個簡單的多,正交攝像機的視口說白了就是一個長方體,效果如下圖所示:


它們的計算公式如下所示:





對於正交攝像機來說,w值不是必須的,也就是說投影矩陣的第四行仍然是(0,0,0,1)。因而,正交攝像機的投影矩陣公式如下所示:


至此,矩陣的推導就全部結束了,3D引擎的底層都是採用矩陣的方式進行了封裝,但是作為開發者一定要清楚其實現原理。不論開發者使用任何引擎,它的底層封裝都需要矩陣計算,為了能夠幫助讀者理解,下面把封裝矩陣的程式碼給讀者展示一下:

class Matrix4
{
public:
    // constructors
    Matrix4();  // init with identity
    Matrix4(const float src[16]);
    Matrix4(float xx, float xy, float xz, float xw,
            float yx, float yy, float yz, float yw,
            float zx, float zy, float zz, float zw,
            float wx, float wy, float wz, float ww);

    void        set(const float src[16]);
    void        set(float xx, float xy, float xz, float xw,
                    float yx, float yy, float yz, float yw,
                    float zx, float zy, float zz, float zw,
                    float wx, float wy, float wz, float ww);
    void        setRow(int index, const float row[4]);
    void        setRow(int index, const Vector4& v);
    void        setRow(int index, const Vector3& v);
    void        setColumn(int index, const float col[4]);
    void        setColumn(int index, const Vector4& v);
    void        setColumn(int index, const Vector3& v);

    const float* get() const;
    const float* getTranspose();                        // return transposed matrix
    float        getDeterminant();

    Matrix4&    identity();
    Matrix4&    transpose();                            // transpose itself and return reference
    Matrix4&    invert();                               // check best inverse method before inverse
    Matrix4&    invertEuclidean();                      // inverse of Euclidean transform matrix
    Matrix4&    invertAffine();                         // inverse of affine transform matrix
    Matrix4&    invertProjective();                     // inverse of projective matrix using partitioning
    Matrix4&    invertGeneral();                        // inverse of generic matrix

    // transform matrix
    Matrix4&    translate(float x, float y, float z);   // translation by (x,y,z)
    Matrix4&    translate(const Vector3& v);            //
    Matrix4&    rotate(float angle, const Vector3& axis); // rotate angle(degree) along the given axix
    Matrix4&    rotate(float angle, float x, float y, float z);
    Matrix4&    rotateX(float angle);                   // rotate on X-axis with degree
    Matrix4&    rotateY(float angle);                   // rotate on Y-axis with degree
    Matrix4&    rotateZ(float angle);                   // rotate on Z-axis with degree
    Matrix4&    scale(float scale);                     // uniform scale
    Matrix4&    scale(float sx, float sy, float sz);    // scale by (sx, sy, sz) on each axis

    // operators
    Matrix4     operator+(const Matrix4& rhs) const;    // add rhs
    Matrix4     operator-(const Matrix4& rhs) const;    // subtract rhs
    Matrix4&    operator+=(const Matrix4& rhs);         // add rhs and update this object
    Matrix4&    operator-=(const Matrix4& rhs);         // subtract rhs and update this object
    Vector4     operator*(const Vector4& rhs) const;    // multiplication: v' = M * v
    Vector3     operator*(const Vector3& rhs) const;    // multiplication: v' = M * v
    Matrix4     operator*(const Matrix4& rhs) const;    // multiplication: M3 = M1 * M2
    Matrix4&    operator*=(const Matrix4& rhs);         // multiplication: M1' = M1 * M2
    bool        operator==(const Matrix4& rhs) const;   // exact compare, no epsilon
    bool        operator!=(const Matrix4& rhs) const;   // exact compare, no epsilon
    float       operator[](int index) const;            // subscript operator v[0], v[1]
    float&      operator[](int index);                  // subscript operator v[0], v[1]

    friend Matrix4 operator-(const Matrix4& m);                     // unary operator (-)
    friend Matrix4 operator*(float scalar, const Matrix4& m);       // pre-multiplication
    friend Vector3 operator*(const Vector3& vec, const Matrix4& m); // pre-multiplication
    friend Vector4 operator*(const Vector4& vec, const Matrix4& m); // pre-multiplication
    friend std::ostream& operator<<(std::ostream& os, const Matrix4& m);

protected:

private:
    float       getCofactor(float m0, float m1, float m2,
                            float m3, float m4, float m5,
                            float m6, float m7, float m8);

    float m[16];
    float tm[16];                                       // transpose m

};

它的原始碼實現如下所示:
Matrix4& Matrix4::transpose()
{
    std::swap(m[1],  m[4]);
    std::swap(m[2],  m[8]);
    std::swap(m[3],  m[12]);
    std::swap(m[6],  m[9]);
    std::swap(m[7],  m[13]);
    std::swap(m[11], m[14]);

    return *this;
}



///////////////////////////////////////////////////////////////////////////////
// inverse 4x4 matrix
///////////////////////////////////////////////////////////////////////////////
Matrix4& Matrix4::invert()
{
    // If the 4th row is [0,0,0,1] then it is affine matrix and
    // it has no projective transformation.
    if(m[12] == 0 && m[13] == 0 && m[14] == 0 && m[15] == 1)
        this->invertAffine();
    else
    {
        this->invertGeneral();
        /*@@ invertProjective() is not optimized (slower than generic one)
        if(fabs(m[0]*m[5] - m[1]*m[4]) > 0.00001f)
            this->invertProjective();   // inverse using matrix partition
        else
            this->invertGeneral();      // generalized inverse
        */
    }

    return *this;
}


Matrix4& Matrix4::invertEuclidean()
{
    // transpose 3x3 rotation matrix part
    // | R^T | 0 |
    // | ----+-- |
    // |  0  | 1 |
    float tmp;
    tmp = m[1];  m[1] = m[4];  m[4] = tmp;
    tmp = m[2];  m[2] = m[8];  m[8] = tmp;
    tmp = m[6];  m[6] = m[9];  m[9] = tmp;

    // compute translation part -R^T * T
    // | 0 | -R^T x |
    // | --+------- |
    // | 0 |   0    |
    float x = m[3];
    float y = m[7];
    float z = m[11];
    m[3]  = -(m[0] * x + m[1] * y + m[2] * z);
    m[7]  = -(m[4] * x + m[5] * y + m[6] * z);
    m[11] = -(m[8] * x + m[9] * y + m[10]* z);

    // last row should be unchanged (0,0,0,1)

    return *this;
}



///////////////////////////////////////////////////////////////////////////////

Matrix4& Matrix4::invertAffine()
{
    // R^-1
    Matrix3 r(m[0],m[1],m[2], m[4],m[5],m[6], m[8],m[9],m[10]);
    r.invert();
    m[0] = r[0];  m[1] = r[1];  m[2] = r[2];
    m[4] = r[3];  m[5] = r[4];  m[6] = r[5];
    m[8] = r[6];  m[9] = r[7];  m[10]= r[8];

    // -R^-1 * T
    float x = m[3];
    float y = m[7];
    float z = m[11];
    m[3]  = -(r[0] * x + r[1] * y + r[2] * z);
    m[7]  = -(r[3] * x + r[4] * y + r[5] * z);
    m[11] = -(r[6] * x + r[7] * y + r[8] * z);

    // last row should be unchanged (0,0,0,1)
    //m[12] = m[13] = m[14] = 0.0f;
    //m[15] = 1.0f;

    return * this;
}



///////////////////////////////////////////////////////////////////////////////
// inverse matrix using matrix partitioning (blockwise inverse)
// It devides a 4x4 matrix into 4 of 2x2 matrices. It works in case of where
// det(A) != 0. If not, use the generic inverse method
// inverse formula.
// M = [ A | B ]    A, B, C, D are 2x2 matrix blocks
//     [ --+-- ]    det(M) = |A| * |D - ((C * A^-1) * B)|
//     [ C | D ]
//
// M^-1 = [ A' | B' ]   A' = A^-1 - (A^-1 * B) * C'
//        [ ---+--- ]   B' = (A^-1 * B) * -D'
//        [ C' | D' ]   C' = -D' * (C * A^-1)
//                      D' = (D - ((C * A^-1) * B))^-1
//
// NOTE: I wrap with () if it it used more than once.
//       The matrix is invertable even if det(A)=0, so must check det(A) before
//       calling this function, and use invertGeneric() instead.
///////////////////////////////////////////////////////////////////////////////
Matrix4& Matrix4::invertProjective()
{
    // partition
    Matrix2 a(m[0], m[1], m[4], m[5]);
    Matrix2 b(m[2], m[3], m[6], m[7]);
    Matrix2 c(m[8], m[9], m[12], m[13]);
    Matrix2 d(m[10], m[11], m[14], m[15]);

    // pre-compute repeated parts
    a.invert();             // A^-1
    Matrix2 ab = a * b;     // A^-1 * B
    Matrix2 ca = c * a;     // C * A^-1
    Matrix2 cab = ca * b;   // C * A^-1 * B
    Matrix2 dcab = d - cab; // D - C * A^-1 * B

    // check determinant if |D - C * A^-1 * B| = 0
    //NOTE: this function assumes det(A) is already checked. if |A|=0 then,
    //      cannot use this function.
    float determinant = dcab[0] * dcab[3] - dcab[1] * dcab[2];
    if(fabs(determinant) <= 0.00001f)
    {
        return identity();
    }

    // compute D' and -D'
    Matrix2 d1 = dcab;      //  (D - C * A^-1 * B)
    d1.invert();            //  (D - C * A^-1 * B)^-1
    Matrix2 d2 = -d1;       // -(D - C * A^-1 * B)^-1

    // compute C'
    Matrix2 c1 = d2 * ca;   // -D' * (C * A^-1)

    // compute B'
    Matrix2 b1 = ab * d2;   // (A^-1 * B) * -D'

    // compute A'
    Matrix2 a1 = a - (ab * c1); // A^-1 - (A^-1 * B) * C'

    // assemble inverse matrix
    m[0] = a1[0];  m[1] = a1[1];  m[2] = b1[0];  m[3] = b1[1];
    m[4] = a1[2];  m[5] = a1[3];  m[6] = b1[1];  m[2] = b1[3];
    m[8] = c1[0];  m[9] = c1[1];  m[10]= d1[0];  m[11]= d1[1];
    m[12]= c1[2];  m[13]= c1[3];  m[14]= d1[2];  m[15]= d1[3];

    return *this;
}



///////////////////////////////////////////////////////////////////////////////
// compute the inverse of a general 4x4 matrix using Cramer's Rule
// If cannot find inverse, return indentity matrix
// M^-1 = adj(M) / det(M)
///////////////////////////////////////////////////////////////////////////////
Matrix4& Matrix4::invertGeneral()
{
    // get cofactors of minor matrices
    float cofactor0 = getCofactor(m[5],m[6],m[7], m[9],m[10],m[11], m[13],m[14],m[15]);
    float cofactor1 = getCofactor(m[4],m[6],m[7], m[8],m[10],m[11], m[12],m[14],m[15]);
    float cofactor2 = getCofactor(m[4],m[5],m[7], m[8],m[9], m[11], m[12],m[13],m[15]);
    float cofactor3 = getCofactor(m[4],m[5],m[6], m[8],m[9], m[10], m[12],m[13],m[14]);

    // get determinant
    float determinant = m[0] * cofactor0 - m[1] * cofactor1 + m[2] * cofactor2 - m[3] * cofactor3;
    if(fabs(determinant) <= 0.00001f)
    {
        return identity();
    }

    // get rest of cofactors for adj(M)
    float cofactor4 = getCofactor(m[1],m[2],m[3], m[9],m[10],m[11], m[13],m[14],m[15]);
    float cofactor5 = getCofactor(m[0],m[2],m[3], m[8],m[10],m[11], m[12],m[14],m[15]);
    float cofactor6 = getCofactor(m[0],m[1],m[3], m[8],m[9], m[11], m[12],m[13],m[15]);
    float cofactor7 = getCofactor(m[0],m[1],m[2], m[8],m[9], m[10], m[12],m[13],m[14]);

    float cofactor8 = getCofactor(m[1],m[2],m[3], m[5],m[6], m[7],  m[13],m[14],m[15]);
    float cofactor9 = getCofactor(m[0],m[2],m[3], m[4],m[6], m[7],  m[12],m[14],m[15]);
    float cofactor10= getCofactor(m[0],m[1],m[3], m[4],m[5], m[7],  m[12],m[13],m[15]);
    float cofactor11= getCofactor(m[0],m[1],m[2], m[4],m[5], m[6],  m[12],m[13],m[14]);

    float cofactor12= getCofactor(m[1],m[2],m[3], m[5],m[6], m[7],  m[9], m[10],m[11]);
    float cofactor13= getCofactor(m[0],m[2],m[3], m[4],m[6], m[7],  m[8], m[10],m[11]);
    float cofactor14= getCofactor(m[0],m[1],m[3], m[4],m[5], m[7],  m[8], m[9], m[11]);
    float cofactor15= getCofactor(m[0],m[1],m[2], m[4],m[5], m[6],  m[8], m[9], m[10]);

    // build inverse matrix = adj(M) / det(M)
    // adjugate of M is the transpose of the cofactor matrix of M
    float invDeterminant = 1.0f / determinant;
    m[0] =  invDeterminant * cofactor0;
    m[1] = -invDeterminant * cofactor4;
    m[2] =  invDeterminant * cofactor8;
    m[3] = -invDeterminant * cofactor12;

    m[4] = -invDeterminant * cofactor1;
    m[5] =  invDeterminant * cofactor5;
    m[6] = -invDeterminant * cofactor9;
    m[7] =  invDeterminant * cofactor13;

    m[8] =  invDeterminant * cofactor2;
    m[9] = -invDeterminant * cofactor6;
    m[10]=  invDeterminant * cofactor10;
    m[11]= -invDeterminant * cofactor14;

    m[12]= -invDeterminant * cofactor3;
    m[13]=  invDeterminant * cofactor7;
    m[14]= -invDeterminant * cofactor11;
    m[15]=  invDeterminant * cofactor15;

    return *this;
}



///////////////////////////////////////////////////////////////////////////////
// return determinant of 4x4 matrix
///////////////////////////////////////////////////////////////////////////////
float Matrix4::getDeterminant()
{
    return m[0] * getCofactor(m[5],m[6],m[7], m[9],m[10],m[11], m[13],m[14],m[15]) -
           m[1] * getCofactor(m[4],m[6],m[7], m[8],m[10],m[11], m[12],m[14],m[15]) +
           m[2] * getCofactor(m[4],m[5],m[7], m[8],m[9], m[11], m[12],m[13],m[15]) -
           m[3] * getCofactor(m[4],m[5],m[6], m[8],m[9], m[10], m[12],m[13],m[14]);
}



///////////////////////////////////////////////////////////////////////////////
// compute cofactor of 3x3 minor matrix without sign
// input params are 9 elements of the minor matrix
// NOTE: The caller must know its sign.
///////////////////////////////////////////////////////////////////////////////
float Matrix4::getCofactor(float m0, float m1, float m2,
                           float m3, float m4, float m5,
                           float m6, float m7, float m8)
{
    return m0 * (m4 * m8 - m5 * m7) -
           m1 * (m3 * m8 - m5 * m6) +
           m2 * (m3 * m7 - m4 * m6);
}



///////////////////////////////////////////////////////////////////////////////
// translate this matrix by (x, y, z)
///////////////////////////////////////////////////////////////////////////////
Matrix4& Matrix4::translate(const Vector3& v)
{
    return translate(v.x, v.y, v.z);
}

Matrix4& Matrix4::translate(float x, float y, float z)
{
    m[0] += m[12]*x;   m[1] += m[13]*x;   m[2] += m[14]*x;   m[3] += m[15]*x;
    m[4] += m[12]*y;   m[5] += m[13]*y;   m[6] += m[14]*y;   m[7] += m[15]*y;
    m[8] += m[12]*z;   m[9] += m[13]*z;   m[10]+= m[14]*z;   m[11]+= m[15]*z;
    return *this;
}



///////////////////////////////////////////////////////////////////////////////
// uniform scale
///////////////////////////////////////////////////////////////////////////////
Matrix4& Matrix4::scale(float s)
{
    return scale(s, s, s);
}

Matrix4& Matrix4::scale(float x, float y, float z)
{
    m[0] = m[0]*x;   m[1] = m[1]*x;   m[2] = m[2]*x;   m[3] = m[3]*x;
    m[4] = m[4]*y;   m[5] = m[5]*y;   m[6] = m[6]*y;   m[7] = m[7]*y;
    m[8] = m[8]*z;   m[9] = m[9]*z;   m[10]= m[10]*z;  m[11]= m[11]*z;
    return *this;
}



///////////////////////////////////////////////////////////////////////////////
// build a rotation matrix with given angle(degree) and rotation axis, then
// multiply it with this object
///////////////////////////////////////////////////////////////////////////////
Matrix4& Matrix4::rotate(float angle, const Vector3& axis)
{
    return rotate(angle, axis.x, axis.y, axis.z);
}

Matrix4& Matrix4::rotate(float angle, float x, float y, float z)
{
    float c = cosf(angle * DEG2RAD);    // cosine
    float s = sinf(angle * DEG2RAD);    // sine
    float xx = x * x;
    float xy = x * y;
    float xz = x * z;
    float yy = y * y;
    float yz = y * z;
    float zz = z * z;

    // build rotation matrix
    Matrix4 m;
    m[0] = xx * (1 - c) + c;
    m[1] = xy * (1 - c) - z * s;
    m[2] = xz * (1 - c) + y * s;
    m[3] = 0;
    m[4] = xy * (1 - c) + z * s;
    m[5] = yy * (1 - c) + c;
    m[6] = yz * (1 - c) - x * s;
    m[7] = 0;
    m[8] = xz * (1 - c) - y * s;
    m[9] = yz * (1 - c) + x * s;
    m[10]= zz * (1 - c) + c;
    m[11]= 0;
    m[12]= 0;
    m[13]= 0;
    m[14]= 0;
    m[15]= 1;

    // multiply it
    *this = m * (*this);

    return *this;
}

Matrix4& Matrix4::rotateX(float angle)
{
    float c = cosf(angle * DEG2RAD);
    float s = sinf(angle * DEG2RAD);
    float m4 = m[4], m5 = m[5], m6 = m[6],  m7 = m[7],
          m8 = m[8], m9 = m[9], m10= m[10], m11= m[11];

    m[4] = m4 * c + m8 *-s;
    m[5] = m5 * c + m9 *-s;
    m[6] = m6 * c + m10*-s;
    m[7] = m7 * c + m11*-s;
    m[8] = m4 * s + m8 * c;
    m[9] = m5 * s + m9 * c;
    m[10]= m6 * s + m10* c;
    m[11]= m7 * s + m11* c;

    return *this;
}

Matrix4& Matrix4::rotateY(float angle)
{
    float c = cosf(angle * DEG2RAD);
    float s = sinf(angle * DEG2RAD);
    float m0 = m[0], m1 = m[1], m2 = m[2],  m3 = m[3],
          m8 = m[8], m9 = m[9], m10= m[10], m11= m[11];

    m[0] = m0 * c + m8 * s;
    m[1] = m1 * c + m9 * s;
    m[2] = m2 * c + m10* s;
    m[3] = m3 * c + m11* s;
    m[8] = m0 *-s + m8 * c;
    m[9] = m1 *-s + m9 * c;
    m[10]= m2 *-s + m10* c;
    m[11]= m3 *-s + m11* c;

    return *this;
}

Matrix4& Matrix4::rotateZ(float angle)
{
    float c = cosf(angle * DEG2RAD);
    float s = sinf(angle * DEG2RAD);
    float m0 = m[0], m1 = m[1], m2 = m[2],  m3 = m[3],
          m4 = m[4], m5 = m[5], m6 = m[6],  m7 = m[7];

    m[0] = m0 * c + m4 *-s;
    m[1] = m1 * c + m5 *-s;
    m[2] = m2 * c + m6 *-s;
    m[3] = m3 * c + m7 *-s;
    m[4] = m0 * s + m4 * c;
    m[5] = m1 * s + m5 * c;
    m[6] = m2 * s + m6 * c;
    m[7] = m3 * s + m7 * c;

    return *this;
}

矩陣實現了,接下來就需要呼叫其介面實現視景體裁剪和投影矩陣了,程式碼如下所示:
///////////////////////////////////////////////////////////////////////////////
// compute 8 vertices of frustum
///////////////////////////////////////////////////////////////////////////////
void ModelGL::computeFrustumVertices(float l, float r, float b, float t, float n, float f)
{
    float ratio;
    float farLeft;
    float farRight;
    float farBottom;
    float farTop;

    // perspective mode
    if(projectionMode == 0)
        ratio     = f / n;
    // orthographic mode
    else
        ratio = 1;
    farLeft   = l * ratio;
    farRight  = r * ratio;
    farBottom = b * ratio;
    farTop    = t * ratio;

    // compute 8 vertices of the frustum
    // near top right
    frustumVertices[0].x = r;
    frustumVertices[0].y = t;
    frustumVertices[0].z = -n;

    // near top left
    frustumVertices[1].x = l;
    frustumVertices[1].y = t;
    frustumVertices[1].z = -n;

    // near bottom left
    frustumVertices[2].x = l;
    frustumVertices[2].y = b;
    frustumVertices[2].z = -n;

    // near bottom right
    frustumVertices[3].x = r;
    frustumVertices[3].y = b;
    frustumVertices[3].z = -n;

    // far top right
    frustumVertices[4].x = farRight;
    frustumVertices[4].y = farTop;
    frustumVertices[4].z = -f;

    // far top left
    frustumVertices[5].x = farLeft;
    frustumVertices[5].y = farTop;
    frustumVertices[5].z = -f;

    // far bottom left
    frustumVertices[6].x = farLeft;
    frustumVertices[6].y = farBottom;
    frustumVertices[6].z = -f;

    // far bottom right
    frustumVertices[7].x = farRight;
    frustumVertices[7].y = farBottom;
    frustumVertices[7].z = -f;

    // compute normals
    frustumNormals[0] = (frustumVertices[5] - frustumVertices[1]).cross(frustumVertices[2] - frustumVertices[1]);
    frustumNormals[0].normalize();

    frustumNormals[1] = (frustumVertices[3] - frustumVertices[0]).cross(frustumVertices[4] - frustumVertices[0]);
    frustumNormals[1].normalize();

    frustumNormals[2] = (frustumVertices[6] - frustumVertices[2]).cross(frustumVertices[3] - frustumVertices[2]);
    frustumNormals[2].normalize();

    frustumNormals[3] = (frustumVertices[4] - frustumVertices[0]).cross(frustumVertices[1] - frustumVertices[0]);
    frustumNormals[3].normalize();

    frustumNormals[4] = (frustumVertices[1] - frustumVertices[0]).cross(frustumVertices[3] - frustumVertices[0]);
    frustumNormals[4].normalize();

    frustumNormals[5] = (frustumVertices[7] - frustumVertices[4]).cross(frustumVertices[5] - frustumVertices[4]);
    frustumNormals[5].normalize();
}



///////////////////////////////////////////////////////////////////////////////
// update projection matrix
///////////////////////////////////////////////////////////////////////////////
void ModelGL::updateProjectionMatrix()
{
    if(projectionMode == 0)         // perspective
    {
        setFrustum(projectionLeft, projectionRight,
                   projectionBottom, projectionTop,
                   projectionNear, projectionFar);
    }
    else if(projectionMode == 1)    // orthographic
    {
        setOrthoFrustum(projectionLeft, projectionRight,
                        projectionBottom, projectionTop,
                        projectionNear, projectionFar);
    }
}
///////////////////////////////////////////////////////////////////////////////
// set a perspective frustum with 6 params similar to glFrustum()
// (left, right, bottom, top, near, far)
// Note: this is for row-major notation. OpenGL needs transpose it
///////////////////////////////////////////////////////////////////////////////
void ModelGL::setFrustum(float l, float r, float b, float t, float n, float f)
{
    matrixProjection.identity();
    matrixProjection[0]  =  2 * n / (r - l);
    matrixProjection[2]  =  (r + l) / (r - l);
    matrixProjection[5]  =  2 * n / (t - b);
    matrixProjection[6]  =  (t + b) / (t - b);
    matrixProjection[10] = -(f + n) / (f - n);
    matrixProjection[11] = -(2 * f * n) / (f - n);
    matrixProjection[14] = -1;
    matrixProjection[15] =  0;
}



///////////////////////////////////////////////////////////////////////////////
// set a symmetric perspective frustum with 4 params similar to gluPerspective
// (vertical field of view, aspect ratio, near, far)
///////////////////////////////////////////////////////////////////////////////
void ModelGL::setFrustum(float fovY, float aspectRatio, float front, float back)
{
    float tangent = tanf(fovY/2 * DEG2RAD);   // tangent of half fovY
    float height = front * tangent;           // half height of near plane
    float width = height * aspectRatio;       // half width of near plane

    // params: left, right, bottom, top, near, far
    setFrustum(-width, width, -height, height, front, back);
}



///////////////////////////////////////////////////////////////////////////////
// set a orthographic frustum with 6 params similar to glOrtho()
// (left, right, bottom, top, near, far)
// Note: this is for row-major notation. OpenGL needs transpose it
///////////////////////////////////////////////////////////////////////////////
void ModelGL::setOrthoFrustum(float l, float r, float b, float t, float n, float f)
{
    matrixProjection.identity();
    matrixProjection[0]  =  2 / (r - l);
    matrixProjection[3]  =  -(r + l) / (r - l);
    matrixProjection[5]  =  2 / (t - b);
    matrixProjection[7]  =  -(t + b) / (t - b);
    matrixProjection[10] = -2 / (f - n);
    matrixProjection[11] = -(f + n) / (f - n);
}

最終實現的透視投影和正交投影效果圖如下所示:



關於投影矩陣的推導和實現就給讀者推導到這裡,後面會把固定流水線的其它矩陣推導也釋出出來,供大家參考。