关于“红宝书”中透视投影矩阵的错误，以及该矩阵的推导

最近在看“红宝书”时发现了其中第165页中关于透视投影矩阵的错误。书中原式如下图所示：

其中 $\mathbf{M}[0,2]$ 和 $\mathbf{M}[1,2]$ 错误出现了除以 2 的操作（应该是宽度），且 $\mathbf{M}[2,3]$ 则漏填了个负号。
接下来我们自己重新推导这个矩阵。

• 首先需要明确的是，在观察空间中摄像机的取景范围将构成一个视锥体，而透视投影矩阵会将视锥体变换成一个立方体（NDC）：可以想象是让视锥体的近平面保持不变，而远平面则缩放到与近平面同等大小，然后映射到 [-1,1] (或者 [0,1] )中。在远平面缩小的这个过程中，在视锥体中的物体也将跟着缩小，从而实现“近大远小”的效果。

• 其次，坐标系的偏手性非常重要，它关乎到 $n,f$ 的符号和大小判定。本文采用的是：观察空间为右手系（换变前），NDC为左手系（换变后）。

• 最后，采用行向量还是列向量，以及映射范围的不同也将导致不一样的结果。本文采用列向量的形式，映射范围为 [-1,1].

  $1.$ 投影点的坐标表示

  我们先观察下图：

由相似三角形可以得出，对于一个点 $\mathbf{P}$ ，它的投影点 ${\mathbf{P}}^{\prime }$ 的 $x$ 和 $y$ 坐标的表达式如下：

$$x=-\frac{d}{z}x$$

$$\begin{array}{}\text{(1.1)}& y=-\frac{d}{z}y\end{array}$$

由于相机望向前方的位置是 $-z$，因此表达式里出现了负号。

$2.$ 视锥体的转换

由于观察空间是右手系，因此看向前方的 $z$ 轴是负值；相机的位置是原点，因此表示近平面的距离记为 $-n$，表示远平面的距离记为 $-f$，它们之间满足 $-f⩽z⩽-n$.在近平面中，分为左右上下，分别记为 $l,r,t,b$，它们分别是围成近平面的边的中点，并且满足 $l⩽x⩽r,b⩽y⩽t$.如下图所示：

它们需要映射到 $[-1,1]$ 中，考虑使用斜截式表达这个映射，如下图所示：

在斜截式方程表达式 $y=kx+b$ 中，$$k=\frac{1-(-1)}{r-l}=\frac{2}{r-l}$$ $$\begin{array}{}\text{(2.1)}& b=(-1)-\frac{2l}{r-l}=-\frac{2l}{r-l}-1\end{array}$$.
因此对于点 $\mathbf{P}$ ，把映射后的【它在近平面的投影点的 $x$ 分量】记为 $x^{\prime }$, 将它应用上面的 $y=kx+b$ 的形式，则有

$$\begin{array}{rl}{x}^{\prime }& =\frac{2}{r-l}x-\frac{2l}{r-l}-1\\ & =(x-l)\frac{2}{r-l}-1\end{array}$$

同理，对于 $y^{\prime }$ 也能得出如下的表达式：

$$y^{\prime }=(y-b)\frac{2}{t-b}-1$$

如下图所示：

以上两个式子就是将 $l,r,t,b$ 映射到 $[-1,1]$ 的线性函数。而由式子（1.1）可以得出

$$x=-\frac{n}{z}x$$

$$y=-\frac{n}{z}y$$

将以上两个变量代入上面 $x^{\prime },y^{\prime }$ 的线性函数，可得

$$\begin{array}{}\text{(2.2)}& x^{\prime }=\frac{2n}{r-l}(-\frac{x}{z})-\frac{r+l}{r-l}\end{array}$$

$$\begin{array}{}\text{(2.3)}& y^{\prime }=\frac{2n}{t-b}(-\frac{y}{z})-\frac{t+b}{t-b}\end{array}$$

到这一步，视锥体已由锥体变为长方体。接着需要将 $z$ 分量同样映射到 $[-1,1]$ 区间，就可以把长方体变为最终齐次除法中的立方体（NDC）。$z$ 分量的计算麻烦一点，因为深度信息的处理涉及到深度插值的处理过程。在大多数的教程中关于 $z$ 分量的推导都是根据与 $x,y$ 的关系一点一点推导出来的，这里采用另外一种思考方法。纵深方向的 $-n,-f$ 分别的映射关系如下：

$$-n\Rightarrow -1$$

$$-f\Rightarrow 1$$

从这里就可以看出，这个映射颠倒了变换前（观察空间）$z$ 轴的正指向，由右手系变成了左手系。
接着，考虑一下$z$ 分量的插值需求。通常对于 $x$ 和 $y$ 分量都是比较直接的，但是对于 $z$ 分量有着特殊的需求。

• 首先，在应用透视投影矩阵后，每个顶点坐标都会被其齐次坐标 $w$ 分量除（齐次除法），对于 $z$ 分量来说，这意味着原始的 $z$ 值会变成 $\frac{z}{w}$ 的形式($w$ 作为分母)。也就是说，$w$ 与 $z$ 成比例。
• 其次，人类视觉系统对深度的感知是非线性的，我们希望远处的物体变化较小（远离摄像机的地方有较低的精度），而近处的物体变化较大（靠近摄像机的地方有更高的精度）。而我们希望在处理 $z$ 分量时也能使用线性插值以获得正确的深度值。

正因为如此，在光栅化阶段顶点的 $z$ 值将会以 $\frac{1}{z}$ 的形式进行插值，随着 $z$ 的增大，$\frac{1}{z}$ 的变化速率减小，这正好符合我们对深度精度中模拟人眼感受的需求。因此，应用斜截式方程可以得出以下的表达式：

$$\begin{array}{}\text{(2.4)}& z^{\prime }=\frac{A}{z}+B\end{array}$$

其中 $A,B$ 是两个待定系数，它们可以联立以下方程组：

$$\{\begin{array}{l}\frac{A}{-n}+B=-1\\ \frac{A}{-f}+B=1\end{array}$$

解方程组，可得

$$A=\frac{2nf}{f-n}$$

$$B=\frac{f+n}{f-n}$$

将其代入以上的式子（2.4）可得关于 $z^{\prime }$ 的线性映射函数：

$$\begin{array}{}\text{(2.5)}& z^{\prime }=-\frac{2nf}{f-n}(-\frac{1}{z})+\frac{f+n}{f-n}\end{array}$$

这里把 $A$ 和 $\frac{1}{z}$ 取负号是为了和式子（2.2）式子（2.3）统一起来。到这一步，长方体就变换成了NDC。观察式子（2.2），式子（2.3）和式子（2.5），它们均含有同样的分母 $-z$. 所以，点 $\mathbf{P}$ 的齐次坐标点为：

$$\mathbf{P}(-x^{\prime }z,-y^{\prime }z,-z^{\prime }z,-z)$$

由式子（2.2），式子（2.3）和式子（2.5），可以得到关于齐次点的表达式如下：

$$\begin{array}{}\text{(2.6)}& \begin{array}{rl}-x^{\prime }z& =\frac{2n}{r-l}x+\frac{r+l}{r-l}z\\ -y^{\prime }z& =\frac{2n}{t-b}y+\frac{t+b}{t-b}z\\ -z^{\prime }z& =-\frac{f+n}{f-n}z-\frac{2nf}{f-n}\end{array}\end{array}$$

$3.$ 提取矩阵因子

现在，可以提取式子（2.6）中关于 $x,y,z$ 的因子了，根据式子的逻辑关系可逐列填出变换矩阵如下：

$${\mathbf{P}}^{\prime }={\mathbf{M}}_{\mathbf{frustum}}\mathbf{P}=\left[\begin{array}{cccc}\frac{2n}{r-l}& 0& \frac{r+l}{r-l}& 0\\ 0& \frac{2n}{t-b}& \frac{t+b}{t-b}& 0\\ 0& 0& -\frac{f+n}{f-n}& -\frac{2nf}{f-n}\\ 0& 0& -1& 0\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\\ 1\end{array}\right]$$

则透视投影矩阵已求。在“红宝书”中，变换前（观察空间）为左手系，$z$ 值随着观察方向是正向累加的，所以 $n,f$ 均为正。代入计算之后 $\mathbf{M}[2,3]$ 的位置仍然为负，因此原文此处有误。而 $\mathbf{M}[0,2]$ 和 $\mathbf{M}[1,2]$ 错误出现了除以 2 的操作则是完全莫须有的错误。

$4.$ GAMES101中关于透视投影矩阵的推导

GAMES101中闫令琪老师引用“虎书”的推导方式：https://www.bilibili.com/video/BV1X7411F744/?spm_id_from=333.788.videopod.episodes&vd_source=57c4e007542d6a99df1a601686225ef6&p=4

其给出的推导方式更直观：首先假设已经存在一个透视投影矩阵 ${\mathbf{M}}_{\mathbf{persp}}$ ，将它设为未知。接着构造一个透视变换到正交的矩阵 ${\mathbf{M}}_{\mathbf{persp}\Rightarrow \mathbf{ortho}}$ ，这一步是将锥体变换为长方体。最后，将这个变换左乘一个正交矩阵 ${\mathbf{M}}_{\mathbf{ortho}}$ ,这一步的目的是将长方体变成正方体，也就是 NDC.表达式为：

$${\mathbf{M}}_{\mathbf{persp}}={\mathbf{M}}_{\mathbf{ortho}}{\mathbf{M}}_{\mathbf{persp}\Rightarrow \mathbf{ortho}}$$

$5.$ 代码实现

代码内容比较简单，因此不注释了。

点击查看代码

public class Vector4
{
    public float X { get; private set; }
    public float Y { get; private set; }
    public float Z { get; private set; }
    public float W { get; private set; }

    public Vector4(float x, float y, float z, float w)
    {
        X = x;
        Y = y;
        Z = z;
        W = w;
    }

    // 矩阵与向量相乘的方法
    public static Vector4 Multiply(float[,] matrix, Vector4 vector)
    {
        float x = matrix[0, 0] * vector.X + matrix[0, 1] * vector.Y + matrix[0, 2] * vector.Z + matrix[0, 3] * vector.W;
        float y = matrix[1, 0] * vector.X + matrix[1, 1] * vector.Y + matrix[1, 2] * vector.Z + matrix[1, 3] * vector.W;
        float z = matrix[2, 0] * vector.X + matrix[2, 1] * vector.Y + matrix[2, 2] * vector.Z + matrix[2, 3] * vector.W;
        float w = matrix[3, 0] * vector.X + matrix[3, 1] * vector.Y + matrix[3, 2] * vector.Z + matrix[3, 3] * vector.W;
        return new Vector4(x, y, z, w);
    }

    public override string ToString()
    {
        return $"({X:F2}, {Y:F2}, {Z:F2}, {W:F2})";
    }
}

public class PerspMatrix
{
    public float L { get; private set; }
    public float R { get; private set; }
    public float T { get; private set; }
    public float B { get; private set; }
    public float N { get; private set; }
    public float F { get; private set; }

    public PerspMatrix(float l, float r, float t, float b, float n, float f)
    {
        if (r == l || t == b || f == n)
            throw new ArgumentException("Invalid frustum parameters: division by zero.");

        L = l;
        R = r;
        T = t;
        B = b;
        N = n;
        F = f;
    }

    public float[,] Frustum()
    {
        float invWidth = 1.0f / (R - L);
        float invHeight = 1.0f / (T - B);
        float invDepth = 1.0f / (F - N);

        float[,] matrix4x4 = new float[4, 4];
        matrix4x4[0, 0] = (float)(2.0 * N * invWidth);
        matrix4x4[0, 1] = 0;
        matrix4x4[0, 2] = (float)((R + L) * invWidth);
        matrix4x4[0, 3] = 0;

        matrix4x4[1, 0] = 0;
        matrix4x4[1, 1] = (float)(2.0 * N * invHeight);
        matrix4x4[1, 2] = (float)((T + B) * invHeight);
        matrix4x4[1, 3] = 0;

        matrix4x4[2, 0] = 0;
        matrix4x4[2, 1] = 0;
        matrix4x4[2, 2] = (float)(-(F + N) * invDepth);
        matrix4x4[2, 3] = (float)(-2.0 * F * N * invDepth);

        matrix4x4[3, 0] = 0;
        matrix4x4[3, 1] = 0;
        matrix4x4[3, 2] = -1;
        matrix4x4[3, 3] = 0;

        return matrix4x4;
    }

    public void PrintMatrix(float[,] matrix)
    {
        Console.WriteLine("Projection Matrix:");
        for (int i = 0; i < 4; i++)
        {
            for (int j = 0; j < 4; j++)
            {
                Console.Write($"{matrix[i, j]:F2}\t");
            }
            Console.WriteLine();
        }
    }

    public Vector4 Apply(Vector4 vector)
    {
        float[,] matrix = Frustum();
        return Vector4.Multiply(matrix, vector);
    }
}