从零实现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矩阵{union{double M[3];struct{double M00, M01, M02;};};} MATRIX1X3, *MATRIX1X3_PTR;typedef struct MATRIX1X4_TYPE // 1×4矩阵{union{double M[4];struct{double M00, M01, M02, M03;};};} MATRIX1X4, *MATRIX1X4_PTR;typedef struct MATRIX3X3_TYPE // 3×3矩阵{union{double M[3][3];struct{double M00, M01, M02;double M10, M11, M12;double M20, M21, M22;};};} MATRIX3X3, *MATRIX3X3_PTR;typedef struct MATRIX4X4_TYPE //4×4矩阵{union{double M[4][4];struct{double M00, M01, M02, M03;double M10, M11, M12, M13;double M20, M21, M22, M23;double M30, M31, M32, M33;};};} MATRIX4X4, *MATRIX4X4_PTR;

 

下面是对原向量和点的结构体的修改：

 typedef struct POINT2D_TYPE // 2D笛卡尔坐标, 2D向量{union{double M[2]; // 数组方式struct // 变量方式{double x;double y;};};} POINT2D, *POINT2D_PTR, VECTOR2D, *VECTOR2D_PTR;typedef struct POINT3D_TYPE // 3D笛卡尔坐标, 3D向量{union{double M[3];struct{double x;double y;double z;};};} POINT3D, *POINT3D_PTR, VECTOR3D, *VECTOR3D_PTR;typedef struct VECTOR4D_TYPE // 4D向量{union{double M[4];struct{double x;double y;double z;double w;};};} VECTOR4D, *VECTOR4D_PTR;

 

 2) 矩阵操作函数实现

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

 

 

 

3. 代码下载

完整项目代码下载：>>点击进入下载页<<