4x4矩阵求逆

最近学习自定义光栅化渲染器，在学习矩阵变换的时候，求矩阵的逆矩阵的时候遇到的一些问题，特记录在此。

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;

// 求 3x3 矩阵的逆矩阵
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;
}

// 求4x4 矩阵的逆矩阵
// m 源矩阵   mi 源矩阵的逆矩阵
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 ) );

    // 先判断行列式是否为0。
    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;
}

如上，是求3x3矩阵和4x4矩阵的逆矩阵，大家都知道矩阵的逆 A-1 = A*/|A|，不知道的可以去看下线代课本。

前面3*3矩阵求法是一般的求法，先求矩阵的行列式|A|，如果不等于0的话，则有逆矩阵。接着求A*（伴随矩阵），这个求法是先求A的代数余子式，再转置即可。

但是在求4x4矩阵的时候，前3行3列求解步骤和上面求3x3的逆是一样的，第四列，我们赋值0,0,0,1。因为在我们这个渲染器中顶点是以行向量的形式存在的，所以，第四列直接赋值，0,0,0,1，第四行才代表的是平移矩阵。如果是列向量的话，则正好相反，第四列代表平移矩阵，第四行赋值0,0,0,1。

言归正传，第四行的前3个参数的求法是不是有点诡异。是的，我们要说的就是这3个数是怎么求出来的。