Differentiable Monte Carlo Ray Tracing

本文是‘Differentiable Monte Carlo Ray Tracing through Edge Sampling’这篇论文的学习总结。论文作者是Tzu-Mao Li(李子懋 mào), Miika Aittala, Frédo Durand 以及Jaakko Lehtinen。李子懋的博士论文（dissertation）获得了SIGGRAPH2020最佳博士论文奖，称其‘为新兴的可微计算机图形学奠定了基础’。而这篇论文(paper)则首次提出了如何解决光路微分的方式(propose a comprehensive solution to the full light transport equation)。插一句题外话，我没有读论文Acknowledgments的习惯，我无意间瞥了这篇论文其中的一段，驱使我读完了整个Acknowledgments，推荐大家读一下该论文的Acknowledgments。

首先，作者提到了该问题的难点在于光路中的遮盖问题，求导时因为被积函数中包含Dirac delta函数，因而无法使用传统的采样方式。之前的方法要么只解决了Primary ray的可见性问题，不支持second ray相关的阴影和间接光，要么只考虑了diffuse的材质， 不支持复杂材质。然后给出了本论文的目标：一个针对全部光线传输求导的完整的解决方案，这样，则可以对场景中所有参数求导并反向传播。这部分分为两种情况，对于被积函数平滑的部分，可以用传统的方式，借助自动微分(automatic differentiation)技术求导，对于不连续的部分，论文中提供了一种全新的方式来捕捉边缘的变化。同时，出于性能的考虑，提供了一种基于视角来优化边缘采样的技术。该论文只讨论静态场景，不考虑时间维度的影响，比如motion blur。

Primary Visibility

第一个解决的是primary ray，也就是直接从相机发射的ray。因为渲染方程通常没有对积分的解析解(closed-form solution)。我们默认所有物体都是三角网构成的mesh，不连续性只存在于三角形的边。这样，一个三角形的边会把空间分为两部分(half plane)，如上图的f_u和f_l。这样，我们用一个Heaviside step function θ来表示（Heaviside step function的导数是Dirac delta function），其中函数α表示对应的三角形的边：

如果一个边的两个端点分别为

，

，则函数α为：

这样，我们可以重写该积分项，显示的表达不连续区域：

这里，f_i函数包含多个Heaviside step functions。这样，当我们对该积分求导时，则有如下的推导：

如上，该积分项分为了两部分，左边的是非连续区域，右侧的是平滑区域。我们用传统的autodiff方式求解右侧积分项，而左侧的有一项δ（Dirac delta function），只有当α函数为0时值为1，我们将该二重积分转为对面积的积分：

这样，当我们想要对场景中某一个参数p（比如相机位置，mesh中的点位置）求导时，我们分别计算积分变化对投影三角形边上某点（x,y）的变化（我对投影的理解是因为primary ray，所以需要针对screen space），然后在计算该点相对p的变化，最终通过chain rule获取我们想要的导数（梯度）：

这里，因为half plane，根据公式（1），我们用f_l和f_u的差值来表示f_i，这里需要两次采样，最终，通过MC采样，对不连续区域求导的求解方式如下：

这里，

是采样边的长度，P(E)是采样的pdf。

至此，我们给出了无偏的一般解，对primary ray中不连续区域求导，但如上图，这种情况下就需要特殊处理了，假设红色和绿色的三角形完全重叠，也就是这个边相对视角完全平行，这时，当我们移动红色三角时，会看到绿色的部分，但在计算时会误认为是空白区域。

Secondary visibility

相比primary ray， secondary ray是一个从2D升级到3D的过程。如上图，我们同样可以用step function θ来表示：

这里，函数α是v_0,v_1和p组成的平面：

同理，在不连续区域，我们对场景中某边积分:

是点m的法线方向。

是边上的一点从点p到点m的投影长度，

等于sinθ，场景表面和边缘表面的无限小的宽度比（面积比）。以下是个人理解，不保证理解是否准确，secondary是primary从2D升级到3D的过程，之前2D下

和3D场景下重合，在secondary中，我们可以认为p点是相机，

是3D场景下的边缘线，而

则是其在相机视角下的2D投影，这里的投影面是上图中褐色的面，黄色的scene surface投影到edge surface时，需要进行面积的补偿（area correction）。n_h所在的个面(edge surface)有很多，个人不太确定是否任选一个就可以，在primary中，相机的朝向是确定的，如果假设在p点处设置一个相机，相机的朝向并没有确定，当然，面(p,v_0,v_1)是最合适的朝向。

我们再进一步，用一个参数t表示边缘的某一点：

，然后用该点来表示表面上的点m，换元后为：

Cameras with Non-linear Projections

在之前的相机都是projective camera，3D空间中的一条线，在2D平面下也是一条线，在非线性投影的相机中，比如鱼眼，这个假设并不成立，如上图。

在这种情况下，我们以相机位置为点p，找到v_0和v_1，基于secondary visibility的方式，绕开这个问题，将非线性的因素转移到相机参数与v_0和v_1之间，利用chain rule最终获取对应参数的导数。

Relation to Reynolds transport theorem

如上的数学推论也可以用流体力学中的雷诺运输定理来解释，以primary visibility为例：

这里，

表示boundary，v是移动的速度，n是边缘点的法线。如果我们对x方向求导，速度是(1,0)，法线方向为(a_y-b_y,b_x-a_x)，这个结果和primary visibility中的推导结果一致。

Importance Sampling the Edges

对边缘的重要性采样。这部分我看的不是很详细。原因是因为对边缘的识别很复杂，对primary，因为相机固定，可以预先找到，我们基于边缘的两个端点在的屏幕距离进行均匀采样，但在secondary中，随着视点的不同，边缘也会不同(view dependent)，这种情况下要高效的采样，性能很差，因此需要优化。在重要性中，要考虑该边缘在视点下是否为边缘（silhouette），和视点距离成反比的边缘长度，材质以及来自该边缘的光照强度（radiance）。

论文中采用了霍夫曼空间（Hough space）的策略，空间中的一条边edge对应霍夫曼空间的一个点（duality），形成一个boundary。这样，在第一层中，基于相机视点下的silhouette，构建一个BVH。第二层则处理其他的边，包括该边的两个端点以及与边的两个面关联的两个平面，在查询时，我们构建一个从shading point p到相机位置的球，然后判断该球是否和边相交。而重要性采样时，我们以边的长度，乘以该边到包围盒中心点的距离的倒数，再乘以材质的权重（Linearly Transformed Cosine Distribution）。

这个采样方式很难对一部分边阻挡小光源的边采样，这些边缘导致具有大梯度的尖锐边界。在论文中采用Next event estimation方式，其中有一个billboard的概念（相对于NEE的ray），个人理解单个billboard有一定的宽度，包含相关的边。然后对所有相交的billboards进行重要性采样，而在采样的billboard中，采样点是该交点投影在边上的最近点。

结果

上图是对应spp下的效果，在实际中，该论文对每个像素采用4spp

上图是和传统的方法进行对比，OpenDR和Neural都是先用栅格化+z buffer的方式得到渲染结果，然后利用渲染结果，对屏幕空间的三角形顶点位置计算其梯度。OpenDR基于亮度恒定假设(a brightness constancy assumption,不懂)进行屏幕像素的平滑，进行过滤的形状则取决于基于三角形ID形成的边缘。Neural会对三角形的边缘进行额外的栅格化渲染，根据边缘两侧的色差计算并累加对用的梯度。当像素间的亮度因为光源而变化时，这两种方式下的计算都不正确。

因为我们的渲染结果可以是自定义的，不仅仅是渲染RGB，可以是用户自定义的渲染结果，比如上图，自然也可以针对该结果计算其梯度。

上图是通过计算梯度，解决可逆的问题。(a)是初始效果，(b)是目标效果，经过迭代不断缩小和(b)的差异，我们最后得到了效果(f)。其中(c)是相机向屏幕方向移动时的像素梯度。(d)是对应于桌子颜色中的红色通道 的梯度，(e)是相对于光源绿色通道的梯度。

在神经网络中，基于梯度的优化方法，我们进行反向传播，让渲染图像更符合我们预期的结果，用于生成对抗样本。上图，原始的场景下，图片分类中，street sign的比率为53%，5次迭代后handrail的比率为26.8%，但其他的可能比率也比较高，在25次迭代后，handrail的比率为23.3%，street sign和traffic light的比率均为3.39%。

不足

论文中提到了该方法的瓶颈在于边缘的采样和对光路的autodiff。同时，我们默认材质的BRDF是连续的。分辨率非常高的纹理也需要进行预先的平滑，减少方差。作者也提到了如何让该算法支持BDPT和photon mapping，也是一个很有意思的方向。