第5章：欧拉视角（下）与大规模渲染

在上一章中，我们深入探讨了欧拉视角下的流体模拟核心算法，包括压力投影、网格求解器等基础技术。本章将继续深化欧拉方法的学习，重点关注高精度输送格式、自由表面追踪，并将视野扩展到物理引擎中的渲染技术。我们将探讨如何高效地可视化大规模物理模拟结果，这对于验证算法正确性和创建视觉特效至关重要。

学习目标

完成本章学习后，你将能够：

掌握高阶输送格式：理解WENO和BFECC等高精度输送方法的数学原理，能够分析其数值耗散特性
实现自由表面追踪：掌握等势面方法的理论基础，理解有符号距离场的演化与重初始化
理解现代渲染技术：从物理引擎角度理解路径追踪、球面追踪等渲染算法
优化大规模可视化：掌握体素渲染、DDA算法等技术，实现高效的物理场可视化
整合物理与渲染：理解物理模拟与渲染管线的交互，优化整体性能

5.1 高级输送格式

在流体模拟中，输送（advection）是最基础也是最具挑战性的操作之一。标准的半拉格朗日方法虽然无条件稳定，但会引入过多的数值耗散。本节介绍两种高精度输送格式：WENO（Weighted Essentially Non-Oscillatory）和BFECC（Back and Forth Error Compensation and Correction）。

5.1.1 WENO格式原理

WENO格式是一类高阶精度的数值格式，特别适合处理包含间断的流场。其核心思想是通过非线性权重组合多个候选模板，自适应地选择光滑区域。

考虑一维标量输送方程： $$\frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x} = 0$$ 对于空间导数的离散，WENO-5格式使用五个点构造三个候选模板：

候选模板：

$S_0$: 使用点 $\{i-2, i-1, i\}$
$S_1$: 使用点 $\{i-1, i, i+1\}$
$S_2$: 使用点 $\{i, i+1, i+2\}$

每个模板可以构造一个三阶精度的插值多项式。以$S_k$表示第k个模板，其对应的数值通量为： $$\hat{f}_k = \sum_{j=0}^{2} c_{kj} f_{i-k+j}$$ 其中系数$c_{kj}$由泰勒展开确定。具体地，对于正速度$u > 0$的情况： $$\begin{aligned} c_{00} &= \frac{1}{3}, \quad c_{01} = \frac{5}{6}, \quad c_{02} = -\frac{1}{6} \\ c_{10} &= -\frac{1}{6}, \quad c_{11} = \frac{5}{6}, \quad c_{12} = \frac{1}{3} \\ c_{20} &= \frac{1}{3}, \quad c_{21} = -\frac{7}{6}, \quad c_{22} = \frac{11}{6} \end{aligned}$$ 非线性权重： WENO的关键在于根据光滑度指标动态调整权重： $$\omega_k = \frac{\alpha_k}{\sum_{l=0}^{2} \alpha_l}, \quad \alpha_k = \frac{d_k}{(\epsilon + \beta_k)^2}$$ 其中：

$d_k$是理想权重（在光滑区域使格式达到五阶精度）：$d_0 = \frac{1}{10}$，$d_1 = \frac{6}{10}$，$d_2 = \frac{3}{10}$
$\beta_k$是光滑度指标，衡量模板$S_k$上解的光滑程度
$\epsilon$是防止除零的小参数（通常取$10^{-6}$）

光滑度指标的计算： $$\beta_k = \sum_{l=1}^{2} \int_{x_{i-1/2}}^{x_{i+1/2}} \Delta x^{2l-1} \left(\frac{\partial^l p_k}{\partial x^l}\right)^2 dx$$ 其中$p_k(x)$是模板$S_k$上的插值多项式。对于均匀网格，光滑度指标有显式表达式： $$\begin{aligned} \beta_0 &= \frac{13}{12}(f_{i-2} - 2f_{i-1} + f_i)^2 + \frac{1}{4}(f_{i-2} - 4f_{i-1} + 3f_i)^2 \\ \beta_1 &= \frac{13}{12}(f_{i-1} - 2f_i + f_{i+1})^2 + \frac{1}{4}(f_{i-1} - f_{i+1})^2 \\ \beta_2 &= \frac{13}{12}(f_i - 2f_{i+1} + f_{i+2})^2 + \frac{1}{4}(3f_i - 4f_{i+1} + f_{i+2})^2 \end{aligned}$$

5.1.2 BFECC方法

BFECC是一种简单而有效的误差补偿方法，通过前后输送来估计和补偿数值误差。其基本思想是利用输送算子的对称性来检测和修正数值误差。

算法步骤：

前向输送：将$\phi^n$输送到$n+1$时刻 $$\phi^{*} = \mathcal{A}(\phi^n, \Delta t)$$
后向输送：将$\phi^{*}$输送回$n$时刻 $$\phi^{**} = \mathcal{A}(\phi^{*}, -\Delta t)$$
误差估计：计算输送误差 $$e = \frac{1}{2}(\phi^n - \phi^{**})$$
误差补偿：修正前向输送结果 $$\tilde{\phi}^{*} = \phi^{*} + e$$
最终输送：使用修正后的值进行输送 $$\phi^{n+1} = \mathcal{A}(\tilde{\phi}^{*}, \Delta t)$$ 其中$\mathcal{A}(\cdot, \Delta t)$表示时间步长为$\Delta t$的输送算子。

数学分析：假设输送算子的截断误差为$O(\Delta x^p)$，则BFECC可以将精度提高到$O(\Delta x^{p+1})$。对于线性输送，可以证明： $$\phi^{n+1} = \phi_{exact}^{n+1} + O(\Delta x^{p+2})$$ 误差估计的理论基础：考虑泰勒展开，前向输送引入的误差为： $$\phi^* = \phi_{exact}^{n+1} + E_f + O(\Delta t^3)$$ 其中$E_f$是前向输送的主要误差项。后向输送得到： $$\phi^{**} = \phi^n - 2E_f + O(\Delta t^3)$$ 因此误差估计$e = \frac{1}{2}(\phi^n - \phi^{**}) = E_f + O(\Delta t^3)$准确捕获了主要误差。

MacCormack格式的关系： BFECC可以看作MacCormack格式的推广。MacCormack格式的预测-校正步骤： $$\begin{aligned} \phi^{*} &= \phi^n - \frac{\Delta t}{\Delta x}(f_{i+1}^n - f_i^n) \quad \text{(预测)} \\ \phi^{n+1} &= \frac{1}{2}[\phi^n + \phi^{*} - \frac{\Delta t}{\Delta x}(f_i^{*} - f_{i-1}^{*})] \quad \text{(校正)} \end{aligned}$$ BFECC通过完整的前后输送提供了更准确的误差估计。

5.1.3 误差分析与稳定性

数值耗散分析：对于波数$k$的正弦波，定义数值耗散因子： $$D(k) = \frac{|\hat{\phi}^{n+1}(k)|}{|\hat{\phi}^n(k)|}$$ 理想情况下$D(k) = 1$。对于不同格式：

一阶迎风：$D(k) = 1 - C(1-C)(1-\cos(k\Delta x))$，其中$C = u\Delta t/\Delta x$是Courant数
WENO-5：在光滑区域接近谱精度，$D(k) \approx 1 - O((k\Delta x)^{10})$
BFECC：显著减少低频耗散，$D(k) \approx 1 - O((k\Delta x)^{2p+2})$

频谱分析：考虑修正波数$k^*$，定义相速度误差： $$\epsilon_\phi = \frac{k^* - k}{k} = \frac{\arg(\hat{\phi}^{n+1}/\hat{\phi}^n)}{Ck\Delta x} - 1$$ 不同格式的相速度误差：

中心差分：$\epsilon_\phi = -\frac{(k\Delta x)^2}{6} + O((k\Delta x)^4)$
WENO-5：$\epsilon_\phi = O((k\Delta x)^4)$（光滑区域）
迎风格式：$\epsilon_\phi = -\frac{1-C}{2}k\Delta x + O((k\Delta x)^2)$

稳定性条件：

显式WENO：CFL条件 $|u|\Delta t / \Delta x \leq 1$
BFECC：继承基础输送算子的稳定性，但需要限制器防止新极值
隐式WENO：无条件稳定，但需要求解非线性系统

von Neumann稳定性分析：对于线性输送方程，增长因子$G = \hat{\phi}^{n+1}/\hat{\phi}^n$满足： $$|G| = |1 - C(1 - e^{-ik\Delta x})|$$ 稳定性要求$|G| \leq 1$，导出CFL条件。

单调性保持：为保证物理量的界限（如密度非负），需要加入限制器： $$\phi^{n+1} = \max(\phi_{min}, \min(\phi_{max}, \phi^{n+1}))$$ 其中$\phi_{min}$和$\phi_{max}$是邻域内的极值。

TVD限制器：总变差递减（TVD）性质保证不产生新的极值： $$TV(\phi^{n+1}) \leq TV(\phi^n)$$ 其中$TV(\phi) = \sum_i |\phi_{i+1} - \phi_i|$。常用的TVD限制器包括：

minmod限制器：$\psi(r) = \max(0, \min(1, r))$
van Leer限制器：$\psi(r) = \frac{r + |r|}{1 + |r|}$
superbee限制器：$\psi(r) = \max(0, \min(2r, 1), \min(r, 2))$

其中$r = (\phi_i - \phi_{i-1})/(\phi_{i+1} - \phi_i)$是连续梯度比。

5.2 等势面方法与自由表面追踪

自由表面流动（如水波、飞溅）是流体模拟中最具视觉冲击力的现象。等势面方法（Level Set Method）提供了一种优雅的数学框架来追踪移动界面。

5.2.1 有符号距离场

等势面方法的核心是用有符号距离函数$\phi(\mathbf{x}, t)$隐式表示界面： $$\phi(\mathbf{x}, t) = \begin{cases} +d(\mathbf{x}, \Gamma) & \text{如果 } \mathbf{x} \text{ 在界面外部} \\ 0 & \text{如果 } \mathbf{x} \text{ 在界面上} \\ -d(\mathbf{x}, \Gamma) & \text{如果 } \mathbf{x} \text{ 在界面内部} \end{cases}$$ 其中$d(\mathbf{x}, \Gamma)$是点$\mathbf{x}$到界面$\Gamma$的最短距离。

几何性质：

法向量：$\mathbf{n} = \frac{\nabla\phi}{|\nabla\phi|}$
曲率：$\kappa = \nabla \cdot \mathbf{n} = \nabla \cdot \left(\frac{\nabla\phi}{|\nabla\phi|}\right)$

对于二维情况，曲率的显式形式为： $$\kappa = \frac{\phi_{xx}\phi_y^2 - 2\phi_x\phi_y\phi_{xy} + \phi_{yy}\phi_x^2}{(\phi_x^2 + \phi_y^2)^{3/2}}$$

5.2.2 等势面演化方程

界面随流体运动的演化由以下偏微分方程描述： $$\frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla\phi = 0$$ 这是一个Hamilton-Jacobi型方程。在数值求解时，空间离散需要特别注意：

HJ-WENO格式：对于$\phi_x$的逼近，使用迎风方向： $$\phi_x \approx \begin{cases} D^{-x}\phi & \text{如果 } u > 0 \\ D^{+x}\phi & \text{如果 } u < 0 \end{cases}$$ 其中$D^{\pm x}$表示单边差分算子，具体使用WENO重构。

时间积分：通常采用TVD Runge-Kutta方法： $$\begin{aligned} \phi^{(1)} &= \phi^n - \Delta t L(\phi^n) \\ \phi^{(2)} &= \frac{3}{4}\phi^n + \frac{1}{4}\phi^{(1)} - \frac{1}{4}\Delta t L(\phi^{(1)}) \\ \phi^{n+1} &= \frac{1}{3}\phi^n + \frac{2}{3}\phi^{(2)} - \frac{2}{3}\Delta t L(\phi^{(2)}) \end{aligned}$$ 其中$L(\phi) = \mathbf{u} \cdot \nabla\phi$。

5.2.3 自由表面边界条件

在自由表面，需要满足动力学和运动学边界条件：

动力学条件（应力平衡）： $$\mathbf{n} \cdot (\boldsymbol{\sigma}_l - \boldsymbol{\sigma}_g) = \sigma\kappa\mathbf{n}$$ 其中$\boldsymbol{\sigma}_l$和$\boldsymbol{\sigma}_g$分别是液体和气体侧的应力张量，$\sigma$是表面张力系数。

对于不可压缩流体，这简化为： $$p_l - p_g = \sigma\kappa - 2\mu\mathbf{n} \cdot \mathbf{D} \cdot \mathbf{n}$$ 其中$\mathbf{D} = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$是应变率张量。

压力边界条件的离散化：在Staggered网格上，界面压力跳跃的处理需要特别注意。使用Ghost Fluid Method（GFM）：

识别界面单元：$\phi_i \cdot \phi_{i+1} < 0$
计算界面位置：$\theta = \frac{\phi_i}{\phi_i - \phi_{i+1}}$
修正压力梯度： $$\frac{\partial p}{\partial x}\bigg|_{i+1/2} = \frac{p_{i+1} - p_i}{\Delta x} + \frac{[p]}{\Delta x}$$ 其中$[p] = \sigma\kappa$是压力跳跃。

表面张力的计算：使用连续表面力（CSF）模型，将表面张力转化为体积力： $$\mathbf{F}_{st} = \sigma\kappa\delta(\phi)\nabla\phi$$ 其中$\delta(\phi)$是Dirac delta函数的光滑近似： $$\delta_\epsilon(\phi) = \begin{cases} \frac{1}{2\epsilon}\left(1 + \cos\left(\frac{\pi\phi}{\epsilon}\right)\right) & |\phi| < \epsilon \\ 0 & \text{otherwise} \end{cases}$$ 通常取$\epsilon = 1.5\Delta x$。

运动学条件（界面随流体运动）： $$\frac{D\phi}{Dt} = \frac{\partial\phi}{\partial t} + \mathbf{u} \cdot \nabla\phi = 0$$ 这正是等势面演化方程。

速度外推：为了在界面附近获得合理的速度场，需要将速度从流体区域外推到空气区域： $$\frac{\partial \mathbf{u}}{\partial \tau} + S(\phi)\mathbf{n} \cdot \nabla\mathbf{u} = 0$$ 其中$\tau$是虚拟时间，通常迭代5-10步即可。

5.2.4 重初始化技术

数值求解过程中，$\phi$会逐渐偏离有符号距离函数。重初始化通过求解以下方程恢复距离性质： $$\frac{\partial\phi}{\partial \tau} + S(\phi_0)(|\nabla\phi| - 1) = 0$$ 其中$\tau$是虚拟时间，$S(\phi_0)$是基于初始值的符号函数： $$S(\phi_0) = \frac{\phi_0}{\sqrt{\phi_0^2 + |\nabla\phi_0|^2\Delta x^2}}$$ HJ-WENO离散化：对于$|\nabla\phi|$的计算，使用Godunov's scheme： $$|\nabla\phi| = \begin{cases} \sqrt{\max(D^{-x}_+\phi, 0)^2 + \min(D^{+x}_+\phi, 0)^2 + \cdots} & S > 0 \\ \sqrt{\min(D^{-x}_-\phi, 0)^2 + \max(D^{+x}_-\phi, 0)^2 + \cdots} & S < 0 \end{cases}$$ 其中$D^{\pm x}$使用WENO重构获得高阶精度。

Sussmann重初始化：为保持界面位置不变，采用约束形式： $$\frac{\partial\phi}{\partial \tau} = S(\phi_0)(1 - |\nabla\phi|) + \lambda H'(\phi_0)$$ 其中$\lambda$通过约束$\int_\Omega H(\phi)d\Omega = \text{const}$确定，$H$是Heaviside函数。

拉格朗日乘子的计算： $$\lambda = -\frac{\int_\Omega S(\phi_0)(1 - |\nabla\phi|)H'(\phi_0)d\Omega}{\int_\Omega [H'(\phi_0)]^2 d\Omega}$$ 快速扫描法：对于稳态Eikonal方程$|\nabla\phi| = 1$，可使用快速扫描法：

固定界面处的$\phi = 0$
按四个（2D）或八个（3D）方向扫描网格
每个点更新：$\phi_{i,j} = \min(s_{i,j}, \phi_{i,j})$

其中$s_{i,j}$由相邻点通过迎风格式计算。具体地，对于2D情况： $$s_{i,j} = \min \begin{cases} \phi_{i-1,j} + \Delta x \\ \phi_{i+1,j} + \Delta x \\ \phi_{i,j-1} + \Delta y \\ \phi_{i,j+1} + \Delta y \\ \text{solve}(\phi_x, \phi_y) \end{cases}$$ 其中$\text{solve}(\phi_x, \phi_y)$求解二次方程： $$\max(D^{-x}\phi, -D^{+x}\phi, 0)^2 + \max(D^{-y}\phi, -D^{+y}\phi, 0)^2 = 1$$ 快速行进法（FMM）：更高效的选择，复杂度$O(N\log N)$：

初始化：将界面点标记为"已知"，邻居标记为"窄带"
循环：从窄带中选择最小值点，更新其邻居
使用最小堆维护窄带

局部重初始化：只在界面附近$|\phi| < \beta\Delta x$的窄带内重初始化，通常$\beta = 5-10$，大大减少计算量。

5.3 路径追踪与球面追踪

渲染是物理引擎不可分割的一部分。本节介绍两种现代渲染技术：路径追踪（Path Tracing）用于物理正确的全局光照，球面追踪（Sphere Tracing）用于隐式表面的高效渲染。

5.3.1 路径追踪基础

路径追踪基于渲染方程，通过蒙特卡洛方法求解光线传输：

渲染方程： $$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) |\cos\theta_i| d\omega_i$$ 其中：

$L_o$：出射辐射度
$L_e$：自发光
$f_r$：双向反射分布函数（BRDF）
$L_i$：入射辐射度
$\theta_i$：入射角

蒙特卡洛估计：使用重要性采样估计积分： $$L_o \approx L_e + \frac{1}{N} \sum_{i=1}^{N} \frac{f_r(\omega_i) L_i(\omega_i) |\cos\theta_i|}{p(\omega_i)}$$ 其中$p(\omega_i)$是采样概率密度函数。

俄罗斯轮盘赌：为避免无限递归，使用俄罗斯轮盘赌终止路径： $$L = \begin{cases} L_e + \frac{f_r L_i |\cos\theta|}{p(\omega)P_{rr}} & \text{以概率 } P_{rr} \\ L_e & \text{以概率 } 1 - P_{rr} \end{cases}$$ 通常选择$P_{rr} = \min(1, \rho_{max})$，其中$\rho_{max}$是最大反射率。

5.3.2 球面追踪算法

对于隐式表面$f(\mathbf{x}) = 0$（如有符号距离场），球面追踪提供了高效的光线求交方法。

基本算法：从点$\mathbf{p}$沿方向$\mathbf{d}$追踪：

t = 0
while t < t_max:
    p = origin + t * direction
    d = SDF(p)  // 有符号距离函数
    if d < epsilon:
        return t  // 找到交点
    t += d  // 安全步进距离
return no_intersection

数学保证：如果$f$是真正的距离函数（满足$|\nabla f| = 1$），则球面追踪保证不会穿过表面。

距离函数组合：

并集：$d_{A \cup B} = \min(d_A, d_B)$
交集：$d_{A \cap B} = \max(d_A, d_B)$
差集：$d_{A \setminus B} = \max(d_A, -d_B)$

软阴影：使用球面追踪可以高效计算软阴影： $$S = \min\left(1, k \cdot \min_{t \in [0,t_{max}]} \frac{d(t)}{t}\right)$$ 其中$k$控制阴影的软硬程度。

5.3.3 运动模糊实现

物理模拟中的运动模糊对于表现快速运动至关重要。

时间采样：对每条光线随机采样时间$t \in [0, 1]$： $$\mathbf{x}(t) = \mathbf{x}_0 + t(\mathbf{x}_1 - \mathbf{x}_0)$$ 对于旋转运动，使用四元数插值： $$q(t) = \text{slerp}(q_0, q_1, t)$$ 变形运动模糊：对于变形物体（如流体），需要插值整个几何： $$\mathbf{v}_i(t) = (1-t)\mathbf{v}_i^0 + t\mathbf{v}_i^1$$ 时空数据结构：为加速，构建4D（3D空间+1D时间）的加速结构：

4D BVH：每个节点存储空间-时间边界框
时间分片：将时间域离散化，每片使用独立的空间结构

采样策略：分层采样减少方差： $$t_i = \frac{i + \xi_i}{N}, \quad \xi_i \in [0, 1]$$ 对于周期性运动（如旋转），使用循环时间映射。

5.4 体素渲染与体积渲染

体素渲染和体积渲染是可视化三维标量场（如密度场、温度场）的核心技术。在物理引擎中，这些技术用于渲染烟雾、云层、医学数据等体积数据。

5.4.1 体积渲染方程

体积渲染基于辐射传输理论，考虑光线穿过参与介质时的吸收、发射和散射。

辐射传输方程：沿光线$s$的辐射度变化： $$\frac{dL(s)}{ds} = -\sigma_t(s)L(s) + \sigma_s(s)L_s(s) + \sigma_a(s)L_e(s)$$ 其中：

$\sigma_t = \sigma_a + \sigma_s$：消光系数（吸收+散射）
$\sigma_a$：吸收系数
$\sigma_s$：散射系数
$L_s$：由于散射进入的辐射度
$L_e$：自发光辐射度

体积渲染积分：忽略散射（或预计算），简化为： $$L(D) = \int_0^D T(s)\sigma_a(s)L_e(s)ds + T(D)L_{bg}$$ 其中透射率： $$T(s) = \exp\left(-\int_0^s \sigma_t(t)dt\right)$$ 离散化形式：将光线离散为$N$个采样点： $$L \approx \sum_{i=0}^{N-1} T_i \alpha_i c_i + T_N L_{bg}$$ 其中：

$\alpha_i = 1 - \exp(-\sigma_t(s_i)\Delta s)$：不透明度
$c_i = L_e(s_i)$：发射颜色
$T_i = \prod_{j=0}^{i-1}(1-\alpha_j)$：累积透射率

5.4.2 体素数据结构

规则体素网格：最简单的体素表示，使用三维数组存储：

voxel[i][j][k] = density/color at position (i,j,k)

内存需求：$O(n^3)$，其中$n$是每个维度的分辨率。

稀疏体素八叉树（SVO）：对于稀疏数据，使用八叉树减少内存：

每个节点代表一个立方体区域
空区域不分配子节点
叶节点存储实际数据

遍历复杂度：$O(\log n)$，内存效率取决于稀疏程度。

层次体素结构：结合多分辨率表示：

低分辨率用于远处或低频区域
高分辨率用于近处或高频区域
支持细节层次（LOD）切换

体素MipMap：预计算多分辨率版本，用于：

反走样
加速远距离渲染
空间-频率分析

下采样时保持物理量守恒： $$\rho_{l+1}(i,j,k) = \frac{1}{8}\sum_{di,dj,dk \in \{0,1\}} \rho_l(2i+di, 2j+dj, 2k+dk)$$

5.4.3 光线积分优化

自适应采样：根据密度梯度调整采样率： $$\Delta s = \Delta s_{base} \cdot \min\left(1, \frac{c}{|\nabla\rho| + \epsilon}\right)$$ 其中$c$控制自适应程度。

提前终止：当累积不透明度接近1时终止： $$\sum_{i=0}^{k} \alpha_i > T_{threshold} \Rightarrow \text{终止}$$ 通常$T_{threshold} = 0.95$或$0.99$。

空间跳跃：使用辅助数据结构跳过空区域：

边界体积层次（BVH）
距离场：存储到最近非空体素的距离
占用图：低分辨率二值网格标记非空区域

预积分传输函数：对于分段线性传输函数，预计算积分表： $$I(s_f, s_b) = \int_{s_f}^{s_b} \tau(\rho(s))ds$$ 查表时使用前后采样点的密度值$(s_f, s_b)$。

5.4.4 高级体积效果

体积光照：考虑光线在体积内的散射，使用相位函数$p(\theta)$描述散射方向分布。

Henyey-Greenstein相位函数： $$p_{HG}(\cos\theta) = \frac{1}{4\pi} \cdot \frac{1-g^2}{(1+g^2-2g\cos\theta)^{3/2}}$$ 其中$g \in [-1,1]$是各向异性参数：

$g > 0$：前向散射
$g < 0$：后向散射
$g = 0$：各向同性散射

多重散射：使用球谐函数（SH）近似散射辐射度： $$L_s(\omega) \approx \sum_{l=0}^{L} \sum_{m=-l}^{l} L_{lm} Y_l^m(\omega)$$ 通常$L=2$（9个系数）足够表示低频散射。

体积阴影：深度阴影图扩展到体积： $$V(s) = \exp\left(-\int_0^{d_{light}(s)} \sigma_t(s')ds'\right)$$ 其中$d_{light}(s)$是从点$s$到光源的距离。

使用级联阴影图或体积阴影图存储不同深度的透射率。

5.5 DDA算法与光线-粒子求交

数字微分分析器（DDA）算法是体素化空间中光线遍历的基础算法。本节探讨DDA的原理、优化及其在粒子系统渲染中的应用。

5.5.1 DDA算法原理

DDA算法通过增量计算高效遍历光线经过的所有体素。

基本思想：给定光线$\mathbf{r}(t) = \mathbf{o} + t\mathbf{d}$，计算光线与每个体素边界的交点，按顺序访问体素。

2D情况分析：考虑光线与垂直/水平网格线的交点：

与垂直线$x = i$的交点：$t_x = (i - o_x) / d_x$
与水平线$y = j$的交点：$t_y = (j - o_y) / d_y$

算法核心：维护下一个交点的$t$值，每步选择最小的$t$前进：

初始化：
tMaxX = (voxel.x + sign(d.x) - o.x) / d.x
tMaxY = (voxel.y + sign(d.y) - o.y) / d.y
tDeltaX = |1.0 / d.x|
tDeltaY = |1.0 / d.y|

遍历循环：
while (在边界内):
    访问当前体素
    if (tMaxX < tMaxY):
        tMaxX += tDeltaX
        voxel.x += sign(d.x)
    else:
        tMaxY += tDeltaY
        voxel.y += sign(d.y)

3D扩展：添加$z$维度，每步选择三个方向中$t$最小的： $$\text{next} = \arg\min(tMaxX, tMaxY, tMaxZ)$$ 数值稳定性：处理方向分量接近零的情况： $$tDelta_i = \begin{cases} |1.0 / d_i| & \text{如果 } |d_i| > \epsilon \\ \infty & \text{否则} \end{cases}$$

5.5.2 DDA优化技术

整数DDA：使用定点数算术避免浮点运算：

// 使用16位小数精度
scale = 65536  // 2^16
ix = int(x * scale)
idx = int(dx * scale)

每步更新：ix += idx

多层次DDA：在稀疏体素结构中，使用层次化DDA：

高层级DDA找到非空区域
低层级DDA遍历具体体素
空间跳跃优化

SIMD加速：并行处理多条光线：

4条光线打包到SIMD寄存器
统一的遍历逻辑
发散处理的掩码操作

光线相干性利用：对于相邻光线，利用增量更新： $$\mathbf{o}_{i+1} = \mathbf{o}_i + \Delta\mathbf{o}$$ $$\mathbf{d}_{i+1} = \mathbf{d}_i + \Delta\mathbf{d}$$ 减少初始化开销。

5.5.3 光线-粒子求交

物理模拟中的粒子系统（SPH、MPM等）需要高效的渲染方法。

球体求交：粒子通常表示为球体，光线-球体求交的二次方程： $$|\mathbf{o} + t\mathbf{d} - \mathbf{c}|^2 = r^2$$ 展开得： $$t^2 + 2t(\mathbf{d} \cdot \mathbf{v}) + |\mathbf{v}|^2 - r^2 = 0$$ 其中$\mathbf{v} = \mathbf{o} - \mathbf{c}$。

判别式：$\Delta = (\mathbf{d} \cdot \mathbf{v})^2 - (|\mathbf{v}|^2 - r^2)$

屏幕空间优化：对于小粒子，使用屏幕空间点精灵：

投影粒子中心到屏幕
计算屏幕空间半径
光栅化为圆形精灵

屏幕空间半径： $$r_{screen} = \frac{f \cdot r_{world}}{z_{eye}}$$ 其中$f$是焦距，$z_{eye}$是眼空间深度。

深度剥离：处理半透明粒子的正确混合：

多遍渲染，每遍剥离最前层
从后向前合成
或使用OIT（顺序无关透明）技术

粒子光照：考虑粒子的次表面散射： $$L_o = L_e + \int_\Omega f_s(\omega_i, \omega_o) L_i \cos\theta_i d\omega_i$$ 对于半透明粒子，使用简化的BSSRDF模型。

5.5.4 大规模粒子渲染

空间哈希加速：使用空间哈希表组织粒子： $$hash(i,j,k) = (i \cdot p_1 \oplus j \cdot p_2 \oplus k \cdot p_3) \mod m$$

其中$p_1, p_2, p_3$是大质数。

层次化表示：远距离粒子聚合为单个代表：

使用八叉树或KD树
远处使用低分辨率表示
支持连续LOD过渡

GPU粒子系统：利用GPU并行性：

粒子数据存储在缓冲区
顶点着色器处理位置/大小
几何着色器生成广告板
片段着色器计算光照

实例化渲染：减少绘制调用：

// 顶点着色器
vec3 position = particlePosition[gl_InstanceID];
float radius = particleRadius[gl_InstanceID];
gl_Position = MVP * vec4(position + vertex * radius, 1.0);

时间相干性：利用帧间相干性：

增量更新空间结构
视锥体剔除缓存
遮挡查询结果重用