Java软光栅渲染器-矩阵

目标

  • 定义并实现3x3矩阵的运算
  • 扩展实现4x4矩阵

源代码:

实现

矩阵的数学定义

在线性代数中,矩阵就是以行和列形式组织的矩形数字块,可以使用二维数组来描述一个矩阵。使用大写加粗的字母来表示矩阵,如:MVP。使用小写的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的转换等。我想把这些内容留到下一节“空间变换”再来继续实现。