目标
- 定义并实现3x3矩阵的运算
- 扩展实现4x4矩阵
源代码：
实现
矩阵的数学定义
在线性代数中，矩阵就是以行和列形式组织的矩形数字块，可以使用二维数组来描述一个矩阵。使用大写加粗的字母来表示矩阵，如：M、V、P。使用小写的m[i][j]来表示矩阵中第i行第j列元素。由于在程序语言中的数组下标一般都是从0开始的，因此实际第1行第1列元素，在程序中被写作m[0][0]
|m[0][0] m[0][1] m[0][2]|
M = |m[1][0] m[1][1] m[1][2]|
|m[2][0] m[2][1] m[2][2]|
行数和列数相同的矩阵称作方阵，我们在3D程序中主要讨论的就是2x2、3x3、4x4方阵。
方阵的对角线元素就是方阵中行号和列号相同的元素。例如，3x3矩阵M的对角线元素为 m[0][0]、m[1][1]、m[2][2]。其他元素均为非对角线元素。
如果所有非对角线元素都为0，那么称这种矩阵为对角矩阵，例如：
|3 0 0 0|
|0 1 0 0|
|0 0 -5 0|
|0 0 0 2|
单位矩阵是一种特殊的对角矩阵。n 维单位矩阵记做In，是n×n矩阵，对角线元素为1，其他元素为0。例如，3x3单位矩阵：
|1 0 0|
I = |0 1 0|
|0 0 1|
Matrix3f
根据矩阵的定义，创建Matrix3f来表示3x3矩阵。
/**
* 3*3矩阵
*
* @author yanmaoyuan
*
*/
public class Matrix3f {
protected float[][] m = new float[3][3];
/**
* 单位矩阵
*/
public final static Matrix3f IDENTITY = new Matrix3f();
/**
* 初始化为单位矩阵。
*/
public Matrix3f() {
loadIdentity();
}
/**
* 使用一个3x3数组来初始化矩阵
* @param mat3
*/
public Matrix3f(float[][] mat3) {
for(int i=0; i<3; i++) {
for(int j=0; j<3; j++) {
m[i][j] = mat3[i][j];
}
}
}
/**
* 使用另一个矩阵来初始化矩阵
* @param mat3
*/
public Matrix3f(Matrix3f mat3) {
for(int i=0; i<3; i++) {
for(int j=0; j<3; j++) {
m[i][j] = mat3.m[i][j];
}
}
}
/**
* 单位矩阵
*/
public void loadIdentity() {
m[0][0] = 1; m[0][1] = 0; m[0][2] = 0;
m[1][0] = 0; m[1][1] = 1; m[1][2] = 0;
m[2][0] = 0; m[2][1] = 0; m[2][2] = 1;
}
}
Matrix4f
由于Matrix4f中的大部分方法，实现方式都与Matrix3f差不多，因此除非有特别的需要，我都不会重写一遍代码了。
转置矩阵
将矩阵内的元素沿着对角线翻折，可以得到它的转置（transpose）矩阵。
|a b c |T |a d g |
|d e f | = |b e h |
|g h i | |c f i |
Matrxi3f中的实现：
/**
* 求转置矩阵
* @return
*/
public Matrix3f transpose() {
float[][] mat3 = new float[3][3];
for(int i=0; i<3; i++) {
for(int j=0; j<3; j++) {
mat3[i][j] = m[j][i];
}
}
return new Matrix3f(mat3);
}
/**
* 求转置矩阵
* @return
*/
public Matrix3f transposeLocal() {
for(int i=0; i<3; i++) {
for(int j=i; j<3; j++) {
if (i == j)
continue;
float tmp = m[i][j];
m[i][j] = m[j][i];
m[j][i] = tmp;
}
}
return this;
}
性能测试1
我有些怀疑这两个方法的性能。对于3x3矩阵来说，实际上只需要交换3组变量的值就可以完成转置，4x4矩阵也只需要交换6组变量。有必要用双重for循环来实现这个功能吗？
为了检验性能，我另写了一个直接交换变量的实现。
/**
* 求转置矩阵
* @return
*/
public Matrix3f transposeLocal2() {
float tmp;
tmp = m[0][1];
m[0][1] = m[1][0];
m[1][0] = tmp;
tmp = m[0][2];
m[0][2] = m[2][0];
m[2][0] = tmp;
tmp = m[1][2];
m[1][2] = m[2][1];
m[2][1] = tmp;
return this;
}
再写一个测试用例：把两个方法分别执行10000次，然后比较执行时间。这个测试用例将执行1000次，以获得比较稳定的数据。
public static void main(String[] args) {
long start, time;
Matrix3f mat3 = new Matrix3f();
for (int j = 0; j<1000; j++) {
// 执行10000次双重for循环
start = System.nanoTime();
for(int i=0; i<10000; i++) {
mat3.transposeLocal();
}
time = System.nanoTime() - start;
System.out.printf("time1:%d ", time);
// 执行10000次直接交换
start = System.nanoTime();
for(int i=0; i<10000; i++) {
mat3.transposeLocal2();
}
time = System.nanoTime() - start;
System.out.printf("time2:%d\n",time);
}
}
测试结果是这样的：
time1:100825 time2:23244
time1:99013 time2:23546
time1:98108 time2:23244
time1:98410 time2:23546
time1:98712 time2:23244
time1:98108 time2:23244
time1:99013 time2:23546
time1:122258 time2:27471
time1:175688 time2:28979
time1:170255 time2:29282
time1:167839 time2:27470
time1:191084 time2:28677
time1:169047 time2:28979
time1:169651 time2:28980
time1:169047 time2:29281
time1:164821 time2:29583
time1:178707 time2:28677
time1:165726 time2:30187
这说明后一种写法的性能更高。
不过你以为这就完了吗？
性能测试2
虽然从矩阵的数学定义上来讲，我们可以用二维数组来定义矩阵。但是二维数组在使用时需要两次寻址，其效率是比不上一维数组的。即使是一维数组，使用时仍然需要进行一次寻址，执行效率比不上普通变量。
因此，我们有另一种方式来实现矩阵：
public Matrix3f {
protected float m00, m01, m02;
protected float m10, m11, m12;
protected float m20, m21, m22;
// 转置矩阵
public Matrix3f transposeLocal3() {
float tmp = m01;
m01 = m10;
m10 = tmp;
tmp = m02;
m02 = m20;
m20 = tmp;
tmp = m12;
m12 = m21;
m21 = tmp;
return this;
}
}
这个方法与方法2相比，谁更快呢？让我们改一下测试用例，把对第3个方法的测试加进去。
public static void main(String[] args) {
long start, time;
Matrix3f mat3 = new Matrix3f();
for (int j = 0; j<1000; j++) {
// 执行10000次双重for循环
start = System.nanoTime();
for(int i=0; i<10000; i++) {
mat3.transposeLocal();
}
time = System.nanoTime() - start;
System.out.printf(" time1:%d ", time);
// 执行10000次直接交换
start = System.nanoTime();
for(int i=0; i<10000; i++) {
mat3.transposeLocal2();
}
time = System.nanoTime() - start;
System.out.printf("time2:%d ",time);
// 执行10000次交换，不用数组
start = System.nanoTime();
for(int i=0; i<10000; i++) {
mat3.transposeLocal3();
}
time = System.nanoTime() - start;
System.out.printf("time3:%d\n",time);
}
}
结果的一部分如下：
time1:93580 time2:20527 time3:2415
time1:124672 time2:29282 time3:3925
time1:107768 time2:21734 time3:2113
time1:85430 time2:20829 time3:2415
time1:96599 time2:26867 time3:4226
time1:96599 time2:22640 time3:4226
time1:95391 time2:26564 time3:3924
time1:145501 time2:31998 time3:4226
time1:143388 time2:28980 time3:3924
time1:144898 time2:28376 time3:3925
time1:173575 time2:40451 time3:5434
time1:147011 time2:27168 time3:3622
time1:131616 time2:27470 time3:3924
time1:131615 time2:27168 time3:3623
time1:132823 time2:27772 time3:3924
time1:133124 time2:27772 time3:3924
time1:134936 time2:35620 time3:3622
time1:150332 time2:24452 time3:3924
time1:94787 time2:20528 time3:2113
time1:91467 time2:23546 time3:2113
time1:91466 time2:20527 time3:2415
time1:91466 time2:20829 time3:2113
相比于使用二维数组，第2种方式提升了4~5倍的性能，而第3种方式提升了40~50倍性能。
不过第三种方式也导致了一些问题，就是代码稍微显得有些冗余。不过我认为这是可以接受。
从下一个方法开始，我会把现有代码改成第三种方式来实现。
性能优化版本
/**
* 3*3矩阵
*
* @author yanmaoyuan
*
*/
public class Matrix3f {
protected float m00, m01, m02;
protected float m10, m11, m12;
protected float m20, m21, m22;
/**
* 零矩阵
*/
public static final Matrix3f ZERO = new Matrix3f(0, 0, 0, 0, 0, 0, 0, 0, 0);
/**
* 单位矩阵
*/
public final static Matrix3f IDENTITY = new Matrix3f();
/**
* 初始化为单位矩阵。
*/
public Matrix3f() {
loadIdentity();
}
/**
* 初始化矩阵
* @param m00
* @param m01
* @param m02
* @param m10
* @param m11
* @param m12
* @param m20
* @param m21
* @param m22
*/
public Matrix3f(float m00, float m01, float m02, float m10, float m11,
float m12, float m20, float m21, float m22) {
this.m00 = m00;
this.m01 = m01;
this.m02 = m02;
this.m10 = m10;
this.m11 = m11;
this.m12 = m12;
this.m20 = m20;
this.m21 = m21;
this.m22 = m22;
}
/**
* 使用另一个矩阵来初始化矩阵
* @param mat3
*/
public Matrix3f(Matrix3f mat3) {
set(mat3);
}
/**
* 复制另一个矩阵的值
* @param matrix
* @return
*/
public Matrix3f set(Matrix3f matrix) {
if (null == matrix) {
loadIdentity();
} else {
m00 = matrix.m00;
m01 = matrix.m01;
m02 = matrix.m02;
m10 = matrix.m10;
m11 = matrix.m11;
m12 = matrix.m12;
m20 = matrix.m20;
m21 = matrix.m21;
m22 = matrix.m22;
}
return this;
}
/**
* 单位矩阵
*/
public void loadIdentity() {
m01 = m02 = m10 = m12 = m20 = m21 = 0;
m00 = m11 = m22 = 1;
}
/**
* 求转置矩阵
* @return
*/
public Matrix3f transpose() {
return new Matrix3f(m00, m10, m20, m01, m11, m21, m02, m12, m22);
}
/**
* 求转置矩阵
* @return
*/
public Matrix3f transposeLocal() {
float tmp = m01;
m01 = m10;
m10 = tmp;
tmp = m02;
m02 = m20;
m20 = tmp;
tmp = m12;
m12 = m21;
m21 = tmp;
return this;
}
}
标量乘法
矩阵M可以和标量k相乘，只要把M中的每个元素都与k相乘即可。
/**
* 标量乘法
* @param scalor
* @return
*/
public Matrix3f multLocal(float scalor) {
m00 *= scalor;
m01 *= scalor;
m02 *= scalor;
m10 *= scalor;
m11 *= scalor;
m12 *= scalor;
m20 *= scalor;
m21 *= scalor;
m22 *= scalor;
return this;
}
矩阵之间的乘法
矩阵乘法的规则很多，网上资料也很多，我就不多写了。下面的代码是照搬 jMonkeyEngine 的源码。
Matrix3f与Matrix3f相乘，结果仍然是一个Matrix3f对象。
/**
* 矩阵乘法
*
* @param mat
* @param product 结果矩阵
* @return 返回一个新的Matrix3f对象，携带计算结果。
*/
public Matrix3f mult(Matrix3f mat, Matrix3f product) {
float temp00, temp01, temp02;
float temp10, temp11, temp12;
float temp20, temp21, temp22;
if (product == null) {
product = new Matrix3f();
}
temp00 = m00 * mat.m00 + m01 * mat.m10 + m02 * mat.m20;
temp01 = m00 * mat.m01 + m01 * mat.m11 + m02 * mat.m21;
temp02 = m00 * mat.m02 + m01 * mat.m12 + m02 * mat.m22;
temp10 = m10 * mat.m00 + m11 * mat.m10 + m12 * mat.m20;
temp11 = m10 * mat.m01 + m11 * mat.m11 + m12 * mat.m21;
temp12 = m10 * mat.m02 + m11 * mat.m12 + m12 * mat.m22;
temp20 = m20 * mat.m00 + m21 * mat.m10 + m22 * mat.m20;
temp21 = m20 * mat.m01 + m21 * mat.m11 + m22 * mat.m21;
temp22 = m20 * mat.m02 + m21 * mat.m12 + m22 * mat.m22;
product.m00 = temp00;
product.m01 = temp01;
product.m02 = temp02;
product.m10 = temp10;
product.m11 = temp11;
product.m12 = temp12;
product.m20 = temp20;
product.m21 = temp21;
product.m22 = temp22;
return product;
}
/**
* 矩阵乘法
* @param mat3
* @return
*/
public Matrix3f mult(Matrix3f mat) {
return mult(mat, null);
}
/**
* 矩阵乘法
* @param mat
* @return
*/
public Matrix3f multLocal(Matrix3f mat) {
return mult(mat, this);
}
矩阵和向量的乘法
Matrix3f与Vector3f相乘，结果是一个Vector3f对象。这个乘法的几何意义，是使用3x3矩阵对三维向量进行线性变换（旋转、缩放）。但是3x3矩阵无法让三维向量平移。
/**
* 计算矩阵与vec的乘积，结果保存在新的Vector3f对象中。
*
* @param vec
* @return
*/
public Vector3f mult(Vector3f vec) {
return mult(vec, null);
}
/**
* 计算矩阵与vec的乘积，结果保存在product对象中。
*
* @param vec
* @param product
* @return
*/
public Vector3f mult(Vector3f vec, Vector3f product) {
if (null == product) {
product = new Vector3f();
}
float x = vec.x;
float y = vec.y;
float z = vec.z;
product.x = m00 * x + m01 * y + m02 * z;
product.y = m10 * x + m11 * y + m12 * z;
product.z = m20 * x + m21 * y + m22 * z;
return product;
}
/**
* 计算矩阵与向量vec的乘积，结果保存在该vec中。
*
* @param vec
* @return
*/
public Vector3f multLocal(Vector3f vec) {
if (vec == null) {
return null;
}
float x = vec.x;
float y = vec.y;
vec.x = m00 * x + m01 * y + m02 * vec.z;
vec.y = m10 * x + m11 * y + m12 * vec.z;
vec.z = m20 * x + m21 * y + m22 * vec.z;
return vec;
}
若把Vector3f“扩展”成一个Vector4f，就可以用它与一个Matrix4f相乘。通过四位空间的“切变”使空间“扭曲”，从而影响原Vector3f的位置。
需要注意的是，将Vector3f扩展为Vector4f时，若第4个分量的值为1，则该向量会计算位移。若第4个分量的值为0，则只会计算旋转和缩放。
注意下面的代码，Vector3f与Matrix4f相乘时，只比Matrix3f多加了第4列的三个分量，从而产生了“位移”。这其实是把 vec3 当做 vec4(vec3, 1.0) 来看待的。
/**
* 4x4矩阵与Vector3f向量相乘，结果保存在一个新的store对象中。
*
* @param vec
* @param store
* @return
*/
public Vector3f mult(Vector3f vec, Vector3f store) {
if (store == null) {
store = new Vector3f();
}
float vx = vec.x, vy = vec.y, vz = vec.z;
store.x = m00 * vx + m01 * vy + m02 * vz + m03;
store.y = m10 * vx + m11 * vy + m12 * vz + m13;
store.z = m20 * vx + m21 * vy + m22 * vz + m23;
return store;
}
/**
* 4x4矩阵与Vector3f向量相乘，结果保存在一个新的Vector3f对象中。
*
* @param vec
* @return
*/
public Vector3f mult(Vector3f vec) {
return mult(vec, null);
}
用Matrix4f于Vector4f的相乘来做对比，就更明显了。当v.w = 1.0时，计算出来的x、y、z就和上面的方法一模一样。
/**
* 矩阵与四位向量相乘，结果返保存在store对象中。
*
* @param vec
* @param store
* @return
*/
public Vector4f mult(Vector4f vec, Vector4f store) {
if (null == vec) {
return null;
}
if (store == null) {
store = new Vector4f();
}
float vx = vec.x, vy = vec.y, vz = vec.z, vw = vec.w;
store.x = m00 * vx + m01 * vy + m02 * vz + m03 * vw;
store.y = m10 * vx + m11 * vy + m12 * vz + m13 * vw;
store.z = m20 * vx + m21 * vy + m22 * vz + m23 * vw;
store.w = m30 * vx + m31 * vy + m32 * vz + m33 * vw;
return store;
}
/**
* 矩阵与四位向量相乘，结果返回一个新的Vector4f对象。
*
* @param vec
* @return
*/
public Vector4f mult(Vector4f vec) {
return mult(vec, null);
}
矩阵的行列式
矩阵的行列式是一个标量，表示物体在经过该矩阵进行空间变换后，体积的变化。若行列式的值为1，说明体积没有变化。
/**
* 计算Matrix3f的行列式
* @return
*/
public float determinant() {
float fCo00 = m11 * m22 - m12 * m21;
float fCo10 = m12 * m20 - m10 * m22;
float fCo20 = m10 * m21 - m11 * m20;
float fDet = m00 * fCo00 + m01 * fCo10 + m02 * fCo20;
return fDet;
}
矩阵的逆
矩阵M的逆记为M^-1，矩阵乘以自身的逆，结果恰好是单位矩阵I。
M * M^-1 = I
从几何意义上来说，矩阵M的作用是让空间变形，而矩阵的逆则可以让空间恢复成原样。并非每一个矩阵都有逆。考虑这样一种变形：一个三维物体被“挤压”成了一个二维平面。
这种变形是无法恢复的，因为压缩后物体表面的“高度”信息完全丢失了，称为维度“坍缩”。将三维物体投影到二维平面很容易，但是根据其投影复原三维物体则几乎不可能。更不用说“坍缩”成一条线段或一个“奇点”了。
对于这类将会导致维度“坍缩”的变形矩阵，其主要特征是行列式的值为0，主要发生在矩阵某一列（或多列）为零向量时，称为奇异矩阵。因此，如果矩阵的行列式为0，这个矩阵是没有逆矩阵的。此外，对于任意可逆矩阵M，当且仅当向量 v 为零向量时，vM=0。
M的“标准伴随矩阵”，记做“adj M”，定义为 M 的代数余子式矩阵的转置矩阵。
/**
* 求Matrix3f的标准伴随矩阵，结果保存为store对象。
*
* @param store
* @return
*/
public Matrix3f adjoint(Matrix3f store) {
if (store == null) {
store = new Matrix3f();
}
store.m00 = m11 * m22 - m12 * m21;
store.m01 = m02 * m21 - m01 * m22;
store.m02 = m01 * m12 - m02 * m11;
store.m10 = m12 * m20 - m10 * m22;
store.m11 = m00 * m22 - m02 * m20;
store.m12 = m02 * m10 - m00 * m12;
store.m20 = m10 * m21 - m11 * m20;
store.m21 = m01 * m20 - m00 * m21;
store.m22 = m00 * m11 - m01 * m10;
return store;
}
/**
* 求Matrix3f的标准伴随矩阵。
*
* @return
*/
public Matrix3f adjoint() {
return adjoint(null);
}
一旦有了伴随矩阵，通过除以M的行列式，就能计算矩阵的逆。
/**
* 求矩阵的逆
*
* @return
*/
public Matrix3f invert() {
return invert(null);
}
/**
* 求矩阵的逆
*
* @return The store
*/
public Matrix3f invert(Matrix3f store) {
if (store == null) {
store = new Matrix3f();
}
float det = determinant();
if (Math.abs(det) <= 0f) {
return store.zero();
}
store.m00 = m11 * m22 - m12 * m21;
store.m01 = m02 * m21 - m01 * m22;
store.m02 = m01 * m12 - m02 * m11;
store.m10 = m12 * m20 - m10 * m22;
store.m11 = m00 * m22 - m02 * m20;
store.m12 = m02 * m10 - m00 * m12;
store.m20 = m10 * m21 - m11 * m20;
store.m21 = m01 * m20 - m00 * m21;
store.m22 = m00 * m11 - m01 * m10;
store.multLocal(1f / det);
return store;
}
/**
* 求矩阵的逆
*
* @return this
*/
public Matrix3f invertLocal() {
float det = determinant();
if (Math.abs(det) <= 0f) {
return zero();
}
float f00 = m11 * m22 - m12 * m21;
float f01 = m02 * m21 - m01 * m22;
float f02 = m01 * m12 - m02 * m11;
float f10 = m12 * m20 - m10 * m22;
float f11 = m00 * m22 - m02 * m20;
float f12 = m02 * m10 - m00 * m12;
float f20 = m10 * m21 - m11 * m20;
float f21 = m01 * m20 - m00 * m21;
float f22 = m00 * m11 - m01 * m10;
m00 = f00;
m01 = f01;
m02 = f02;
m10 = f10;
m11 = f11;
m12 = f12;
m20 = f20;
m21 = f21;
m22 = f22;
multLocal(1f / det);
return this;
}
总结
本文是实现了矩阵的基本性质，但是还有很多运算都还没有实现。
这些运算主要是不同类型的空间变换。比如：绕轴旋转矩阵、缩放矩阵、平移矩阵等。还有Quaternion与Matrix3f的互相转换，Matrix4f与Matrix3f的转换等。我想把这些内容留到下一节“空间变换”再来继续实现。