本文是论文‘Mitsuba 2: A Retargetable Forward and Inverse Renderer’的读后感（review）。
Mitsuba和pbrt可以说是最主流的两款PBR开源渲染引擎，两者各有利弊，我对pbrt较为熟悉，pbrt的好处是有PBR book辅助，如果你能读懂PBR book，则可以帮助你很好的理解代码。Mitsuba的好处是很多论文的算法都有正确的实现，功能相比PBR更为丰富，如果你希望验证某个算法，或者基于该算法做一些自己的扩展，Mitsuba是一个很好的实验选择。
Mitsuba 2（M2）基于Mitsuba 0.6进行了升级，增强了并行的能力（SIMD和GPU），同时支持了自动求微分（automatic differentiation，AD），具备了对渲染求逆（微分渲染）的能力。这篇论文是对M2的详细介绍，作者是Wenzel教授领导的EPFL RGL团队成员MERLIN NIMIER-DAVID，DELIO VICINI和TIZIAN ZELTNER，此处假设可以点赞。
下面是对该论文的重点概括以及个人解读。包括简介，相关工作，Mitsuba设计思路，应用案例以及总结。
图(a),分析光学实验中的偏振光渲染，(b),coherent MCMC算法，用一组光路来代替单个光路，增强并行能力，并借助M2并行能力提高性能，(c)，inverse rendering的应用，利用梯度优化算法不断调整前方玻璃板的每个grid的折射率，进而模拟出一幅图画，(d) inverse rendering的应用，根据效果图反推participating media中单位voxel对应的density。
简介
客套了几句PBR的重要性后，首先，针对并行的需求，谈到了Vectorization的重要性和难点，其次则是基于梯度的优化技术在机器学习和视觉领域的显著发展，differentiable rendering在理论上已经解决了不连续性这个难题，因此，理论上我们可以进行AD，实现逆渲染。但实际中，如何在CUDA中实现光线追踪本身和AD还是一个难点。这正是M2要解决的问题，用户调用M2提供的一组标准组件编写自己的应用，根据用户的需求和场景，最终会系统的将这些抽象接口编译为对应的执行程序。个人理解，这里的编译包括两步，第一是M2自己的预编译，针对用户的需求进行优化，比如应用场景是CPU，则将抽象接口对应倒AVX的实现，如果是GPU，则对应CUDA的APU来实现，是否需要对某参数求导，这类逻辑的优化；其次则是程序的编译，比如C 或者CUDA将代码转为机器码，还有一些内容，可以借助C 的template特性，在编译期执行，减少了运行期间的计算量，这类优化介于两者之间，属于Mitsuba实现的功能，但会在编译期间由编译器来实现。
总体上，M2的主要贡献是：
- 一个可组合类型的通用框架，一个渲染系统重定位适合一系列不同任务的具体实现，比如上图中不同的应用，有forward rendering的，也有reverse rendering的。
- 一个 lazy just-in-time (JIT)编译器，能够将算术运算和控制流生成为适合GPU执行的kernel
- 一个graph-based的方案，用于计算forward或reverse-mode的AD
- 一个简化算法，用于优化这个graph，减少AD对内存的消耗
相关工作
Coherent and vectorized rendering
这里，谈到了MC方法，比如光线追踪path tracing，因为采样的随机性，算法很难在并行中有很好的性能，coherent rendering技术则试图解决该问题，大概的思路是物以类聚，减少一个kernel中divergence的情况，基于材质或ray的相似性进行分类，提高并行计算的效率。这里，问题的关键是如何在原本incoherent的算法中发现coherent的部分。
Automatic differentiation（AD）
之前在《differentiable rendering小结》中，也介绍过AD分为forward和reverse两种，取决于input参数和output参数的多少，input参数多的，可以用reverse，而output参数多的则可以采用forward的方式。pbr渲染中，output就是image，比如1024*768像素，每个像素对应output的一个参数，看上去很多，但input可能的参数是相机位置，不同位置的object对应材质的不同属性，灯光等，每一次intersection都会产生新的input参数，因此，通常会比output多很多，这也是可微分渲染通常采用reverse AD的方式。
Domain Specific Languages（DSL）
M2中支持更多的转换，这里我不太确定这里提到的DSL具体指什么方面的转换，比如支持偏振光，光路求导（可微分渲染）这些功能，集成到M2中，另一个理解是用户可以编码，然后根据具体的场景编译成CPU或CUDA支持的代码。如上是个人的理解和猜测，论文并没有过多介绍这部分。
Differentiable rendering
M2提供了基于GPU的，从光线传输算法中高效求导的应用，因此可以实现可微分渲染方面的应用。同时，对光路求导并不仅仅用于inverse rendering，在复杂的场景中可以利用导数的特点来优化MCMC算法的采样。
Template metaprogramming（TMP）
M2采用的模板元编程，支持C 17标准，通过Variadic，Compile-time conditionals（比如constexpr），Type computation（比如decltype）等方式，把计算尽可能的放在编译期，进而提高运行时的效率。
Expression templates（ET）
M2尝试使用expression templates技术（具体可以阅读C templates第二版的章节），将计算尽可能的延后，避免创建过多的临时变量，实现效率的优化，在Eigen中大量采用这种技术。但发现在渲染复杂场景下，ET的效果并不理想，因此最后放弃了这种技术。
系统设计
该系统由两部分组成，基础的enoki，这个是一个模板库，负责实现vectorization，JIT编译以及程序在不同端的转换，这部分不涉及到渲染。第二部分则是基于Mitsuba 0.6的渲染引擎，提供了forward和inverse的渲染管线。
Static arrays
M2中通过C template，比如Array<Value, Size>，可以固定数组的长度，这种方式的好处是可以预知数组的长度，方便在编译期的并行化改造，并不会导致运行时的消耗。M2在兼容Intel的CPU中支持SSE4.2，AVX，AVX2以及AVX512。
在GPU中就有点复杂了，如下代码：
1 using Float = GPUArray<float>;
2 using UInt64 = GPUArray<uint64_t>;
3 using Vector3f = Array<Float, 3>;
4 PCG32<UInt64> rng(arange<UInt64>(1000000));
5 Vector3f v(rng.next_float(), rng.next_float(), rng.next_float());
6 size_t inside = count(norm(v) < 1.f);
上述代码如果转化为CUDA的代码，比如4-6行对应三个操作，可以用三个kernels实现。这种逐行生成GPU代码的方式简单直接，但kernel之间需要内存的交换（GPU-CPU-GPU），导致了效率的降低。对此，M2中采用JIT的思路，尽可能延后计算，使用NV的Parallel Thread Expression（PTX）构建符合single static assignment（SSA）的程序，这里对变量采用了引用计数，如果某变量只是局部的，不会作为返回值，这样，当函数执行后，其引用计数为0。对于这一类的方法，可以合并在一个kernel中一起执行，这些变量只需要保存在寄存器中，而不需要用global memory。这种方式可以提高性能，降低内存的消耗。
Relation to existing frameworks
如上的思路和tensorflow和pytorch中的思路类似。在TF和Pytorch中有两种模式。第一种是eager mode，直接计算。这种方式在计算密集的操作中具有较好的性能，比如卷积，矩阵计算等。第二种操作模式需要预先规范完整的graph以生成单个优化的GPU kernel。深度学习中这个graph通常不复杂，而在M2中，渲染的流程非常复杂，graph需要通过虚函数实时分成不同的应用分支，从而动态的构建这个graph，在GPU中，最终通过JIT技术，融合多个指令（通常超过100k的指令），创建一个大的kernel。
Autodiff backend
M2中支持forward和reverse两种方式，同时，可以使用DiffArray<gpuarray<float>>,该方式是autodiff和GPU的组合，借助JIT技术，实现在GPU中的AD，测试发现，无论是float还是double，该方式相比Redner（另一款可微分渲染）下的AD性能要高。</gpuarray<<>
Graph simplification
如图，因为pbr中input参数很多，复杂场景下会消耗较多内存，出现内存崩溃的问题，所以M2尝试在graph中去掉一些节点，以增加计算量的方式来减少内存的消耗。
Custom data structures
用户可以自定义数据结构，如下结构体用于记录射线和物体相交时的相关属性：
template <typename Point3f> struct SurfaceInteraction {
using Float = value_t<Point3f>;
using Vector3f = Vector<Float, 3>;
using Frame3f = Frame<Vector3f>;
using UInt32 = uint32_array_t<Float>;
using Shape3f = replace_scalar_t<Float, const Shape<Float> *>;
Float t; // ray distance
Point3f p; // position
Vector3f wi; // incident direction
Frame3f sh_frame; // shading coordinate frame
UInt32 prim_id; // intersected primitive (e.g. triangle ID)
Shape3f shape; // pointer to Shape<..> instance
/// ...
};
自定义结构体在CPU和GPU中都支持，但M2还会进行Structure of arrays（SOA）改造。
Horizontal operations
这类操作，比如count会极大的阻碍GPU的性能，比如kernel1 kernel2，其中kernel2属于horizontal，此时，GPU会在执行kernel2时做一次flush，这是一个同步操作，需要等待所有的kernel1计算全部完成，性能上会有损失。M2中提供了一种缺省选择，比如any_or<true>，此时，如果需要生成GPU的代码，则返回值为缺省值（true），在CPU则会判断真实值。这种方式则可以避免GPU中不必要的horizontal操作，减少对异步计算的干预。
Scatters and gathers
简单说，scatters是有一个点，然后推断其他N点的值，gathers则相反，有N个点，推断一个点的值。M2还提供了atomic scatter-add，用来查询纹理或体（可以理解为3D纹理）。这里，还提到了在求微分时，两者之间互逆，比如c=a*b,求导结果是dc=da*b a*db，在forward传播中，这是一个scatters，而在reverse传播中，则是一个gathers。
Method dispatch
这就是一个物以类聚的过程，在CPU中，通过虚函数指针来进行分类，而在GPU中，专门有一个kernel来执行分类（并行的基数排序），在渲染中，比如有n1个射线命中了某个BSDF材质，有n2个命中了另一个BSDF材质，这个变量会更新，属于runtime kernel parameter，M2不需要重新编译kernel就可以更新这些变量（个人理解就是设置一个全局变量或结构体来动态更新这些变量就可以实现）。
Mathematical support library
M2提供了一些数据函数，比如复数，矩阵，四元数等，以及求矩阵的秩，逆等操作，这些方法也都支持AD，JIT以及vectorization的能力。
Language bindings
封装成python接口
应用案例
Coherent MCMC sampling
该算法的思路是用一组光线的贡献值来取代一个光路的贡献值，本身其实是一种平滑手段，在光线追踪中，这增强了local exploration。同时，因为都是相邻的光路，增强了coherent，可以更好的利用并行计算提高性能。
Caustic design
论文中主要有二种方式，surface displacement和gradient-index optics。如上图，当平行光射入玻璃板时，根据当前的折射率，得到投影后的效果，然后和目标效果进行对比，利用梯度调整折射率，通过迭代的方式，投影效果会逐步接近目标效果。
如图，通过inverse rendering，基于梯度调整单位grid对应的折射率，迭代的方式最终达到期望的效果(v)。(a)是简单的玻璃板，调整折射率后达到最终的效果；(b), (c)是有颜色的玻璃板，在这个情况下调整折射率。(e)-(d)属于gradient-index optics，根据Eikonal公式计算光路的弯曲，最终在不同的侧面达到对应的效果。
Heterogeneous Participating Media
如图，基于VRE，我们利用梯度调整density，最终接近期望效果。这里，需要注意的是，在participating media的边界处是不连续的，无法求导，而我们在Primary Sample Space中无法识别这种不连续性，因为丢失了path space信息。所以，我们在对VRE求导时，先基于path space求导，然后在把对应的结果对应到PSS。
另一个应用属于subsurface，次表面散射。比如3D打印，因为材质半透明的特性，看到的颜色是多次折射后的综合效果，不仅仅取决于当前点对应的颜色，因此打印后的效果会有失真的问题。
论文中对问题进行简化，认为σ_s σ_t是一个常量，而σ_s/σ_t则会变化。直观可以理解为假设3D打印的材料的extinction（该参数决定光线在材质中的衰减率）相同，只有颜色(albedo)这个参数会变化。然后基于次表面散射的公式，通过inverse rendering，最终确定这个面板每一个grid对应的颜色值。实验中，目标图片为128*128像素，3D打印的材料是256*256*64(grid)，采样64spp。
总结
论文首先给出了渲染引擎的目标，准确，快速，可微分。在光线追踪下，这种高度专业化的实现并不容易理解和维护。这也是M2致力解决的问题。M2通过TMP，JIT，AD等技术，减少了运行时的计算量，发挥硬件的并行能力，并通过梯度来解决inverse problem。
论文也给出了目前M2的两点不足，第一是，因为运用template，不容易调试。文章提到了C 20的concepts，或许可以更好的定位缺陷。其实，我也有类似的问题，如果有机会，我一定要问问为什么非要在Linux下开发，VS不香吗？第二个问题涉及到可微分渲染，目前需要较多的内存，并且在解决不连续区域时性能不高这两个问题，在另外几篇论文中提到了一些可行的解决思路。
我个人读完后的总结是，第一，GPU是大势所趋。相比Mitsuba 0.6，M2扩展了并行能力。无独有偶，今年的HPC会议上，PBR book的作者之一Matt Phar也谈到了将pbrt从CPU扩展到GPU的大概思路。传统的GPU ray tracing也有了很多发展，随着硬件的不断加强，如何利用硬件来改善光线传输算法，甚至是纯光线传输的渲染，这一领域有很多有意义的事情可以尝试。第二，我对C 语法的学习停滞在十年前，读这篇论文时有很多涉及到C 和CUDA方面的知识盲区，这些年，随着CPU Cache硬件方面的优化，SIMD等并行技术的推广，还有C 标准的快速发展，真的有必要系统化的重学C 了，有一些好的新特性都可以有选择的学习，重新定义自己的编码风格。第三，一个人的精力和能力都是有限的，真的需要一个团队来互相配合，互相学习，感觉，在图形学这个领域，需要Geometry和渲染两方面的专家，还需要一个编程专家，如果有条件，最好还能有仿真，视觉方面的专家，这样的一个团队就如同路飞的草帽军团，是可遇而不可求的一种事情。火影中君麻吕看到沼泽中的一朵花，像极了他自己，于是愤怒道“你为什么要在这里绽放，又没有人会看到，这样的生命有什么意义？”这时，大蛇丸出现在他的身后，说：“生命本身没有意义，但是会遇到很多有意义的事情，比如你遇到了这朵花，而我遇见了你。”这是我最喜欢的火影片段，就如同喜欢的歌不愿意随意分享。毕竟能看到这的人不多，就多说两句了，不知道有几个人能忍到这里：）
