[图形学]smallpt代码详解（上）

一、简介

本文介绍了著名的99行代码实现全局光照的光线跟踪代码smallpt。
包括对smallpt的功能介绍、编译运行介绍，和对代码的详细解释。希望能够帮助读者更进一步的理解光线跟踪。

二、smallpt介绍

1.smallpt是什么

smallpt(small Path Tracing) 是一个全局光照渲染器。smallpt只有99行C++代码组成，代码开源，并使用 无偏蒙特卡罗路径跟踪 (unbiased Monte Carlo path tracing)，可以渲染得到以下场景结果：

无偏蒙特卡罗路径跟踪 unbiased Monte Carlo path tracing：
无偏：指可以确保计算出的平均值（期望结果）与实际光照结果一致，不存在系统性的偏差。算法中，随着采样次数增加，最终图像会逐渐逼近真实的物理模拟结果。

渲染结果

2.smallpt支持的功能（特点）

1).基于无偏蒙特卡洛路径跟踪的全局光照（global illumination）；
2).只有99行（每行最大77列字符）代码的开源C++代码；
3).支持使用OpenMP进行并行加速；
4).使用漫反射产生软阴影（soft shadows）；
5).支持镜面反射（specular）、漫反射（diffuse）和玻璃材质（glass）的BRDFs；
6).使用在2*2子像素上的重要性采样+超采样的抗锯齿算法（antialiasing）；
7).支持球面与光线求交计算；
8).实现了对类康奈尔盒（modified Cornell box）场景的渲染；
9).实现了在漫反射材质上半球面上的余弦重要性采样；
10).使用俄罗斯轮盘赌用于确定是否结束光线的传播；
11).使用俄罗斯轮盘赌和光线分裂由于计算支持镜面反射和折射的玻璃材质；
12).编译后的可执行文件小于 4K；

在本文对代码的详细讲解部分也涉及了对上述部分功能特点的讲解。

3.如何运行smallpt

1).smallpt源代码

smallpt.cpp：

#include <math.h>   // smallpt, a Path Tracer by Kevin Beason, 2008 
#include <stdlib.h> // Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt 
#include <stdio.h>  //        Remove "-fopenmp" for g++ version < 4.2 
struct Vec {        // Usage: time ./smallpt 5000 && xv image.ppm 
  double x, y, z;                  // position, also color (r,g,b) 
  Vec(double x_=0, double y_=0, double z_=0){ x=x_; y=y_; z=z_; } 
  Vec operator+(const Vec &b) const { return Vec(x+b.x,y+b.y,z+b.z); } 
  Vec operator-(const Vec &b) const { return Vec(x-b.x,y-b.y,z-b.z); } 
  Vec operator*(double b) const { return Vec(x*b,y*b,z*b); } 
  Vec mult(const Vec &b) const { return Vec(x*b.x,y*b.y,z*b.z); } 
  Vec& norm(){ return *this = *this * (1/sqrt(x*x+y*y+z*z)); } 
  double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross: 
  Vec operator%(Vec&b){return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);} 
}; 
struct Ray { Vec o, d; Ray(Vec o_, Vec d_) : o(o_), d(d_) {} }; 
enum Refl_t { DIFF, SPEC, REFR };  // material types, used in radiance() 
struct Sphere { 
  double rad;       // radius 
  Vec p, e, c;      // position, emission, color 
  Refl_t refl;      // reflection type (DIFFuse, SPECular, REFRactive) 
  Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_): 
    rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {} 
  double intersect(const Ray &r) const { // returns distance, 0 if nohit 
    Vec op = p-r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0 
    double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad; 
    if (det<0) return 0; else det=sqrt(det); 
    return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0); 
  } 
}; 
Sphere spheres[] = {//Scene: radius, position, emission, color, material 
  Sphere(1e5, Vec( 1e5+1,40.8,81.6), Vec(),Vec(.75,.25,.25),DIFF),//Left 
  Sphere(1e5, Vec(-1e5+99,40.8,81.6),Vec(),Vec(.25,.25,.75),DIFF),//Rght 
  Sphere(1e5, Vec(50,40.8, 1e5),     Vec(),Vec(.75,.75,.75),DIFF),//Back 
  Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(),Vec(),           DIFF),//Frnt 
  Sphere(1e5, Vec(50, 1e5, 81.6),    Vec(),Vec(.75,.75,.75),DIFF),//Botm 
  Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(),Vec(.75,.75,.75),DIFF),//Top 
  Sphere(16.5,Vec(27,16.5,47),       Vec(),Vec(1,1,1)*.999, SPEC),//Mirr 
  Sphere(16.5,Vec(73,16.5,78),       Vec(),Vec(1,1,1)*.999, REFR),//Glas 
  Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12),  Vec(), DIFF) //Lite 
}; 
inline double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; } 
inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); } 
inline bool intersect(const Ray &r, double &t, int &id){ 
  double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20; 
  for(int i=int(n);i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;} 
  return t<inf; 
} 
Vec radiance(const Ray &r, int depth, unsigned short *Xi){ 
  double t;                               // distance to intersection 
  int id=0;                               // id of intersected object 
  if (!intersect(r, t, id)) return Vec(); // if miss, return black 
  const Sphere &obj = spheres[id];        // the hit object 
  Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c; 
  double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl 
  if (++depth>5) if (erand48(Xi)<p) f=f*(1/p); else return obj.e; //R.R. 
  if (obj.refl == DIFF){                  // Ideal DIFFUSE reflection 
    double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2); 
    Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u; 
    Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm(); 
    return obj.e + f.mult(radiance(Ray(x,d),depth,Xi)); 
  } else if (obj.refl == SPEC)            // Ideal SPECULAR reflection 
    return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi)); 
  Ray reflRay(x, r.d-n*2*n.dot(r.d));     // Ideal dielectric REFRACTION 
  bool into = n.dot(nl)>0;                // Ray from outside going in? 
  double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t; 
  if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)    // Total internal reflection 
    return obj.e + f.mult(radiance(reflRay,depth,Xi)); 
  Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm(); 
  double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n)); 
  double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P); 
  return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ?   // Russian roulette 
    radiance(reflRay,depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) : 
    radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr); 
} 
int main(int argc, char *argv[]){ 
  int w=1024, h=768, samps = argc==2 ? atoi(argv[1])/4 : 1; // # samples 
  Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir 
  Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h]; 
#pragma omp parallel for schedule(dynamic, 1) private(r)       // OpenMP 
  for (int y=0; y<h; y++){                       // Loop over image rows 
    fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps*4,100.*y/(h-1)); 
    for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w; x++)   // Loop cols 
      for (int sy=0, i=(h-y-1)*w+x; sy<2; sy++)     // 2x2 subpixel rows 
        for (int sx=0; sx<2; sx++, r=Vec()){        // 2x2 subpixel cols 
          for (int s=0; s<samps; s++){ 
            double r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1); 
            double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2); 
            Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) + 
                    cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d; 
            r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps); 
          } // Camera rays are pushed ^^^^^ forward to start in interior 
          c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25; 
        } 
  } 
  FILE *f = fopen("image.ppm", "w");         // Write image to PPM file. 
  fprintf(f, "P3\n%d %d\n%d\n", w, h, 255); 
  for (int i=0; i<w*h; i++) 
    fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z)); 
}

2).使用g++直接编译运行

可以使用g++编译器直接编译smallpt，使用以下命令编译smallpt：

g++ -O3 -fopenmp smallpt.cpp -o smallpt

编译成功后同目录下会生成一个smallpt可执行文件，然后运行以下命令可以得到最终的渲染结果：

./smallpt

也可以运行以下命令指定渲染中的每个像素上的采样光线数spp（sample per pixel），例如：

./smallpt 128

运行结束后会在smallpt可执行文件的同目录下生成一个名为image.ppm的图片文件，该文件即为渲染结果。

3).使用CMake编译运行

也可以使用CMake编译smallpt，CMakeLists.txt代码如下：

cmake_minimum_required(VERSION 3.5.0)
project(smallpt VERSION 0.1.0 LANGUAGES C CXX)

# 搜索 OpenMP 库
find_package(OpenMP REQUIRED)
add_executable(smallpt smallpt.cpp)
# 链接 OpenMP 库
if(OpenMP_CXX_FOUND)
    target_link_libraries(smallpt PUBLIC OpenMP::OpenMP_CXX)
endif()

将CMakeList.txt放到smallpt.cpp同目录下，然后在该目录下运行以下命令编译、运行smallpt：

mkdir build
cd build
cmake ..
make 
./smallpt

运行结束后会在smallpt可执行文件的同目录下生成一个名为image.ppm的图片文件，该文件即为渲染结果。

三、smallpt代码详解

0.代码结构

本文将smallpt的代码分为自定义数据结构（第4到第22行）、自定义函数（第23到47行）、递归路径跟踪（第48到第74行）和main函数（第75到99行）四部分。
其中：

1).自定义数据结构（第4到第22行）包含代码中用于表示坐标(x,y,z)和颜色值(r,g,b)的Vec类（尽管代码中使用的是struct结构体，但为了便于理解本文依旧将其称为class类），用于描述光线的Ray类，用于描述物体表面材质类型（BRDF，Bidirectional Reflectance Distribution Function）的Refl_t枚举结构，和用于描述球面对象的Sphere类；
2).自定义函数（第23到47行）包含用于防止数据溢出、避免产生非法值的 clamp() 函数、用于类型转换的 toInt() 函数，和用于求光线与球面物体相交的 intersect() 函数；
3).递归路径跟踪（第48到第74行）部分包含路径跟踪算法中关键的，基于递归的场景渲染函数radiance() ；
4).main函数（第75到99行）里包含对场景的相机camera的描述、在每个像素上进行光线采样、使用多重采样进行抗锯齿（antialiasing）和将渲染结果保存为.ppm文件部分功能的代码；

下面详细介绍各部分的功能和代码细节。

1.自定义数据结构部分（第4到第22行）

1).Vec类

Vec(double x_=0, double y_=0, double z_=0){ x=x_; y=y_; z=z_; } 
Vec operator+(const Vec &b) const { return Vec(x+b.x,y+b.y,z+b.z); } 
Vec operator-(const Vec &b) const { return Vec(x-b.x,y-b.y,z-b.z); } 
Vec operator*(double b) const { return Vec(x*b,y*b,z*b); } 
Vec mult(const Vec &b) const { return Vec(x*b.x,y*b.y,z*b.z); } 
Vec& norm(){ return *this = *this * (1/sqrt(x*x+y*y+z*z)); } 
double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross: 
Vec operator%(Vec&b){return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);} 
};

Vec类是一个三维向量，既可以用于描述3D空间中点的位置(x,y,z)，也可以用于描述颜色(r,g,b)。
同时代码中给出了几种用于三维向量的几种运算符号，包括：

+，两向量相加，两向量的各分量分别相加：
$x_{1},y_{1},z_{1})+(x_{2},y_{2},z_{2})=(x_{1}+x_{2},y_{1}+y_{2},z_{1}+z_{2})$
-，两向量相减，两向量的各分量分别相减：
$x_{1},y_{1},z_{1})-(x_{2},y_{2},z_{2})=(x_{1}-x_{2},y_{1}-y_{2},z_{1}-z_{2})$
*，向量乘以标量，向量的各分量分别乘以标量：
$x_{1},y_{1},z_{1})*c=(c*x_{1},c*y_{1},c*z_{1})$
mult，两向量相乘，两向量的各分量分别相乘：
$x_{1},y_{1},z_{1})\ mult\ (x_{2},y_{2},z_{2})=(x_{1}*x_{2},y_{1}*y_{2},z_{1}*z_{2})$
norm()，将向量进行归一化：
$norm((x_{1},y_{1},z_{1})) = (x_{1},y_{1},z_{1})/length((x_{1},y_{1},z_{1}))$
dot，两向量进行点乘：
$x_{1},y_{1},z_{1})\ dot\ (x_{2},y_{2},z_{2})=x_{1}*x_{2}+y_{1}*y_{2}+z_{1}*z_{2}$
%，两向量进行叉乘：
$x_{1},y_{1},z_{1})\ \% (x_{2},y_{2},z_{2})=(y_{1}*z_{2}-z_{1}*y_{2},z_{1}*x_{2}-x_{1}*z_{2},x_{1}*y_{2}-y_{1}*x_{2})$

2).Ray类

struct Ray { Vec o, d; Ray(Vec o_, Vec d_) : o(o_), d(d_) {} };

Ray类用于描述一根光线（射线），包含两个变量：

Vec o：光线的起始点位置；
Vec d：光线的方向；

一条光线可以用以下公式表示：
$\tag{1}$
其中 $t$ 表示光线的长度参数， $O$ 即光线的起始点位置， $D$ 为光线的方向。

3).Refl_t枚举结构

enum Refl_t { DIFF, SPEC, REFR };  // material types, used in radiance()

Refl_t枚举类型 用于描述场景中物体表面的材质类型（BRDF），smallpt中包含三类材质：

DIFF，diffuse，漫反射材质。漫反射材质表面可以向表面法向所在的半球范围内均匀的反射入射光。
SPEC，specular，镜面反射材质。镜面反射材质表面可以认为是理想的镜面，只向理想反射光方向反射入射光。
REFR，refractor，折射材质。折射材质表面会令入射光发生反射和折射两种行为。若发生反射，入射光会从理想的反射光方向反射出去。若发生折射，光线会进入物体内部，折射过程中光线的传播路线遵循折射定律（斯涅尔定律，Snell’Law）：
$n_{a}*sin(\theta_{a}) = n_{b}*sin(\theta_{b}) \tag{2}$
其中：
$n_{a}$ 和 $n_{b}$ 是两个介质的折射率；
$\theta_{a}$ 和 $\theta_{b}$ 分别是入射角和折射角，如下图所示：

上图中 $D$ 为入射光线， $N$ 为点 $P$ 的法向， $B$ 为点 $P$ 处垂直法向 $N$ 的单位向量， $T$ 为折射光线。
折射材质表面发生反射和折射的比例，反射率 $R$ ，与入射角度 $\theta_{a}$ 、折射角 $\theta_{t}$ 、两种材质的折射率 $n_{a}$ 和 $n_{b}$ 相关，可以根据菲涅尔方程（Fresnel）计算：
$\frac{1}{2} \left( \frac{(n_a \cos \theta_a - n_b \cos \theta_b)^2}{(n_a \cos \theta_a + n_b \cos \theta_b)^2} + \frac{(n_b \cos \theta_a - n_a \cos \theta_b)^2}{(n_b \cos \theta_a + n_a \cos \theta_b)^2} \right) \tag{3}$
反射比率 $R$ 表示的是，入射光线中有 $R*100\%$ 的光被反射，有 $(1-R)*100\%$ 的光被折射。
但是由于菲涅尔方程需要大量计算，在实际的计算机图形学中常使用 Schlick近似来简化反射率的计算，Schlick近似公式如下：
$n=\frac{n_a}{n_b} \\ R_{0}=\frac{(n-1)^2}{(n+1)^2} \\ R_{r}(\theta_{a}) = R_{o} + (1-R_{o})(1-cos(\theta_{a}))^{5} \tag{4}$
其中 $R_{0}$ 为光线垂直入射（即入射角 $\theta_{a}=0$ ）时的反射率， $R_{r}(\theta_{a})$ 为入射角为 $\theta_{a}$ 时的反射率。

4).Sphere类

struct Sphere { 
double rad;       // radius 
Vec p, e, c;      // position, emission, color 
Refl_t refl;      // reflection type (DIFFuse, SPECular, REFRactive) 
Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_): 
    rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {} 
double intersect(const Ray &r) const { // returns distance, 0 if nohit 
    Vec op = p-r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0 
    double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad; 
    if (det<0) return 0; else det=sqrt(det); 
    return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0); 
} 
};

Sphere类 用于描述球面物体，包含五个变量和一个用于计算光线与该球面进行相交检测的函数intersect()：

double rad：球面的半径；
Vec p：球面球心的坐标位置；
Vec e：球面的自发光颜色；
Vec c：球面的反射颜色；
Refl_t refl：球面的材质类型；
double intersect(const Ray&r) const：计算光线r是否与该球面相交，若相交返回交点处到光线原点的距离。

下面详细讲解intersect()函数的原理：
前面讲过，一条光线可以用以下公式表示：
$\tag{5}$
类似的，一个球面可以用以下公式表示：
$r^2=0 \tag{6}$
其中 $P$ 为球面上的点集， $C$ 为球心坐标， $r$ 为球半径。
将公式(5)代入公式(6)可以得到：
$r^2=0 \\ (O+tD-C)*(O+tD-C)-r^2=0 \\ (D*D)t^2+2D*(O-C)t+(O-C)*(O-C)-r^2=0 \tag{7}$
公式(7)是一个自变量为 $t$ 的一元二次函数，公式(7)中各项可以写成：
$a'*t^2+b'*t+c'=0 \\ a'=(D*D)=1\\ b'=2D(O-C)=-2D(C-O) \\ c'=(O-C)*(O-C)-r^2=(C-O)*(C-O)-r^2 \tag{8}$
因此：
$t=\frac{-b' \pm \sqrt{(b'^2-4a'c')}}{2a'}=\frac{-b' \pm \sqrt{(b'^2-4c')}}{2}$
如果光线 $P (t)$ 与球面没有交点，那么公式(7)无实数解；如果有一个交点则公式(7)有一个实数解；如果有两个交点则有两个实数解。

因此intersect()函数即进行求解未知量 $t$ ：

// op 为公式(8)中的 (C-O)
Vec op = p-r.o; 
// t 为待求的自变量
// b 为公式(8)中的 D*(C-O)=-b'/2
// det 为公式(8)中的：
// det = (-b'/2)*(-b'/2)-{(C-O)*(C-O)-r^2} = (b'*b'/4-c')
double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad; 
// 若det<0 说明公式(8)无解，光线与球面没有交点
if (det<0) return 0; else det=sqrt(det); 
// 若 det>=0 说明公式(8)有解，光线与球面存在交点

// 令 det = sqrt(det) = sqrt(b'*b'/4-c')=sqrt(b'*b'-4c')/2
// 此时我们需要求距离光线原点最近的交点P(t)，并且要保证 t>0 即在光线的正方向处与球面相交

// 如果 t = b-det > 0，即
//     t = -b'/2 - sqrt(b'*b'-4c')/2
//       = (-b'-sqrt(b'*b'-4c'))/2 > 0
// 说明 t=b-det 是光线与球面相交的第一个合法交点

// 如果 t=b-det < 0，
// 说明 t=b-det 时的交点在光线的负方向处，因此合法交点为 
// t = b+det 
//   = -b'/2 + sqrt(b'*b'-4c')/2 
//   = (-b'-sqrt(b'*b'-4c'))/2
return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0);