第7章：光线追踪基础

光线追踪是计算机图形学中实现真实感渲染的核心技术之一。不同于光栅化的前向渲染方式，光线追踪采用逆向追踪光线的方式，从相机出发追踪光线在场景中的传播路径。本章将深入探讨光线追踪的基本原理、加速结构的设计与实现，以及光线与物体相交算法的优化技术。我们将重点关注算法的数学基础和性能优化策略，为后续的全局光照和高级渲染技术打下坚实基础。

7.1 光线追踪基本原理

7.1.1 光线的数学表示

光线可以用参数方程表示： $$\mathbf{r}(t) = \mathbf{o} + t\mathbf{d}$$ 其中：

• $\mathbf{o}$ 是光线起点（origin）
• $\mathbf{d}$ 是归一化的方向向量（direction）
• $t \geq 0$ 是参数，表示沿光线方向的距离

扩展表示：为了数值稳定性，实际应用中常使用： $$\mathbf{r}(t) = \mathbf{o} + t\mathbf{d}, \quad t \in [t_{\min}, t_{\max}]$$ 其中 $t_{\min} > 0$ 避免自相交，$t_{\max}$ 限制光线长度。典型值：$t_{\min} = 10^{-4}$ 到 $10^{-3}$，取决于场景规模。

光线微分：用于抗锯齿和纹理过滤，追踪光线的微分信息： $$\frac{\partial \mathbf{r}}{\partial x} = \frac{\partial \mathbf{o}}{\partial x} + t\frac{\partial \mathbf{d}}{\partial x}$$ 这对于计算纹理footprint至关重要。考虑相邻像素的光线： $$\begin{aligned} \frac{\partial \mathbf{d}}{\partial x} &\approx \mathbf{d}_{i+1,j} - \mathbf{d}_{i,j} \ \frac{\partial \mathbf{d}}{\partial y} &\approx \mathbf{d}_{i,j+1} - \mathbf{d}_{i,j} \end{aligned}$$ 光线锥（Ray Cone）表示： 为了更准确的LOD选择和过滤，可以将光线扩展为圆锥： $$r(t) = r_0 + t \cdot \tan(\theta)$$ 其中 $r_0$ 是像素footprint在近平面的半径，$\theta$ 是锥角，通常： $$\theta \approx \frac{1}{\text{image_width}} \cdot \text{fov}$$ 光线的齐次坐标表示： 在某些情况下（如投影变换），使用Plücker坐标更方便： $$\mathbf{L} = (\mathbf{d}, \mathbf{o} \times \mathbf{d})$$ 这种表示在处理光线-光线距离和相交测试时特别有用。

7.1.2 基本光线追踪算法

光线追踪的核心步骤：

1. 主光线生成（Primary Ray Generation） 对于图像平面上的每个像素 $(i, j)$，生成从相机位置出发的主光线： $$\mathbf{d}_{i,j} = \text{normalize}(\mathbf{p}_{i,j} - \mathbf{e})$$ 其中 $\mathbf{e}$ 是相机位置，$\mathbf{p}_{i,j}$ 是像素在世界坐标系中的位置。

子像素采样：为了抗锯齿，在像素内部进行多次采样： $$\mathbf{p}_{i,j}^{(k)} = \mathbf{p}_{i,j} + \xi_x \Delta_x \mathbf{right} + \xi_y \Delta_y \mathbf{up}$$ 其中 $\xi_x, \xi_y \in [-0.5, 0.5]$ 是采样偏移。

采样模式：

• 规则采样：$\xi_x = \frac{k_x + 0.5}{n_x} - 0.5$
• 随机采样：$\xi_x \sim U(-0.5, 0.5)$
• 分层采样（Stratified）：将像素分为 $n \times n$ 子格，每格内随机采样
• 低差异序列：使用Halton或Sobol序列

重要性采样优化： 根据前帧或预览pass的信息，在高频区域增加采样： $$p(\xi) \propto |\nabla I(x, y)|^2 + \epsilon$$

2. 光线-场景相交测试 对每条光线，找到最近的相交点： $$t_{\text{hit}} = \min_{k \in \text{objects}} t_k$$ 其中 $t_k$ 是光线与第 $k$ 个物体的相交参数。

相交信息记录： struct HitRecord { float t; // 相交参数 vec3 point; // 相交点 vec3 normal; // 表面法线 vec2 uv; // 纹理坐标 Material* mat; // 材质指针 float pdf_area; // 面积概率密度 vec3 dpdu, dpdv; // 表面参数化导数 vec3 tangent; // 切线向量 float curvature; // 曲率信息 }

浮点误差处理： 相交点的精确计算考虑误差界： $$\mathbf{p}_{\text{hit}} = \mathbf{o} + t_{\text{hit}}\mathbf{d} \pm \epsilon_{\text{hit}}$$ 其中 $\epsilon_{\text{hit}} = \gamma_3(|\mathbf{o}| + t_{\text{hit}}|\mathbf{d}|)$， $\gamma_n = \frac{n\epsilon_{\text{machine}}}{1 - n\epsilon_{\text{machine}}}$

3. 着色计算 在相交点处计算颜色，考虑：

• 直接光照（使用阴影光线）
• 反射（递归追踪反射光线）
• 折射（对透明物体）
• 环境光照（环境贴图采样）

光线相干性利用：

• 光线包（Ray Packets）：同时追踪多条相邻光线
• 光线排序：按方向或起点聚类以改善缓存局部性
• 光线重排序缓冲：延迟处理以批量化相似操作

7.1.3 递归光线追踪

递归光线追踪通过追踪次级光线来模拟复杂的光照效果： $$L_o(\mathbf{x}, \omega_o) = L_e(\mathbf{x}, \omega_o) + \int_{\Omega} f_r(\mathbf{x}, \omega_i, \omega_o) L_i(\mathbf{x}, \omega_i) (\omega_i \cdot \mathbf{n}) d\omega_i$$ Whitted-style光线追踪简化为： $$L_o = L_e + k_d L_{\text{direct}} + k_s L_{\text{reflect}} + k_t L_{\text{refract}}$$ 其中：

• $L_{\text{direct}}$：直接光照（Lambert + Phong）
• $L_{\text{reflect}}$：镜面反射贡献
• $L_{\text{refract}}$：折射贡献
• $k_d, k_s, k_t$：材质参数，满足能量守恒 $k_d + k_s + k_t \leq 1$

直接光照计算的细节：

对于点光源： $$L_{\text{direct}} = \sum_{i} \frac{I_i}{r_i^2} \cdot V_i \cdot \max(0, \mathbf{n} \cdot \mathbf{l}_i) \cdot \left(\frac{k_d}{\pi} + k_s \frac{(\mathbf{h}_i \cdot \mathbf{n})^n}{(\mathbf{n} \cdot \mathbf{l}_i)}\right)$$ 其中：

• $I_i$：光源强度
• $r_i$：到光源距离
• $V_i$：可见性函数（0或1）
• $\mathbf{h}_i = \text{normalize}(\mathbf{l}_i + \mathbf{v})$：半向量
• $n$：Phong指数

反射方向计算： $$\mathbf{r}_{\text{reflect}} = \mathbf{d} - 2(\mathbf{d} \cdot \mathbf{n})\mathbf{n}$$ 改进的反射模型： 考虑微表面理论，可以添加粗糙度扰动： $$\mathbf{r}'_{\text{reflect}} = \text{normalize}(\mathbf{r}_{\text{reflect}} + \alpha \mathbf{\xi})$$ 其中 $\mathbf{\xi}$ 是随机扰动向量，$\alpha$ 是粗糙度参数。

扰动向量的生成： 使用切线空间：

1. 构建正交基：${\mathbf{n}, \mathbf{t}, \mathbf{b}}$
2. 在半球内采样：$\mathbf{\xi} = \sqrt{1-u^2}\cos(2\pi v)\mathbf{t} + \sqrt{1-u^2}\sin(2\pi v)\mathbf{b} + u\mathbf{n}$
3. 其中 $u, v \in [0,1]$ 是随机数

折射方向计算（Snell定律）： $$\mathbf{r}_{\text{refract}} = \frac{\eta_i}{\eta_t}\mathbf{d} + \left(\frac{\eta_i}{\eta_t}(\mathbf{n} \cdot \mathbf{d}) - \sqrt{1 - \sin^2\theta_t}\right)\mathbf{n}$$ 其中 $\sin^2\theta_t = \left(\frac{\eta_i}{\eta_t}\right)^2(1 - (\mathbf{n} \cdot \mathbf{d})^2)$

稳定的折射计算：

令 k = 1 - η²(1 - (n·d)²)
if k < 0: 全内反射
else: t = η*d + (η*(n·d) - √k)*n

全内反射判断：当 $\sin^2\theta_t > 1$ 时发生全内反射。

临界角：$\theta_c = \arcsin(\eta_t/\eta_i)$（当 $\eta_i > \eta_t$）

Fresnel方程： 计算反射和折射的能量分配： $$F_r = \frac{1}{2}\left[\left(\frac{\eta_i\cos\theta_i - \eta_t\cos\theta_t}{\eta_i\cos\theta_i + \eta_t\cos\theta_t}\right)^2 + \left(\frac{\eta_t\cos\theta_i - \eta_i\cos\theta_t}{\eta_t\cos\theta_i + \eta_i\cos\theta_t}\right)^2\right]$$ Schlick近似（计算效率更高）： $$F_r \approx F_0 + (1 - F_0)(1 - \cos\theta_i)^5$$ 其中 $F_0 = \left(\frac{\eta_i - \eta_t}{\eta_i + \eta_t}\right)^2$

金属材质的Fresnel： 复数折射率 $\tilde{\eta} = \eta + i\kappa$ $$F_r = \frac{(\eta - 1)^2 + \kappa^2}{(\eta + 1)^2 + \kappa^2}$$ 色散效果： 不同波长具有不同折射率： $$\eta(\lambda) = A + \frac{B}{\lambda^2} + \frac{C}{\lambda^4}$$ （Cauchy方程）

典型值（玻璃）：

• 红光（700nm）：$\eta_r \approx 1.513$
• 绿光（550nm）：$\eta_g \approx 1.517$
• 蓝光（450nm）：$\eta_b \approx 1.522$

实现色散：

1. 将白光分解为RGB分量
2. 对每个分量使用不同折射率
3. 分别追踪三条折射光线
4. 合成最终颜色

递归深度控制：

简单深度限制：

if (depth > MAX_DEPTH) return background_color;

基于贡献的终止：

if (throughput < THRESHOLD) return black;

俄罗斯轮盘赌终止策略：

float p = max(color.r, color.g, color.b);
if (random() > p) return black;
return color / p;  // 能量补偿

改进的终止策略： 考虑路径贡献和计算成本的平衡： $$p_{\text{continue}} = \min\left(1, \frac{\text{throughput} \cdot \text{max_radiance}}{\text{current_radiance}}\right)$$ 自适应深度策略：

• 镜面反射：允许更深递归
• 漫反射：较早终止
• 基于材质重要性调整深度限制

光线树优化： 记录光线贡献路径，避免重复计算：

struct RayNode {
    Ray ray;
    float contribution;
    Material* mat;
    RayNode* children[2];  // reflect, refract
}

多次反射间的相互作用：

• 焦散（Caustics）：光线经过折射/反射聚焦
• 间接光照：多次漫反射
• 色彩渗透：相邻表面的颜色互相影响

这些效果在Whitted模型中是近似的，更准确的模拟需要路径追踪或光子映射。

7.1.4 相机模型与光线生成

透视投影相机的光线生成：

给定视场角 $\text{fov}$、宽高比 $\text{aspect}$，像素 $(i, j)$ 对应的光线方向： $$\begin{aligned} u &= \frac{2i - \text{width}}{\text{width}} \cdot \text{aspect} \cdot \tan(\text{fov}/2) \ v &= \frac{2j - \text{height}}{\text{height}} \cdot \tan(\text{fov}/2) \ \mathbf{d} &= \text{normalize}(u\mathbf{right} + v\mathbf{up} - \mathbf{forward}) \end{aligned}$$ 改进的相机模型考虑：

• 传感器偏移：模拟tilt-shift效果
• 非方形像素：像素宽高比校正
• 畸变模型：径向和切向畸变

薄透镜相机模型（景深效果）：

1. 计算焦平面上的目标点： $$\mathbf{p}_{\text{focus}} = \mathbf{o} + \frac{\text{focus_distance}}{|\mathbf{d} \cdot \mathbf{forward}|} \mathbf{d}$$

2. 在透镜上采样： $$\mathbf{o}' = \mathbf{o} + r(\cos\theta \mathbf{right} + \sin\theta \mathbf{up})$$ 其中 $r \leq \text{aperture}/2$，采样分布可以是：

• 均匀圆盘：$r = \sqrt{\xi_1} \cdot \text{aperture}/2$, $\theta = 2\pi\xi_2$
• 多边形光圈：模拟真实光圈叶片
• 自定义形状：实现散景（bokeh）效果

3. 新的光线方向： $$\mathbf{d}' = \text{normalize}(\mathbf{p}_{\text{focus}} - \mathbf{o}')$$ 景深的物理参数：

• 弥散圆直径：$c = \frac{A|S_2 - S_1|}{S_2}$
• 景深范围：$\text{DOF} = \frac{2Nc(S - f)^2}{f^2 - N^2c^2}$ 其中 $N$ 是光圈数，$S$ 是对焦距离，$f$ 是焦距

鱼眼相机： $$\begin{aligned} r &= \sqrt{u^2 + v^2} \ \theta &= r \cdot \text{fov} / 2 \ \phi &= \arctan2(v, u) \ \mathbf{d} &= \sin\theta\cos\phi \mathbf{right} + \sin\theta\sin\phi \mathbf{up} - \cos\theta \mathbf{forward} \end{aligned}$$ 其他投影模型：

• 等距投影：$\theta = r$
• 等立体角投影：$\theta = 2\sin^{-1}(r/2)$
• 正交投影：$\theta = \sin^{-1}(r)$
• 立体投影：$\theta = 2\tan^{-1}(r/2)$

全景相机：

• 等矩形投影：$(u, v) \rightarrow (\theta, \phi)$ 直接映射
• 立方体贴图：6个透视投影相机
• 双抛物面映射：前后两个抛物面投影

运动模糊： 在快门时间内对相机位置和方向进行积分： $$\mathbf{o}(t) = \mathbf{o}_0 + t \cdot \mathbf{v}$$ $$\mathbf{d}(t) = \mathbf{R}(t) \cdot \mathbf{d}_0$$ 其中 $t \in [0, \text{shutter_time}]$

7.2 加速结构（BVH、KD-Tree）

7.2.1 空间数据结构的必要性

朴素的光线追踪需要测试每条光线与场景中所有物体的相交，复杂度为 $O(N)$。对于包含百万级三角形的场景，这是不可接受的。空间加速结构可以将平均复杂度降低到 $O(\log N)$。

理论分析：

• 无加速结构：每条光线需要 $N$ 次相交测试
• 理想加速结构：$O(\log N)$ 次节点访问 + $O(1)$ 次图元测试
• 实际性能：取决于场景分布和构建质量

性能模型： 设光线追踪总时间 $T = T_{\text{build}} + R \cdot T_{\text{trace}}$，其中：

• $T_{\text{build}}$：构建时间
• $R$：光线数量
• $T_{\text{trace}}$：单条光线追踪时间

对于单条光线： $$T_{\text{trace}} = T_{\text{traverse}} + N_{\text{leaf}} \cdot T_{\text{intersect}}$$ 其中 $N_{\text{leaf}}$ 是访问的叶节点数。

空间连贯性原理：

• 光线连贯性：相邻光线倾向于相交相似的物体集合
• 几何连贯性：空间相近的物体应该在加速结构中相近
• 时间连贯性：连续帧之间的相交模式相似

常见加速结构分类：

1. 空间分割：KD-Tree、Octree、BSP Tree - 优点：无重叠，内存效率高 - 缺点：物体可能跨越多个节点

2. 物体分割：BVH、R-Tree - 优点：每个物体只属于一个叶节点 - 缺点：包围盒可能重叠

3. 混合方法：SBVH、Dual-Split BVH - 结合空间和物体分割的优点 - 更复杂的构建算法

加速结构选择准则：

• 静态场景：KD-Tree（遍历效率最高）
• 动态场景：BVH（易于更新）
• GPU实现：BVH（更规则的内存访问）
• 内存受限：压缩BVH或隐式Grid

7.2.2 层次包围盒（BVH）

BVH是一种物体层次结构，每个节点包含其所有子节点的轴对齐包围盒（AABB）。

BVH节点结构：

struct BVHNode {
    AABB bounds;          // 包围盒
    union {
        struct {          // 内部节点
            BVHNode* left;
            BVHNode* right;
        };
        struct {          // 叶节点
            int primOffset;   // 图元起始索引
            int primCount;    // 图元数量
        };
    };
}

BVH构建算法：

1. 自顶向下构建 ``` function BuildBVH(primitives, start, end): if (end - start <= leaf_threshold): return CreateLeaf(primitives[start:end])

   axis = ChooseSplitAxis(primitives[start:end]) mid = Partition(primitives, start, end, axis)

   left = BuildBVH(primitives, start, mid) right = BuildBVH(primitives, mid, end)

   return CreateNode(left, right) ```

2. SAH（Surface Area Heuristic） 最优分割的代价函数： $$C = C_{\text{trav}} + \frac{A_L}{A} \cdot N_L \cdot C_{\text{isect}} + \frac{A_R}{A} \cdot N_R \cdot C_{\text{isect}}$$ 其中：

• $C_{\text{trav}}$ 是遍历节点的代价（典型值：1.0）
• $C_{\text{isect}}$ 是相交测试的代价（典型值：4.0）
• $A_L, A_R$ 是左右子节点的表面积
• $N_L, N_R$ 是左右子节点包含的图元数

完整SAH实现考虑： $$C_{\text{split}} = C_{\text{trav}} + P_L \cdot C_L + P_R \cdot C_R$$ 其中概率 $P_L = \frac{A_L}{A}$，$P_R = \frac{A_R}{A}$

表面积计算（AABB）： $$A = 2(dx \cdot dy + dy \cdot dz + dz \cdot dx)$$ 其中 $dx = \text{max}_x - \text{min}_x$ 等。

3. 分割策略比较：

中点分割（Middle Split）：

• 简单快速：$p = (\text{min} + \text{max}) / 2$
• 可能产生不平衡的树

等数分割（Equal Count）：

• 保证平衡：每边 $N/2$ 个图元
• 可能产生重叠严重的包围盒

SAH分割：

• 最小化期望遍历代价
• 计算开销大，但质量最高

混合策略： if (图元数 < 阈值) 使用中点分割 else 使用SAH分割

4. 高效SAH实现：

桶排序方法（Binned SAH）： ``` 将包围盒范围分为K个桶（典型K=32） for each 图元: 计算质心 分配到对应桶

for each 可能的分割位置（K-1个）: 计算左右两边的图元数和包围盒 评估SAH代价 选择最小代价的分割 ```

复杂度：$O(KN)$ 而非 $O(N^2)$

5. BVH遍历

递归遍历： ``` function Intersect(node, ray, tmin, tmax): if (IsLeaf(node)): return IntersectPrimitives(node.primitives, ray)

   t1 = IntersectAABB(node.left.bbox, ray)
   t2 = IntersectAABB(node.right.bbox, ray)

   if (t1.hit && t2.hit):
       first = (t1.tmin < t2.tmin) ? node.left : node.right
       second = (t1.tmin < t2.tmin) ? node.right : node.left

       hit1 = Intersect(first, ray, tmin, tmax)
       if (hit1.t < second.tmin) return hit1

       hit2 = Intersect(second, ray, tmin, tmax)
       return Closer(hit1, hit2)

```

栈式遍历（更适合GPU）： stack[0] = root while (stack not empty): node = stack.pop() if (IntersectAABB(node.bbox, ray)): if (IsLeaf(node)): ProcessPrimitives(node) else: stack.push(node.farChild) stack.push(node.nearChild)

优化的遍历顺序： 基于光线方向符号预计算子节点访问顺序： dirIsNeg[3] = {ray.dir.x < 0, ray.dir.y < 0, ray.dir.z < 0} // 使用dirIsNeg[node.axis]决定先访问哪个子节点

6. BVH质量评估：

SAH代价： $$\text{SAH}(T) = \sum_{\text{leaves}} P(\text{leaf}) \cdot N(\text{leaf})$$ EPO（Expected Primary Operations）： 平均每条光线的期望操作数

树平衡度： $$\text{Balance} = \frac{\text{avg_depth}}{\log_2 N}$$

7. 高级BVH技术：

SBVH（Spatial Split BVH）：

• 允许空间分割，不仅是物体分割
• 可以减少包围盒重叠
• 代价：可能增加图元引用数

实现策略： if (SAH代价 > 阈值 && 重叠率 > 阈值): 尝试空间分割 选择最优的物体或空间分割

压缩BVH：

• 量化包围盒坐标（16位或8位）
• 节点合并（将多个节点打包）
• 典型压缩率：50-75%

量化公式： $$\text{quantized} = \text{round}\left(\frac{\text{value} - \text{parent_min}}{\text{parent_max} - \text{parent_min}} \cdot (2^{bits} - 1)\right)$$ Wide BVH：

• 每个节点有4个或8个子节点
• 更好的SIMD利用率
• 减少内存访问次数

8. 动态场景的BVH更新：

重拟合（Refit）： 自底向上更新包围盒 function Refit(node): if (IsLeaf(node)): node.bounds = ComputeBounds(node.primitives) else: Refit(node.left) Refit(node.right) node.bounds = Union(node.left.bounds, node.right.bounds)

局部重建：

• 标记"脏"节点
• 只重建变化较大的子树
• 使用启发式判断是否需要重建

双缓冲策略：

• 维护两个BVH
• 异步更新非活动BVH
• 快速切换

9. 内存布局优化：

深度优先布局：

• 改善缓存局部性
• 适合递归遍历

van Emde Boas布局：

• 优化缓存性能
• 适合任意遍历模式

节点压缩布局： struct CompactBVHNode { float bounds_min[3]; float bounds_max[3]; uint32_t offset; // 子节点偏移或图元偏移 uint8_t nPrims; // 0表示内部节点 uint8_t axis; // 分割轴 uint16_t pad; // 对齐 } // 32字节，正好两个缓存行

7.2.3 KD-Tree

KD-Tree特点：

• 空间分割而非物体分割
• 分割平面与坐标轴对齐
• 可能需要处理跨越分割平面的三角形

构建策略：

1. 分割位置选择：

中位数分割： $$p_{\text{split}} = \text{median}({p_i \cdot \mathbf{axis} : i \in \text{primitives}})$$ SAH分割： 最小化代价函数： $$C(p) = C_{\text{trav}} + \frac{V_L(p)}{V} N_L(p) C_{\text{isect}} + \frac{V_R(p)}{V} N_R(p) C_{\text{isect}}$$ 空盒优化： $$C_{\text{empty}} = 0.8 \cdot C_{\text{trav}}$$

2. 精确SAH计算： for each candidate position p: leftCount = CountPrimitivesLeft(p) rightCount = CountPrimitivesRight(p) leftVolume = ComputeVolume(min, p) rightVolume = ComputeVolume(p, max) cost = EvaluateSAH(leftCount, rightCount, leftVolume, rightVolume)

3. KD-Tree遍历算法： ``` function TraverseKDTree(ray, tmin, tmax, node): if (IsLeaf(node)): return IntersectPrimitives(node.primitives, ray)

   axis = node.splitAxis t_split = (node.splitPos - ray.origin[axis]) / ray.dir[axis]

   nearNode = (ray.origin[axis] < node.splitPos) ? node.left : node.right farNode = (ray.origin[axis] < node.splitPos) ? node.right : node.left

   if (t_split > tmax || t_split < 0): return TraverseKDTree(ray, tmin, tmax, nearNode) if (t_split < tmin): return TraverseKDTree(ray, tmin, tmax, farNode)

   hit = TraverseKDTree(ray, tmin, t_split, nearNode) if (hit.valid) return hit

   return TraverseKDTree(ray, t_split, tmax, farNode) ```

遍历优化：

• Mailboxing：避免重复相交测试 struct Mailbox { int rayID; float t; }

• 早期终止：当找到相交点后提前退出

• Rope技术：存储邻居指针加速遍历

7.2.4 加速结构比较

| 特性 | BVH | KD-Tree | Octree |

性能模型： 给定场景包含 $N$ 个图元，光线数量为 $R$：

1. 构建时间： - BVH: $T_{\text{build}} = O(N\log N)$ - KD-Tree: $T_{\text{build}} = O(N\log^2 N)$（SAH）

2. 遍历时间： - 期望深度：$D = O(\log N)$ - 每条光线：$T_{\text{ray}} = D \cdot C_{\text{node}} + L \cdot C_{\text{leaf}}$ - 总时间：$T_{\text{total}} = R \cdot T_{\text{ray}}$

3. 内存访问模式： - BVH：更好的空间局部性 - KD-Tree：更深的树，更多cache miss

混合加速结构：

顶层：BVH（场景级别）
  ├── 中层：KD-Tree（物体级别）
  └── 底层：Grid（密集三角形）

7.3 光线-物体相交算法优化

7.3.1 光线-三角形相交

光线-三角形相交是光线追踪中最频繁的操作，其性能直接影响渲染速度。

Möller-Trumbore算法：

给定三角形顶点 $\mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2$，光线 $\mathbf{r}(t) = \mathbf{o} + t\mathbf{d}$： $$\begin{aligned} \mathbf{e}_1 &= \mathbf{v}_1 - \mathbf{v}_0 \ \mathbf{e}_2 &= \mathbf{v}_2 - \mathbf{v}_0 \ \mathbf{h} &= \mathbf{d} \times \mathbf{e}_2 \ a &= \mathbf{e}_1 \cdot \mathbf{h} \end{aligned}$$ 如果 $|a| < \epsilon$，光线与三角形平行。否则： $$\begin{aligned} f &= 1/a \ \mathbf{s} &= \mathbf{o} - \mathbf{v}_0 \ u &= f(\mathbf{s} \cdot \mathbf{h}) \end{aligned}$$ 如果 $u < 0$ 或 $u > 1$，无相交。继续： $$\begin{aligned} \mathbf{q} &= \mathbf{s} \times \mathbf{e}_1 \ v &= f(\mathbf{d} \cdot \mathbf{q}) \end{aligned}$$ 如果 $v < 0$ 或 $u + v > 1$，无相交。否则： $$t = f(\mathbf{e}_2 \cdot \mathbf{q})$$ 算法复杂度分析：

• 1次除法，27次乘法，17次加减法
• 2次叉积，3次点积
• 分支预测友好（早期退出）

Baldwin-Weber算法（避免除法）： 使用齐次坐标避免除法： $$\begin{aligned} \mathbf{o}' &= \mathbf{o} - \mathbf{v}_0 \ \mathbf{q}_1 &= \mathbf{d} \times \mathbf{e}_2 \ \mathbf{q}_2 &= \mathbf{o}' \times \mathbf{e}_1 \ \text{det} &= \mathbf{q}_1 \cdot \mathbf{e}_1 \ u &= \mathbf{q}_1 \cdot \mathbf{o}' \ v &= \mathbf{q}_2 \cdot \mathbf{d} \ t &= \mathbf{q}_2 \cdot \mathbf{e}_2 \end{aligned}$$ 相交条件（避免除法）：

if (det > 0):
    return (u >= 0) && (v >= 0) && (u + v <= det) && (t >= 0)
else:
    return (u <= 0) && (v <= 0) && (u + v >= det) && (t <= 0)

Woop算法（预计算优化）： 预计算变换矩阵，将三角形变换到单位三角形： $$\mathbf{M} = \begin{bmatrix} \mathbf{v}_0 - \mathbf{v}_2 & \mathbf{v}_1 - \mathbf{v}_2 & \mathbf{v}_2 \end{bmatrix}^{-1}$$ 光线变换： $$\begin{aligned} \mathbf{o}' &= \mathbf{M} \cdot \mathbf{o} \ \mathbf{d}' &= \mathbf{M} \cdot \mathbf{d} \end{aligned}$$ 相交测试简化为： $$t = -\frac{o'_z}{d'_z}, \quad u = o'_x + t \cdot d'_x, \quad v = o'_y + t \cdot d'_y$$ 相交条件：$u \geq 0, v \geq 0, u + v \leq 1$

SIMD优化： 同时测试4个三角形：

__m128 e1x = _mm_sub_ps(v1x, v0x);
__m128 e1y = _mm_sub_ps(v1y, v0y);
__m128 e1z = _mm_sub_ps(v1z, v0z);
// ... 继续向量化计算

精度改进： 使用Shewchuk的精确谓词算法处理边界情况：

• 精确的方向测试
• 鲁棒的重心坐标计算
• 避免数值不一致

7.3.2 光线-包围盒相交

Slab方法：

对于轴对齐包围盒（AABB），计算光线与每对平行平面的相交： $$\begin{aligned} t_{\text{min}} &= \max(t_{x,\text{min}}, t_{y,\text{min}}, t_{z,\text{min}}) \ t_{\text{max}} &= \min(t_{x,\text{max}}, t_{y,\text{max}}, t_{z,\text{max}}) \end{aligned}$$ 相交条件：$t_{\text{min}} \leq t_{\text{max}}$ 且 $t_{\text{max}} \geq 0$

优化实现：

invD = 1.0 / ray.dir  // 预计算
t0 = (box.min - ray.origin) * invD
t1 = (box.max - ray.origin) * invD

// 处理负方向
tmin = min(t0, t1)
tmax = max(t0, t1)

// 计算最终区间
tmin_final = max(tmin.x, tmin.y, tmin.z)
tmax_final = min(tmax.x, tmax.y, tmax.z)

return tmin_final <= tmax_final && tmax_final >= 0

Williams等人的优化： 消除min/max中的分支：

sign[3] = {ray.dir.x < 0, ray.dir.y < 0, ray.dir.z < 0}
tmin = (box[sign[0]].x - ray.origin.x) * invD.x
tmax = (box[1-sign[0]].x - ray.origin.x) * invD.x

SSE优化版本：

__m128 tmin = _mm_mul_ps(_mm_sub_ps(box_min, origin), inv_dir);
__m128 tmax = _mm_mul_ps(_mm_sub_ps(box_max, origin), inv_dir);
__m128 t0 = _mm_min_ps(tmin, tmax);
__m128 t1 = _mm_max_ps(tmin, tmax);

鲁棒性改进： 处理光线方向分量为0的情况：

if (abs(ray.dir.x) < epsilon):
    if (ray.origin.x < box.min.x || ray.origin.x > box.max.x):
        return false
    tx_min = -infinity
    tx_max = +infinity

7.3.3 光线-球体相交

球体方程：$||\mathbf{p} - \mathbf{c}||^2 = r^2$

代入光线方程得到二次方程： $$at^2 + bt + c = 0$$ 其中： $$\begin{aligned} a &= \mathbf{d} \cdot \mathbf{d} = 1 \text{ (归一化方向)} \ b &= 2\mathbf{d} \cdot (\mathbf{o} - \mathbf{c}) \ c &= ||\mathbf{o} - \mathbf{c}||^2 - r^2 \end{aligned}$$ 判别式 $\Delta = b^2 - 4ac$，相交参数： $$t = \frac{-b \pm \sqrt{\Delta}}{2a}$$ 几何优化方法：

L = center - ray.origin
tca = dot(L, ray.dir)
if (tca < 0) return false  // 球在光线后方

d2 = dot(L, L) - tca * tca
if (d2 > radius2) return false  // 光线错过球体

thc = sqrt(radius2 - d2)
t0 = tca - thc
t1 = tca + thc

数值稳定性改进： 使用以下形式避免相减相消： $$t = \frac{-b - \text{sign}(b)\sqrt{\Delta}}{2a}$$ 椭球体相交： 使用变换矩阵 $\mathbf{M}$ 将椭球变换为单位球： $$\begin{aligned} \mathbf{o}' &= \mathbf{M}(\mathbf{o} - \mathbf{c}) \ \mathbf{d}' &= \mathbf{M}\mathbf{d} \end{aligned}$$ 然后使用球体相交算法。

7.3.4 其他几何体相交

光线-圆柱相交： 无限圆柱：$(x - c_x)^2 + (z - c_z)^2 = r^2$

相交计算： $$a = d_x^2 + d_z^2$$ $$b = 2(d_x(o_x - c_x) + d_z(o_z - c_z))$$ $$c = (o_x - c_x)^2 + (o_z - c_z)^2 - r^2$$ 有限圆柱需要额外检查端盖。

光线-圆锥相交： 圆锥方程：$x^2 + z^2 = (r(y/h))^2$

二次方程系数： $$a = d_x^2 + d_z^2 - (r/h)^2 d_y^2$$ $$b = 2(o_x d_x + o_z d_z - (r/h)^2 o_y d_y)$$ $$c = o_x^2 + o_z^2 - (r/h)^2 o_y^2$$ 光线-环面（Torus）相交： 四次方程，使用数值方法求解： $$(x^2 + y^2 + z^2 + R^2 - r^2)^2 = 4R^2(x^2 + z^2)$$ 可以使用Sturm序列或Ferrari方法求解。

7.3.5 相交测试优化策略

1. 早期拒绝（Early Rejection） - 使用简单包围体进行预测试 - 利用空间连贯性 - 背面剔除（对闭合物体）

2. SIMD并行化 - 同时测试多条光线（Ray Packets） - 批量处理三角形 - 使用SOA（Structure of Arrays）布局

3. 缓存优化 - 数据结构对齐（避免false sharing） - 热数据分离（将常用数据打包） - 预取（prefetch）关键数据

4. 精度考虑 - 使用稳定的数值算法 - 处理自相交问题（shadow acne） - 考虑舍入误差的传播

5. 分支优化 - 减少条件分支（使用位运算） - 使用模板特化消除运行时分支 - Profile-guided optimization

6. 间隔算术（Interval Arithmetic） 用于保守估计和鲁棒性： struct Interval { float min, max; Interval operator*(const Interval& b) { float products[4] = {min*b.min, min*b.max, max*b.min, max*b.max}; return {min(products), max(products)}; } }

7. 自适应精度 - 远处物体使用低精度 - 重要区域使用高精度 - 基于误差估计动态调整

8. 光线差分（Ray Differentials） 追踪光线的微分信息用于：

• LOD选择
• 纹理过滤
• 自适应细分 $$\frac{\partial \mathbf{p}}{\partial u} = t \frac{\partial \mathbf{d}}{\partial u} + \frac{\partial t}{\partial u} \mathbf{d}$$

本章小结

本章介绍了光线追踪的核心概念和关键技术：

核心概念：

• 光线的参数表示：$\mathbf{r}(t) = \mathbf{o} + t\mathbf{d}$
• 递归光线追踪实现全局光照效果
• 空间加速结构将复杂度从 $O(N)$ 降至 $O(\log N)$

关键算法：

• Möller-Trumbore 三角形相交算法
• SAH启发式用于构建最优BVH
• Slab方法用于包围盒相交测试

性能优化：

• 使用BVH或KD-Tree加速结构
• SIMD并行化相交测试
• 缓存友好的数据布局

重要公式：

• SAH代价函数：$C = C_{\text{trav}} + \sum_i \frac{A_i}{A} N_i C_{\text{isect}}$
• 透视投影光线生成
• 各种几何体的相交测试公式

光线追踪为后续的全局光照、路径追踪等高级渲染技术奠定了基础。掌握本章内容是理解现代渲染技术的关键。

练习题

基础题

练习7.1 推导正交投影相机的光线生成公式。

练习7.2 给定光线 $\mathbf{o} = (0, 0, 0)$, $\mathbf{d} = (1, 0, 0)$ 和三角形顶点 $\mathbf{v}_0 = (2, -1, -1)$, $\mathbf{v}_1 = (2, 1, -1)$, $\mathbf{v}_2 = (2, 0, 1)$，计算相交点和重心坐标。

练习7.3 证明BVH中使用SAH构建的树，期望遍历代价是最优的。

练习7.4 设计一个算法，检测光线是否穿过由多个三角形组成的封闭网格。

挑战题

练习7.5 推导四叉BVH（4个子节点）相对于二叉BVH的理论优势和劣势。在什么情况下四叉BVH更优？

练习7.6 设计并分析一种自适应的光线追踪算法，能够根据场景复杂度动态调整采样密度。

练习7.7 分析光线追踪中的数值精度问题，并提出一套完整的解决方案。

练习7.8 探讨如何将光线追踪扩展到非欧几何空间（如球面几何或双曲几何）。

常见陷阱与错误

1. 数值精度问题

问题：自相交导致的阴影痤疮（Shadow Acne）

错误：shadow_ray.origin = hit_point
正确：shadow_ray.origin = hit_point + ε * normal

问题：光线方向未归一化导致t值含义错误

检查：assert(|ray.direction| ≈ 1.0)

2. 性能陷阱

问题：BVH构建时未考虑空间局部性

• 使用Morton编码改善缓存性能
• 考虑节点大小与缓存行对齐

问题：过早优化导致代码复杂

• 先实现正确的算法
• 使用性能分析工具定位瓶颈

3. 算法错误

问题：BVH遍历时错误的节点访问顺序

错误：总是先访问左子节点
正确：先访问距离较近的子节点

问题：忽略了背面剔除的边界情况

• 透明物体需要双面测试
• 法线翻转的处理

4. 边界情况

问题：退化三角形（共线顶点）导致除零错误

• 预处理移除退化几何
• 相交测试中添加检查

问题：极小或极大的场景规模

• 使用相对误差而非绝对误差
• 考虑分层LOD

最佳实践检查清单

设计阶段

• [ ] 选择合适的加速结构（BVH vs KD-Tree）
• [ ] 确定精度要求和误差容限
• [ ] 规划内存布局和缓存优化策略
• [ ] 考虑目标平台（CPU/GPU）特性

实现阶段

• [ ] 光线方向始终保持归一化
• [ ] 正确处理数值精度问题
• [ ] 实现鲁棒的相交测试算法
• [ ] 添加适当的debug可视化

优化阶段

• [ ] 使用性能分析确定瓶颈
• [ ] 实现SIMD加速的相交测试
• [ ] 优化内存访问模式
• [ ] 考虑光线包（Ray Packet）技术

验证阶段

• [ ] 测试各种退化情况
• [ ] 验证数值稳定性
• [ ] 比较不同场景下的性能
• [ ] 确保结果的确定性（用于调试）

维护阶段

• [ ] 文档化关键算法选择
• [ ] 提供性能调优参数
• [ ] 支持渐进式渲染
• [ ] 保持代码的可扩展性