谈一谈3D编程中的矩阵
常见约定行主(row major)和列主(column major)行主列主3D编程中的例子行向量左乘矩阵和列向量右乘矩阵左乘右乘为什么要约定左乘还是右乘？对矩阵构造的影响对矩阵串接的影响坐标系变换矩阵的构造向量空间和向量空间的基坐标空间的变换矩阵思考推导结论验证看透变换矩阵的小黑盒向量变换矩阵逆矩阵计算正交矩阵的逆矩阵矩阵乘积的逆矩阵一般矩阵的逆矩阵计算矩阵变换法线通用的正确变换方法：使用逆转置矩阵推导正交矩阵可直接使用原变换矩阵3D场景中物体的法线变换常见约定
谈论矩阵之前，要先明确一下使用的约定。约定不同，用法迥异。
行主(row major)和列主(column major)
图形API或shader语言中常常把矩阵存储为一维数组，这就带来了一个问题，按什么顺序将一维数组中的元素填入到矩阵中，以及要访问矩阵的元素应该指定什么样的下标。
行主
如果遍历一维数组的元素，然后依次填入矩阵的每一行，这就是行主。以CG语言为例：
float3x3 mat = float3x3(1.1, 1.2, 1.3, 2.1, 2.2, 2.3, 3.1, 3.2, 3.3);
以上的代码构建出的矩阵是这样的：
[1.11.21.32.12.22.33.13.23.3]\left[\begin{matrix} 1.1 & 1.2 & 1.3 \\ 2.1 & 2.2 & 2.3 \\ 3.1 & 3.2 & 3.3 \end{matrix}\right] ⎣⎡1.12.13.11.22.23.21.32.33.3⎦⎤
而如果要访问第1行第2列的元素，应该使用mat[0][1], 下标指定时先行后列（索引从0开始）。
列主
按顺序填入元素时，是按列填入的，则是列主。以glsl语言为例：
mat3 mat = mat3(1.1, 2.1, 3.1, 1.2, 2.2, 3.2, 1.3, 2.3, 3.3);
以上代码构建出的矩阵同样是：
[1.11.21.32.12.22.33.13.23.3]\left[\begin{matrix} 1.1 & 1.2 & 1.3 \\ 2.1 & 2.2 & 2.3 \\ 3.1 & 3.2 & 3.3 \end{matrix}\right] ⎣⎡1.12.13.11.22.23.21.32.33.3⎦⎤
但是一维数组（或者说参数列表）中的顺序和之前CG的完全不同，每个元素都是依次从第一列开始填入，填满后填第二列。而在glsl中，要访问第1行第2列的元素，则应该使用mat[1][0],下标指定时先列后行。
3D编程中的例子
不光是shader中需要分清楚行主列主的约定，在3D程序代码中使用矩阵时，也是要注意这个约定，无论是使用3D引擎的数学库的时候，还是自己写3D程序时自行撰写Matrix类。以笔者的开源项目mini3d.js中的Matrix4类为例。这个矩阵类使用列主的约定，其中设置位移的方法：
setTranslate(x,y,z){let e = this.elements;e[0] = 1; e[4] = 0; e[8] = 0; e[12] = x;e[1] = 0; e[5] = 1; e[9] = 0; e[13] = y;e[2] = 0; e[6] = 0; e[10] = 1; e[14] = z;e[3] = 0; e[7] = 0; e[11] = 0; e[15] = 1;return this;}
这是一个4x4齐次变换矩阵，因为是列主约定，位置向量被设置到一维数组的12,13,14三个下标位置中。在后期从世界矩阵获取世界坐标时，遵循这一约定，从12,13,14三个下标位置取出世界坐标。不过如果数学库封装的足够好，对于一般的应用，倒是不会直接去获取元素，但总有一些特殊的需求，需要自己构建矩阵，此时必须知道所使用引擎数学库的约定。
行向量左乘矩阵和列向量右乘矩阵
矩阵和向量相乘时，也就是使用矩阵变换向量时，按照向量所在的位置，可以分为左乘和右乘两种约定。
左乘
向量在左边，比较符合从左向右的阅读顺序，例如：
float3 vec2 = mul(vec1, someMat3);
如果vec1是[1,0,0], someMat3是上面的3x3矩阵，左乘就是这样写：
[1,0,0]∗[1.11.21.32.12.22.33.13.23.3][1, 0, 0]* \left[\begin{matrix} 1.1 & 1.2 & 1.3 \\ 2.1 & 2.2 & 2.3 \\ 3.1 & 3.2 & 3.3 \end{matrix}\right] [1,0,0]∗⎣⎡1.12.13.11.22.23.21.32.33.3⎦⎤
可以看到，左乘时，向量要写成行向量的形式，才能满足矩阵乘法的规则。
右乘
向量在右边，例如：
float3 vec2 = mul(someMat3, vec1);
还是上面的向量和矩阵，右乘就是这样写：
[1.11.21.32.12.22.33.13.23.3]∗[100]\left[\begin{matrix} 1.1 & 1.2 & 1.3 \\ 2.1 & 2.2 & 2.3 \\ 3.1 & 3.2 & 3.3 \end{matrix}\right] * \left[\begin{matrix} 1 \\ 0 \\ 0 \end{matrix}\right] ⎣⎡1.12.13.11.22.23.21.32.33.3⎦⎤∗⎣⎡100⎦⎤
可以看到，右乘时，向量要写成列向量的形式。
为什么要约定左乘还是右乘？
在CG/HLSL/GLSL的shader中，向量和矩阵之间相乘既可以写成左乘也可以写成右乘。而一般的3D引擎会提供一个函数去用矩阵变换向量，或者对于一些C++引擎，会使用运算符重载，一般而言你没有机会在调用时放错向量的位置。但是是左乘还是右乘，决定了矩阵应该如何构造和串接。
对矩阵构造的影响
以上面的mini3d.js中的位移矩阵为例，这个矩阵约定是要使用向量右乘，因此位移矢量<Tx,Ty,Tz><T_x, T_y, T_z><Tx,Ty,Tz>位于矩阵的第四列，这个矩阵乘以列向量的时候是这样的：
[100Tx010Ty001Tz0001]∗[2341]=[2+Tx3+Ty4+Tz1]\left[\begin{matrix} 1 & 0 & 0 & T_x \\ 0 & 1 & 0 & T_y \\ 0 & 0 & 1 & T_z \\ 0 & 0 & 0 & 1 \end{matrix}\right] * \left[\begin{matrix} 2 \\ 3 \\ 4 \\ 1 \end{matrix}\right] = \left[\begin{matrix} 2+T_x \\ 3+T_y \\ 4+T_z \\ 1 \end{matrix}\right] ⎣⎢⎢⎡100001000010TxTyTz1⎦⎥⎥⎤∗⎣⎢⎢⎡2341⎦⎥⎥⎤=⎣⎢⎢⎡2+Tx3+Ty4+Tz1⎦⎥⎥⎤
这样正确的将齐次坐标<2,3,4,1><2,3,4,1><2,3,4,1>偏移了<Tx,Ty,Tz><T_x,T_y,T_z><Tx,Ty,Tz>。
如果这个地方将<Tx,Ty,Tz><T_x,T_y,T_z><Tx,Ty,Tz>放在矩阵的最后一行，结果就完全不对了。
有兴趣的读者可以自己乘一下试试。
那么这儿可能有个误解，就是列主矩阵用右乘，行主矩阵用左乘。其实行主列主和左乘右乘并没有关系，不需要强行绑定。例如在CG中，矩阵是行主构造的，但是可以随便你使用左乘还是右乘，并没有限制。行主列主只是决定了矩阵如何存储在内存中，而左乘右乘才会决定矩阵的样子，因为要适应向量的位置。不过有意思的是，列主加右乘的约定和行主加左乘的约定，对于矩阵的一维数组存储来说是会得到相同的一维数组。这是因为同样的变换，使用左乘约定和右乘约定的矩阵是互为转置的。而同一个一维数组按行主和按列主构建的矩阵也是互为转置的。
对矩阵串接的影响
使用矩阵往往需要把变换矩阵相乘，即串接起来，然后再使用相乘的结果去变换向量，这样大大减少了矩阵向量乘法的次数。矩阵相乘是不可交换顺序的，比如我们要把位移T，旋转R和缩放S三个变换矩阵相乘。如何决定矩阵相乘的顺序呢？这就要看向量的位置了。靠近向量的矩阵先起作用。例如，在左乘约定下：
P∗S∗R∗TP * S * R * T P∗S∗R∗T
对于点向量P先缩放，再旋转，再平移。
而在右乘约定下，同样是先缩放，再旋转，最后平移，矩阵的串接顺序是这样的：
T∗R∗S∗PT * R * S * P T∗R∗S∗P
因此矩阵串接的顺序取决与向量的位置。
坐标系变换矩阵的构造
3D编程中坐标系的变换特别常见，比如我们常用的物体空间到世界空间变换，世界空间到物体空间变换。那么我们如何轻松构造出坐标系变换矩阵，或者说给你一个变换矩阵，能从中得出什么信息呢？先介绍一点数学知识。
向量空间和向量空间的基
线性无关和线性相关：对于一组向量[e1,e2,...,en][e_1,e_2, ... , e_n][e1,e2,...,en],如果不存在一组实数a1,a2,...ana_1,a_2,...a_na1,a2,...an,其中至少有一个aia_iai不为零，使得下式成立：
a1e1+a2e2+...+anen=0a_1e_1+a_2e_2+ ... + a_ne_n = 0a1e1+a2e2+...+anen=0
则称这组向量线性无关，否则称这组向量线性相关
什么意思呢？通俗的说，一组向量如果是线性无关的，那么其中任何一个向量不能用剩下的向量通过线性运算（线性运算是加法和数量乘法）得到。比如3D坐标系的x,y,z三个轴向量，就是线性无关的。一个n维向量空间可以由n个线性无关的向量组成的向量集合进行线性组合来生成，这n个线性无关的向量组合成为该n维空间的基。
还是以3D坐标系的x,y,z三个轴为例，它们是线性无关的，且可以组合得到3D空间下的任意向量，因此x,y,z轴是3D坐标系的基。正交基：如果一组基向量两两之间点积为0，则该基为正交基
3D坐标系的x,y,z就是一组正交基。
坐标空间的变换矩阵
假设我们有坐标系A和坐标系B，已知坐标系A的三个坐标轴在坐标系B中的向量为XA,YA,ZAX_A,Y_A,Z_AXA,YA,ZA，坐标系A的原点在坐标系B中的坐标为O=<Ox,Oy,Oz>O=<O_x,O_y,O_z>O=<Ox,Oy,Oz>求坐标系A到坐标系B的变换矩阵。注意：使用列向量右乘约定。
思考
坐标之所以能被定义出来，就是因为有个参考，对于坐标系A，其所有向量是根据A的三个轴定义出来的，其实就是上面向量空间的概念，三个轴是基向量，所有的向量都能使用基向量线性表出。而点可以看成从原点出发的向量的端点，这也是一个线性运算，因此三个轴再加上原点就能线性表出坐标系中的所有的点。当我们在B坐标系中表示A坐标系的三个轴和原点时，可以认为B是A的父坐标空间。坐标系A的三个轴本身也是坐标系B中的向量，是由坐标系B的三个轴线性表出的。这其实具有传递性。子空间的向量由子空间的轴表示，而子空间的轴作为父空间中的一个向量，又由父空间的轴表示。
推导
坐标系A到坐标系B的变换矩阵的作用，是将坐标系A中的向量和点在坐标系B中表示。我们先看向量。根据上面的数学知识，向量可以使用向量空间的基线性表出，设向量为PA=<a,b,c>P_A=<a,b,c>PA=<a,b,c>是坐标系A中的向量，则：
PA=aX+bY+cZP_A = aX + bY+cZ PA=aX+bY+cZ
其中X,Y,ZX,Y,ZX,Y,Z是A坐标系的x,y,z三个轴，即<1,0,0>,<0,1,0>,<0,0,1><1,0,0>,<0,1,0>,<0,0,1><1,0,0>,<0,1,0>,<0,0,1>。
坐标系A的三个坐标轴在坐标系B中的向量为XA,YA,ZAX_A,Y_A,Z_AXA,YA,ZA，将其代入上式，那么PPP在坐标系B中的表示PBP_BPB，可以用下式得到：
PB=aXA+bYA+cZAP_B = aX_A+ bY_A + cZ_A PB=aXA+bYA+cZA
我们按分量展开，并且将a,b,ca,b,ca,b,c写到右边：
{PBx=XAxa+YAxb+ZAxcPBy=XAya+YAyb+ZAycPBz=XAza+YAzb+ZAzc\left\{ \begin{aligned} P_{B_x} = X_{A_x}a + Y_{A_x}b + Z_{A_x}c \\ P_{B_y} = X_{A_y}a + Y_{A_y}b + Z_{A_y}c \\ P_{B_z} = X_{A_z}a + Y_{A_z}b + Z_{A_z}c \end{aligned} \right. ⎩⎪⎨⎪⎧PBx=XAxa+YAxb+ZAxcPBy=XAya+YAyb+ZAycPBz=XAza+YAzb+ZAzc
上式的右边可以写成向量的点积：
{PBx=<XAx,YAx,ZAx>∗<a,b,c>PBy=<XAy,YAy,ZAy>∗<a,b,c>PBz=<XAz,YAz,ZAz>∗<a,b,c>\left\{ \begin{aligned} P_{B_x} = <X_{A_x}, Y_{A_x}, Z_{A_x}> *<a,b,c> \\ P_{B_y} = <X_{A_y}, Y_{A_y}, Z_{A_y}> * <a,b,c> \\ P_{B_z} = <X_{A_z}, Y_{A_z}, Z_{A_z}> * <a,b,c> \end{aligned} \right. ⎩⎪⎨⎪⎧PBx=<XAx,YAx,ZAx>∗<a,b,c>PBy=<XAy,YAy,ZAy>∗<a,b,c>PBz=<XAz,YAz,ZAz>∗<a,b,c>
因此可以转换成矩阵：
[PBxPByPBz]=[XAxYAxZAxXAyYAyZAyXAzYAzZAz][abc]\left[\begin{matrix} P_{B_x} \\ P_{B_y} \\ P_{B_z} \end{matrix}\right] = \left[\begin{matrix} X_{A_x} & Y_{A_x} & Z_{A_x} \\ X_{A_y} & Y_{A_y} & Z_{A_y} \\ X_{A_z} & Y_{A_z} & Z_{A_z} \end{matrix}\right] \left[\begin{matrix} a \\ b \\ c \end{matrix}\right] ⎣⎡PBxPByPBz⎦⎤=⎣⎡XAxXAyXAzYAxYAyYAzZAxZAyZAz⎦⎤⎣⎡abc⎦⎤
那么我们得到了，列向量右乘约定下，将A中的向量转换到B的矩阵。这个矩阵是将A的在B中表示的三个坐标轴按列填入矩阵得到。
为了变换点，需要将原点加入，这需要使用齐次坐标和4X4齐次变换矩阵，具体就不推导了，直接看结论。
结论
将XA,YA,ZAX_A, Y_A, Z_AXA,YA,ZA分别填入矩阵的前三列，将Ox,Oy,OzO_x,O_y,O_zOx,Oy,Oz填入矩阵的第4列，得到的4x4矩阵MA−>BM_{A->B}MA−>B就是从坐标系A变换到坐标系B的矩阵。使用这个矩阵能将坐标系A中的点和向量变换到坐标系B中，也即将相对于坐标系A的坐标值变换成相对于坐标系B的坐标值。
MA−>B=[XAxYAxZAxOxXAyYAyZAyOyXAzYAzZAzOz0001]M_{A->B} = \left[\begin{matrix} X_{A_x} & Y_{A_x} & Z_{A_x} & O_x \\ X_{A_y} & Y_{A_y} & Z_{A_y} & O_y \\ X_{A_z} & Y_{A_z} & Z_{A_z} & O_z \\ 0 & 0 & 0 & 1 \end{matrix}\right] MA−>B=⎣⎢⎢⎡XAxXAyXAz0YAxYAyYAz0ZAxZAyZAz0OxOyOz1⎦⎥⎥⎤
验证
将矩阵MA−>BM_{A->B}MA−>B乘以单位向量<1,0,0,0><1,0,0,0><1,0,0,0>，即将坐标系A中的X轴变换到坐标系B中。
MA−>B∗[1000]=[XAxYAxZAxOxXAyYAyZAyOyXAzYAzZAzOz0001]∗[1000]=[XAxXAyXAz0]M_{A->B} * \left[\begin{matrix} 1 \\ 0 \\ 0 \\ 0 \end{matrix}\right] = \left[\begin{matrix} X_{A_x} & Y_{A_x} & Z_{A_x} & O_x \\ X_{A_y} & Y_{A_y} & Z_{A_y} & O_y \\ X_{A_z} & Y_{A_z} & Z_{A_z} & O_z \\ 0 & 0 & 0 & 1 \end{matrix}\right]* \left[\begin{matrix} 1 \\ 0 \\ 0 \\ 0 \end{matrix}\right]= \left[\begin{matrix} X_{A_x} \\ X_{A_y} \\ X_{A_z} \\ 0 \end{matrix}\right] MA−>B∗⎣⎢⎢⎡1000⎦⎥⎥⎤=⎣⎢⎢⎡XAxXAyXAz0YAxYAyYAz0ZAxZAyZAz0OxOyOz1⎦⎥⎥⎤∗⎣⎢⎢⎡1000⎦⎥⎥⎤=⎣⎢⎢⎡XAxXAyXAz0⎦⎥⎥⎤
得到的向量正是A坐标系的X轴在B坐标系中表示。
同样如果我们用该矩阵变换<0,1,0><0,1,0><0,1,0>和<0,0,1><0,0,1><0,0,1>，会分别得到A坐标系的Y轴和Z轴。而如果我们用该矩阵变换A坐标系中的原点<0,0,0,1><0,0,0,1><0,0,0,1>，得到的就是A坐标系的原点在B中的坐标点<Ox,Oy,Oz,1><O_x,O_y,O_z,1><Ox,Oy,Oz,1>。
看透变换矩阵的小黑盒
知道了上面的结论之后，给你一个从A空间到B空间的变换矩阵，提取它的前三列并归一化（因为可能存在缩放）就得到了A空间的坐标轴在B空间中的向量。提取它的第四列，就得到了A空间的原点在B空间中的坐标。如此，矩阵不再是黑盒子。
向量变换矩阵
在3D中，我们经常需要变换各种向量，如光线方向，视角方向，法线，切线等等。在前面的推导过程中，我们得到了变换向量的坐标空间3x3变换矩阵，只要将空间A在空间B中的三个轴按列填入矩阵，就得到了从A到B的变换矩阵。例如，下面的glsl shader代码构建了一个从切线空间到世界空间的变换矩阵。
//TBN向量按列放入矩阵，构造出 TangentToWorld矩阵//注意：glsl矩阵是列主的mat3 tangentToWorld = mat3(worldTangent, worldBinormal, worldNormal);
填入矩阵的三列是世界空间下的TBN向量，因此该矩阵是从切线空间到世界空间的向量变换矩阵。
那么我们如果想得到世界空间到切线空间的向量变换矩阵，应该怎么办呢？按照上面的思路，只要找到世界空间的坐标轴在切线空间中的表示，然而世界空间是切线空间的父空间，所以除非首先有世界到切线空间的变换矩阵，否则无法得到世界空间的坐标轴在切线空间中的表示，这显然陷入了死循环。但是我们可以换个思路，只要求得tangentToWorld 矩阵的逆矩阵就可以了。而tangentToWorld 矩阵是一个正交矩阵，因此其逆矩阵等于其转置。也就是说，将世界空间的TBN向量按行放入矩阵，构造出worldToTangent矩阵：
//将TBN向量按行放入矩阵，构造出worldToTangent矩阵//注意glsl中mat3是列主的mat3 worldToTangent = mat3(worldTangent.x, worldBinormal.x, worldNormal.x,worldTangent.y, worldBinormal.y, worldNormal.y, worldTangent.z, worldBinormal.z, worldNormal.z);
逆矩阵计算
上面讨论到了逆矩阵，本节讨论一下3D引擎和shader中常见的逆矩阵计算方法。
正交矩阵的逆矩阵
正交矩阵的逆矩阵等于其转置矩阵。什么样的矩阵是正交矩阵呢？
首先，根据正交矩阵的定义，正交矩阵MMM的和它的转置矩阵MTM^TMT的乘积是单位阵：MMT=IMM^T = IMMT=I。那么计算这个等式是否成立就可以判断了。不过我们一般不会去判断某矩阵是否是正交矩阵。因为根据矩阵的构造方式，我们就可以知道它是否是正交矩阵。正交矩阵一般是由一组标准正交基构造。例如上面的TBN矩阵，它的三个轴T,B,NT,B,NT,B,N是一组标准正交基。再比如旋转矩阵，因为坐标轴只进行了旋转，三个轴仍然是互相正交的，且长度为1，所以也是正交矩阵。
矩阵乘积的逆矩阵
根据如下性质：
(AB)−1=B−1A−1(AB)^{-1} = B^{-1}A^{-1} (AB)−1=B−1A−1
可计算矩阵乘积的逆矩阵，当A,B是正交矩阵或可根据含义推导出逆矩阵时，这个方法是可用的。例如我们推导camera的view矩阵时，view矩阵其实是camera的世界矩阵的逆矩阵，而camera本身只有旋转和平移，因此其世界矩阵是：
Mcamera=TcameraMUVNM_{camera} = T_{camera}M_{UVN} Mcamera=TcameraMUVN
其中MUVNM_{UVN}MUVN是使用Camera本地坐标系的三个轴UVN在世界空间的表示构造的矩阵（和我们上面的方法一样），这是一个正交矩阵。而平移矩阵TcameraT_{camera}Tcamera的逆矩阵可以简单的把位移值取反得到。
因此可得：
Mview=Mcamera−1=(TcameraMUVN)−1=MUVN−1Tcamera−1=[UxUyUz−dot(U,camPos)VxVyVz−dot(V,camPos)NxNyNz−dot(N,camPos)0001]M_{view} = M_{camera}^{-1} = (T_{camera}M_{UVN})^{-1} = M_{UVN}^{-1} T_{camera}^{-1} \\ = \left[\begin{matrix} U_x & U_y & U_z& -dot(U,camPos) \\ V_x & V_y & V_z& -dot(V,camPos) \\ N_x & N_y & N_z& -dot(N,camPos) \\ 0 & 0 & 0 & 1 \end{matrix}\right] Mview=Mcamera−1=(TcameraMUVN)−1=MUVN−1Tcamera−1=⎣⎢⎢⎡UxVxNx0UyVyNy0UzVzNz0−dot(U,camPos)−dot(V,camPos)−dot(N,camPos)1⎦⎥⎥⎤
一般矩阵的逆矩阵计算
上面计算视图矩阵的方法虽然不错，但是也只是写写简单demo时能用。在开发引擎或渲染框架时，camera可能挂在某个节点下面，也就是存在一系列的父子坐标系嵌套，这种情况下直接计算出上面UVN矩阵需要的U,V,NU,V,NU,V,N向量并不方便。以N向量为例：
N=worldEyePos−worldTargetN = worldEyePos - worldTarget N=worldEyePos−worldTarget
需要世界空间中的camera坐标和目标点的坐标。虽然我们可以先计算出camera的物体到世界矩阵，然后从中获取到世界坐标。但是还有一个问题是，我们往往需要通过旋转来控制camera，一般会使用到四元数。所以最后的结果往往是我们需要一个通用的逆矩阵计算方法。一般会使用标准伴随阵去计算。具体不讨论了，给出mini3d.js中的方法，有详细的注释：
/*** Calculate the inverse matrix of source matrix, and set to this.* @param {Matrix4} source The source matrix.* @returns this*/setInverseOf(source){let s = source.elements;let d = this.elements;let inv = new Float32Array(16);//使用标准伴随阵法计算逆矩阵：//标准伴随阵 = 方阵的代数余子式组成的矩阵的转置矩阵//逆矩阵 = 标准伴随阵/方阵的行列式//计算代数余子式并转置后放入inv矩阵中inv[0] = s[5]*s[10]*s[15] - s[5] *s[11]*s[14] - s[9] *s[6]*s[15]+ s[9]*s[7] *s[14] + s[13]*s[6] *s[11] - s[13]*s[7]*s[10];inv[4] = - s[4]*s[10]*s[15] + s[4] *s[11]*s[14] + s[8] *s[6]*s[15]- s[8]*s[7] *s[14] - s[12]*s[6] *s[11] + s[12]*s[7]*s[10];inv[8] = s[4]*s[9] *s[15] - s[4] *s[11]*s[13] - s[8] *s[5]*s[15]+ s[8]*s[7] *s[13] + s[12]*s[5] *s[11] - s[12]*s[7]*s[9];inv[12] = - s[4]*s[9] *s[14] + s[4] *s[10]*s[13] + s[8] *s[5]*s[14]- s[8]*s[6] *s[13] - s[12]*s[5] *s[10] + s[12]*s[6]*s[9];inv[1] = - s[1]*s[10]*s[15] + s[1] *s[11]*s[14] + s[9] *s[2]*s[15]- s[9]*s[3] *s[14] - s[13]*s[2] *s[11] + s[13]*s[3]*s[10];inv[5] = s[0]*s[10]*s[15] - s[0] *s[11]*s[14] - s[8] *s[2]*s[15]+ s[8]*s[3] *s[14] + s[12]*s[2] *s[11] - s[12]*s[3]*s[10];inv[9] = - s[0]*s[9] *s[15] + s[0] *s[11]*s[13] + s[8] *s[1]*s[15]- s[8]*s[3] *s[13] - s[12]*s[1] *s[11] + s[12]*s[3]*s[9];inv[13] = s[0]*s[9] *s[14] - s[0] *s[10]*s[13] - s[8] *s[1]*s[14]+ s[8]*s[2] *s[13] + s[12]*s[1] *s[10] - s[12]*s[2]*s[9];inv[2] = s[1]*s[6]*s[15] - s[1] *s[7]*s[14] - s[5] *s[2]*s[15]+ s[5]*s[3]*s[14] + s[13]*s[2]*s[7] - s[13]*s[3]*s[6];inv[6] = - s[0]*s[6]*s[15] + s[0] *s[7]*s[14] + s[4] *s[2]*s[15]- s[4]*s[3]*s[14] - s[12]*s[2]*s[7] + s[12]*s[3]*s[6];inv[10] = s[0]*s[5]*s[15] - s[0] *s[7]*s[13] - s[4] *s[1]*s[15]+ s[4]*s[3]*s[13] + s[12]*s[1]*s[7] - s[12]*s[3]*s[5];inv[14] = - s[0]*s[5]*s[14] + s[0] *s[6]*s[13] + s[4] *s[1]*s[14]- s[4]*s[2]*s[13] - s[12]*s[1]*s[6] + s[12]*s[2]*s[5];inv[3] = - s[1]*s[6]*s[11] + s[1]*s[7]*s[10] + s[5]*s[2]*s[11]- s[5]*s[3]*s[10] - s[9]*s[2]*s[7] + s[9]*s[3]*s[6];inv[7] = s[0]*s[6]*s[11] - s[0]*s[7]*s[10] - s[4]*s[2]*s[11]+ s[4]*s[3]*s[10] + s[8]*s[2]*s[7] - s[8]*s[3]*s[6];inv[11] = - s[0]*s[5]*s[11] + s[0]*s[7]*s[9] + s[4]*s[1]*s[11]- s[4]*s[3]*s[9] - s[8]*s[1]*s[7] + s[8]*s[3]*s[5];inv[15] = s[0]*s[5]*s[10] - s[0]*s[6]*s[9] - s[4]*s[1]*s[10]+ s[4]*s[2]*s[9] + s[8]*s[1]*s[6] - s[8]*s[2]*s[5];//计算行列式，选择方阵的第一列，对该列中的四个元素S[0],s[1],s[2],s[3]分别乘以对应的代数余子式，然后相加let det = s[0]*inv[0] + s[1]*inv[4] + s[2]*inv[8] + s[3]*inv[12];//注：选择任意一行，例如第一行，也是可以的//let det = s[0]*inv[0] + s[4]*inv[1] + s[8]*inv[2] + s[12]*inv[3];if(det===0){return this;}det = 1 / det;for(let i=0; i<16; ++i){d[i] = inv[i] * det;}return this;}
矩阵变换法线
通用的正确变换方法：使用逆转置矩阵
很多书上都讲了，使用一个非正交矩阵去变换法线时，变换后的法向量最终不垂直于变换后的表面。变换矩阵需要使用原空间变换矩阵的逆转置矩阵。
推导
因为切线方向和法线方向是垂直的，所以同一顶点的切向量TTT和法向量NNN必须满足等式N⋅T=0N \cdot T = 0N⋅T=0。而变换后的切线T′T'T′和法线N′N'N′仍然满足T′⋅N′=0T' \cdot N'=0T′⋅N′=0。设变换矩阵为MMM，则有T′=MTT'=MTT′=MT。设变换法向量NNN使用的矩阵为GGG，那么有：
N′⋅T′=(GN)⋅(MT)=0N' \cdot T' = (GN) \cdot (MT) = 0 N′⋅T′=(GN)⋅(MT)=0
将(GN)⋅(MT)(GN)\cdot(MT)(GN)⋅(MT)这个点积写成矩阵乘法，则有：
(GN)⋅(MT)=(GN)T(MT)=NTGTMT=0(GN)\cdot(MT) = (GN)^T(MT) = N^TG^TMT=0 (GN)⋅(MT)=(GN)T(MT)=NTGTMT=0
由于NTT=0N^TT=0NTT=0，则当GTM=IG^TM=IGTM=I时上式成立，因此可以得到
G=(M−1)TG = (M^{-1})^T G=(M−1)T
即，使用变换点的矩阵的逆矩阵的转置矩阵，也就是逆转置矩阵，可正确变换法线。
正交矩阵可直接使用原变换矩阵
因为正交矩阵的逆矩阵就等于其转置矩阵，所以其逆转置矩阵就是其自身。因此如果确定矩阵是正交矩阵，就可以直接变换法线。例如我们在世界空间中计算法线贴图时，计算得到的切线空间到世界空间变换矩阵就是一个正交矩阵，因此可以直接用其变换法线贴图中解压出来的切线空间的法线，将其变换到世界空间，在世界空间中进行光照计算。
3D场景中物体的法线变换
在shader中我们经常需要将物体的法线变换到世界空间，这样就需要一个物体空间到世界空间的法线变换矩阵，这个矩阵是物体空间到世界空间变换矩阵的逆转置矩阵，也就是世界空间到物体空间变换矩阵的转置矩阵。由于从物体空间到世界空间的变换比较复杂，可能是很多层坐标系变换的叠加，且可能包含非等比缩放，所以一般游戏引擎都会计算出物体空间到世界空间变换矩阵的逆矩阵（使用上面的标准伴随阵方法），得到世界空间到物体空间的变换矩阵。理论上，将这个矩阵转置一下，就可以作为法线变换矩阵传入shader中。不过可以通过反转矩阵向量乘法的顺序，省掉转置矩阵的计算，也这样可以省掉法线矩阵的计算和uniform的传入，直接使用世界空间到物体空间的变换矩阵
v_worldPos = u_object2World*a_Position;v_worldNormal = normalize(a_Normal * mat3(u_world2Object));
在上面的glsl代码中，我们看到计算世界空间的位置v_worldPos时，使用的是向量右乘矩阵u_object2World的方法，这说明我们在使用右乘约定。但是下面计算世界空间法线的时候，我们使用了向量左乘u_world2Object矩阵，这并不是我们破坏了约定，而是因为这样等价于向量右乘该矩阵的转置矩阵。这样我们就可以直接使用世界空间到物体空间变换矩阵u_world2Object来变换法线。这个矩阵还有其他作用，因此这么使用一专多能，还省去了转置计算法线矩阵并传入。算是一个小技巧。
