從零實現3D映像引擎：(7)矩陣函數庫

1. 數學分析

1) 矩陣到底是什麼，用它來幹嘛？

 千萬不要把矩陣想複雜了，說白了，他和算盤是一類東西，矩陣根本就沒有實際意義，他就是個數學工具而已，就好像敲算盤，珠子上上下下的按照一定的規律去撥弄，就能得到結果一樣，矩陣是數學家們發明的一個工具，這個工具也有一系列規則，通過這些規則也可以計算出一些結果，只不過他不像算盤是用來做數字加減乘除用的，他是用來求解方程組的。學習完了矩陣的運算規則，就可以發現，所有對向量、點座標的變換都可以通過矩陣運算來完成，這也是為什麼3D映像的運算與矩陣息息相關了。

 

比如方程組：

x + 2y = 3

4x + 5y = 6

就可以用三個矩陣來表示他們。

 

係數矩陣：A =

| 1  2 |

| 4  5 |

 

變數矩陣：B =

| x |

| y |

 

常量矩陣：C =

| 3 |

| 6 |

 

所以上面的方程組也可以寫成是：A * B = C

所以解方程組的問題也變成了求B。後面我會介紹到矩陣的逆，就會說到如何用A和C矩陣來求B。

 

 

2) 單位矩陣

定義：矩陣的主對角線所有元素都是1，其餘都是0。用I表示。

如：

| 1  0 |

| 0  1 |

或：

| 1  0  0 |

| 0  1  0 |

| 0  0  1 |

單位矩陣的重要用途是因為他滿足下面的公式：

A * I = I * A = A

 

 

3) 零矩陣

非常簡單：

| 0  0 |

| 0  0 |

所有元素都是零。

 

 

4) 矩陣加法和減法

兩個矩陣的相應元素相加或相減即可。

 

 

5) 矩陣的轉置

將矩陣的行和列元素交換，就是轉置，比如A =

| 0  1  2 |

| 3  4  5 |

的轉置為 A' =

| 0  3 |

| 1  4 |

| 2  5 |

 

 

6) 矩陣與標量相乘

每個元素都乘以標量即可。

 

 

7) 矩陣與矩陣相乘

並不是任意兩個矩陣都可以相乘的，矩陣A和矩陣B相乘，只有它們的內維數相等，才有意義。

比如，A是m×nj矩陣，那麼B必須是n×r矩陣。

運演算法則是：A*B=C，那麼Cij = A中第i行的向量與B中第j列的向量的點積。

比如A =

| 0  1 |

| 2  3 |

 

B =

| 0  1  2 |

| 3  4  5 |

 

那麼C =

| (0,1).(0,3)  (0,1).(1,4)  (0,1).(2,5) |

| (2,3).(0,3)  (2,3).(1,4)  (2,3).(2,5) |  =

 

| (0*0+1*3)  (0*1+1*4)  (0*2+1*5) |

| (2*0+3*3)  (2*1+3*4)  (2*2+3*5) |  =

 

|  3   4  5  |

|  9 14 19 |

 

8) 矩陣運算定律

1. A + B = B + A

2. A + (B + C) = (A + B) + C

3. A * (B * C) = (A * B)  * C

4. A * (B + C) = A * B + A * C

5. k * (A + B) = k * A+ k * B

6. (A + B) * C = A * C + B * C

7. A * I = I * A = A

 

注意：矩陣一般不滿足乘法交換律,即 A * B 一般不等於 B * A

 

 

9) 矩陣的行列式

行列式來源：每一個方形矩陣可以和一個成為矩陣的行列式的實數相對應，這個數值將可以告訴我們矩陣是否是奇異的。

矩陣A的行列式表示為det(A)

在實際應用中，我們需要通過求行列式來求矩陣的逆，所以行列式的意義非常重大。

具體如何通過上面的來源來推匯出求行列式的公式，在任何一本線性代數數上都有，我就省略了。公式如下：

對於n×n矩陣A，

det(A) = a11  (當n=1)

det(A) = a11 * A11 + a12 * A12 + ... + a1n * A1n  (當n>1時)

其中A1j = (-1)的1+j次方 * detA(M1j)  (其中j = 1, ... , n)

為第一行元素的餘子式

 

說點具體的吧，2×2矩陣

| a  b |

| c  d |

的行列式為det(A) = a*d - b*c

 

3×3矩陣

| a00  a01  a02 |

| a10  a11  a12 |

| a20  a21  a22 |

的行列式為det(A) = a00*a11*a22 + a01*a12*a20 + a02*a10*a21 - a02*a11*a20 - a01*a10*a22 - a00*a12*a21

這個公式是使用克萊姆法則求出的，為什麼上面我介紹餘子式，因為使用克萊姆法則來求3×3矩陣的行列式要進行9次乘法和5次加法。但如果我們使用餘子式來求行列式，則可以減少兩次乘法：

det(A) = a00*(a11*a22 - a21*a12) - a01*(a10*a22 - a20*a12) - a02*(a10*a21 - a20*a11)

 

 

10) 矩陣的逆

定義：對於矩陣A，有矩陣A-1，使得A * A-1 = I，則A-1為矩陣A的逆。其實這個A-1就是A的乘法逆元。

數學上還有這種說法：如果一個n×n矩陣，如果不存在乘法逆元，則稱這個矩陣是奇異的。

他的作用十分重要，比如上面說的方程組：

A * B = C，求矩陣B

那麼通過將方程兩邊同時乘以A-1，則：

(A-1 * A) * B = A-1 * C

B = A-1 * C，這樣就可以求得B了。

 

 

2. 代碼實現

1) 我提供了1×3、1×4、3×3、4×4矩陣，因為對於2D點或向量，可以使用齊次座標在3×3矩陣中存放，同理，對於3D點或向量可以使用齊次座標在4×4矩陣中存放，所以只需要寫上面這4種矩陣的運算函數，就可以滿足所有的3D矩陣變換了。這次我把以前寫的向量的結構體也做了修改，因為向量、點，通常可以表示為一個1×3或1×4矩陣，所以我把他們的記憶體結構修改為和矩陣的相同，以便直接強制轉換。

下面是矩陣結構體的定義：

typedef struct MATRIX1X3_TYPE // 1×3矩陣<br />{<br />union<br />{<br />double M[3];<br />struct<br />{<br />double M00, M01, M02;<br />};<br />};<br />} MATRIX1X3, *MATRIX1X3_PTR;</p><p>typedef struct MATRIX1X4_TYPE // 1×4矩陣<br />{<br />union<br />{<br />double M[4];<br />struct<br />{<br />double M00, M01, M02, M03;<br />};<br />};<br />} MATRIX1X4, *MATRIX1X4_PTR;</p><p>typedef struct MATRIX3X3_TYPE // 3×3矩陣<br />{<br />union<br />{<br />double M[3][3];<br />struct<br />{<br />double M00, M01, M02;<br />double M10, M11, M12;<br />double M20, M21, M22;<br />};<br />};<br />} MATRIX3X3, *MATRIX3X3_PTR;</p><p>typedef struct MATRIX4X4_TYPE //4×4矩陣<br />{<br />union<br />{<br />double M[4][4];<br />struct<br />{<br />double M00, M01, M02, M03;<br />double M10, M11, M12, M13;<br />double M20, M21, M22, M23;<br />double M30, M31, M32, M33;<br />};<br />};<br />} MATRIX4X4, *MATRIX4X4_PTR;

 

下面是對原向量和點的結構體的修改：

 typedef struct POINT2D_TYPE // 2D笛卡爾座標, 2D向量<br />{<br />union<br />{<br />double M[2]; // 數組方式<br />struct // 變數方式<br />{<br />double x;<br />double y;<br />};<br />};<br />} POINT2D, *POINT2D_PTR, VECTOR2D, *VECTOR2D_PTR;</p><p>typedef struct POINT3D_TYPE // 3D笛卡爾座標, 3D向量<br />{<br />union<br />{<br />double M[3];<br />struct<br />{<br />double x;<br />double y;<br />double z;<br />};<br />};<br />} POINT3D, *POINT3D_PTR, VECTOR3D, *VECTOR3D_PTR;</p><p>typedef struct VECTOR4D_TYPE // 4D向量<br />{<br />union<br />{<br />double M[4];<br />struct<br />{<br />double x;<br />double y;<br />double z;<br />double w;<br />};<br />};<br />} VECTOR4D, *VECTOR4D_PTR;

 

 2) 矩陣操作函數實現

void _CPPYIN_Math::MatrixCreate(MATRIX3X3_PTR pm,<br />double m00, double m01, double m02,<br />double m10, double m11, double m12,<br />double m20, double m21, double m22)<br />{<br />pm->M00 = m00;<br />pm->M01 = m01;<br />pm->M02 = m02;<br />pm->M10 = m10;<br />pm->M11 = m11;<br />pm->M12 = m12;<br />pm->M20 = m20;<br />pm->M21 = m21;<br />pm->M22 = m22;<br />}</p><p>void _CPPYIN_Math::MatrixCreate(MATRIX4X4_PTR pm,<br />double m00, double m01, double m02, double m03,<br />double m10, double m11, double m12, double m13,<br />double m20, double m21, double m22, double m23,<br />double m30, double m31, double m32, double m33)<br />{<br />pm->M00 = m00;<br />pm->M01 = m01;<br />pm->M02 = m02;<br />pm->M03 = m03;<br />pm->M10 = m10;<br />pm->M11 = m11;<br />pm->M12 = m12;<br />pm->M13 = m13;<br />pm->M20 = m20;<br />pm->M21 = m21;<br />pm->M22 = m22;<br />pm->M23 = m23;<br />pm->M30 = m30;<br />pm->M31 = m31;<br />pm->M32 = m32;<br />pm->M33 = m33;<br />}</p><p>void _CPPYIN_Math::MatrixAdd(MATRIX3X3_PTR ma, MATRIX3X3_PTR mb, MATRIX3X3_PTR msum)<br />{<br />for (int row = 0; row < 3; ++row)<br /> {<br />for (int col = 0; col < 3; ++col)<br /> {<br />msum->M[row][col] = ma->M[row][col] + mb->M[row][col];<br />}<br /> }<br />}</p><p>void _CPPYIN_Math::MatrixAdd(MATRIX4X4_PTR ma, MATRIX4X4_PTR mb, MATRIX4X4_PTR msum)<br />{<br />for (int row = 0; row < 4; ++row)<br /> {<br />for (int col = 0; col < 4; ++col)<br /> {<br />msum->M[row][col] = ma->M[row][col] + mb->M[row][col];<br />}<br /> }<br />}</p><p>void _CPPYIN_Math::MatrixMul(MATRIX1X3_PTR ma, MATRIX3X3_PTR mb, MATRIX1X3_PTR mprod)<br />{<br />for (int col = 0; col < 3; ++col)<br />{<br />double sum = 0;<br />for (int index = 0; index < 3; ++index)<br />{<br />sum += (ma->M[index] * mb->M[index][col]);<br />}<br /> mprod->M[col] = sum;<br />}<br />}</p><p>void _CPPYIN_Math::MatrixMul(MATRIX1X4_PTR ma, MATRIX4X4_PTR mb, MATRIX1X4_PTR mprod)<br />{<br />for (int col = 0; col < 4; ++col)<br />{<br />double sum = 0;<br />for (int index = 0; index < 4; ++index)<br />{<br />sum += (ma->M[index] * mb->M[index][col]);<br />}<br /> mprod->M[col] = sum;<br />}<br />}</p><p>void _CPPYIN_Math::MatrixMul(MATRIX3X3_PTR ma, MATRIX3X3_PTR mb, MATRIX3X3_PTR mprod)<br />{<br />for (int row = 0; row < 3; ++row)<br />{<br />for (int col = 0; col < 3; ++col)<br />{<br />double sum = 0.0;<br />for (int index = 0; index < 3; ++index)<br />{<br />sum += ma->M[row][index] * mb->M[index][col];<br />}<br />mprod->M[row][col] = sum;<br />}<br />}<br />}</p><p>void _CPPYIN_Math::MatrixMul(MATRIX4X4_PTR ma, MATRIX4X4_PTR mb, MATRIX4X4_PTR mprod)<br />{<br />for (int row = 0; row < 4; ++row)<br />{<br />for (int col = 0; col < 4; ++col)<br />{<br />double sum = 0.0;<br />for (int index = 0; index < 4; ++index)<br />{<br />sum += ma->M[row][index] * mb->M[index][col];<br />}<br />mprod->M[row][col] = sum;<br />}<br />}<br />}</p><p>double _CPPYIN_Math::MatrixDet(MATRIX3X3_PTR m)<br />{<br />return (<br />m->M00 * (m->M11 * m->M22 - m->M21 * m->M12) - m->M01 * (m->M10 * m->M22 - m->M20 * m->M12) + m->M02 * (m->M10 * m->M21 - m->M20 * m->M11)<br />);<br />}</p><p>int _CPPYIN_Math::MatrixInverse(MATRIX3X3_PTR m, MATRIX3X3_PTR mi)<br />{<br />double det = m->M00*(m->M11*m->M22 - m->M21*m->M12) -<br /> m->M01*(m->M10*m->M22 - m->M20*m->M12) +<br /> m->M02*(m->M10*m->M21 - m->M20*m->M11);</p><p>if (abs(det) < EPSILON)<br />return 0;</p><p>double det_inv = 1.0/det;</p><p>mi->M00 = det_inv*(m->M11*m->M22 - m->M21*m->M12);<br />mi->M10 = -det_inv*(m->M10*m->M22 - m->M20*m->M12);<br />mi->M20 = det_inv*(m->M10*m->M21 - m->M20*m->M11);</p><p>mi->M01 = -det_inv*(m->M01*m->M22 - m->M21*m->M02);<br />mi->M11 = det_inv*(m->M00*m->M22 - m->M20*m->M02);<br />mi->M21 = -det_inv*(m->M00*m->M21 - m->M20*m->M01);</p><p>mi->M02 = det_inv*(m->M01*m->M12 - m->M11*m->M02);<br />mi->M12 = -det_inv*(m->M00*m->M12 - m->M10*m->M02);<br />mi->M22 = det_inv*(m->M00*m->M11 - m->M10*m->M01);</p><p>return 1;<br />}</p><p>int _CPPYIN_Math::MatrixInverse(MATRIX4X4_PTR m, MATRIX4X4_PTR mi)<br />{<br />double det = ( m->M00 * ( m->M11 * m->M22 - m->M12 * m->M21 ) -<br /> m->M01 * ( m->M10 * m->M22 - m->M12 * m->M20 ) +<br /> m->M02 * ( m->M10 * m->M21 - m->M11 * m->M20 ) );</p><p>if (abs(det) < EPSILON)<br /> return 0;</p><p>double det_inv = 1.0 / det;</p><p>mi->M00 = det_inv * ( m->M11 * m->M22 - m->M12 * m->M21 );<br />mi->M01 = -det_inv * ( m->M01 * m->M22 - m->M02 * m->M21 );<br />mi->M02 = det_inv * ( m->M01 * m->M12 - m->M02 * m->M11 );<br />mi->M03 = 0.0;</p><p>mi->M10 = -det_inv * ( m->M10 * m->M22 - m->M12 * m->M20 );<br />mi->M11 = det_inv * ( m->M00 * m->M22 - m->M02 * m->M20 );<br />mi->M12 = -det_inv * ( m->M00 * m->M12 - m->M02 * m->M10 );<br />mi->M13 = 0.0;</p><p>mi->M20 = det_inv * ( m->M10 * m->M21 - m->M11 * m->M20 );<br />mi->M21 = -det_inv * ( m->M00 * m->M21 - m->M01 * m->M20 );<br />mi->M22 = det_inv * ( m->M00 * m->M11 - m->M01 * m->M10 );<br />mi->M23 = 0.0;</p><p>mi->M30 = -( m->M30 * mi->M00 + m->M31 * mi->M10 + m->M32 * mi->M20 );<br />mi->M31 = -( m->M30 * mi->M01 + m->M31 * mi->M11 + m->M32 * mi->M21 );<br />mi->M32 = -( m->M30 * mi->M02 + m->M31 * mi->M12 + m->M32 * mi->M22 );<br />mi->M33 = 1.0;</p><p>return 1;<br />}

 

 

 

3. 代碼下載

完整項目代碼下載：>>點擊進入下載頁<<