蒙特卡洛(Monte Carlo)方法的介绍和应用

蒙特卡洛(Monte Carlo)方法

在渲染中，我们经常听到术语“蒙特卡洛”（通常缩写为MC）。但是这是什么意思？实际上，它所指的是一个非常简单的想法，蒙特卡洛方法指的是一系列统计方法，这些方法本质上用于查找事物的解决方案，例如计算函数的期望值，或者对由于没有封闭形式而无法进行分析积分的函数进行积分。我们可以用该原理来解决不同的问题，并且针对这些问题中都可以关联不同的技术或算法。所有这些算法的共同点是它们使用随机（或随机）采样。

MC方法都具有使用随机抽取的样本来计算给定问题的解决方案的概念。这些问题通常分为两个主要类别：

模拟：使用蒙特卡洛进行近似模拟，如果我们要计算从A点到B点所花费的时间，在某些情况下，例如在旅途中下雨或下雪的可能性，交通阻塞的可能性，加油等我们将不得不停止的操作。您可以在模拟开始时设置这些条件，并运行模拟1,000次或者更多以获取估计的时间。与往常一样，运行或试验的次数（此处为1,000）越多，我们的的估计就越准确。比如我们想求国家人口高度，从总人口（样本）中随机选择的N个成年人的身高除以数字N（样本数量）。从本质上讲，这就是我们所说的蒙特卡洛近似，其公式为：

我们得到的值则是我们全国平均身高的一个近似解，当样本N越多时得到的解越精确。但要注意的是我们这里的样本必须满足均匀随机。

积分：这听着像数学家们使用的技术，但我们的着色技术会使用到它，例如我们的基于图像的照明(IBL)中对图像的采样得到最后的光线就使用了蒙特卡洛积分技术。随着积分尺寸的增加，此方法可能会变的非常昂贵，因为MC积分的收敛速度非常的慢，但在许多情况，它依然能带来非常好的效果。我们的积分公式为：

我们积分中除了选择均匀分布之外，还可以选择其他有助于减少方差的分布，来更快的收敛。

蒙特卡洛方法应用

MC方法有个非常经典的方法，即用来求解不规则图形的面积：

如上图所示，我们如果想求解黑色区域的面积，我们如果采用均匀随机的撒点，如果点击中黑色区域，则我们为其计数，最后只要将击中数除以总数乘上四边形区域面积就是我们的不规则图形面积：

但注意的是，这里我们计算的面积永远是一个近似值，当发射粒子数量越多，则越接近正确答案。而且该例子我们的采样是要均匀随机的。但是均匀采样并不是MC的绝对条件，比如重要性采样(非均匀采样)，该方法可以减少方差。因此我们可以使用一些特殊的非均匀采样可以让我们有更快的收敛速度。

到这里我们可以慢慢明白MC方法指的是一个概念，即利用随机采样去逼近真实解，我们利用模拟中子运动来全面了解蒙特卡洛方法如何在实际中应用。

模拟光子传输

我们利用程序来模拟光子的传输，其中会在许多方面涉及到蒙特卡洛方法。了解该方法之后，我们再看路径追踪就能很快明白其原理了。

当光子与由某种材料制成的物体相互作用时（假设该物体具有一定的厚度），该光子传播时会有三种情况。如果光子与构成该材料的原子相互作用，则它可以被吸收（光子的能量传递到原子上，然后原子以热的形式返回到其环境），或者可以被散射（它的行进方向被改变）。最后一种情况则是光子根本不与原子相互作用，在这种情况下，光子一直沿相同的方向传播，直到穿过物体。从物体（或层）的一侧传递到另一侧的现象称为透射。这些术语的说明可在下图中找到。

请注意，尽管光子被散射，其方向也会改变。如果光子在原子周围的任何随机方向上散射，则我们称之为各向同性散射（isotropic，即散射在所有方向上都是相等的）。但是有时，光子会在以光子入射方向为中心的给定圆锥方向内发生偏转。在这种情况下，我们说的是各向异性散射（anisotropic，即光子并不是在所有方向上都被散射，而是在以光子入射行进方向为中心的方向）。在各向异性散射的情况下，可以使用我们所谓的相位或角度函数来计算散射光子的新方向。在本章中，我们主要只处理各向同性散射。如下图直观地看到这两种散射模式。

我们的光子若从入射表明离开，则这种光子现象称为反射或者漫反射，我们用Rd表示。但注意，光子不一定从同一点进入或者离开表面，这是由于光子在物体内部发生了散射。若光子穿透物体，则该现象称为透射，可以用T表示，详细图示如下：

在物理学中，可以使用两个量来定义材料吸收和散射光子的属性。它们简称为吸收系数和散射系数。它们通常用希腊字母σ（sigma）表示：为吸收系数，为散射系数。他们测量单位长度的材料（单位长度可以是毫米，厘米等）吸收或散射光子的可能性。换句话说，吸收系数例如表示光子通过材料板时被吸收的可能性。请注意，这些值可以大于1。即使概念非常相似，也不要将它们与概率混淆（概率永远不能大于1）。在物理学中，光子通过材料平板时吸收或散射的速率由一个非常著名的定律-比尔-朗伯定律（Beer-Lambert law）定义。该定律的方程可以采用不同的形式。在本中，我们将使用以下公式：

其中，=+。 该术语称为消光系数。 正如下图所看到的，距离越大，传输率越小，并且传输率随距离降低的速率随指数衰减而变化。 透射率随距离降低的原因是，光子沿途被吸收。 因此，当然，它们越过距离要越过另一侧表面的距离，就越可能吸收它们：

如果我们将称为从入射点到平板中第i个光子被散射或吸收的点的距离，如图：

一些光子在它们进入平板的位置分散很短的距离，而另一些光子则在与原子相互作用之前传播更远的距离。 实际上，该距离（称为自由路径长度）是随机的。 随机变量的概念对于本课程和蒙特卡洛方法的概念至关重要。 光子将被吸收或散射（作为距离的函数）的概率可以通过以下公式计算：

该公式也是我们的概率密度函数，其函数图：

表示了距离X越大的地方被消光(吸收和散射)的概率越低，这个很好理解，当距离(X)越长时，光在之前被吸收或者散射的概率越高。就比如我们在商场买东西，排队越靠后买到东西的概率越低是一样的。我们还需要知道光子被吸收的概率为/，而其散射的概率为/。

我们代码模拟光子经过平板时的步骤如下：

步骤1：计算光子的新位置。沿着方向v运动距离为x。
步骤2：光子是否仍在穿过平板？这可以通过计算光子的位置并将其与平板的边界进行比较轻松地进行测试。
如果光子在平板边界之外，则停止模拟。
如果光子仍在穿过平板，则继续。
步骤3：在此阶段，光子可以被吸收或散射。计算被吸收或散射的光子概率（使用吸收和散射系数）。
如果吸收了光子，则停止模拟（光子消失）
如果光子被散射，则更新光子的行进方向（v），然后返回到步骤1。
该流程图为：

由于我们需要得到每次运动的距离x，所以我们要用上述PDF函数求反函数得出：

则是我们的概率，这里正好可以使用满足[0,1]的均匀随机数，这一步很好的应用了我们的蒙特卡洛方法。但是我们获得x之后每次运动需要判别是否离开平板，我们可以使用简单的相似三角形原理，如下图：

我们可以利用速度获得最后要运动的距离s：

，为距离平板顶部的距离。

我们把速度标准化则有：

但当v方向朝下时，我们的公式为：

当获得s后我们更新位置，这就十分简单了，代码为：

P.x = P.x + x * V.x; 
P.y = P.y + y * V.y; 
P.z = P.z + z * V.z; 

我们的v如何获得会在后面提到。每次运动完我们需要知道它下一次是被吸收还是散射。如果光子被吸收，它将“死亡”。 然后，我们可以更新计数器以跟踪吸收的光子数量，打破循环并继续前进到下一个光子。 我们检查光子是否被吸收的方法是简单地比较均匀分布的随机变量是否低于光子被吸收的概率/。但是，大多数MC模拟使用不同的方法来解决此问题。 实际上，逐个模拟光子是非常繁琐的。 相反，我们可以考虑的是，我们实际上遵循的是一包光子，它们遵循的是同一条路径。 当这些光子通过介质传播时，其中一些会被吸收。 我们跟踪已吸收了多少光子的方法是为光子包分配权重。 该权重最初设置为1。光子被吸收的可能性（/）也可以看作是例如100个光子在与光子相互作用时被吸收的部分（以百分比表示）。 因此，在每次散射事件中，我们要做的就是以该比率降低数据包的重量。

上述所说的代码如下：

void MCSimulation(){int nphotons=10000;//粒子总数int count=0;//没被吸收的粒子(反射和透射的总和)double sigma_a=1,sigma_s=2;//(sigma_a吸收系数，sigma_s散射系数)double sigma_t=sigma_a+sigma_s;double d=0.5,g=0.75;//d为平板厚度，g为散射系数const int m=10;//轮盘参数double Rd=0,Tt=0;//反射和穿透positions.clear();indices.clear();for(int n=0;n<nphotons;n++){double w=1;//包权重double x=0,y=0,z=0,mux=0,muy=0,muz=1;//xyz为粒子坐标，mux,muy,muz为粒子运动单位方向while (w!=0) {double s=-log(drand48())/sigma_t;double distToBoundary=0;//到边界的距离if (muz>0) {distToBoundary=(d-z)/muz;}else if(muz<0){distToBoundary=-z/muz;}if (s>distToBoundary) {//离开平板边界时判断反射还是透射if(muz>0)Tt+=w;elseRd+=w;//记录离开平板的点x+=distToBoundary*mux;y+=distToBoundary*muy;z+=distToBoundary*muz;positions.push_back(x);positions.push_back(y);positions.push_back(z);positions.push_back(w);count++;indices.push_back(count);break;}//若粒子没有离开平板则继续在平板内运动x+=s*mux;y+=s*muy;z+=s*muz;double dw=sigma_a/sigma_t;w-=dw;w=std::max(0.0, w);//轮盘模式if (w<0.001) {if (drand48()>1.0/m) {break;}elsew*=m;}//旋转spin(mux, muy, muz, g);}}
}

positions和indices为OpenGL渲染中传给VBO的缓存，为顶点缓存和索引缓存。我们设定粒子包总数，这里我们每次跟踪粒子时，把它看成一个包，所以在循环中开始的w包权重设为1。在粒子运动时，首先用PDF函数获得移动距离s。然后利用三角相似获得到边界的距离，当s超出到边界的距离时，根据移动方向判断为反射还是透射，并记录粒子的位置。若s没有超出边界，则继续运动。我们根据计算的s和v移动相应的距离然后判断是否被吸收，但是这里我们使用的是光子包，所以我们每次迭代后递减权重w，直到w小于一个阈值(这里我们设为0.001)，则表示所有的光子被吸收了。但是这里还多使用了一个轮盘模式，利用轮盘参数m，在w足够小时，如果直接判定吸收会使结果不准确。因此避免过早的“杀死”这些光子，我们采用轮盘模式。轮盘模式依旧使用到了我们MC方法的原理，利用随机数，当随机数大于1/m时我们才选择丢弃粒子，否则权重w乘上m让粒子继续运动。

最后，如果光子或光子包未被杀死，则进行散射。从代码的角度来看，这意味着我们需要为此数据包计算一个新的方向（更新向量v）。正如我们前面提到的，组成平板的材料可以是各向同性的（光子随后将在所有可能的方向上散射）或各向异性的（光子主要是“分散在以光子传播方向为中心的圆锥形方向内）。可以通过采样所谓的相位函数来计算新方向。计算机图形学中最流行的相位函数之一称为Henyey-Greenstein散射相位函数。您可以在“体积渲染”课程中找到有关此主题的详细说明。此函数仅采用一个通常表示为g的参数。当g等于1时，光子沿与光子入射方向相同的方向散射。值为0表示介质是各向同性的。中间值表示光子大部分散布在以光子传播方向为中心的圆锥形方向内。实际上，相位函数表示光子还可以沿与入射方向相反的光锥居中散射。因此，光子实际上是向前和向后散射的，但是前向散射与后向散射的量取决于参数g。当g大于0时，前向散射占主导，而对于g小于0的情况，后向散射占主导。说明如下图：

散射函数返回作为偏转角的角度θ的余弦，为了对Henyey-Greenstein函数进行采样，我们将使用以下公式：

该公式得到。如上图，我们的方向需要获得两个角和，我们还需要来使向量v进行旋转。旋转公式为：

该公式得到旋转后单位向量v的各个值。该公式是在不等于-1或1的情况下，当时：

当时：

到这里我们的v向量也清楚该如何更新了，我们列出全部代码：

//计算costheta
double getCosFromG(float g){if (g==0) {return 1.0-2.0*drand48();}else{double fT=((1.0-g*g)/(1.0-g+2.0*g*drand48()));return (1.0+g*g-fT*fT)/(2.0*g);}
}//计算旋转量
void spin(double &mu_x,double &mu_y,double &mu_z,const double& g){double costheta=getCosFromG(g);double sintheta=sqrt(1.0-costheta*costheta);double phi=2*PI*drand48();double sinphi=sin(phi);double cosphi=cos(phi);if (mu_z == 1.0) {mu_x = sintheta * cosphi;mu_y = sintheta * sinphi;mu_z = costheta;}else if (mu_z == -1.0) {mu_x = sintheta * cosphi;mu_y = -sintheta * sinphi;mu_z = -costheta;}else {double denom = sqrt(1.0 - mu_z * mu_z);double muzcosphi = mu_z * cosphi;double ux = sintheta * (mu_x * muzcosphi - mu_y * sinphi) / denom + mu_x * costheta;double uy = sintheta * (mu_y * muzcosphi + mu_x * sinphi) / denom + mu_y * costheta;double uz = -denom * sintheta * cosphi + mu_z * costheta;mu_x = ux;mu_y = uy;mu_z = uz;}
}void MCSimulation(){int nphotons=10000;//粒子总数int count=0;//没被吸收的粒子(反射和透射的总和)double sigma_a=1,sigma_s=2;//(sigma_a吸收系数，sigma_s散射系数)double sigma_t=sigma_a+sigma_s;double d=0.5,g=0.75;//d为平板厚度，g为散射系数const int m=10;//轮盘参数double Rd=0,Tt=0;//反射和穿透positions.clear();indices.clear();for(int n=0;n<nphotons;n++){double w=1;//包权重double x=0,y=0,z=0,mux=0,muy=0,muz=1;//xyz为粒子坐标，mux,muy,muz为粒子运动单位方向while (w!=0) {double s=-log(drand48())/sigma_t;double distToBoundary=0;//到边界的距离if (muz>0) {distToBoundary=(d-z)/muz;}else if(muz<0){distToBoundary=-z/muz;}if (s>distToBoundary) {//离开平板边界时判断反射还是透射if(muz>0)Tt+=w;elseRd+=w;//记录离开平板的点x+=distToBoundary*mux;y+=distToBoundary*muy;z+=distToBoundary*muz;positions.push_back(x);positions.push_back(y);positions.push_back(z);positions.push_back(w);count++;indices.push_back(count);break;}//若粒子没有离开平板则继续在平板内运动x+=s*mux;y+=s*muy;z+=s*muz;//轮盘模式double dw=sigma_a/sigma_t;w-=dw;w=std::max(0.0, w);if (w<0.001) {if (drand48()>1.0/m) {break;}elsew*=m;}//旋转spin(mux, muy, muz, g);}}
}

以上代码在OpenGL中简单实现后的效果为：

在模拟玩10000个光子包撞击平板后，我们反射粒子和透射粒子可以清晰的在上图中反应出来。蒙特卡洛方法在该模拟中多处的使用，也让我们到时候理解蒙特卡洛路径追踪有了很好的概念。