基于Metropolis光线跟踪的组合滤波器

2016-11-01吴熙徐庆卜红娟王征

计算机应用 2016年9期

吴熙　徐庆　卜红娟　王征

摘要：

蒙特卡罗方法是计算全局光照的基础，目前已经有很多基于蒙特卡罗的全局光照算法，但大多数算法在渲染时间上都有一定局限性。在蒙特卡罗方法基础上，结合Metropolis光线跟踪算法和组合滤波器，提出一种新的全局光照算法。该算法分为两个部分，首先使用多组不同尺度的滤波器对图像进行处理，然后将多组滤波器处理后的结果组合成最终的结果。该算法使用相对均方根误差作为选择滤波尺度的依据，在采样和重建过程中自适应地为每个像素选择合适的滤波器，以最大化降低误差，得到更好的重建结果。实验结果表明，该算法相对于传统Metropolis算法在效率和图像质量上都有较大提高。

关键词：

Metropolis算法；蒙特卡罗方法；光线跟踪；组合滤波器；全局光照算法

中图分类号：

TP391.41

文献标志码：A

Abstract：

The Monte Carlo method is the basis of calculating global illumination. Many Monte Carlobased global illumination algorithms have been proposed. However， most of them have some limitations in terms of rendering time. Based on the Monte Carlo method， a new global illumination algorithm was proposed， combining the Metropolis ray tracing algorithm with an integrated filter. The algorithm is composed of two parts. In the first part， multiple sets of filters with different scales were used to smooth the image； in the second part， filtered images were combined into the final result. Relative Mean Squared Error （RMSE） was used as a basis for the selection of filtering scale， and an appropriate filter was adaptively selected for each pixel during the process of sampling and reconstruction， aiming to reduce the errors to a minimum degree and gain better reconstruction results. Experimental results show that the proposed method outperforms many traditional Metropolis algorithms in terms of both efficiency and image quality.

英文关键词Key words：

Metropolis algorithm； Monte Carlo method； ray tracing； integrated filter； global illumination algorithm

0引言

目前Monte Claro渲染热门研究方向大致分为两个。第一个方向是通过寻找更加有效的光径来提高图像合成质量。Kelemen等[1]提出一种新的突变策略来提高Metropolis效率，将路径突变用一个空间里面的随机数来代替，能很好地提高Metropolis算法效率。Lehtinen等[2]提出一种基于图像梯度的Metropolis算法，直接计算图像梯度，通过求解Poisson方程直接从图像梯度重建图像。Hachisuka等[3]通過将多重重要性采样和Metropolis算法相结合，提出一种多路复用Metropolis算法，用视点路径的长短来计算多重重要性的权值，增加路径利用率，提高合成图像效率。

第二个方向是对图像进行自适应渲染滤波。Rousselle等[4]提出一种自适应的NonLocal Means算法，使用双缓存的方法，在合成图像过程中生成两个采样数目相同的缓存，并分别计算滤波系数，交叉处理，去除噪声对真实像素的影响。Kaplanyan等[5]提出一种局部自适应的Photon Mapping算法，该算法通过平衡渲染过程中的噪声和偏差来使渲染结果错误率最小。Delbracio等[6]提出一种基于光线颜色直方图来促进图像渲染的方法，该方法通过统计每个像素采样点的颜色分布来判断像素之间是否共享光线。Liu等[7]提出通过梯度来对图像进行聚类分析，然后通过多项式拟合图像的方式进行自适应渲染。Zwicker等[8]对自适应渲染方法作了一个详尽的概述。

为进一步提高蒙特卡罗渲染效率，本文将这两个方向相结合，提出一种新的蒙特卡罗渲染方法。该算法相比传统渲染算法，不是通过暴力增加采样点来提高图像质量，而是通过寻找更有效光径和更适宜的滤波尺度来提高渲染质量。

1Metropolis原理

假设定义一个状态空间Ω和一个非负函数f：Ω→R+，Metropolis算法通过f产生一系列采样点xi。当选定起始点x0，下一个采样点x1将通过x0突变产生，之后每一个采样点xi都是由它的前一个采样点xi-1通过突变产生。这种每个采样点（也称为状态）xi只与它的前一个采样点xi-1有关的随机过程叫作马尔可夫链。每次突变都可能被接收或者是被拒绝，但无

论从哪一个采样点x0开始，最后得到的采样点xi都是服从和函数f成比例的分布，这个分布就叫作平稳分布。

为得到服从平稳分布的采样点，需要使用特定方式来接受或拒绝突变。Metropolis算法采用细节平衡的方法来保证马尔可夫过程能够趋向于平稳分布。

1.1细节平衡

假设当前时刻状态为xi，通过式（1）可以得到状态xi+1。在生成状态xi+1时，首先根据预先制定的突变策略得到临时状态x′。这里突变函数定义为M（|），该函数表示从状态转移到状态的概率，函数M（|）被称为临时转移函数。在状态转移过程中，每次突变得到的临时状态x′都可能被接收或拒绝，将取决于接收概率r（|）。生成新状态表达式：

xi+1=x′，概率为r（x′|xi）

xi，概率为1-r（x′|xi）（1）

如何定义接收概率r（|）是马尔可夫链快速达到平稳分布的关键。假设已经达到平稳分布，那么必须定义转移概率T，使已经达到的平稳分布保持。考虑状态空间中两个状态和，从宏观上来看，它们达到平稳分布。然而在微观上，它们一直进行着相互转移，状态之间的相互转移一直没有停止，它们之间处于一个动态平衡，即从状态到状态的采样点和从状态到状态的采样点一定相等，这就称之为细节平衡[9]，即：

f（）M（|）r（|）=f（）M（|）r（|）（2）

r（|）=min1，f（）M（|）f（）M（|）（3）

这样，就能够使得r（|）和r（|）中较大的被置为1，也就表示两个状态之间的转移，一定会有一个方向被接收，而另一个方向的转移有一定概率会被接收。这样设定接收概率，既可以保证细节平衡也可以保证马尔可夫过程以较快的速度达到平稳分布。

1.2Metropolis渲染

在图像渲染过程中，渲染方程[10-11]描述能量从光源到达人眼的过程，对它求解是图像合成的关键问题。在物体表面上某点x沿方向α出射的辐射亮度的定义如式（4）所示：

L（x→α）=Le（x→α）+

∫Ωxfr（x， β α）L（x ← β）cos β dωβ（4）

这里Ωx是点x所有半球方向的集合，Le（x→α）是点x自身所发出的辐射亮度， fr是双向散射分布函数（Bidirectional Scattering Distribution Function， BSDF），L（x←β）是点x从β方向得到的入射辐射度。

式（4）是渲染方程最基本形式，由它可得到多种变形，比如路径积分形式[12]等，它们最终目的都是方便求取出射辐射亮度。Metropolis光线传输[9，11]将式（4）改写成蒙特卡罗积分形式：

Ij≈CN∑ihj（xi）L*（xi）L（xi）（5）

其中：Ij是像素j处的辐射亮度，C是使用蒙特卡罗方法求得L*的积分值，N是通过马尔可夫链产生的随机路径个数，hj是滤波权值系数，L*是图像的贡献函数，马尔可夫链生成的路径就是和L*成比例的随机数。

2组合滤波器

本文提出的方法中，组合滤波器指在渲染过程中，使用多个尺度不同但类型相同的滤波器对图像进行处理，处理后会得到多组图像滤波结果，然后根据图像之间的特征来选取不同尺度的滤波器进行组合得到最后结果。因为在蒙特卡罗渲染中，图像不同区域噪声水平不同，例如场景边界、阴影等比较复杂的区域渲染时会有较多噪声，而一些平滑区域则噪声较少。所以不同区域需要不同尺度的滤波器以取得最好效果。使用组合滤波器可以在噪声水平不同的区域使用不同尺度滤波器。在这个方法中，有三个关键问题需要考虑：首先是滤波器种类的选取，其次为滤波器之间尺度差异的选取，最后是组合结果时依据的标准。本文算法基本流程如圖1所示。

2.1组合滤波

本节将详细介绍组合滤波器原理。滤波器的选择上，本文选择较为简单的高斯滤波器，因为高斯滤波器能够快速地对图像进行滤波，在渲染过程中能够满足快速合成图像的需求。本文提出的算法选取m组高斯滤波器，定义它们尺度为ai（i=1，2，…，m），它们满足a1

组合滤波器关键点是如何选取组合图像滤波结果的依据。图像渲染过程中，渲染后处理在很多情况下都选取方差作为滤波器的输入信息，Rousselle等[13]使用均方根误差（Mean Squared Error， MSE）作为选择滤波器的依据，Li等[14]选择MSE的一个无偏估计量（Steins Unbiased Risk Estimator， SURE）作为选择双边滤波器的依据。本文提出的方法中使用相对均方根误差（Relative Mean Squared Error， RMSE）作为选择滤波器的依据。

在渲染过程中，因为无法直接知道渲染的真实结果，也就无法直接计算MSE，所以本文使用RMSE选择滤波器。在渲染过程中，可以估算偏差和方差，并以此估算RMSE。在计算RMSE时有三个关键点需要注意。首先，滤波尺度从小到大对应着图像的结果是越来越平滑；其次，对于大多数像素，随着滤波尺度从小到大变化其偏差呈单调递增趋势，而方差则是呈单调递减趋势；最后，偏差和方差之间需要一个权衡系数使二者达到一个相对平衡。RMSE计算公式如下：

RMSE（ai→ai+1）=MSE[ai+1]-MSE[ai]=

ΔVar+ΔBias=ξ*（Var[ai+1]-Var[ai]）+Bias[ai+1]2-Bias[ai]2 （6）

其中：RMSE表示相对均方根误差，ai是滤波器的尺度，MSE指在各自尺度下的均方根误差。本文所定义的RMSE是在不知道图像真实值得情况下估算的结果。为了能够消去真实值，将RMSE定义为两个不同尺度的MSE之差。RMSE定义分为方差和偏差两部分，本文提出参数ξ来均衡方差和偏差之间的差距。

2.2参数计算

在未知图像真实值得情况下估算图像偏差是本文算法很重要的模块。本文用ai和ai+1分别表示小尺度和大尺度的滤波器，pi和pi+1分别表示使用小尺度和大尺度滤波器滤波后的像素值。文献[15]表明滤波尺度和偏差之间的关系：

偏差可以通过式（10）求得，那么求取RMSE就只需解出ΔVar部分即可。对于各个尺度下的方差求取，本文参考Rousselle等[13]中的做法。本文和Rousselle方法的不同之处在于计算RMSE所需的ΔVar需要在两个相邻尺度下的方差之差上乘上权衡系数ξ，系数ξ的选取直接影响滤波效果。因为RMSE是由偏差和方差两部分组成，并且当滤波尺度由小变大时，由于方差和偏差变化的单调性，通常情况下会有ΔVar<0，ΔBias>0，所以RMSE的值就会在0附近变化。因为方差和偏差往往计算出来后不在同一个数量级，所以需要一个系数ξ来平衡它们。通过大量实验得出，方差和偏差之间相差的数量级约为10000，因此ξ取值为10000。为了使该算法适应不同的场景，在不同的场景中都能得到最佳的效果，针对不同的场景还需要对ξ进行细微调整。为实现算法方便引入一个新的系数n用于微调系数ξ，在具体实现时，直接将两个系数相乘即可。通过实验分析，n的取值在1～10的大多数场景都能取得较好的渲染效果。

本文通过RMSE来选择所需要的滤波器尺度。ΔVar和ΔBias分别表示同一像素较大尺度和较小尺度之间的方差和偏差变化。通常在尺度较小时方差变化会大于偏差变化，所以RMSE為负数。本文选择的滤波器是使RMSE绝对值最小时的滤波器，这样能够很好地平衡方差和偏差的变化得到最优的结果。在实际操作中，随着滤波器的尺度由小变大，当RMSE>0时，就选择较小尺度的滤波器，否则就继续遍历。

3实验结果与分析

本文算法基于PBRTV2系统[16-17]实现，所有结果都是在Intel Core i52400 CPU@ 3.10GHz，内存4GB上实现。本算法易实现，本文选用8组高斯滤波器进行组合。为选择合适的系数ξ，进行了大量实验。由实验得出，大多数图像首先需要一个改变方差数量级的系数10000，才能够细微调整方差和偏差之间的平衡系数n，经过实验这个系数可以选择1～10。

图2是平衡系数n为5时得到的Toasters场景每个像素分布16个采样点的尺度选择图像。图中由黑到白分别对应着滤波尺度1～8。由图中可以看出，在三角面中心区域，图像比较平滑，该区域图像的像素值较为相似。所以选择的滤波尺度往往比较大，这样可以使图像更加平滑，去除渲染所带来的噪声；而在三角面边缘区域图像像素值变化较大，所以选择的滤波尺度往往比较小，这样有利于保持图像的特征，使图像不至于变得过于模糊。

图3中给出cornellmlt和Toasters两个场景的渲染结果。

图3（a）是cornellmlt场景渲染结果，图3（b）是toasters场景渲染结果。Lehtinen等[2]和Hachisuka等[3]所提出的对Metropolis方法的改进侧重点和本文有所不同，所以仅选择将本文提出的算法和经典Metropolis方法[1]以及只用单个的高斯滤波进行对比。cornellmlt场景渲染生成的图像大小是250×250，平衡系数n=4，每个像素分布128个采样点，首先计算场景中的直接光照，采用双向路径追踪渲染图像。合成图像大概需要30s，合成速度很快，在可以接受范围内。Toasters场景中，图片的原始大小是400×400，每个像素分布16个采样点，合成图像所用时间大约需要10s。图3中，左边是原始Metropolis所合成的图像，中间是采用单个高斯滤波器所合成的图像，右边是本文算法合成图像结果。总的来说，本文算法能够根据图像特征来自适应选择滤波器尺度，提高渲染效率。

从图3的结果可以直观看出，本文提出的组合滤波器算法渲染结果优于经典的Metropolis算法和单个的高斯滤波器算法。为更加进一步证明本文的方法确实有效，又对PBRT官方网站[17]上的多个经典场景进行了实验。以每个像素2048个采样点所生成的图像作为参考图像，计算图像的峰值信噪比（Peak SignaltoNoise Ratio， PSNR）和MSE，各个场景的渲染结果计算的结果如表1所示。通过表1中MSE和PSNR的比较，可以看出本文的组合滤波器相比原始的Metropolis算法和单一的高斯滤波效果有明显的提高。

4结语

本文基于蒙特卡罗方法提出一种新的渲染方法，该算法将组合滤波器和Metropolis算法相结合，大大提高了蒙特卡罗渲染效率，在较短时间内便可以得到高质量的真实感图像。后续研究中，我们的主要工作是利用更少的采样点，寻找更有效路径来得到高质量的合成图像；同时，将进一步改进本算法，寻找自适应的方法来设置该算法中需要的平衡系数，并且尝试将本算法应用到动画的合成过程中，以提高其渲染的效率。

参考文献：

[1]

KELEMEN C， SZIRMAYKALOS L， ANTAL G， et al. A simple and robust mutation strategy for the metropolis light transport algorithm [J]. Computer Graphics Forum， 2002， 21（3）： 531-540.

[2]

LEHTINEN J， KARRAS T， LAINE S， et al. Gradientdomain metropolis light transport [J]. ACM Transactions on Graphics， 2013， 32（4）： Article No. 95.

[3]

HACHISUKA T， KAPLANYAN A S， DACHSBACHER C. Multiplexed metropolis light transport [J]. ACM Transactions on Graphics， 2014， 33（4）： Article No. 100.

[4]

ROUSSELLE F， KNAUS C， ZWICHER M. Adaptive rendering with nonlocal means filtering [J]. ACM Transactions on Graphics， 2012， 31（6）： Article No. 195.

[5]

KAPLANYAN A S， DACHSBACHER C. Adaptive progressive photon mapping [J]. ACM Transactions on Graphics， 2013， 32（2）： Article No. 16.

[6]

DELBRACIO M， MUS P， BUADES A， et al. Boosting Monte Carlo rendering by ray histogram fusion [J]. ACM Transactions on Graphics， 2014， 33（1）： Article No. 8.

[7]

LIU X D， ZHENG C W. Adaptive cluster rendering via regression analysis [J]. Visual Computer， 2015， 31（1）： 105-114.

[8]

ZWICKER M， JAROSZ W， LEHTINEN J， et al. Recent advances in adaptive sampling and reconstruction for Monte Carlo rendering [J]. Computer Graphics Forum， 2015， 34（2）： 667-681.

[9]

鮑世强.基于双向路径追踪的 Metropolis 全局光照算法的研究[D].天津：天津大学，2006. （BAO S Q. The study of bidirectional path tracing based metropolis global illumination method [D]. Tianjin： Tianjin University， 2006.）

[10]

KAIJYA J T. The rendering equation [J]. ACM SIGGRAPH Computer Graphics， 1986， 20（4）： 143-150.

[11]

VEACH E， GUIBAS L J. Metropolis light transport [C] //Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques. New York： ACM Press/AddisonWesley Publishing Company， 1997： 65-76.

[12]

PHARR M， HUMPHREYS G. Physically Based Rendering： From Theory to Implementation [M]. San Francisco， CA： Morgan Kaufmann， 2004： 652-660.

[13]

ROUSSELLE F， KNAUS C， ZWICKER M. Adaptive sampling and reconstruction using greedy error minimization [J]. ACM Transactions on Graphics， 2011， 30（6）： Article No. 159.

[14]

LI T M， WU Y T， CHUANG Y Y. SUREbased optimization for adaptive sampling and reconstruction [J]. ACM Transactions on Graphics， 2012， 31（6）： Article No. 194.

[15]

SILVERMAN B W. Density Estimation for Statistics and Data Analysis [M]. Boca Raton， FL： CRC press， 1986：1-22.

[16]

Implementation of pbrtv2 [S/OL]. [20151205]. https：//github.com/mmp/pbrtv2/wiki.

[17]

Scenes and data for rendering [S/OL]. [20151205]. http：//www.pbrt.org.