理解投影矩阵

理解投影矩阵

本文结合基础的显卡渲染流程，解释了如何推导透视矩阵，本文不涉及正交矩阵的计算。

引言

如下图1所示，在一个3D物体的渲染过程中，需要进行多次的坐标空间的转换。其中前3次的坐标空间转换用到的变换矩阵分辨是World Matrix（又名Model Matrix）, View Matrix, Projection Matrix。这三次坐标空间的转换通常发生在vertex shader中，以unity shader为例，在vertex shader我们通常通过调用unity提供的UnityObjectToClipPos函数完成这三次转换。而后面的两次转换，即Clip Space-> NDC Space-> Screen Space是在从vertex shader到fragment shader的裁剪插值过程中由硬件自动完成的。更准确地说，是在先进行了裁剪插值过程，然后进行了两次坐标空间转换的过程。 进行裁剪插值和坐标变换的过程是一个不可配置的过程，也就是说硬件会强制进行上述两次的坐标空间的转换。以unity shader为例，编程人员需要使用SV_POSITION标记vertex shader的一个float4类型的返回值来告诉硬件这个变量就是需要进行两次坐标空间转换的变量，请在裁剪插值时对其进行处理。所以在fragment shader中读到的被SV_POSITION标记的变量已经是在Screen Space中了。

理解Model Matrix和View Matrix相对于理解Projection Matrix的更简单一些，这主要是因为Model Space-> World Space-> View Space的转换都是发生在3D坐标空间之间，而Project Matrix本意是要完成从3D坐标空间到2D坐标空间的映射，即从View Space-> Clip Space。但如龙书《3d_game_programming_with_DirectX11.pdf》5.6.3 所描述，矩阵这样一个转换形式不能够完整地表达出一个3D坐标空间到2D坐标空间的变换。所以整个投影变换被拆成了两部分，一部分由Projection Matrix完成View Space到Clip Space的转换（Clip Space还是一个3D坐标空间），另一个部分由硬件完成Clip Space到NDC Space的转换（NDC Space是一个2D坐标空间了）。所以计算Projection Matrix需要考虑到和其后续硬件工作的约定，这一定程度上增加了对Projection Matrix的理解难度。为了更好地理解Projection Matrix，下文将从数学原理入手依次说明投影变换到底想要做什么工作、为什么需要在坐标变换过程中引入硬件及硬件到底做了什么工作、和投影变换相关的一些应用问题。

Directx下投影变换的数学原理

通常情况下，使用View Space中的一个截头椎体(Frustum)来描述用做投影变换的“相机”。DirectX和OpenGL在相机的朝向上略有不同：DirectX中所有3D坐标空间中都是用左手坐标系，DirectX中相机位于View Space坐标原点，看向z轴正方向；OpenGL中View Space使用的是右手坐标系，其余都是左手坐标系，相机位于坐标原点，看向z轴的负方向。这一小节中接下来提到的所有坐标变换都是在DirectX环境进行的。常用描述Frustum的方式有两种，我们这里只讨论如下图2所示的其中一种，即：近平面n，远平面f，沿着x轴负方向看过去的垂直夹角α，宽高比r。需要注意的是，①在View Space中，近平面和远平面是和xy-平面平行的，所以n、f是近平面和远平面到坐标原点的距离；②宽高比r是投影平面宽w和高度h的比值，也可以说是Frustum的宽高比，因为这两个比值是一样的，为了后续描述方便这里通一称为投影平面宽高比；③投影平面到坐标原点的距离为d，投影平面并不是Frustum的一部分，但在后续的计算中有着重要的作用。这样，通过四元组(n,f,α,r)就能够唯一确定一个Frustum，也就能够进行View Space到投屏平面的映射的了。这里之所以没有说是到NDC Space的映射是因为在此处探讨数学原理的时候，我们只关心一个3D物体是怎么样被映射到一个2D平面上的。 至于NDC Space、Screen Space等探讨完数学原理之后其概念会顺势变得清晰可见了。

设A(x,y,z)是存在于Frustum中的一点，A’(x’,y’,z)是点A在投影平面上的一个投影点。实际应用中，所有上述坐标空间都是用的是齐次坐标系(Homegeneous Coordinates),其还有第四个分量w，用来表示一个坐标到底是向量还是点。如果w=0，那么表示该坐标是一个向量，如果w=1，表示该坐标是一个点，也就是说A,A’完整坐标是(x,y,z,1),(x’,y’,z’,1)。更多关于齐次坐标系的介绍请参照龙书3.2.1。这一时间，暂时只关注坐标分量x，y和x’，y’的映射关系。通过下图3这样一个Frustum的顶视图以及三角形相似性很容易得到公式1
$$\frac{x^{\prime }}{d}=\frac{x}{z}\Rightarrow x^{\prime }=\frac{xd}{z}公式1$$

并且根据简单的三角形几何相关知识，有如下公式成立
$$\frac{h}{2d}=\tan \left(\frac{\alpha }{2}\right)\Rightarrow d=\frac{h}{\text{2}\tan \left(\frac{\alpha }{2}\right)}公式2$$
此外按照设定有如下r和w、h的关系成立
$$r=\frac{w}{h}\Rightarrow w=rh公式3$$
结合公式1和公式3可以得到点A’的x轴坐标如下公式4所示。并且通过图3可以知道，最终x’的取值范围为是$\left[-\frac{w}{2},\frac{w}{2}\right]$
$$x^{\prime }=\frac{xd}{z}=\frac{xh}{2z\tan \left(\frac{\alpha }{2}\right)}\in \left[-\frac{w}{2},\frac{w}{2}\right]公式4$$
同理通过下图4可以知道点A’在y轴的坐标如下公式5所示，并且可以知道y’的取值范围是$\left[-\frac{h}{2},\frac{h}{2}\right]$
$$y^{\prime }=\frac{yd}{z}=\frac{yh}{2z\tan \left(\frac{\alpha }{2}\right)}\in \left[-\frac{h}{2},\frac{h}{2}\right]公式5$$

如果单从数学的角度来说，似乎已经推导出来了A和A’坐标分量x、y的关系，但是公式4和公式5中x’和y’的取值范围分别是$\left[-\frac{w}{2},\frac{w}{2}\right]$、$\left[-\frac{h}{2},\frac{h}{2}\right]$。这也就是说在进行后续的向显示屏进行映射的时候还需要知道投影平面的宽w和高h。因为只有知道了投影平面的宽w和高h，才能通过比例关系知道点A’最终在显示屏上面的坐标。由于后续的映射工作都是在硬件上面完成，传递投影平面宽w和高h增加了硬件的设计成本。所以目前的做法是将A’做一次坐标映射（线性变换），使x’和y’最终的取值范围都是[-1,1]，这样的话，就无需告诉硬件投影平面的信息也能够计算出最终A在显示屏上面的坐标了。通过公式4和公式5可以知道，只需要将公式4里面的x’除以$\frac{w}{2}$,公式5里面的y’除以$\frac{h}{2}$就能够完成上述坐标映射。那么结合公式3可以得到最终A’的坐标和A的坐标关系如下公式6、7所示。需要注意的是，由于又进行了一次坐标变换此时A’就不是投影平面上面的一点了，我们约定此时A’所在的坐标空间叫做NDC Space。
$$\begin{gathered} x^{\prime }=\frac{\frac{xh}{2z\tan \left(\frac{\alpha }{2}\right)}}{\frac{w}{2}}=\frac{xh}{zw\tan \left(\frac{\alpha }{2}\right)}=\frac{x}{rz\tan \left(\frac{\alpha }{2}\right)}\in \left[-\text{1,}1\right]\text{公式6} \\ {y}^{\prime }=\frac{\frac{yh}{2z\tan \left(\frac{\alpha }{2}\right)}}{\frac{h}{2}}=\frac{y}{z\tan \left(\frac{\alpha }{2}\right)}\in \left[-\text{1,}1\right]\text{公式7} \end{gathered}$$
仅仅知道了A点x，y坐标是如何映射还远远不够，因为Frustum中可能有多个点同时映射到投影平面上面的同一个点。试想一条View Space中的从原点发出的穿过Frustum的射线，射线上面的所有点都会映射到投影平面上面的同一位置，这就需要知道映射到投影平面同一位置的所有Frustum点中，哪一个是在“最前面”的。最容易想到是存储View Space中的点的z值，这样的直接存储z值或者z经过线性变换的buffer叫做W-Buffer（大概因为z值一般存储在float4的w字段，所以叫w-buffer）。一些比较老的硬件支持W-Buffer，比较新的硬件很多都不再支持W-Buffer转而使用存储$\frac{1}{z}$或者$\frac{1}{z}$线性变换值的方式，即Z-Buffer。这主要是因为W-Buffer在View Space中是线性的，但是在Screen Space中不是线性的，而Z-Buffer在View Space中不是线性的，而在Screen Space中是线性的，即下表
|Buffer类型|View Space|Screen Space|
|-----|
|W-Buffer|线性的|非线性的|
|Z-Buffer|非线性的|线性的|

为什么$\frac{1}{z}$在Screen Space中是线性的？线性的意义又是什么？

为了弄懂这个问题，设有如下图5所示的在View Space中的一个平行于z轴的线段PQ，线段PQ到Z的距离为$y_0$，$A(y_0,z)$是PQ上面一点，$A^{\prime }(y^{\prime },d)$是点A’在投影平面上的投影点。根据三角形相似性很容易知道
$$\frac{d}{z}=\frac{y^{\prime }}{y_0}\Rightarrow y^{\prime }=\frac{1}{z}d{y}_0公式8$$

其中$d,y_0$都是常量，所以根据公式8能够看出$y^{\prime }$和$\frac{1}{z}$线性相关，而投影平面到Screen Space只是进行了缩放这样的线性变换，所以最终$\frac{1}{z}$在Screen Space中是线性的。那么在Screen Space中是线性的有什么意义呢？主要的意义还是能够简化硬件设计，节省成本。由于光栅化是在Screen Space中进行的，在进行光栅化的时候，在$\frac{1}{z}$是线性的这样的条件下，硬件只需简单的进行一次线性插值就能够片元三角形顶点的$\frac{1}{z}$值求得三角形内部一点的$\frac{1}{z}$的值了。此外，在Shadow Map和Post Processing中一个具有Screen Space线性的$\frac{1}{z}$也具有重大意义（方便运算？待考证）。
不管是Z-Buffer还是W-Buffer，用来存储深度信息的Buffer都被称作深度缓冲区(Depth-Buffer)，深度缓冲区中存储的值叫做深度值。不同的硬件API设定的深度缓冲区的取值范围不同，譬如DirectX要求的范围是[0,1],即近平面n处的深度值是0，远平面f处的深度值是1；而OpenGL设定的深度缓冲区的取值范围是[-1,1]。在选用DirectX作为探讨环境的情况下，因为取值范围不同是不能够直接将$\frac{1}{z}\in \left[\frac{1}{f},\frac{1}{n}\right]$写入深度缓冲区的。为了适应DirectX关于深度缓冲区取值范围的约定，根据$\frac{1}{z}$在Screen Space中是线性的属性，只需要对$\frac{1}{z}$做一次线性变换即可。那么最终设定点A坐标z到深度值的映射关系为$f\left(\frac{1}{z}\right)=a+b\frac{1}{z}$，并且根据以上讨论已知,当z为n时$f\left(n\right)=0$,当z为f时$f\left(f\right)=1$。解二元一次方程组可以最终求得
$$\begin{gathered} a=\frac{f}{f-n}公式9 \\ b=\frac{-nf}{f-n}公式10 \\ f\left(z\right)=\frac{f}{f-n}+\frac{-nf}{f-n}\frac{1}{z}公式11 \end{gathered}$$

矩阵形式。及硬件做了什么

上一小节中我们推导出了View Space中的点A到NDC Space是的点A’的坐标映射关系。但是因为View Space到NDC Space之间的转换不是一个线性转换，所以找不到有效的矩阵表达形式（只有线性变换才有矩阵表达形式，参照龙书3.1 - 3.4）。通过观察公式6、7、11可以看到，如果将转换因子提出一个公约数z，剩下的部分是一个线性的表达式！这也是目前实际应用的做法，即将公式6、7、11拆解成两个部分：①线性变换的部分：可以转换一个矩阵表达形式，约定这个矩阵叫做投影矩阵(Projection Matrix，简称P)。；②非线性变换的部分：非线性变换的部分就是除以z，并且非线性的部分会有硬件自动帮我们做，这个除以z的流程叫做齐次除法(homogeneous divide)*。
点A(x,y,z,1)这样一个行向量右乘矩阵P之后会得到一个传递给硬件的行向量$A_p$，我们需要在$A_p$中记录下来A的z值，以便硬件进行齐次除法。这里用到的小技巧是让P的第3行第4列的元素（记为P.m34）等于1，这样$A_p$的w分量就等于A.z。至此可以写出矩阵P的完整性形式如下
$$\left[\begin{array}{cccc}\frac{1}{r\tan \left(\frac{\alpha }{2}\right)}& 0& 0& 0\\ 0& \frac{1}{\tan \left(\frac{\alpha }{2}\right)}& 0& 0\\ 0& 0& \frac{f}{f-n}& 1\\ 0& 0& \frac{-nf}{f-n}& 0\end{array}\right]$$
现在再回过头来看一下软硬件是如何配合完成View Space中的点A到Screen Space中的点$A_s$的变换的。如下图6所示，点A通过右乘投影矩阵P被转换到Clip Space中得到点$A_p$。P.m34的值为1，故而A的z值被复制到了$A_p$的第四个分量w上。此时$A_p$的四个坐标分量应该满足条件$-z⩽x_p⩽z,-z⩽y_p⩽z,0⩽z_p⩽z,0⩽z$。如果 $A_p$不满足这样的条件，就说明A不再Frustum中，那么$A_p$就会被裁剪掉（如果有必要会有新的恰好在裁剪边缘的点生成）。然后会对$A_p$做齐次除法，即对$A_p$的每一个分量都除以它的第四个分量，数值上来看就是A的z值，进而得到在NDC Space中的点A’(x’,y’,z’,1)。点A’的前三个分量的取值范围是$x^{\prime }\in \left[-\text{1,}1\right],y^{\prime }\in \left[-\text{1,}1\right],z^{\prime }\in \left[0,1\right]$。随后一个被称作View Port Transform的线性变换将NDC Space中的点A’转换成Screen Space中的点$A_s$，然后硬件会对片元进行光栅化（也就是线性插值）操作。这些转换中，只有View Space到Clip Space的转换是通过软件（也就是Shader代码）控制的，后续的Clip Space-> NDC Space-> Screen Space都是由硬件自动完成的，不可配置的过程。

了解了硬件工作原理就能够更好的理解投影矩阵P。譬如P.m34只能是1么？通过上文可以知道，P.m34可以是任何大于0的值。回想，DirectX在View Space中使用的是左手坐标系，所以Frustum中的所有点的z值都是大于0的，需要在图6的裁剪阶段满足$0⩽A_p.w$，就需要P.m34满足大于0的条件。此外还需要P的其他元素都乘上m34，这样在图6中的齐次除法阶段正好将m34约去得到正确的x’,y’,和z’。刚刚描述的m34（> 0）不为1的P如下所示
$$\left[\begin{array}{cccc}\frac{m_{34}}{r\tan \left(\frac{\alpha }{2}\right)}& 0& 0& 0\\ 0& \frac{m_{34}}{\tan \left(\frac{\alpha }{2}\right)}& 0& 0\\ 0& 0& \frac{f{m}_{34}}{f-n}& m_{34}\\ 0& 0& \frac{-nf{m}_{34}}{f-n}& 0\end{array}\right]$$

OpenGL投影矩阵推导

在OpenGL环境下，可以直接照搬上文中DirectX的推导的过程。不过需要注意的，OpenGL环境下View Space使用的左手坐标系，并且OpenGL环境下使用坐标列向量左乘投影矩阵。在OpenGL环境下投影矩阵如下，同上理$m_{34}$是一个小于0的事情，通常情况下取 -1。
$$\left[\begin{array}{cccc}\frac{m_{34}}{r\tan \left(\frac{\alpha }{2}\right)}& 0& 0& 0\\ 0& \frac{34}{\tan \left(\frac{\alpha }{2}\right)}& 0& 0\\ 0& 0& -m_{34}\frac{f+n}{f-n}& -2m_{34}\frac{nf}{f-n}\\ 0& 0& m_{34}& 0\end{array}\right]$$