计算相机射线

# 计算摄像机光线 来源：https://momentsingraphics.de/CameraRays.html 发布时间：2026-06-20 我们正处于光线追踪与光栅化并存的过渡期。典型的实时渲染管线仍然使用光栅化（通常结合延迟着色）处理主可见性，然后通过光线追踪处理阴影、反射或全局光照。不过，我们正在接近这样一个阶段：至少在某些平台上，可以认真考虑抛弃光栅化，转而使用光线追踪处理主可见性。这样一来，就需要计算摄像机光线——由光线原点、光线方向和光线长度描述——并且要与光栅化方式保持一致。在光栅化中，通常使用世界到裁剪空间变换矩阵（也称为视图-投影矩阵）来定义摄像机。本文将推导如何基于这样的矩阵计算摄像机光线。 目标是得到一种同时适用于透视投影、正交投影以及该矩阵可能表示的其他投影的方法。直观的方法存在数值抵消问题，我将介绍一种更可靠的替代方案。总体而言，这不是一个难题：对于任何给定的摄像机模型（例如已知视场角的透视投影），很容易想出针对性的临时解决方案。但我认为，拥有一个基于现成变换矩阵、无需针对每种摄像机模型额外调试的解决方案很有价值。如果你不关心推导过程，可以直接复制代码列表 4 中的着色器代码（此处以及本博文中所有代码均已发布到公共领域）。 ## 裁剪空间中的摄像机光线 光栅化管线与裁剪空间严重依赖于齐次坐标，因此我们先回顾一下这个概念。如果一个点的三维笛卡尔坐标为 \((x', y', z')^\mathsf{T}\)，我们可以通过附加第四坐标 1 来获得齐次坐标：\((x', y', z', 1)^\mathsf{T}\)。这仍然描述一个三维点，但现在我们获得了用任意非零因子 \(w \neq 0\) 缩放其坐标的自由度，即 \[(x, y, z, w)^\mathsf{T} = (wx', wy', wz', w)^\mathsf{T}.\] 无论 \(w\) 如何选择，这些坐标仍然描述同一个点。我们可以通过除以第四分量 \(w\) 来恢复非齐次坐标（即去齐次化）： \[\frac{1}{w}(x, y, z, w)^\mathsf{T} = \left(\frac{x}{w}, \frac{y}{w}, \frac{z}{w}, 1\right)^\mathsf{T} = (x', y', z', 1)^\mathsf{T}.\] 齐次坐标使许多公式更简洁。例如，下面会看到，我们也可以写出三维空间中平面的齐次坐标，然后通过点积来检查点是否在平面上。它们还允许我们用 \(4 \times 4\) 矩阵表达平移。此外，它们对于透视投影的光栅化也很有用。透视投影在某个时刻必须进行除法运算。使用齐次坐标，这个除法发生在光栅化之前，并且就是上面提到的去齐次化。 与这种光栅化设计相伴的是裁剪空间的概念。对于屏幕空间，我们使用的坐标系中，坐标 \(x'_c\) 和 \(y'_c\) 在摄像机视锥体范围内从 -1 到 1（下标 \(c\) 表示裁剪空间）。在齐次坐标中，这些边界转化为 \(-w_c \leq x_c \leq w_c\) 和 \(-w_c \leq y_c \leq w_c\)。此外，我们根据裁剪空间 z 坐标定义近裁剪平面和远裁剪平面。远裁剪平面位于 \(z'_c = 1\)，即 \(z_c \leq w_c\)。近裁剪平面的定义因 API 而异：对于 Direct3D，不等式为 \(0 \leq z_c\)。对于 OpenGL，默认近裁剪平面位于 \(-w_c \leq z_c\)，但通过扩展 GL_ARB_clip_control 可配置为 Direct3D 的 \(0 \leq z_c\) 行为。该扩展已进入 OpenGL 4.5 的核心功能。对于 Vulkan，行为同样可配置。为处理这些差异，我使用变量 \(z'_n\)，对于 Direct3D 约定取 0，对于旧的 OpenGL 默认值取 -1，即无论哪种方式，近裁剪平面为 \(z'_n w_c \leq z_c\)。 典型渲染器会使用许多其他坐标系，如摄像机空间和每个物体的物体空间，但我们这里只关心世界空间，因为我们要在世界空间中获取光线。在光线追踪的上下文中，世界空间的定义很容易：它就是顶层加速结构所使用的空间。对于光栅化，我们需要准备世界到裁剪空间矩阵 \(M_{w,c} \in \mathbb{R}^{4 \times 4}\)，将点从世界空间变换到裁剪空间。如果 \(p_w = (x_w, y_w, z_w, w_w)^\mathsf{T}\) 表示点的世界空间坐标，对应的裁剪空间坐标为 \[p_c := (x_c, y_c, z_c, w_c)^\mathsf{T} := M_{w,c} p_w.\] 我在这里使用 OpenGL 约定，这与常见的线性代数约定一致。在 Direct3D 约定中，矩阵和向量会转置，乘法顺序也会颠倒（不知道是谁认为这是个好主意，也不知道为什么）。更多内容请参见下文。 现在，如果我们有屏幕空间坐标 \(x'_c, y'_c\) 在 -1 到 1 范围内，我们可以很容易地确定该像素在近裁剪平面和远裁剪平面上的点，以其在裁剪空间中的齐次坐标表示。它们分别是 \[n_c := (x'_c, y'_c, z'_n, 1)^\mathsf{T} \quad \text{和} \quad f_c := (x'_c, y'_c, 1, 1)^\mathsf{T}.\] 我们的目标是将此转换为世界空间中的光线原点和光线方向。听起来很简单，对吧？ ## 世界空间光线原点 就光线原点而言，确实很简单。我们取近裁剪平面上的点，将其变换到世界空间，然后去齐次化。为此，我们需要裁剪到世界空间变换矩阵 \(M_{c,w} := M_{w,c}^{-1}\)。那么所需的点就是 \(n_w := M_{c,w} n_c\)。由于光线追踪 API 不处理齐次坐标，我们需要对其进行去齐次化，即除以 w 坐标。代码列表 1 提供了 GLSL 实现。 **代码列表 1：** 返回给定裁剪到世界空间变换的摄像机在近裁剪平面上的一个点。当 `clip_space_xy` 为 \((-1, -1)\) 时，该点位于视口左上角；为 \((1, -1)\) 时位于右上角；为 \((1, 1)\) 时位于右下角（至少在默认的 OpenGL、Vulkan 和 Direct3D 行为下如此）。 ```glsl vec3 get_camera_ray_origin( vec2 clip_space_xy, mat4 clip_to_world) { // 使用旧的 OpenGL 约定时，z_n = -1 float z_n = 0.0; vec4 n_c = vec4(clip_space_xy, z_n, 1.0); vec4 n_w = clip_to_world * n_c; return n_w.xyz / n_w.w; } ``` ## 光线方向的错误方法 要描述摄像机光线，除了光线原点，还需要光线方向。直观的策略如下：使用上面的代码片段，其中 `z_n = 1` 计算远平面上的点。然后取远、近平面点的差，并对该向量进行归一化。从数学角度看，这种方法没有问题。它给出正确结果，而且计算效率也很高。尽管如此，我强烈建议不要使用它。实际上，我甚至不会给出其代码列表，以免有人无意中复制。 相反，我们直接看看这种方法的结果，如图 1 所示。我喜欢好的图形故障，但这并不是我们想要的结果。一个缓慢、连续的摄像机运动变成了抖动的混乱。发生了什么？ **图 1：** 一个摄像机沿直线缓慢移动的视频。摄像机距离原点相当远，并注视着以原点为中心的场景（即 Lumberyard Bistro）。由于光线方向向量计算中的舍入误差，摄像机运动出现抖动。 观察裁剪到世界空间矩阵 \(M_{c,w}\) 的最后两列，大致如下： ``` -5995.15332031 5994.86083984 11787.78613281 -11786.94140625 -4032.25268555 4031.86547852 -50.01052094 50.01062393 ``` 注意这两列几乎完全相同，只是符号相反。当我们计算远平面上的点时，计算 \(M_{c,w} f_c\)，而 \(f_c\) 的最后两个元素是 1。因此，在计算此矩阵-向量乘积时，我们是在对裁剪到世界空间矩阵的最后两列求和。由于它们几乎相等但符号相反，这种加法会导致灾难性抵消：这些大浮点数具有较大的绝对舍入误差，它们微小的差值将具有同样大的绝对舍入误差。于是远平面上的点精度降低。这种舍入误差的具体表现取决于具体矩阵元素，而矩阵元素又取决于摄像机位置。因此图 1 的结果不光滑。使用双精度算术计算裁剪到世界空间矩阵（在传递给着色器之前强制转换为浮点数）会减小这些伪影的幅度，但不会消除它们。可能还有其他问题，但抵消似乎是主要罪魁祸首。 ## 光线方向的正确方法 那么让我们尝试不同的方法。我们不通过指定直线上的两个点 \(n_c, f_c\) 来描述摄像机光线，而是定义两个平面，使得直线是这两个平面的交线。我们使用齐次坐标来描述这两个平面： \[G_c := (1, 0, 0, -x'_c)^\mathsf{T} \quad \text{和} \quad H_c := (0, 1, 0, -y'_c).\] 根据定义，裁剪空间中的点 \(p_c \in \mathbb{R}^4\) 在平面 \(G_c\) 上当且仅当 \(G_c^\mathsf{T} p_c = 0\)。在这个等式中，我们将行向量（由于转置）与列向量相乘，得到一个点积。展开后得到 \[G_c^\mathsf{T} p_c = x_c - x'_c w_c = 0,\] 等价于 \[\frac{x_c}{w_c} = x'_c.\] 因此，该平面包含了所有非齐次裁剪空间 x 坐标为 \(x'_c\) 的点。特别地，它包含了我们摄像机光线上的所有点。类似地，我们可以推导出 \[H_c^\mathsf{T} p_c = 0 \quad \Leftrightarrow \quad \frac{y_c}{w_c} = y'_c.\] 所以摄像机光线上的所有点 \(p_c\) 都在这两个平面的交线上，即满足 \[G_c^\mathsf{T} p_c = H_c^\mathsf{T} p_c = 0.\] 下一步推导是将这些平面变换到世界空间。注意 \[G_c^\mathsf{T} p_c = 0 \quad \Leftrightarrow \quad G_c^\mathsf{T} M_{w,c} p_w = 0 \quad \Leftrightarrow \quad (M_{w,c}^\mathsf{T} G_c)^\mathsf{T} p_w = 0.\] 因此，\(G_w := M_{w,c}^\mathsf{T} G_c\) 保存了平面 \(G_c\) 的世界空间坐标。注意我们使用了世界到裁剪空间变换来从裁剪空间变换到世界空间。这是因为变换平面需要使用逆的转置。我们刚才看到的就是这个规则的推导。 我们用同样的方式变换另一个平面：\(H_w := M_{w,c}^\mathsf{T} H_c\)。现在，我们需要两个平面交线的方向向量 \(d \in \mathbb{R}^3\)。为此，注意到平面 \(G_w, H_w\) 的前三个齐次坐标就是它们未归一化的世界空间法向量。我们将这些法向量记为 \(u, v \in \mathbb{R}^3\)。那么方向向量必须与这两个法向量都正交。因此，我们只需取叉积：\(d := u \times v\)。这就给出了一个实用的算法来计算光线方向，避免了前面方法的数值问题。代码列表 2 用 GLSL 实现了该算法。 **代码列表 2：** 返回给定世界到裁剪空间变换的摄像机在指定像素处的归一化光线方向。该函数中的部分计算可以移到常量的预计算中，见下文。`clip_space_xy` 的含义与代码列表 1 相同。 ```glsl vec3 get_camera_ray_direction_cross( vec2 clip_space_xy, mat4 world_to_clip) { vec3 mx = transpose(world_to_clip)[0].xyz; vec3 my = transpose(world_to_clip)[1].xyz; vec3 mw = transpose(world_to_clip)[3].xyz; float x = clip_space_xy.x; float y = clip_space_xy.y; vec3 u = fma(vec3(-x), mw, mx); vec3 v = fma(vec3(-y), mw, my); vec3 d = cross(u, v); return normalize(d); } ``` 这个推导有一个注意事项：我们没有说明 \(d\) 的符号，所以我们的光线可能指向后方。至少，反转 `world_to_clip` 的符号不会改变 `d` 的符号，但其他操作（如交换 `mx` 和 `my` 对应的行）会改变。总之，只要摄像机不产生镜像图像，上述代码将给出预期的符号。如果你想要无论如何都正确的结果，请继续阅读，我们将回到这个问题。 ## 光线方向的更快解法 代码列表 2 中仍有微小的优化空间。计算 `d` 需要 \(3 + 3 + 6 = 12\) 个融合乘加 (FMA) 指令。通过将部分计算移到着色器传入的常量中，可以减少到 6 个 FMA。如果我们拥有 10 万亿次浮点运算的 GPU 吞吐量（现代标准下并不算高），在 4k 分辨率下每个像素省一次，并且我们确实受限于计算吞吐量，那么每帧节省的时间为 \[\frac{3840 \cdot 2160 \cdot 6}{10\,\mathrm{THz}} = 0.005\,\mathrm{ms}.\] 每帧节省 5 微秒确实可以忽略，因此坚持使用代码列表 2 完全没问题。但也许你需要运行千...