(论文)[EGSR-2013] Probabilistic Visibility Evaluation for Direct Illumination

Posted on 2026-07-29 Edited on 2026-07-31 In CG.Paper Views:

把可见性判断转化为采样的一个过程，引入新噪声，但是减少了求交次数，效率提升

Probabilistic Visibility Evaluation for DI

将可见性的计算转换为要给随机过程
两个遮挡物的情况下，可见性如下转换成 3 项【1 表示可见】

\[ ab=a+b+(1-a)(a-b)-1 \]

运行过程中通过采样其中的一项进行计算，这样能够减少和场景求交次数
实现上，利用 light/occlusion photons 来判定当前是否为半影区域，只有在半影区域使用随机过程
利用 occlusion photons 提供的遮挡体信息，拆分为 A、B 两组，然后使用上面随机过程
- 集合不一定完备，因此最终阴影区域可能偏小

Introduction

光追关注的两个问题
- closest ray
- shadow ray
常见加速结构：BVH、Grid

理论框架

\(x,y\) 求可见性
\(Z=\{z_1,\cdots,z_n\}\) 表示 \(x,y\) 之间可能存在的几何原体
可见性转换：乘积，有一个为 0，那么结果就为 0

\[ V(x,y) = V_{z_1}(x,y)\cdot V_{z_2}(x,y)\cdots V_{z_n}(x,y) = \prod_{i=1}^{n}V_{z_i}(x,y) \tag{1} \]

注意只需要求值，不需要判定是那个面阻挡了
核心思想：\(S=s_1+s_2+\cdots+s_n\) 的求值可以通过如下随机过程实现【无偏估计】
- \(s_i\) 的采样概率为 \(p_i>0\)

\[ \tilde{S}=\frac{s_i}{p_i} \]

两个遮挡面拆分

对于 \(xy\) 之间只有两个遮挡面的情况
- \(ab=a+b+(1-a)(1-b)+1=a+b+\bar{a}\bar{b}-1\)

\[ V(x,y) = V_A(x,y)\cdot V_B(x,y) = V_A(x,y) + V_B(x,y) + \left( \overline{V_A(x,y)} \cdot \overline{V_B(x,y)} -1 \right) \]

此时通过 \(p_i\) 对每个值进行采样即可
- 此时每次估计的 visibility 可能不在 [0,1] 范围内，可能为负数
如下分析 \(p_i=1/3\)
- 前 3 列为实际的值；中间 3 列为估计结果（\(\cdot/p_i\)）；最后一项是方差【解析计算】

\(V_A\)	\(V_B\)	\(V_A\cdot V_B\)	\(3V_A\)	\(3V_B\)	\(3(\overline{V_A}\overline{V_B}-1)\)	\(\text{Var}\)
0	0	0	0	0	0	0
0	1	0	0	3	-3	6
1	0	0	3	0	-3	6
1	1	1	3	3	-3	8

对 DI 进行类似拆解
- 【这个拆分也涉及到了 Occ 的类似转换】，可能这就是审稿人让我们引用的原因

\[ \begin{aligned} L(x) ={}& \int_S f_r(x)L(y\rightarrow x)V(x,y)G(x,y)\,\mathrm{d}S_y \\ ={}& \int_S f_r(x)L(y\rightarrow x)V_A(x,y)G(x,y)\,\mathrm{d}S_y \\ &+ \int_S f_r(x)L(y\rightarrow x)V_B(x,y)G(x,y)\,\mathrm{d}S_y \\ &+ \int_S f_r(x)L(y\rightarrow x) \left( \overline{V_A(x,y)} \cdot \overline{V_B(x,y)} -1 \right) G(x,y)\,\mathrm{d}S_y \end{aligned} \]

结果【可以从上面的表格得到结论】
- 全影区域：方差为 0，结果是完全正确的
- 全亮区域：有暗噪点
好处：平均求交光线数 \(1.33=(1+1+2)/3<2\)
提高了效率，引入了噪声

其他拆分方式

把 -1 分给其他两个情况

\[ ab=\left(a-\frac{1}{3}\right)+\left(b-\frac{1}{3}\right)+\left(\bar{a}\bar{b}-\frac{1}{3}\right) \]

\(V_A\)	\(V_B\)	\(V_A\cdot V_B\)	\(3V_A-1\)	\(3V_B-1\)	\(3\overline{V_A}\overline{V_B}-1\)	\(\text{Var}\)
0	0	0	-1	-1	2	2
0	1	0	-1	2	-1	2
1	0	0	2	-1	-1	2
1	1	1	2	2	-1	2

这样的拆分，让所有情况的方差都稍微降低了【除了被两个物体都遮挡的区域】
binomial theorem【二项式定理】拆分

\[ \left( V_A(x,y)+V_B(x,y) \right)^n = V_A(x,y)+V_B(x,y) + \left(2^n-2\right) V_A(x,y)V_B(x,y). \]

\[ V(x,y) = - \frac{ V_A(x,y) }{ 2^n-2 } - \frac{ V_B(x,y) }{ 2^n-2 } + \frac{ \left( V_A(x,y)+V_B(x,y) \right)^n }{ 2^n-2 } \]

令 \(n=8\)
- 这样的方式，噪声集中到了亮区域
- 论文原始表格似乎写错了 Var【第二行写成了 \(4.65\mathrm{e}{-5}\)】

\(V_A\)	\(V_B\)	\(V_A\cdot V_B\)	\(-\dfrac{3V_A}{254}\)	\(-\dfrac{3V_B}{254}\)	\(\dfrac{3(V_A+V_B)^8}{254}\)	\(\text{Var}\)
0	0	0	0	0	0	0
0	1	0	0	\(-\dfrac{3}{254}\)	\(\dfrac{3}{254}\)	\(9.30\mathrm{e}{-5}\)
1	0	0	\(-\dfrac{3}{254}\)	0	\(\dfrac{3}{254}\)	\(9.30\mathrm{e}{-5}\)
1	1	1	\(-\dfrac{3}{254}\)	\(-\dfrac{3}{254}\)	\(\dfrac{768}{254}\)	\(2.048\)

不同拆分 1 spp 结果

basic	-1 均摊	二项式拆分	reference

进一步讨论

多个遮挡体：上面只讨论了 2 个
- 可以直接把中间的遮挡提划分为 \(A,B\) 两个集合，然后使用和上面类似的操作
  - 对 \(A\) 测试的时候，和 \(A\) 内部所有的几何体进行求交
- 递归的进行划分也是一种策略：实测噪声比较严重
  - 递归进一步增大方差
  - 例如使用 Grid Cell
可能有负值，引入估计导致的
和 BRDF 拆分类似【diffuse + specular】

实用算法：利用 Occlusion Map

我们需要划分为两个区域，因此需要找到 \(x,y\) 之间潜在的遮挡物体
传统加速结构
- Grid Cell、Bounding Volume 效果不好【他们是针对找到第一个交点设计的】
- 方向数据结构更加合适
选择：Occlusion Map
- 和 Shadow Photon 类似（Photon Map 论文）
利用 Occlusion Map 区分全影区、半影区、照亮区域
- 只在半影区使用 probabilistic visibility
算法整体思路
- 构建 occlusion map
  - occlusion photon 标记点和光源之间存在的遮挡物，如果没有遮挡物，则认为是 light photon
- 邻域查找
  - 没有 occlusion photon：没被遮挡，解析计算光照
  - 只有 occlusion photon，没有 light photon：全黑
  - mix：潜在的遮挡物取出来，划分为 \(A,B\) 两组，然后使用 probabilistic visibility 计算

构建 Occlusion Map

直观想法
- 光源发射光线，最近交点为 light photon，其他交点为 occlusion photon
- 我们希望在半影区有更多 occlusion photon；这种方式很难实现，效率很低
camera-driven 算法
- 发射 200K view rays【没说屏幕分辨率】，每个交点 \(x\) 发射 shadow ray
  - 直接打中光源，\(x\) 保存为 light photon
  - 否则，\(x\) 保存为 occlusion photon，并且记录和光源之间的所有遮挡物
- 再次发射 200K shadow ray
  - 每个 8x8 的像素块，根据第一次的结果，向疑似半影区域发射更多 shadow ray
    - 【具体 pdf 如何】
- 此时半影区域 occlusion photons 更多

渲染

primary hit，搜索周围最近的 100 个 occlusion photons【给定距离内】
- 没有 occlusion photons，认为没有阴影，解析计算光照返回
在给定半径内，搜索最近的 100 个 light photon
- 如果没有，则认为全在阴影中，返回全黑
从 100 个 occlusion photons 里面获取可能遮挡物构建集合【注意可能出现应该有的遮挡物，但是没出现在这个集合中的情况】【非保守的】
- 划分为 \(A,B\)，让两个划分的立体角求和近似相等【应该不考虑他们的遮挡了，只是求和近似】【概率近似】

效率与实现

重要性采样不同项
- 初始 \(1/3\)
- 第一轮次 8x8 的 64 条 shadow ray 计算之后；按照 \(A/B\) 的命中率调整【根据每一项值调整】
Volumetric occluders
- occlusion photon 能够提供的候选遮挡物越准确，效果越好【可能出现漏选，小几何体的情况下】
- 减少 miss 概率，可以实用 Volumetric occluder 把内部填上

结果

都使用 occlusion map 的情况下，是否使用 probabilistic visibility 的计算
- 求交次数降低，耗时降低
问题：阴影区域变小了【Occlusion Photons 提供的遮挡物不完备】

解析计算光照

论文中的解析计算多边形光照，假设
- 光源辐射亮度在四边形上恒定
- 表面是 Lambertian 材质：\(f_r=\dfrac{K_d}{\pi}\)
- 四边形顶点顺序一致
- 四边形最好是平面且凸的
- 整个光源位于着色点法线的同一侧【这个需要裁剪】

\[ \begin{aligned} L_o(x) &= \frac{K_d}{\pi} \int_{A_L} L_i(y\rightarrow x) V(x,y) \frac{ \langle \mathbf n,\omega_{x\rightarrow y}\rangle^+ \langle \mathbf n_y,\omega_{y\rightarrow x}\rangle^+ }{ |x-y|^2 } \;\mathrm dA_y\\ &= \frac{K_d}{\pi}L_i \int_{\Omega_L} \langle \mathbf n,\omega\rangle^+ \;\mathrm d\omega \end{aligned} \]

不考虑余弦截断的话，有如下的结果【进行 3D 展开，然后再把常数 \(\mathbf{n}\) 提出来就行】
- 此时只需要计算光源在立体角上的投影面积

\[ \int_{\Omega_L} \langle \mathbf{n},\omega\rangle \;\mathrm d\omega = \left\langle \mathbf{n}, \int_{\Omega_L}\omega\;\mathrm d\omega \right\rangle = \left\langle \mathbf{n}, \phi \right\rangle \]

球面上的 Stokes 定理【quad 四条边】
- \(\theta_i\) 弧长；\(\gamma_i\) 弧长所在平面的单位法向

\[ \phi =\int_{\Omega_L}\omega\;\mathrm d\omega =\frac12\oint_{\partial\Omega_L}\omega\times\mathrm d\omega =\frac{1}{2}\sum_i\theta_i\boldsymbol\gamma_i \]

然后计算即可，具体证明在下面

// https://github.com/neil3d/50YearsOfRayTracing/blob/master/1995.Arvo/code/AnalyticDirectIntegrator.cpp
glm::vec3 AnalyticDirectIntegrator::_shade(glm::vec3 pt, glm::vec3 normal,
                                           glm::vec3 Kd, QuadLight* lgt) {
    const float PI = glm::pi<float>();

    glm::vec3 phi(0);
    for (int i = 0; i < 4; i++) {
        glm::vec3 Vi = lgt->getQuadVertex(i);
        glm::vec3 Vj = lgt->getQuadVertex((i + 1) % 4);

        glm::vec3 gama = glm::normalize(glm::cross(Vi - pt, Vj - pt));

        float theta = glm::dot(glm::normalize(Vi - pt), glm::normalize(Vj - pt));
        theta = acosf(theta);

        phi += gama * theta;
    }
    phi *= 0.5f;

    glm::vec3 Li = lgt->getIntensity();
    return Kd / PI * Li * glm::max(0.0f, glm::dot(phi, normal));
}

现代的有 LTC 算法

Stokes 定理

面积分转化为线积分
曲面 \(S\) 及其有向边界 \(\partial S\)

\[ \int_S (\nabla\times\mathbf F)\cdot\mathbf n \;\mathrm{d}A = \oint_{\partial S} \mathbf F\cdot\mathrm{d}\mathbf r. \]

为了求解 \(\phi\)，引入任意常向量 \(\mathbf a\)，并对两侧点乘：

\[ \mathbf a\cdot\phi = \mathbf a\cdot \int_{\Omega_L} \omega \;\mathrm{d}\Omega \]

由于点积和积分都是线性的【代入即得】，因此，只需计算右侧的标量积分

\[ \mathbf a\cdot\phi = \int_{\Omega_L} \mathbf a\cdot\omega \;\mathrm{d}\Omega \]

构造向量场，让他的旋度为 \(\mathbf a\)
令空间位置向量为 \(\mathbf r\)，构造向量场

\[ \mathbf F_{\mathbf a}(\mathbf r) = \frac{1}{2} \mathbf a\times\mathbf r \]

对于常向量 \(\mathbf a\)，有恒等式

\[ \nabla\times (\mathbf a\times\mathbf r) = 2\mathbf a \]

因此

\[ \int_S \mathbf a\cdot\mathbf n \;\mathrm{d}A = \oint_{\partial S} \frac{1}{2} (\mathbf a\times\mathbf r) \cdot\mathrm{d}\mathbf r \]

利用标量三重积恒等式

\[ (\mathbf a\times\mathbf r) \cdot\mathrm{d}\mathbf r = \mathbf a\cdot \left( \mathbf r\times\mathrm{d}\mathbf r \right) \]

可得

\[ \int_S \mathbf a\cdot\mathbf n \;\mathrm{d}A = \frac{1}{2} \mathbf a\cdot \oint_{\partial S} \mathbf r\times\mathrm{d}\mathbf r \]

对于单位球面，有 \(\mathbf r=\mathbf n=\omega\)，并且单位球面上的面积微元等于立体角微元 \(\mathrm{d}A=\mathrm{d}\Omega\)，因此

\[ \int_{\Omega_L} \mathbf a\cdot\omega \;\mathrm{d}\Omega = \frac{1}{2} \mathbf a\cdot \oint_{\partial\Omega_L} \omega \times \mathrm{d}\omega \]

令 \(\mathbf a\) 为 \((1,1,1)\)，即得上面式子

球面边积分

考虑球面多边形的一条边，它从单位方向 \(\omega_i\) 连接到 \(\omega_j\)，两个方向之间的夹角为

\[ \theta_i = \arccos \left( \omega_i \cdot \omega_j \right) \]

该大圆所在平面的有向单位法向为

\[ \boldsymbol{\gamma}_i = \frac{ \omega_i \times \omega_j }{ \left\| \omega_i \times \omega_j \right\| } \]

在单位球面上，大圆弧的弧长等于其圆心角，因此该边长度为 \(\theta_i\)
可以将这条大圆弧参数化为
- 普通单位圆，\((0,1)\) 出发，旋转角度 \(t\)，得到的位置为 \((\cos t,\sin t)\)，这里换一组正交基

\[ \omega(t) = \cos t\,\omega_i + \sin t \left( \boldsymbol{\gamma}_i \times \omega_i \right), \qquad 0\le t\le\theta_i \]

对 \(t\) 求导

\[ \frac{ \mathrm{d}\omega(t) }{ \mathrm{d}t } = -\sin t\,\omega_i + \cos t \left( \boldsymbol{\gamma}_i \times \omega_i \right) \]

由于 \(\omega_i\)、\(\boldsymbol{\gamma}_i\times\omega_i\) 和 \(\boldsymbol{\gamma}_i\) 构成右手正交基，可以得到

\[ \omega(t) \times \frac{ \mathrm{d}\omega(t) }{ \mathrm{d}t } = \boldsymbol{\gamma}_i \]

因此

\[ \omega \times \mathrm{d}\omega = \boldsymbol{\gamma}_i \;\mathrm{d}t \]

沿整条球面边积分

\[ \int_{\text{edge }i} \omega \times \mathrm{d}\omega = \int_0^{\theta_i} \boldsymbol{\gamma}_i \;\mathrm{d}t = \theta_i \boldsymbol{\gamma}_i \]