(转)从零实现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. 代码下载

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

转自:http://blog.csdn.net/cppyin/archive/2011/02/09/6175733.aspx