声明：这篇文章主要是关于SIG23年这篇黑狗毛论文的个人笔记。应该也没啥阅读门槛，因为我自己也没有入门图形学。公式有但是不多。公式tag跟原论文一样。所有内容都是我瞎捣鼓的，公式如果有理解错误请多多指正，感谢感谢！

原文：https://zhuanlan.zhihu.com/p/809636731

关键词：图形学入门、离线渲染、基于波动光学渲染、毛发渲染

Mengqi Xia, Bruce Walter, Christophe Hery, Olivier Maury, Eric Michielssen, and Steve Marschner, “A Practical Wave Optics Reflection Model for Hair and Fur,” ACM Transactions on Graphics (TOG), vol. 42, no. 4, article 39, pp. 1-15, Jul. 2023.

1. Related Work

• 基于Ray的纤维模型 在人类毛发的研究中，Marschner[2003]的模型广泛应用于业界，通过分析介电圆柱和圆锥中的光线路径，将散射分为R、TT和TRT。Zinke[2009]加入了漫反射成分。Sadeghi[2010]提出了艺术家控制参数化方法。d’Eon[2014]和Huang[2022]提出了非可分离的表征方法，也就是说方位角和纵向角度之间有耦合效应，不可简单分离。Chiang[2016] 进一步优化了该模型，使其适用于生产级渲染。 除了人类毛发，Khungurn和Marschner[2017]探讨了椭圆形毛发的建模。Yan[2015, 2017]研究了带有内部髓质的动物毛发。Aliaga[2017]将模型推广至有更加复杂截面的纺织纤维。
• 基于波动光学的纤维模型 Linder[2014]解析了具有完全圆形截面的圆柱纤维，并且探讨了散射行为。Xia[2020]通过在具有任意截面的圆柱上进行二维波动模拟，展示了与几何光学模型相比的若干重要差异。然而这个模型是规则的圆形，并没有譬如毛鳞片等微观几何结构。Benamira和Pattanaik[2021]提出了更加快速的混合模型，这个模型仅仅在前向散射使用波动光学求解，其余部份则是依赖几何光学。
• 基于物理光学的平面模型模拟 物理光学平面模型使用物理光学近似来模拟光线在接近平面、可以表示为高度场的粗糙表面上的散射行为，包括高斯随机、周期性、预计算和划痕表面。它们结合了Kirchhoff标量衍射理论和路径追踪方法来处理散射与反射，并通过各种模型（如Beckmann-Kirchhoff和Harvey-Shack）计算光在不同粗糙表面上的衍射。尽管这些模型在平面表面上有效，但纤维几何的封闭性和复杂的光线交互需要更复杂的处理方法。
• 计算电磁学、Spickle Effect 与 程式化噪声 上一篇文章已经提及了，这里略过。

https://zhuanlan.zhihu.com/p/776529221

2. Background Research

概述

毛发建模是基于毛发的扫描电子显微镜（SEM）图像。可以准确还原毛发的微观结构，例如毛鳞片、粗糙度。

接下来，用「WAVE SIMULATION WITH 3D FIBER MICROGEOMETRY」计算粗糙纤维表面的反射和衍射。

在这个基础上，引入了散斑理论来分析散射模式的统计特性。并且用噪声来描述这些散斑，大幅优化模型。

然后通过与实际测量数据的对比，推导出合理的纤维参数（例如尺寸、表皮角度和表面粗糙度）。最后集成进渲染系统中。

3. OVERVIEW

3.1 Fiber scattering models

这个我翻译成纤维散射模型。这个东西描述的是单根纤维发生相互作用。这里的关键核心是BCSDF（双向曲线散射分布函数）。有别于在表面反射和折射中常用的双向散射分布函数（BSDF），BCSDF是专门为曲线形状的纤维设计的。下面这个公式讲的是，给定某波长的光照射在一根纤维上时，它从某个方向射入后，会从另一个方向反射或透射出来。
$$
L_r(\omega_r, \lambda) = \int L_i(\omega_i, \lambda) S(\omega_i, \omega_r, \lambda) \cos \theta_i d\omega_i \tag{1}
$$
公式左边表示给定波长 $\lambda$ ，出射方向 $\omega_r$ 下的辐射亮度。$ L_i(\omega_i, \lambda)$ 是入射方向 $\omega_i$ 的辐射亮度。 $S(\omega_i, \omega_r, \lambda)$ 是双向曲线散射分布函数，描述了光线如何被纤维”打散”。$\cos \theta_i$ 是为了考虑入射角度的影响，如果光线以一个很平的角度照射到纤维上，它的影响会比垂直照射时小。

把上面这个公式写成球坐标，并且把每一种不同的光线与毛发纤维的交互当作不同的模式，然后把不同的散射项累加起来：
$$
S(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) = \sum_{p=0}^{\infty} S_p(\theta_i, \theta_r, \phi_i, \phi_r, \lambda)
\tag{2}
$$
第一散射项 $S_0$ 描述表面反射，即以前经常说的直接反射项 $R$ 。这一项是一般代表从光滑纤维的反射或者粗糙纤维表面反射的统计平均值。论文这里能更加精确地计算这一项。

回想以前的渲染方法（比如Marschner[2003]），通常是将每一个散射模式 $S_p$ 分解为两个独立的函数：longitudinal function $M_p$ 和anazimuthal function $N_p$ ，这样的方法已经被批判过是不准确的，因此应该避免使用下面这种可分离的近似方法。
$$
S_p(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) = M_p(\theta_i, \theta_r)N_p(\theta_i, \phi_i, \phi_r, \lambda)
\tag{3}
$$
因此，XIA[2023]采样了多个粗糙纤维的散射参数并且取均值，记为 $S_0,avg$ 。 $f(\theta_h, \phi_h, \lambda)$ 是噪声分量，通过两个半程向量角度和波长来表示，用来修正当前的特定纤维实例和均值的偏差。
$$
S_{0,\text{sim}}(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) \approx S_{0,\text{avg}}(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) f(\theta_h, \phi_h, \lambda)
\tag{4}
$$
当前，完整的纤维散射模型如下。第一项表示纤维表面的反射的散射模式，结合了3D波动模拟。后面的求和项是其他高阶的散射模式之和。
$$
S_{\text{prac}}(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) = S_{0,\text{prac}}(\theta_i, \theta_r, \phi_i, \phi_r, \lambda) + \sum_{p=1}^{\infty} S_p(\theta_i, \theta_r, \phi_i, \phi_r, \lambda)
\tag{5}
$$
实际上，最终的散射公式更为简洁，以适应有更复杂几何结构的表面。

3.2 Speckle theory

散斑理论（Speckle Theory）描述光与粗糙表面相互作用时产生的随机光强分布现象。Goodman[2007]公式中的 $A$ 表示所有相位向量的叠加，即结果相位向量（Phasor）：
$$
\mathbf{A} = \frac{1}{\sqrt{N}} \sum_{n=1}^{N} a_n = \frac{1}{\sqrt{N}} \sum_{n=1}^{N} a_n e^{i\phi_n}
$$
上式为散斑强度的综合表达，换而言之，用来表示散射后的光线强度。产生该现象的原因是光线会在许多微小表面之间相互反射、折射并相互干涉。光线传播发生的相位差导致明暗不均的图案。通过该公式，可以统计性地计算出光线如何在纤维表面相互干涉，得到散斑图案。

4. WAVE SIMULATION WITH 3D FIBER MICROGEOMETRY

在波动光学模拟中，光看作一种电磁波。由相互垂直的磁场和电场组成。光线与毛发的相互作用（如散射）可以转化为分析电磁场如何受到物体的影响。

4.1 Wave optics

首先需要明白，当电磁场随着时间周期变化时，称该场为时偕场（Time-Harmonic Fields）。在时谐场中，电场和磁场可以通过复数表示为相位向量（phasors）。巧妙的是，电场和磁场本身就是相互垂直。
$$
E_{\text{inst}} = \Re(E e^{j\omega t}), \quad H_{\text{inst}} = \Re(H e^{j\omega t})
\tag{6}
$$
分别是电场和磁场，只不过是分别取了实数部分。复数部分包含了场的振幅和相位信息。

下面这两组方程是麦克斯韦方程组在时谐场中的形式：
$$
\nabla \times \mathbf{E} = -\mathbf{M} – j\omega \mu \mathbf{H}
\
\nabla \times \mathbf{H} = \mathbf{J} + j\omega \varepsilon \mathbf{E}
\tag{7}
$$
其中， $\varepsilon$ 是介电常数（影响电场）， $\mu$ 是磁导率（影响磁场）。这两项描述了材料怎么影响电磁场的传播。

当物体（如纤维）被一束入射波照射时，入射电场和磁场分别用 $\mathbf{E}_i$ 和 $\mathbf{H}_i$ 表示。但是物体的存在改变了这些场，使得我们观察到的是总场。
$$
\mathbf{E}_1 = \mathbf{E}_i + \mathbf{E}_s, \quad \mathbf{H}_1 = \mathbf{H}_i + \mathbf{H}_s
\tag{8}
$$
简而言之，总场 = 入射场 + 散射场。

光能量随着散射场 $\mathbf{E}_s$ 和 $\mathbf{H}_s$ 向外传播，因此计算纤维的散射函数是关键。怎么算呢？

全波模拟是最为准确的。进行全波模拟需要将物体离散化为网格，并且需要高分辨率（一般来说，每个波长至少10个网格单元），这导致需要处理百万级别的网格单元。在上一篇文章也有提及。

也就是说，即便是模拟一个仅仅几十微米长的短纤维段，也需要处理数百万个网格单元。

因此，采用 物理光学近似（PO）。在PO中，物体表面的电流和磁流（分别记作 $J$ 和 $M$ ）可以看作是散射场的次级源头。电磁流产生了二次辐射，形成散射波。PO假设物体表面只会发生单次反射，忽略多次反射和复杂的衍射效应。得到表面的电流和磁流之后，计算他们产生的散射波。通过从表面电流和磁流中推导出这些远场波的性质。

光一个PO还不够，还得揉个八叉树算法进去加速远场计算。

4.2 Physical Optics Approximation

具体地，PO做了两个简化假设：单次散射与局部平面假设。这个方法也足够通用。

如上图，某物体表面采样了多个点，近似为平面。计算切平面电流和磁场，进而产生一个散射场。通过这个散射场，计算这些电流和磁流产生的远场散射波。与此同时，使用Octree结构划分为更小的体素。每个叶节点的计算结果会聚合至父节点，得到总的辐射贡献。

Surface current calculation.

表面电流、磁流的计算是PO模拟中非常关键的一步。对于每一个采样点，存储其法向量 $n(r{\prime})$ 和面积信息。接下来再每个小平面计算入射波和表面电流的相互作用。

根据入射波的方向 $e_i$ 和表面法向 $n(r{\prime})$ ，将入射电场分解为平行极化和垂直极化分量。分解的原因是，平行极化和垂直极化在菲涅尔方程中，是完全不同的表达。

平行极化分量：光的电场平行于入射光线的反射面。
垂直极化分量：光的电场垂直于入射光线的反射面。

因此，得到入射场 $E_i$ 被分解为平行极化的分量 $E_i^p$ 和垂直极化的分量 $E_i^s$ ，长这个样子： $E_i = E_i^p + E_i^s$ 。

反射场 $E_r$ 被表示为平行极化和垂直极化的分量之和，入射场旁边的系数表示菲涅尔方程中的反射：
$$
E_r = E_r^p + E_r^s = F^p E_i^p + F^s E_i^s
\tag{9-1}
$$
然后，总电场 $E_1$ 是入射场和反射场的总和：
$$
E_1 = E_i + E_r
\tag{9-2}
$$
根据中学物理，我们知道电会生磁，磁也会生电。因此有如下表达：
$$
M = -n \times E_1,
J = n \times H_1
\tag{10}
$$
因此，根据入射场和反射场就可以计算出表面的感应电流和磁流。得到电流和磁流后，计算远场的散射波。

在理论上，入射场可以是任意的。但是这里作者用了Gaussian-windowed平面波。这种波的振幅遵循正态分布，好算。

总结一下，这里将入射光分解为两个分量，用菲涅尔方程算反射场。进而将反射场与入射场相加，得到总场。通过电/磁场与电/磁流的关系，计算表面的感应电/磁流。这样就可以模拟纤维的散射了。

也就是说，计算纤维表面的电流和磁流的确可以得到光的散射行为。

Far-field radiation in 3D

上一节通过表面电流得到表面电流和磁流。本节用这个信息来计算远场散射。

基于惠更斯原理（Huygens’s Principle）将原本的散射问题转化为辐射问题。

惠更斯原理指出，每个波前的电都可以看作新的次级波源。一环扣一环的感觉。

毛发纤维表面的电流 $J$ 和 磁流 $M$ 被视为光线的次级光源，然后重新辐射出来，得到散射场。

这里就有点难了，涉及到电磁学中的矩量法（Method of Moments, MoM）。读者可以深入研究一下Gibson[2021]这本书，学会了一定要好好教我。但是没关系，我从闫老师2021年的paper中找了一张图，保证你能看懂。

图中的红色小球表示的次级辐射源，各自向外辐射，类似于表面电流和磁流产生的次级辐射波。
次级辐射源向各个方向发射等强度的波，这与本文散射公式中每个表面点对所有方向都有散射贡献的思想一致。
图中观察了距离光源远场区域（距离为 $r$ ）的一小片区域（ $\delta \mathbf{r}$ ）。在远场区域中，分析波在不同位置 $\mathbf{r}_1$ 和 $\mathbf{r}_2$ 的行为，分别沿着 $\hat{r}_1$ 和 $\hat{r}_2$ 两个方向观察。这与散射公式中，在远场处计算散射电场的行为非常相似。
图中右侧的光波会在远场叠加形成复杂的干涉条纹，这类似于散射公式中的积分项，将次级辐射源的电磁波在远处叠加。

用公式则表示为：
$$
E_s(\mathbf{r}) = j \omega \mu_0 \frac{e^{-jk_0 R}}{4\pi R} \hat{r} \times \int_\Gamma \left[ \hat{r} \times \mathbf{J}(r{\prime}) + \frac{1}{Z_0} \mathbf{M}(r{\prime}) \right] e^{jk_0 r{\prime} \cdot \hat{r}} d\mathbf{r{\prime}}
\tag{11}
$$
结合图来理解，公式描述的是散射电场 $E_s(\mathbf{r})$ 在远处（远场）某一点 $\mathbf{r}$ 处的表现，与距离呈 $1/R$ 关系衰减。这个公式是一个二维复数平面，也就是把麦克斯韦方程组的解看作时谐电磁波的形式。

看这个上式结构，形式上跟 $\mathbf{E}(\mathbf{r}, t) = \mathbf{E}_0 e^{j(\omega t – k \cdot \mathbf{r})}$ 很像，叫电场的常见时谐解。

对应起来。 $ j \omega \mu_0$ 这虚数项描述磁场对电场的影响，与电磁场的频率 $\omega$ 和真空中的磁导率 $\mu_0$ 相关。光波的频率越高，电场越强。 $e^{-jk_0 R}$ 是相位因子，表示光波在传播中的相位变化。波数 $k_0 = \frac{2\pi}{\lambda}$ ， $\lambda$ 是电磁波的波长。表示波传播到距离 R 的地方时，电场的相位会发生变化。 $\hat{r}$ 是从散射物体指向观察点的单位向量，表示波的传播方向。叉乘 $\times$ 操作符表示计算的是电场的方向。确保计算出的电场和波传播的方向一致。

积分项 $ \int_\Gamma \left[ \hat{r} \times \mathbf{J}(r{\prime}) + \frac{1}{Z_0} \mathbf{M}(r{\prime}) \right] e^{jk_0 r{\prime} \cdot \hat{r}} d\mathbf{r{\prime}} $ 是该公式的核心。对物体表面 $\Gamma$ 进行积分，得到毛发表面每个点对散射电场的贡献。每个点的表面电流加表面磁流之和，再点乘相位因子。进一步的， $\mathbf{J}(r{\prime})$ 是表面电流密度， $\hat{r} \times$ 确保了电流产生的散射电场与波的传播方向 $\hat{r}$ 正交。 $Z_0$ 是自由空间阻抗（free space impedance），具体数值 $Z_0 = \sqrt{\frac{\mu_0}{\varepsilon_0}} \approx 377 \, \Omega$ ，个人理解该常数是真空中电场和磁场的比例关系，电场和磁场的抗阻不同因此需要统一尺度才能线性相加，即将磁流规范化到和电流类似的形式。指数项也是相位因子，不做过多解释。

细心的读者可能发现，为何表面电流密度 $\mathbf{J}(r{\prime})$ 有叉乘 $\hat{r} \times \mathbf{J}(r{\prime})$ ，而表面磁流密度 $\mathbf{M}(r{\prime})$ 没有叉乘。

根据麦克斯韦方程组，电场和磁场是互相正交的。电流 $\mathbf{J}$ 会产生磁场，变化的磁场再反过来产生电场。电场和磁场是相互耦合的。但是，我们重新看看麦克斯韦方程组的下面两条。
$$
\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}
\
\nabla \times \mathbf{B} = \mu_0 \mathbf{J} + \mu_0 \varepsilon_0 \frac{\partial \mathbf{E}}{\partial t}
$$
磁场的生成是直接通过位移电场来实现的，方向已经和电场正交。而电场的产生是磁场的变化率，方向性需要调整。

换一个角度来理解，积分项左边描述的是电磁场强度随距离衰减，也考虑了相位变化和方向性。积分项描述每个表面点的电磁流与相位贡献。

总结，公式(11)是一个精确的表达，描述物体表面电磁场如何在距离 $r$ 处散射。得到了电场的公式，磁场就很容易计算了。

当散射距离 $R$ 足够大时，远场的简化表达式如下，可以直接忽略近场效应。换句话说，在远场中，传播特性已经趋于稳定，因此可以把积分项简化为与散射方向 $\hat{r}$ 相关的项，即 $E_s^{\text{far}}(\hat{r})$ 。
$$
E_s(r) = \frac{e^{-jk_0 R}}{R} E_s^{\text{far}}(\hat{r}), \quad H_s(r) = \frac{e^{-jk_0 R}}{R} H_s^{\text{far}}(\hat{r})
\tag{12}
$$
此外，若直接计算每个点的贡献，时间复杂度为 $O(MN)$ ，即离散点数$M$ $*$ 散射方向数$N$ 。此处引入八叉树，通过空间划分减少计算表面的离散点数，时间复杂度降到了 $O(M+log(M)N)$ 。具体实现方式在4.3。

4.3 Multilevel fast Physical Optics

用多层快速物理光学（Multilevel Fast Physical Optics, MFPO）加速远场散射计算。

Chew[2001]提出的多层快速多极子算法（MLFMA），用于加速求解电磁场散射问题。这个算法首先为例如毛发表面构建一个八叉树结构，每个叶节点表示一个采样点。构建完Octree后计算每个叶节点的表面电流、磁流。接着从叶节点出发逐层向上累积。从原本的 $O(MN)$ 降低到 $O(M + \log(M)N)$ 。

接着作者介绍了八叉树加速算法的三个关键部分。首先是远场散射核（far-field scattering kernel）。
$$
e^{jk_0 r{\prime} \cdot \hat{r}} = e^{jk_0 (r{\prime} – c_L) \cdot \hat{r}} \prod_{i=1}^{L} e^{jk_0 (c_i – c_{i-1}) \cdot \hat{r}}
\tag{13}
$$
最终目标是计算物体表面每个点的电磁波对某个远场观察区域（距离为 $R$ ，方向为 $\hat{r}$ ）的散射贡献。具体是求散射电场 $E_s(r)$ 和磁场 $H_s(r)$ ，即从毛发表面上的每个点 $r{\prime}$ 发出的电磁波在远场中的贡献。公式等号左边是某表面点 $r{\prime}$ 到远场观察方向 $\hat{r}$ 的相位变化。最终的时间复杂度为 $O(MN)$ 。作者把表面分为不同的区域，每个区域指定一个参考中心点，公式中的 $c_0, c_1, …, c_L$ 是八叉树中不同层次的节点中心。

光是这样也没办法减少计算量。因此需要将距离较近的表面点，由于它们的相位变化差别很小，通过八叉树的高层节点合并这些点的贡献。

还是用闫老师这幅图来解释。对于远场的某个区域 $\delta \mathbf{r} $ ，首先还是会精确计算所有采样点对参考点 $\mathbf{r}$ 的贡献。但是该区域的其他点则用八叉树的父节点进行近似计算。也就是说，旁边的 $\mathbf{r}_1$ 和 $\mathbf{r}_2$ 就不再需要考虑这么多的采样点了，而是直接用八叉树得到的总贡献。

作者定义了一个方向集，八叉树每个父节点都存储了关于不同方向集的累积贡献数据。因此，父节点不仅包含空间上的信息，还包含多个散射方向上的累积信息。最终在根节点得到完整的360度散射场分布。

接下来，从树的第二层开始，依次向上合并。上采样具体方法是对子节点的方向集中的散射贡献做正向FFT，接着zero-padding扩充频域数据最后逆向FFT转换为空间域。最终，使得父节点可以在不同方向上得到更加精确的散射信息。

• Performance

八叉树加速效果很好。三层树效果最佳。越复杂的纤维优化效果越好。

• Fiber microgeometry and scattering patterns

毛发纤维横截面一般不是完美的圆形，而是椭圆形。通过椭圆的主半径 $r_1$ 和次半径 $r_2$ 来定义纤维的几何参数。为了模拟纤维表面的微观粗糙度，作者在椭圆柱体的表面上叠加了一个高斯随机高度场（Gaussian random height field），模拟真实的纤维表面。进一步地，加入了角质层倾斜（cuticle tilt），模拟薄片在纤维表面倾斜排列。

通过对比传统的基于光线的毛发模型，波动光学模拟发现，除了XIA[2020]预测到的显著的向前散射（forward-scattering）现象，还观察到了复杂的波长依赖颗粒图案（wavelength-dependent granular patterns）。在转换为RGB颜色时，会生成丰富的色彩效果。

研究指出了一些从模拟中观察到的规律：

• 无论它们的实际位置如何，具有相同几何参数（如半径、粗糙度、角质层倾斜等）的纤维，在散射行为上会生成相似的颗粒图案。
• 如果纤维的几何参数不同，则它们会生成不同统计特性的颗粒图案，即散射模式明显不同。
• 散射斑点的位置依赖于光线的入射角度，偏移的方向跟随half vector方向。
• 斑点（speckle）图案的大小随光波波长的增加而增大。这一现象与Goodman[2007]结果一致。

5. A PRACTICAL FIBER SCATTERING MODEL

实用纤维散射模型（A PRACTICAL FIBER SCATTERING MODEL）旨在解决纤维的微观几何变化和复杂散射行为。

以前的研究一般是用Lut存散射分布函数，这种方法空间消耗巨大，需要同时记录纵向和方位角散射分布。

现在，作者提出了一种基于小波噪声表示的紧凑纤维散射模型，可表示的几何复杂度提高了，散射效果更好了。

简单的说，作者想要紧凑的统计散斑现象，使其可以通过均值、方差、自相关函数（ACF）等统计量来表示。因此作者借助了散斑理论（Speckle Theory）来描述光随机干涉所产生的图样。

5.1 Speckle statistics

这里作者上来就提到完全发展的散斑 (Fully Developed Speckle)这个概念。以下是个人理解。在光线照射到一个粗糙的表面（比如纤维/毛发表面）时，表面上的每个小区域都会散射光线。由于表面的微小不规则性，散射的光线之间会相互干涉，产生一个复杂的光强分布。这种分布表现为一系列亮点和暗点，我们称之为散斑。表面的微小特征（比如粗糙度）在整个照射区域（相对于光线的波长来说）变得足够不规则，导致散射光线在各个点上的相位和强度是随机的，这就是所谓的完全发展的散斑。

这个时候，可以用Goodman[2007]的complex Gaussian distribution来描述Fully Developed Speckle。即，电磁场的实部 $\mathcal{R}$ 和虚部 $\mathcal{I}$ 在空间上服从复高斯分布（complex Gaussian distribution）。
$$
p_{\mathcal{R},\mathcal{I}}(\mathcal{R}, \mathcal{I}) = \frac{1}{2 \pi \sigma^2} \exp\left( – \frac{\mathcal{R}^2 + \mathcal{I}^2}{2\sigma^2} \right)
\tag{14}
$$
场的实部和虚部是独立且正态分布的。均值为零，方差相同。

电磁场光强 $I$ 和该光强分布的服从指数分布（exponential distribution）的概率密度函数：
$$
I = \mathcal{R}^2 + \mathcal{I}^2 \ p_I(I) = \frac{1}{2\sigma^2} \exp\left( -\frac{I}{2\sigma^2} \right)
\tag{15}
$$
这些公式没啥可讲的，总之散斑场就是很随机哒！

作者让光强遵循指数分布，确保在单个方向的统计特性。其次，通过研究两个点之间的光强统计关系，衡量两个点的ensemble average。这里用autocorrelation function (ACF)。
$$
C(I_{p_1}, I_{p_2}) = \frac{\overline{(I_{p_1} – \overline{I_{p_1}})(I_{p_2} – \overline{I_{p_2}})}}{\sigma(I_{p_1}) \sigma(I_{p_2})}
\tag{16}
$$
自相关函数的值介于 -1 和 1 之间，当值接近 1 时，说明两个光强的散射行为非常相近。这是高效再现纤维散射场中颗粒结构的关键。

5.2 Wavelet noise representation of the speckles

作者引入小波噪声（Wavelet noise）来表示散斑的噪声分量 $f(\theta_h, \phi_h, \lambda)$ 。具体公式如下：
$$
f(\mathbf{x}) = \sum_{b=0}^{n-1} w_b(\mathbf{x}) I\left(2^b g_{\lambda}(\mathbf{x})\right)
\tag{17}
$$
公式的核心思想是将散斑的光强分解为不同频率层次的噪声，并且加权组合。

通过调整不同频率带的权重，生成自相关函数 $ C_f(\mathbf{x}_1, \mathbf{x}_2)$ 接近与目标自相关函数的最终噪声。

根据维纳-辛钦定理（Wiener-Khinchin theorem），自相关函数可以通过傅里叶变换（Fourier Transform）来计算。具体公式如下：

$$
C_f(\mathbf{x}_1, \mathbf{x}2) = \mathcal{F} \left( \mathcal{F}^{-1} \left( \sum{b=0}^{n-1} w_b I_b \right)^2 \right)
\tag{18}
$$
这表示我们可以通过计算小波噪声的傅里叶变换，来获得自相关函数。自相关函数和功率谱密度函数是一对傅里叶变换，太Amazing了。

这里原论文给了一个通过频率带加权求和近似计算ACF的证明。

省流：通过对噪声的各个频率带进行加权，可以近似得到整个噪声的自相关函数，而不用处理不同频率带之间的交互。

通过最小二乘法找到非负权重 $v_b$ ，以使得每个频率带的自相关函数 $C_b(\mathbf{x}1, \mathbf{x}_2)$ 的加权和可以逼近目标自相关函数 $C_t(\mathbf{x}_1, \mathbf{x}_2)$ 。

$$ C_t(\mathbf{x}_1, \mathbf{x}2) \approx \sum{b=0}^{n-1} v_b C_b(\mathbf{x}_1, \mathbf{x}_2) \tag{19} $$ 计算完权重后，要确保能量守恒。也就是调整噪声函数 $f(\mathbf{x})$ 的期望值 $\mathbb{E}[f(\mathbf{x})] = 1$ 。 $$ \begin{aligned} & \mathrm{E}\left[S{\text {avg }}\left(\theta_i, \theta_r, \phi_r, \phi_r, \lambda\right) f\left(\theta_h, \phi_h, \lambda\right)\right] \\\
& =S_{\text {avg }}\left(\theta_i, \theta_r, \phi_r, \phi_r, \lambda\right) \mathrm{E}[f(\mathbf{x})] \
& \approx S_{\text {avg }}\left(\theta_i, \theta_r, \phi_r, \phi_r, \lambda\right)
\end{aligned}
\tag{20}
$$

虽然效果好，但是有局限性。比如在grazing incidence angles的地方（光线几乎平行表面的情况），拟合准确性下降，导致散射准确度下降。Degeneration in the Forward Direction问题，当光线处于正前方方向（光线方向与表面法线一致）时，半向量方向会退化。

6. Validation

6.1 Wave simulation validation

计算了x-y平面中的散射强度，并通过3600个方位角 $\phi_r$ 进行计算，最终将这些角度平均到360个方向上。

首先对比Mie散射。在grazing angles不如Mie散射。物体的半径相对于波长较大时，PO近似更为精确。当物体曲率较小，就不准。

然后对比BEM。模拟一个小椭球体的波散射。BEM用了三个小时，PO花了两秒。

最后对比2D BEM。用一维高斯高度场（1D Gaussian height fields）包裹在圆形和椭圆形的横截面上。PO完胜。

6.2 Measurement

使用了波长为 633nm 的 HeNe 激光器，光束的光斑大小为 0.7mm（沿着头发的长度方向）× 3mm（垂直于头发方向）。激光通过一个小孔照射在人类头发样本的一个小区域。

总之就是效果好！

6.3 Noise representation validation

一个字，好！

7. Rendering

整合进PBRT-v3。原本的基于光线的模型记作 $S_{\text{ray}}$ ，而衍射模型记作 $S_{\text{diffract}}$ 。

对于衍射，这里近似为单缝衍射，而缝的宽度等于圆柱体（纤维）的直径。
$$
f_{\text{diffract}}(\theta_i, \phi_d, a) = a \cos \theta_i \cdot \text{sinc}^2(a \cos \theta_i \sin \phi_d)
\tag{22}
$$
这种衍射模型被用来和纵向函数结合，得到一个完整的双向曲面散射分布函数（BCSDF）。将衍射因子以 $50 \times 50 \times 200$ 的表格形式进行预计算。并使用同样大小的表格进行重要性采样。

消光截面可以简单理解为光与物体相互作用的“有效面积”。 纤维的消光截面会比它的实际几何截面大。消光截面会接近于几何截面的两倍。光在纤维上既发生反射也发生衍射，因此我们需要把总能量分配给这两种现象。根据经验，可以把一半的能量用于衍射，另一半用于反射。
$$
S_{\text{diffract}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) = \frac{1}{2} \left[S_{\text{ray}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) + f_{\text{diffract}}(\theta_i, \phi_d, D/\lambda)\right]
\tag{23}
$$
通过重要性采样，公平地考虑光的反射和衍射。

接下来，最终的描述光散射现象的渲染公式，形式整得挺好：
$$
S_{0,\text{prac}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) = S_{0,\text{avg}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) f(\theta_h, \phi_h, \lambda)\ \\
= \frac{1}{2} \left[ S_{\text{ray}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) + f_{\text{diffract}}(\theta_i, \phi_d, D/\lambda) \right] f(\theta_h, \phi_h, \lambda)\
\tag{24}
$$
反射+衍射+噪声。噪声函数 $f(\theta_h, \phi_h, \lambda)$ 在最后，使得散射的光线不在集中在少数方向。

重点看这个部分：
$$
S_{0,\text{prac}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) = S_{0,\text{avg}}(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) f(\theta_h, \phi_h, \lambda)
$$
$S_{0,\text{avg}}$ 表示毛发表面平均反射/折射和衍射的行为，这是通过预处理（BSDF表格）计算得到的。

接下来讲讲如何从PO提取BCSDF。

通过麦克斯韦方程得到散射的电场和磁场（$E_s^{\text{far}}$ 和 $H_s^{\text{far}}$）。

接着用坡印廷矢量（Poynting vector）计算能量流。相当于计算光的强度。

$$
\langle S \rangle = \frac{1}{2} \text{Re}(E \times H^) \tag{25}
$$

为了将模拟结果用于渲染，将散射强度与入射功率关联起来。散射强度的计算公式，其中 $R^2$ 是远场球面面积：

$$
I_s(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) = \langle S(\mathbf{r}) \rangle \cdot \hat{n} R^2 \tag{26}
$$

散射功率 $P_s$ 和吸收功率 $P_a$ 分别通过积分散射光和吸收光的强度来计算。

$$
P_s = \int_{\Omega} I_s(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) \, d\omega \tag{27}
$$

对于吸收功率 $P_a$ ，通过表面上的坡印廷矢量积分计算。

$$
\begin{aligned}
P_a & =\int_{\Gamma} \frac{1}{2} \operatorname{Re}\left(\mathbf{E}1 \times \mathbf{H}_1^{}\right) \cdot \hat{\mathbf{n}}_1(A) d A \ & = \\
\int{\Gamma} \frac{1}{2} \operatorname{Re}\left(\mathbf{J}^{} \times \mathbf{M}\right) \cdot \hat{\mathbf{n}}_1(s) d s
\end{aligned}
\tag{28}
$$

最终，BCSDF搞出来了。

$$
S(\theta_i, \phi_i, \theta_r, \phi_r, \lambda) = \frac{I_s(\theta_i, \phi_i, \theta_r, \phi_r, \lambda)}{|P_a – P_s|}
\tag{30}
$$

最后用了一个5D的表格。

1D 维度描述波长。
2D 描述入射光的方向。
2D 描述出射光的方向。

可以通过散射的记忆效应减少存储需求。

渲染时，对于每个入射光方向，查询相邻的表格数据，应用相应的角度偏移（记忆效应中的角度偏移量）。

最终使用的表格大小为 $25 \times 32 \times 72 \times 180 \times 360$ ，非常大，约15GB的内存需求。

8. Result

与XIA[2020]相比，新模型更加鲜艳。色彩光斑（colorful glints）更好。体现多种波长对毛发产生的光斑。

9. DISCUSSION AND CONCLUSION

第一篇能够实际应用的 3D 波动光学纤维散射模型。

能够生成光学散射中常见的细微光斑（speckle patterns）。

虽然 3D 模拟的精度很高，但它的内存占用和计算成本非常大。提出了噪声实用模型，通过捕捉纤维散射光斑的统计特性（如自相关函数），减少计算量。

目前该模型主要用于第一阶反射模式的模拟。未来可以扩展到更高阶的散射模式。

高阶模式时可能需要进一步开发新技术来预测或测量高阶模式的统计特性。