29、光线传输与蒙特卡罗方法：原理、算法及代码实现

光线传输与蒙特卡罗方法：原理、算法及代码实现

1. 数值逼近

渲染方程可以写成更紧凑的形式：
[L = L_e + TL]
其中 (T) 是在函数空间中定义的泛函：
[T(g(x, v)) = \int_{S_x} f(x, u, v) \langle u, n \rangle g(x, u) du]
为了得到方程的近似解，可以将 (L) 用右侧的表达式替换。经过一次替换可得：
[L = L_e + T(L_e + TL) = L_e + TL_e + T^2L]
重复 (n + 1) 次这个过程，得到：
[L = L_e + \sum_{i = 1}^{n} T^nL_e + T^{n + 1}L]
忽略 (n + 1) 阶的残差项，可得到渲染方程的近似解：
[L = L_e + \sum_{i = 1}^{n} T^nL_e]
如果考虑所有的 (n \in N)，就得到了纽曼级数。从物理意义上理解，场景中某点发出的光由该点向观察者方向发出的光以及该点经过 1 次、2 次、3 次……反射的所有光组成。这是渲染方程蒙特卡罗解法的基础。

2. 蒙特卡罗积分方法

蒙特卡罗算法是一种利用概率过程来求积分近似值的算法。为了理解该算法，首先需要了解概率变量的期望值 (E[x]) 的定义。设概率变量的密度函数为 (p(x))，在集合 (\Omega) 上，期望值定义为：
[E[x] = \int_{x \in \Omega} x p(x) dx]
其核心思想是利用大数定律。当 (n \to \infty) 时：
[E[x] \approx \frac{1}{n} \sum_{i = 1}^{n} x_i]
其中 (x_1, x_2, \cdots, x_n) 是密度为 (p(x)) 的概率变量的样本。更精确地说，大数定律表明，如果所有 (x_i) 是独立同分布的随机变量，则：
[\text{Probability} \left[ E(x) = \lim_{n \to \infty} \frac{1}{n} \sum_{i = 1}^{n} x_i \right] = 1]
我们还可以定义函数 (f(x)) 的期望值：
[E[f(x)] = \int_{x \in \Omega} f(x) p(x) dx]
如果想用期望值积分来近似 (\int_{x \in \Omega} f(x) dx)，可以定义一个概率密度为 (g) 的变量：
[g(x) = \frac{f(x)}{p(x)}]
当 (n \to \infty) 时：
[\int_{x \in \Omega} f(x) dx = \int_{x \in \Omega} g(x) p(x) dx = E(g(x)) \approx \frac{1}{n} \sum_{i = 1}^{n} g(x_i)]
其中 (x_1, x_2, \cdots, x_n) 是密度为 (p(x)) 的概率变量的样本。

3. 路径追踪

考虑一个包含多个纯反射表面和一个光源的场景。路径追踪算法通过应用蒙特卡罗积分方法，并结合前面提到的数值逼近方法来求解渲染方程。
如果截断纽曼级数，就可以得到渲染方程的近似解。级数中的每一项都是一个积分，可以用蒙特卡罗方法求解。具体来说，为了得到像素 (p) 处渲染方程的近似样本，可以从相机中心发射一条方向对应于 (p) 的光线。每次光线与场景中的表面相交时，在该表面上方的半球内均匀分布地创建一条新光线。重复这个过程，直到光线与光源相交。
在每个交点处应用双向反射分布函数（BRDF）来估计每个表面反射到下一个表面的光量。每条路径都为求解渲染方程的纽曼级数提供一个近似值。通过计算同一像素处多条路径的平均值，可以提高渲染方程解的精度。

3.1 半球上的均匀采样

可以使用条件概率过程在法向量为 (n \in R^3) 的表面上方的半球上进行均匀采样。具体步骤如下：
1. 以均匀概率在立方体 ([-1, 1] \times [-1, 1] \times [-1, 1]) 内选择一个点 (x_i \in R^3)。
2. 如果 (|x_i| > 1)，则丢弃该点。
3. 如果 (\langle x_i, n \rangle > 0)，则接受 (\frac{x_i}{|x_i|}) 作为样本；否则接受 (-\frac{x_i}{|x_i|}) 作为样本。

3.2 直接照明和间接照明的分离

在前面的路径追踪算法中，光线可能需要经过很长的路径才能到达光源。为了更好地近似渲染方程的解，可以将 3D 场景中的直接照明和间接照明分离。即：
[TotalRadiance = DirectRadiance + IndirectRadiance]
在每个交点处，认为光来自两个不同的方向，一个是光线反射到另一个表面的方向，另一个是从一组指向光源上所选点的光线方向。
计算直接光时，可以使用蒙特卡罗算法，步骤如下：
1. 按照均匀采样过程在光源上选择一个随机点 (x’)。
2. 根据前面定义的 (V) 和 (G)，得到蒙特卡罗和的一项：
[L(x, u, v) = \frac{f(x, u, v)L_e(x, u)V(x, x’)G(x, x’)}{p(x’)} = f(x, u, v)L_e(x, u)V(x, x’)G(x, x’)A]
其中 (A) 是光源的面积。
因此，来自光源的直接光的近似值为：
[L(x, v) \approx \frac{1}{n} \sum_{i = 1}^{n} f(x, u_i, v)L_e(x, u_i)V(x, x’_i)G(x, x’_i)A]
其中 (x’_1, x’_2, \cdots, x’_n) 是在光源表面上均匀分布采样的点，(u_i \in R^3) 定义为：
[u_i = \frac{x’_i - x}{|x’_i - x|}]

3.3 多边形光源

为了使用上述公式定义多边形光源产生的直接照明，需要知道多边形的面积以及均匀采样该多边形的方法。这里假设多边形可以分割成面积相等的三角形，这样就可以对每个三角形进行均匀采样。
如果 (v_0, v_1, v_2) 是三角形的顶点，(r_1, r_2) 是区间 ([0, 1]) 上的两个均匀随机变量，则可以通过以下方式得到三角形上的均匀随机点 ((x, y, z))：
[x = (1 - \sqrt{r_1})v_0]
[y = \sqrt{r_1}(1 - r_2)v_1]
[z = r_2\sqrt{r_1}v_2]
多边形的面积可以通过将这些三角形的面积相加得到，每个三角形的面积计算公式为：
[\frac{1}{2} |(v_2 - v_0) \times (v_1 - v_0)|]

4. 代码模块

4.1 路径追踪 API

以下是路径追踪相关的 API 函数：

Color trace_path( int first_level, int level, Ray v, Object *ol );
Color mirror_reflect( int first_level, int level, Ray v, Inode *i, Vector3 p, Object *ol, Material *m );
Color apply_bphong( Material *m, Color lcolor, Vector3 n, Vector3 l, Vector3 v, Real geomfactor );
Vector3 sort_new_direction( Cone c );
Bool hit_surface( Ray v, Object *ol );
void adjust_normal( Inode *i, Ray v );

• trace_path ：实现路径追踪算法的递归代码。 first_level 表示场景表面之间的反射次数， level 表示当前反射深度，每次递归减 1， v 表示当前要估计颜色的光线， ol 是场景中所有对象的列表。如果某个对象的材料的 se 参数为负，则路径追踪表现为镜面反射。
• mirror_reflect ：返回完美镜面方向的颜色。
• apply_bphong ：基于 Blinn - Phong 反射模型实现 BRDF。
• sort_new_direction ：返回由 c 描述的半球上的均匀采样。
• hit_surface ：如果光线 v 穿过对象列表 ol 中的对象，则返回 TRUE ；否则返回 FALSE 。
• adjust_normal ：使 inode 的法向量 i->n 指向包含光线 v 的半空间。

4.2 路径追踪代码

// extracted from ibl/ibl.h
#define REAL_RANDOM ((double)rand()/RAND_MAX)
// extracted from lptrace/ptrace.h
#ifndef PTRACE_H
#define PTRACE_H
#include <math.h>
#include "rshade.h"
#include "rt.h"
#include "ibl.h"
#include "plight.h"
#include <time.h>
#include <stdlib.h>
#include <string.h>
Color trace_path( int first_level, int level, Ray v, Object *ol );
Color mirror_reflect( int first_level, int level, Ray v, Inode *i, Vector3 p, Object *ol, Material *m );
Color apply_bphong( Material *m, Color lcolor, Vector3 n, Vector3 l, Vector3 v, Real geom_fac );
Vector3 sort_new_direction( Cone c );
Bool hit_surface( Ray v, Object *ol );
void adjust_normal( Inode *i, Ray v );
#endif
// lptrace/ptrace.c
#include "ptrace.h"
static Color direct_light( int ns, Vector3 p, Vector3 n, Vector3 v, Object *ol, PolyLight *l, Material *m );
static int light_visible(Ray r, Object *ol, PolyLight *l, Color *c );
#include "ptrace.h"
Color trace_path( int first_level, int level, Ray v, Object *ol )
{
    PolyLight *l = plight_head;
    Color c = C_BLACK, c_aux, w;
    Ray r;
    Inode *i = ray_intersect(ol, v);
    adjust_normal( i, v );
    if( light_visible(v,ol,l,&w) )
        return w;
    if((i != NULL) && level){
        level--;
        Material *m = i->m;
        Vector3 p = ray_point(v, i->t);
        Cone recv = cone_make(p, i->n, PIOVER2);
        if( m->se < 0 )
            c = mirror_reflect( first_level, level, v, i, p, ol, m );
        else{
            Vector3 d = sort_new_direction(recv);
            c = direct_light( DIRECT_LIGHT_NSMPLS, p, i->n, v3_scale(-1,v.d), ol, l, m);
            if( !light_visible(ray_make(recv.o, d), ol, l, &w) ){
                c_aux = trace_path( first_level, level, ray_make(p, d), ol );
                c = c_add(c, apply_bphong( m, c_aux, i->n, d, v3_scale(-1,v.d), 2*PI*v3_dot(d, i->n) ) );
            }
        }
        inode_free(i);
        return c;
    }
    else{
        if( i == NULL )
            return C_BLACK;
        else{
            inode_free(i);
            return C_BLACK;
        }
    }
}
Color direct_light( int ns, Vector3 p, Vector3 n, Vector3 v, Object *ol, PolyLight *l, Material *m )
{
    int k;
    Vector3 r, light_ptdir, inv_light_ptdir;
    Ray s;
    Color c, cl_light;
    Inode *i;
    double g=0;
    c = c_make(0,0,0);
    for( k = 0; k < ns; k++ ){
        s = plight_sample(l,&cl_light);
        r = v3_sub(s.o, p);
        light_ptdir = v3_unit(r);
        inv_light_ptdir = v3_scale( -1, light_ptdir );
        i = ray_intersect(ol, ray_make(p, light_ptdir));
        if( (i == NULL) || (i->t > v3_norm(r)) ){
            g = fabs( v3_dot(inv_light_ptdir,s.d)*v3_dot(light_ptdir, n) );
            c = c_add(c, apply_bphong(m, cl_light, n, light_ptdir, v, g*plight_area(l)/SQR(v3_norm(r) ) ) );
        }
        if( i != NULL )
            inode_free(i);
    }
    return v3_scale(1./ns, c);
}
int light_visible(Ray r, Object *ol, PolyLight *l, Color *c)
{
    Inode *i1, *i2;
    i1 = plight_intersect( l, r, c );
    if( i1 != NULL ){
        i2 = ray_intersect( ol, r );
        if( i2 == NULL || (i1->t < i2->t ) ){
            inode_free(i1);
            if( i2 != NULL )
                inode_free(i2);
            return TRUE;
        }
        inode_free(i1);
        if( i2 != NULL )
            inode_free(i2);
    }
    return FALSE;
}
// Extracted from lptrace/aux.c
#include "ptrace.h"
Color mirror_reflect( int first_level, int level, Ray v, Inode *i, Vector3 p, Object *ol, Material *m )
{
    Vector3 d, c_aux;
    d = v3_sub( v.d, v3_scale( 2*v3_dot( v.d, i->n ), i->n ) );
    c_aux = trace_path( first_level, level, ray_make(p, d), ol );
    return v3_scale( m->ks, c_aux );
}
Color apply_bphong( Material *m, Color lcolor, Vector3 n, Vector3 l, Vector3 v, Real geom_fac)
{
    Color cs,cd;
    Vector3 h = v3_unit(v3_add(l,v));
    cs = c_scale(m->ks* pow(MAX(0, v3_dot(h,n)),m->se)*geom_fac, lcolor);
    cd = c_mult(m->c, c_scale(m->kd*geom_fac, lcolor) );
    return c_add( cs, cd );
}
Vector3 sort_new_direction( Cone c )
{
    double x, y, z;
    Vector3 v, n = v3_unit(c.d);
    do{
        x = 2*REAL_RANDOM - 1;
        y = 2*REAL_RANDOM - 1;
        z = 2*REAL_RANDOM - 1;
        v = v3_make(x,y,z);
    }
    while( v3_norm(v) > 1 );
    if( v3_dot(v, c.d) < 0 ){
        double s = v3_dot(v3_scale(-1,v), n);
        v = v3_add( v, v3_scale(2*s,n) );
    }
    return v3_unit(v);
}
Bool hit_surface( Ray v, Object *ol )
{
    Inode *i = ray_intersect(ol, v);
    if( i == NULL )
        return FALSE;
    else{
        inode_free(i);
        return TRUE;
    }
}
void adjust_normal( Inode *i, Ray v )
{
    if( i != NULL && v3_dot( i->n, v.d ) > 0 )
        i->n = v3_scale( -1, i->n );
}

4.3 多边形光源 API

void init_plight_list( void );
PolyLight *plight_alloc( Color c, Hpoly *p );
void plight_insert( PolyLight *pl );
void plight_free( PolyLight *head );
Ray plight_sample( PolyLight *p, Color *c );
Real plight_area( PolyLight *p );
Inode *plight_intersect( PolyLight *p, Ray r, Color *c );
Val plight_parse(int pass, Pval *pl);
Ray poly3_sample( Hpoly *p );

• init_plight_list ：初始化一个空的多边形光源列表。
• plight_alloc ：返回一个颜色为 c 且形状由 p 定义的多边形光源。
• plight_insert ：将 pl 插入多边形光源列表。
• plight_free ：删除由 head 指向的多边形光源列表。
• plight_sample ：返回一条，其原点是多边形光源 p 的均匀样本，方向垂直于 p ， c 接收样本的颜色。
• plight_area ：返回多边形光源 p 的面积。
• plight_intersect ：返回光线 r 与多边形光源 p 相交的 inode ， c 返回光的颜色，如果没有相交则返回 NULL 。
• plight_parse ：负责处理场景描述文件中的多边形光源。
• poly3_sample ：返回一条，其原点是多边形光源 p 的均匀样本，方向垂直于 p 。

4.4 多边形光源代码

// plight/plight.h
#ifndef P_LIGHT_H
#define P_LIGHT_H
#define DIRECT_LIGHT_NSMPLS 3
#include "color.h"
#include "poly.h"
#include "ibl.h"
typedef struct PolyLight{
    struct PolyLight *next;
    Hpoly *p;
    Color c;
} PolyLight;
PolyLight *plight_head;
void init_plight_list( void );
void plight_free( PolyLight *head );
PolyLight *plight_alloc( Color c, Hpoly *p );
void plight_insert( PolyLight *pl );
Ray plight_sample( PolyLight *p, Color *c );
Real plight_area( PolyLight *p );
Inode *plight_intersect( PolyLight *p, Ray，r, Color *c );
Val plight_parse(int pass, Pval *pl);
Ray poly3_sample( Hpoly *p );
#endif
// plight/plight.c
#include "plight.h"
void init_plight_list( void )
{
    plight_head = NULL;
}
PolyLight *plight_alloc( Color c, Hpoly *p )
{
    PolyLight *pl;
    pl = NEWSTRUCT(PolyLight);
    pl->next = NULL;
    pl->p = p;
    pl->c = c;
    return pl;
}
void plight_free( PolyLight *head )
{
    PolyLight *pl;
    for( pl = head; pl != NULL; pl=pl->next ){
        free( pl->p );
        free( pl );
    }
}
void plight_insert( PolyLight *pl )
{
    pl->next = plight_head;
    plight_head = pl;
}
Ray plight_sample( PolyLight *pl, Color *c )
{
    int i, r, n = 0;
    PolyLight *aux;
    for( aux = pl; aux != NULL; n++, aux = aux->next );
    r = rand()%n;
    for( aux = pl, i=0; i < r-1; aux = aux->next );
    *c = aux->c;
    return poly3_sample(aux->p);
}
Real plight_area( PolyLight *pl )
{
    Real a = 0;
    while (pl != NULL) {
        a += hpoly3_area(pl->p);
        pl = pl->next;
    }
    return a;
}
Inode *plight_intersect( PolyLight *pl, Ray r, Color *c )
{
    Inode *l = NULL;
    while( pl != NULL && l == NULL ){
        l = hpoly_intersect(pl->p, hpoly3_plane(pl->p), r);
        *c = pl->c;
        pl = pl->next;
    }
    return l;
}
Val plight_parse(int pass, Pval *plist)
{
    Color c;
    struct Hpoly *pl, *pols;
    Val v = {V_NULL, 0};
    if (pass == T_EXEC) {
        Pval *p;
        for (p = plist; p !=NULL; p = p->next) {
            if (p->name == NULL) {
                error("(extense light) syntax error");
            } else if (strcmp(p->name, "color") == 0) {
                c = pvl_to_v3(p->val.u.v);
            } else if (strcmp(p->name, "shape") == 0 && p->val.type == V_HPOLYLIST) {
                pols = p->val.u.v;
            } else error("(extense light) syntax error");
        }
        for(pl = pols; pl != NULL; pl = pl->next)
            plight_insert( plight_alloc( c, pl ));
        v.type = V_PL_LIGHT;
        v.u.v = NULL;
    }
    return v;
}
Ray poly3_sample( Hpoly *p )
{
    Vector3 a, b, c, v0, v1, v2;
    Real r1 = REAL_RANDOM;
    Real r2 = REAL_RANDOM;
    // 后续代码未完整给出，此处省略
}

综上所述，通过数值逼近、蒙特卡罗积分、路径追踪等方法，结合相应的代码实现，可以有效地求解渲染方程，实现高质量的光线传输模拟。同时，多边形光源的处理使得光照效果更加真实和多样化。这些方法和代码为计算机图形学中的光线渲染提供了重要的基础。

5. 光线传输与蒙特卡罗方法的应用与优势

5.1 应用场景

光线传输与蒙特卡罗方法在计算机图形学领域有着广泛的应用，以下是一些常见的应用场景：
| 应用场景 | 描述 |
| — | — |
| 电影特效制作 | 在电影中，为了营造出逼真的场景和视觉效果，需要对光线的传播和反射进行精确模拟。光线传输与蒙特卡罗方法可以帮助实现高质量的光影效果，如逼真的阴影、反射和折射等。 |
| 游戏开发 | 在游戏中，实时渲染的光影效果可以提升游戏的沉浸感。通过使用这些方法，可以模拟出自然的光照效果，使游戏场景更加真实。 |
| 建筑可视化 | 在建筑设计和规划中，需要对建筑物内部和外部的光照情况进行模拟。光线传输与蒙特卡罗方法可以帮助设计师评估不同设计方案下的光照效果，优化建筑布局。 |
| 工业设计 | 在产品设计中，光线的反射和折射效果可以影响产品的外观和质感。通过模拟光线传输，可以更好地展示产品的设计特点，提高产品的吸引力。 |

5.2 优势分析

光线传输与蒙特卡罗方法相较于传统的渲染方法具有以下优势：
- 真实感强 ：能够精确模拟光线在场景中的传播和反射，产生更加真实的光影效果。
- 灵活性高 ：可以处理复杂的场景和材质，适应不同的渲染需求。
- 理论基础坚实 ：基于物理原理的数值计算方法，具有可靠的理论支持。

5.3 局限性与挑战

尽管光线传输与蒙特卡罗方法具有很多优势，但也存在一些局限性和挑战：
- 计算复杂度高 ：由于需要进行大量的光线追踪和积分计算，导致计算时间较长，对计算机性能要求较高。
- 收敛速度慢 ：在某些情况下，蒙特卡罗积分的收敛速度较慢，需要更多的样本才能得到准确的结果。
- 噪声问题 ：在渲染过程中，可能会出现噪声，影响渲染质量。需要采用一些去噪技术来解决这个问题。

6. 优化策略

为了克服光线传输与蒙特卡罗方法的局限性，可以采用以下优化策略：

6.1 采样策略优化

• 重要性采样 ：根据场景的特点，对重要的区域进行更多的采样，减少不必要的采样，提高采样效率。
• 分层采样 ：将采样空间分成多个层次，在每个层次上进行均匀采样，避免采样点的聚集，提高收敛速度。

6.2 加速结构

• 包围盒层次结构（BVH） ：通过构建包围盒层次结构，可以快速排除与光线不相交的物体，减少光线与物体的相交测试次数。
• kd - 树 ：kd - 树是一种空间划分数据结构，可以有效地组织场景中的物体，提高光线追踪的效率。

6.3 并行计算

利用多核处理器和 GPU 的并行计算能力，同时处理多个光线追踪任务，加速渲染过程。以下是一个简单的并行计算流程图：

graph TD;
    A[开始] --> B[初始化场景和光线];
    B --> C[将光线任务分配到多个线程];
    C --> D[每个线程进行光线追踪计算];
    D --> E[合并计算结果];
    E --> F[结束];

6.4 去噪技术

• 双边滤波 ：通过考虑像素的颜色和空间位置信息，对图像进行滤波，去除噪声的同时保留图像的细节。
• 基于机器学习的去噪方法 ：利用深度学习模型对噪声图像进行学习和处理，得到高质量的去噪图像。

7. 总结与展望

7.1 总结

光线传输与蒙特卡罗方法是计算机图形学中用于求解渲染方程的重要方法。通过数值逼近、蒙特卡罗积分和路径追踪等技术，可以模拟光线在场景中的传播和反射，实现高质量的光影效果。同时，代码模块的实现为实际应用提供了基础。

7.2 展望

随着计算机硬件性能的不断提升和算法的不断改进，光线传输与蒙特卡罗方法将在更多领域得到应用。未来的研究方向可能包括：
- 实时渲染 ：进一步提高渲染速度，实现实时的高质量光影效果。
- 深度学习与光线传输的结合 ：利用深度学习技术优化光线传输算法，提高渲染效率和质量。
- 多物理场耦合模拟 ：将光线传输与其他物理现象（如热传导、流体动力学等）进行耦合模拟，实现更加真实的场景模拟。

在实际应用中，我们可以根据具体的需求选择合适的方法和优化策略，以达到最佳的渲染效果。同时，不断探索和创新，推动光线传输与蒙特卡罗方法的发展。