最近在看“红宝书”时发现了其中第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 \]

\[y=-\frac{d}{z}y \tag{1.1} \]

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

$2.$ 视锥体的转换

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

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

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

\[\begin{align*} x^{\prime}&=\frac{2}{r-l}x-\frac{2l}{r-l}-1 \\ &=(x-l)\frac{2}{r-l}-1 \end{align*} \]

同理，对于 $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}$ 的线性函数，可得

\[x^{\prime}=\frac{2n}{r-l}(-\frac{x}{z})-\frac{r+l}{r-l} \tag{2.2} \]

\[y^{\prime}=\frac{2n}{t-b}(-\frac{y}{z})-\frac{t+b}{t-b} \tag{2.3} \]

到这一步，视锥体已由锥体变为长方体。接着需要将 $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}$ 的变化速率减小，这正好符合我们对深度精度中模拟人眼感受的需求。因此，应用斜截式方程可以得出以下的表达式：

\[z^{\prime} = \frac{A}{z} + B \tag{2.4} \]

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

\[ \begin{cases} \frac{A}{-n} + B = -1 \\ \frac{A}{-f} + B = 1 \end{cases} \]

解方程组，可得

\[A=\frac{2nf}{f-n} \]

\[B=\frac{f+n}{f-n} \]

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

\[z^{\prime}=-\frac{2nf}{f-n}(-\frac{1}{z})+\frac{f+n}{f-n} \tag{2.5} \]

这里把 $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{align*} -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{align*} \tag{2.6}\]

$3.$ 提取矩阵因子

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

\[\mathbf{P^{\prime}}=\mathbf{M_{frustum}}\mathbf{P}=\begin{bmatrix} \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{bmatrix} \begin{bmatrix} x\\y\\z\\1 \end{bmatrix} \]

则透视投影矩阵已求。在“红宝书”中，变换前（观察空间）为左手系，$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_{persp}}$ ，将它设为未知。接着构造一个透视变换到正交的矩阵 $\mathbf{M_{persp \Rightarrow ortho}}$ ，这一步是将锥体变换为长方体。最后，将这个变换左乘一个正交矩阵 $\mathbf{M_{ortho}}$ ,这一步的目的是将长方体变成正方体，也就是 NDC.表达式为：

\[\mathbf{M_{persp}}=\mathbf{M_{ortho}} \mathbf{M_{persp \Rightarrow 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);
    }
}

posted on 2025-01-22 00:28 小花猫喵喵阅读(64) 评论(0) 收藏举报

刷新页面返回顶部

\(1.\) 投影点的坐标表示

\(2.\) 视锥体的转换

\(3.\) 提取矩阵因子

\(4.\) GAMES101中关于透视投影矩阵的推导

\(5.\) 代码实现

公告