蒙特卡洛方法详解：从基本原理到光子传输模拟

创作时间:

作者:

@小白创作中心

蒙特卡洛方法详解：从基本原理到光子传输模拟

引用

CSDN

https://blog.csdn.net/qq_39300235/article/details/104565901

蒙特卡洛方法是一种通过随机抽样来解决问题的统计方法，广泛应用于计算函数期望值、积分计算等领域。本文将详细介绍蒙特卡洛方法的基本原理，并通过模拟光子传输的实际应用，帮助读者深入理解这一重要技术。

蒙特卡洛方法简介

在计算机图形学中，我们经常听到“蒙特卡洛”（Monte Carlo，简称MC）这个术语。蒙特卡洛方法本质上是一系列统计方法，用于解决各种问题，如计算函数的期望值或对无法进行解析积分的函数进行数值积分。所有这些方法的核心思想都是使用随机采样。

蒙特卡洛方法主要应用于两大类问题：

模拟：通过多次随机模拟来估计结果。例如，计算从A点到B点的时间时，可以考虑各种不确定因素（如天气、交通等），通过多次模拟获得更准确的估计。
积分：在许多情况下，特别是处理复杂的函数时，蒙特卡洛积分提供了一种有效的数值解法。例如，在基于图像的光照（IBL）中，就需要用到蒙特卡洛积分技术。

蒙特卡洛方法的应用示例

一个经典的蒙特卡洛方法应用是求解不规则图形的面积。具体做法是：在一个已知面积的矩形区域内随机撒点，统计落在目标区域内的点数，然后用这个比例乘以矩形的面积，得到目标区域的近似面积。这种方法虽然简单，但非常有效，且随着样本数量的增加，结果会越来越精确。

模拟光子传输

接下来，我们将通过模拟光子在介质中的传输过程，深入理解蒙特卡洛方法的实际应用。当光子与物质相互作用时，主要有三种可能的结果：被吸收、散射或透射。这些过程可以用吸收系数（σ_a）和散射系数（σ_s）来描述，它们表示单位长度内光子被吸收或散射的概率。

比尔-朗伯定律描述了光子在介质中传输时的衰减规律。根据该定律，光子被吸收或散射的概率可以通过以下公式计算：

其中，σ_t = σ_a + σ_s 是消光系数。这个概率密度函数（PDF）表明，随着距离的增加，光子被消光（吸收或散射）的概率逐渐降低。

在模拟中，我们需要计算光子每次运动的距离x。这可以通过对上述PDF函数求反函数来实现：

为了判断光子是否离开平板，我们可以使用相似三角形原理：

当光子在平板内部运动时，需要判断其是被吸收还是散射。这可以通过比较随机数与光子被吸收的概率（σ_a / σ_t）来实现。为了提高效率，通常采用“光子包”的概念，即每次模拟一个包含多个光子的数据包，通过调整数据包的权重来跟踪光子的吸收情况。

代码实现

下面是一个完整的蒙特卡洛模拟光子传输的代码实现：

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;
                else
                    Rd+=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;
                }
                else
                    w*=m;
            }
            //旋转
            spin(mux, muy, muz, g);
            
        }
    }
}

在上述代码中，我们使用了轮盘模式（roulette）来处理光子包权重过低的情况，以避免过早终止模拟。此外，我们还实现了Henyey-Greenstein相位函数来计算光子的散射方向：

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;
    }
}

通过OpenGL渲染，我们可以直观地看到模拟结果：

在模拟了10000个光子包撞击平板后，反射粒子和透射粒子的分布清晰可见。这个例子充分展示了蒙特卡洛方法在实际应用中的强大能力。