SmallPT —— 99 行代码光线追踪解析

光线追踪（Ray tracing）是三维计算机图形学中的特殊渲染算法：根据光路可逆原理，对每一个像素，沿着入射光线逆向追踪若干次反射、折射，进而计算此光线的颜色，把场景渲染出来。

理论上，光线追踪算法可以完整地模拟物理世界中的光照，分毫不差地计算出每个像素的颜色，但是这样做的算力消耗趋近于正无穷。所以实际上所有光线追踪算法都包含了一些近似优化的逻辑，以将运算开销控制在可接受的数量级内。

Rhythm & Hues Studios 公司的程序员 Kevin Beason 曾于 2010 年编写过一个名为 SmallPT 的 C++ 程序，仅包含 99 行代码，即实现了最简单的光线追踪效果。此程序可视为光追算法的可运行最小集，是初学者学习和理解光追原理的极佳材料，其运行结果如下图所示。作为一个门外汉，我花了好几个晚上研究这 99 行代码，并在这个极好的 PPT 的帮助下，总算基本弄明白了其运行的原理。不妨记录下来：

Small PT 程序运行结果Small PT 程序运行结果

渲染方程

还是从最基础的渲染方程开始：

解释：

  1. 此方程描述的问题是：从物体表面上的一点 处，射入到观察者眼中的某条光线的强度，是如何确定的。
  2. 为射入观察者眼中的光线的颜色，即需要求取的值； 为出射方向。
  3. 为物体表面在点 向观察者方向自发射的光线的颜色（灯）。
  4. 表示环境入射到点 的光的颜色； 代表入射方向（为计算方便，取真实入射方向的反方向，不影响代表关系，后面也简称为入射方向）。
  5. 表示在给定 和 时，由 方向的入射单位光强产生的 方向的出射光的强度，此函数与表面的性质有关，又称表面的 BRDF 函数。
  6. 表示 与表面法线的夹角。
  7. 整个积分项表示：对半球（不透明材质）或全球（透明介质）内的所有入射方向 进行积分，得到的 方向的出射光强度。

依照光线追踪的原理，我们沿着 方向追踪到一处交点 ，就需要进行一次积分操作。在计算机程序中，积分是用求和模拟的，求和的次数越多（自变量的间隔越小），结果就越准确。假设每次积分都要进行 次求和操作，那么当追踪的光线遇到第一个交点时，会发散成 条光线；同时追踪这 条光线，到下一个交点，每条光线又会发散成 条光线……随着追踪深度的增加，计算开销的量级将按照指数级上升。

蒙特卡洛方法

蒙特卡洛方法（Monte Carlo method）是一种使用概率理论（通过大量随机数采样）进行数值计算以求取积分的方法。一个常用的有助理解的例子是：对「如何计算圆的面积」这个问题（圆的面积公式求取其实也是一个积分问题），蒙特卡洛的解法是「撒豆」：在包含圆的已知面积为 的矩形内随机采样（撒豆） 次，统计豆在圆内的次数为 ，则圆的面积为 。

撒豆问题不仅可以解圆的面积，还可以解任意形状，甚至不规则形状的面积求取。

更具体地，蒙特卡洛方法可以表述为：

  1. 左侧表示函数 在区间 的积分，即下图中的部分阴影部分面积。
  2. 此积分的值，可以这样求取：随机在 区间取值，采样计算 ，然后计算所有样本的均值并乘以区间的长度。当取样数量 越大，最后的值就越接近真实的积分值。

有同学可能会问（我也曾有此困惑），为什么不直接均等分采样，而要随机采样呢？其实，对于计算一维的函数 ，区别确实不大，但对二维甚至更高维度的函数（此时需要求取重积分）如 ，情况就不一样了。

此时，如果采取均等分采样，就需要选择在 上均等分为若干份，选择在 上均等分为若干份（其实撒豆问题已经是二维积分了，只不过积分函数是最简单的二值函数）。随着维度的增加，我们采样的数量也更难以控制，甚至还会出现维度不确定的情况。相比之下，蒙特卡洛方法的随机采样，可以轻易地控制或调整采样的次数：取两次（或更多次数的）随机数，可以视为单次随机行为，其背后包含的某种「随机性」是一致的。

对于光线追踪而言，我们会追踪若干次反射或折射（这个次数又称深度），每一次反射或折射都需要进行一次采样，相当于增加了一个维度。如果每次采样的次数过多，随着深度的增加，总采样次数很快就会不可接受，而如果每次的采样次数过少，（直觉告诉我）那么第一次采样对后续的采样将造成比较大的偏差（这背后应该有更完整的数学解释）。而且，在光线追踪算法中，有时候是否继续进行采样取决于具体的采样值，这就让等分采样变得更加困难。实际上，这份光追算法的实现，在每次反射或折射时只取一个随机样本，而通过增加总总采样次数 来保证最后的结果满足期望。

光线追踪算法

使用蒙特卡洛方法，对每一次反射或折射不再进行积分，而是随机选取一条可能的反射或折射光线进行追踪。然后在开始追踪的源头处，重复多次追踪操作以求取期望，这就是 SmallPT 光追的算法。简述一下具体步骤：

  1. 从每个像素 处发出一条射线 ，其方向与入射到相机并产生该像素的光线相反。
  2. 求取此射线照射到的物体表面的点 ，即与场景中物体的交点；如有多个交点，取距离相机最近的那个。如果未求到交点，则返回背景色。
  3. 射向 的光线强度，为 本身发射光强度 ，加上反射或折射环境光的强度。
    1. ：如果 处是光源，则 为光源的强度；否则， 为 0。
    2. ：根据 处根据表面性质，按概率随机取一次反射的射线 ，然后重复 2 的步骤，分别递归地求取 ，，， ……等等，直到满足一些特定条件停止递归。

以上便是单个像素的光追的算法。对单个像素，重复大量的次数求取平均值，作为此像素的颜色。

对每个像素完成以上步骤，就渲染出了整幅图像。

代码摘录

完整的代码摘抄如下：

#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/8, h=768/8, samps = argc==2 ? atoi(argv[1])/4 : 30;
  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];
  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,(unsigned short)(y*y*y)}; x<w; x++)
      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));
}

本文的主要部分，就是对这段代码的解析：

矢量运算

定义 Vec 结构体描述三维矢量，重载运算符来实现矢量运算。

struct Vec {
  double x, y, z;
  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; }
  Vec operator%(Vec&b){return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}
};

矢量运算包括（假设 ；）：

  1. 加法 operator+：。
  2. 减法 operator-：。
  3. 乘法 mult()：。这种矢量分量直接相乘的运算通常用于颜色的计算，在三维空间中并无特殊的物理意义。
  4. 点乘 dot()：，其中 为 与 的夹角。点乘的结果是一个标量，物理含义是 在 方向上的投影。
  5. 叉乘 operator%：，叉乘的结果是一个矢量，物理意义是垂直于 和 的矢量（按 顺序右手螺旋），其长度为 和 构成的平行四边形面积。
  6. 矢量乘以标量 operator*：。
  7. 归一化 norm()：，其物理意义是保持矢量方向不变，将其长度缩放为单位长度 1。

射线

通过定义射线的出射点 和出射方向 来定义射线：

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

定义场景

如前文 SmallPT 的渲染效果图所示，整个场景由一个房间，一个光源和两个球体组成。为了简单，SmallPT 将场景全部用球体来表示，房间的墙壁就用直径极大（使观者感受不到墙壁是弯曲的）的球体来表示。

球体对象的表示

  1. 首先枚举出三种表面特性：散射面——用于墙壁，镜面——用于左侧的镜面小球，折射面——用于右侧的玻璃小球。
  2. 定义球体，内容包括：
    1. 数值类型成员：半径 ；
    2. 矢量类型成员：圆心位置 ，表面发光强度（颜色），表面散射颜色 ；
    3. 表面特性 ；
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_) {}
  // ...
};

射线与球体相交

Sphere 上定义了 intersect 方法，以求取给定的射线与球体是不是相交。如果相交则返回射线原点与交点的距离。

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

注意，下面为了清楚，我们使用大写字母为矢量，小写字母为标量来进行推导。我们需要求解的是射线的原点 与交点 的距离 ，而条件是：

  1. 点在射线上，即 。
  2. 点在球面上，即 与球心 的距离等于球的半径 。即 ，也就是 。

将 (1) 中关于 的表达式代入 (2)，得到：

展开得到：

Undefined control sequence \*

这已经是一个一元二次方程，自变量 是未知的。所有矢量经过平方或点积后都成为了标量，可以直接根据一元二次方程 的求根公式：

计算出根 。

  1. 一元二次方程可能没有实根的，此时射线与球体不相交。这对应着代码中 det<0 的情况，直接返回 0。
  2. 一元二次方程可能有两个实根，此时射线（所在的直线）与球体有两个交点。我们需要取的是最小的正根。如果另一个根 （对应交点为 ）大于 ，表明是左图的情况；如果 小于等于 0，表明是右图的情况，这条光线是在介质内的折射光线。

使用球体表示场景

代码定义了一个 spheres 数组来表达整个场景，每一个球体分别是：

  1. 左侧墙体，红色；墙体的半径都很大，观者感知不到曲率；墙体的反射类型都是散射，但是颜色有所不同。
  2. 右侧墙体，蓝色。
  3. 后面墙体，灰色。
  4. 前面墙体，黑色。
  5. 底部墙体，灰色。
  6. 顶部墙体，灰色。
  7. 左侧球体，镜面材质。
  8. 右侧球体，透明玻璃材质，光线可能折射进入玻璃中继续传播。
  9. 顶部的光源，发射出强度为 12 的白光。
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
};

可以看出，坐标轴的原点位于盒子的左面、后面和底面的交点附近。宽度 100 左右，高度 80 左右，深度 170 左右。

准备渲染

在真正开始渲染之前，还需要做一些准备工作。

确定相机和输出图像

先跳过 clamp，toInt，radiance 这几个函数，直接来看 main 函数。

int main(int argc, char *argv[]){
  int w=1024, h=768, samps = argc==2 ? atoi(argv[1])/4 : 30;
  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];
  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,(unsigned short)(y*y*y)}; x<w; x++)
      // ... 暂时先省略
  }
  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));
}

首先，定义需要渲染图像的宽度（1024），高度（768），对每个像素的采样次数（命令行参数传入，或 120 次）。注意，变量 samps 的值是采样次数的 1/4，在这个例子中，最终渲染在图像上的每一个像素，又会拆分成 4 个子像素进行追踪和采样，而 samps 的值是子像素的采样次数。

然后，使用射线类的变量 cam 来表示相机，射线的原点 cam.o 即是相机的位置，而射线的方向 cam.d 是相机朝向的方向。由于相机本身的位置只是一个点，所以需要把渲染出的图像视为在屏幕前方（亦或是后方，参考一下眼球和相机的构造）不远处的一块矩形幕布，抵达相机的光线穿透幕布时会在与幕布的交点 处留下颜色，这样就在幕布上投影了。只要幕布的分辨率和宽高比不变，其距离相机的远近与最终结果是无关的，所以程序直接设定幕布与相机的距离是 1，设定幕布本身的高度为 0.5135（配合像素尺寸的宽高比，其实这里已经暗含了相机的水平和垂直视场角的信息）。然后，求取长度等于幕布真实宽度 0.5135*w/h 的水平矢量 cx（平行于 X 轴），求取长度等于幕布真实高度 0.5135 的垂直矢量 cy（通过 cx 与 cam.d 叉乘而来）。

再接着，生成长度为 w*h 的 Vec 类型数组 c 备用。这个数组的作用是存储渲染的结果。

最后，遍历宽度和高度，对每一个像素进行渲染，并将渲染得到值写入到数组 c 中，并将其输出为 PPM 格式的图片，以便打开查看渲染效果。

PPM 是一种基于文本的图片格式，能够帮助你避免考虑图片文件编码问题，轻易地按像素输出图片。比如以下文本内容就是一张 4x4 的图片。

P3
# feep.ppm
4 4
15
 0  0  0    0  0  0    0  0  0   15  0 15
 0  0  0    0 15  7    0  0  0    0  0  0
 0  0  0    0  0  0    0 15  7    0  0  0
15  0 15    0  0  0    0  0  0    0  0  0

总结一下各个变量的含义：

  1. w 和 h：宽度和高度像素数。
  2. samps：每个像素采样次数的 1/4。
  3. cam：相机描述。
  4. cx 和 cy：长度与幕布真实尺寸相同的水平、垂直矢量。
  5. c：按像素存储渲染结果的数组。

子像素采样

接下来，我们看一下如何计算每个像素的值的（即之前省略的部分）。对每一个像素，都将其拆分为 2x2 一共 4 个子像素。一个当然的想法是，以子像素的中心点 作为射线的出发点开始追踪，但是由于每个子像素都会进行 samps 次光追取均值，所以这里就引入了一次随机过程：以子像素的中心点，在一个像素大小的范围内进行一次随机采样，并以采样到的点 为射线的出发点，进行一次光追。对每个子像素完成 samps 次光追，依次以 ， ……为出发点，最后求均值。

而且，为了使渲染的质量更高一些（达到相同渲染质量所需的采样次数更小一些），这里对随机过程进行了优化：对随机数进行了一次 tent 滤波。此滤波函数将 从 区间映射到了 ，当 为随机的值时， 的分布更向中点值 集中（从图中可以看出）。

Missing \left or extra \right

然后，拿到经过滤波后的随机数 ，，计算 在整个画布中的归一化坐标（假设成为 和 ，如 nx=(sx+0.5+dx)/2+x)/w-0.5），再乘以之前求得的画布真实尺寸矢量 和 ，加上相机本身的方向矢量 ，就得到了最终从相机原点射向 的矢量 ，即代码中的 d。

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

最后，调用 radiance，沿着穿过此子像素的射线，在 140 倍远（让光追的开始点进到盒子里面，而不会直接被朝前的表面挡掉）的地方开始进行光线追踪算法。进行 samps 次，取均值，再基于 4 个子像素再取一次均值，就得到了这个像素的颜色。

erand48() 函数是 C 标准库提供的随机函数。传入长度为 3 的数组作为随机种子，返回 [0, 1] 区间的双精度浮点数，同时改写传入的随机种子，便于下一次调用 erand() 时传入（以获取一个不同的结果）。

总结一下各个变量的含义：

  • i：为当前像素在 c 中的索引。
  • r：通过光追计算出的子像素的颜色。
  • sx 和 sy：子像素索引，为 0 或 1。
  • r1 和 r2：滤波用自变量，在 [0, 2] 区间随机取得。
  • dx 和 dy：滤波后的随机自变量，分布在 [-1, 1] 区间内。
  • d ：从相机原点指向子像素位置的矢量。

光线追踪

讲了这么多，终于到正餐了。函数 radiance() 即进行了一次光线追踪。

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
      // ...
    } else if (obj.refl == SPEC)            // Ideal SPECULAR reflection
      // ...
    // ...                                   // Ideal dielectric REFRACTION
}

函数 radiance() 接收三个参数，射线 r，深度 depth 和随机种子 Xi。此函数是一个递归函数，在计算射线颜色时候，遇到反射或折射，会计算出下一次反射或折射的射线，然后递归调用自己求算下一次射线的颜色。深度 depth 标记了反射/折射的次数。

下面看 radiance() 函数的逻辑：

求取交点和递归结束条件

首先，对场景中的所有球体对象，尝试求取射线与之交点。函数 intersect() 负责做这件事：遍历 spheres 数组，依次调用 sphere 对象上的 intersect() 方法，并保留最小的正值 t 和对应球体在 spheres 数组中的索引 i。

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

如果没有任何对象与之相交，直接返回黑色（递归的结束条件 1）。

如果检测到有球体与射线相交，首先检查当前深度是否大于了 5（之前已经经过了 5 次反射或折射），如果是，那么遵循一个概率来随机地决定采取以下哪一种方案：

  1. 放弃反射光的贡献，直接返回物体的出射光——虽然大部分时候是 0（递归的结束条件 2）。
  2. 继续追踪光线。

这个概率就是当前物体表面的颜色 RGB 分量中最大的那个分量，归一化后的值。换言之，物体越暗（的颜色值越低），越容易被筛选到第 1 种方案（因为物体暗，环境光反射贡献的值也会小），而物体越亮，越容易被筛选到第二种方案。值得注意的是，当选择第二种方案时，需要把物体的颜色按照概率调高一些 f=f*(1/p)，这是因为那些运气不够好，被筛选到 1 方案光线，其反射分量完全地消失了，这会导致最终计算出的颜色的期望（也就是说，即使采样无穷次）与真实值产生偏差。在运气好的 2 方案中调整颜色，可以消除这个偏差。

按照概率进行随机选择是无处不在的，只要采样的次数足够多，那么遵照概率选择计算得到的均值就能够逼近真实值。

接着，如果深度小于等于 5，那么不管物体的颜色值如何，都会正常地递归地进行光线追踪。在此之前，还需要预先求取一些变量的值，比如法线 nl 等，便于后续使用。

总结一下各变量的意义：

  • t：交点距离相机原点的距离。
  • id：相交的球体在 spheres 数组中的索引。
  • obj：相交的球体 Sphere 对象。
  • x：交点的位置。
  • n：球体在交点处的归一化的法线（从球心射向表面）。
  • nl：与反射/折射上下文契合的归一化的法线：如果射线在球体外部，nl 与 n 相同；如果射线在球体内部（折射），nl 与 n 相反。

散射（漫反射）

追踪光线，当代表光线的射线与球体相交时，根据球体表面的特性，进行不同的操作（代码中省略的部分）。

如果物体表面是漫反射，将随机选取一个角度取归一化的反射光线 ：

  1. 确定互相正交的三个归一化矢量 ， 和 ，其中 为反射面的法线。
  2. 在区间 中随机确定一个值 ，并使 作为反射光线 在法线 方向上投影的长度，使 也就是 成为反射光线在表面上投影的长度。
  3. 在区间 中随机确定一个角度值 ，使之成为反射光点 在反射表面（即 平面）的投影与 轴的交角。
  4. 极坐标公式转笛卡尔坐标， 计算出在半球面的随机矢量，再归一化得到反射矢量 。

最后，递归地调用 radiance() 函数，将反射光线作为参数传入，传入已自增过一次的深度和随机种子。其结果乘以表面的颜色，再加上球体表面的发射光颜色，作为结果返回。这就是此光追算法对散射表面的处理。

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

总结一下各变量的含义：

  • u，v，w：互相正交的三个归一化矢量，w 即为法向量。
  • r1：反射光线在反射表面与 u 的交角。
  • r2 与 r2s：r2s 为反射光线在反射表面投影的长度，r2 为此值的平方。
  • d：随机求取的反射方向。

镜面反射

镜面反射比较简单，因为对于确定的入射光线，反射的方向也是确定的，因此没有随机过程，直接取镜面反射的方向继续进行光追。

镜面反射的求法比较简单，不过也记录一下吧：如下图对称的两条射线 和 ，而法线为 （法线是归一化的）。此时有 ，求解 。在代码中， 就是 -r.d，带入得到 r.d-n*2*n.dot(r.d)，即反射光线的。由于 r.d 也已经是归一化的了，所以求取得到的结果也一定是归一化的。

代码如下：

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

折射

折射是比较复杂的。不同的物体有不同的折射率，空气的折射率为 1，玻璃的折射率为 1.5。根据折射定理：光线从折射率为 的介质折射到折射率为 的介质，折射角满足（此定理可以从更一般的麦克斯韦电磁方程组推导出来，厉害厉害）：

全反射

首先需要考虑的是全反射现象，即光线从光密介质（玻璃）向光疏介质（空气）传播时（上图右半部分）， 大于临界角（使得 为 90° 的 ），以至于不存在 满足上述折射定理时，所有光都被反射的现象，即：

两边平方并作简单变换，得到：

由于 nnt 为 ，而 ddn 为 ，代入上式，得到全反射的条件即是 1-nnt*nnt*(1-ddn*ddn))<0。满足全反射条件时，镜面反射的方式进行处理。

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

总结一下各变量的含义：

  • reflRay：镜面反射方向，与上一节中的求法完全一致。
  • into：由空气（试图）进入玻璃则为 true，由玻璃（试图）进入空气则为 false。
  • nc 和 nt：空气和玻璃的折射率。
  • nnt 为折射率比值，当追踪的射线是从空气进入玻璃时，nnt 为 1.0/1.5，当射线是从玻璃进入空气时，nnt 为 1.5。
  • ddn：即 的值，由 和 点乘而来。

正常折射

所谓正常折射，即是光线一部分镜面反射，一部分折射进入另一种介质。此时，我们需要：

  1. 计算折射光线矢量 。
  2. 计算有多少比例的光线被反射，有多少比例的光线被折射。
  3. 根据条件，选择以随机或非随机的形式，递归地进行光线追踪。

首先计算折射光线矢量 ：

首先，假设入射光线在表面平面的投影矢量（并归一化）为 已知，折射角 也已知，我们可以得到 为：

其中 就是入射光线平面上的投影矢量（因为入射和折射在同一个平面内）归一化得到，而此投影矢量可以使用 减去 在法线上的投影 （注意这里点积是个负数）来得出，因此：

同时 的正弦和余弦值也可以通过折射定律求得：

全部代入，可得：

按照上述公式，将变量代入，即可计算出折射光线归一化矢量 tdir：

Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();

接下来，计算有多少比例的光线被反射，有多少比例的光线被折射。根据电动力学的一些理论，入射光线（普通的非偏振光）被反射/折射的比例，与入射角 ，两种介质的折射率都有关系。当入射角等于 0（即垂直入射到表面）时，被反射的比例 为：

而当入射角为 时的反射比例 为：

最后，基于当前光追的深度：

  1. 如果深度小于等于 2，那么分别递归地追踪反射光和折射光，将两者对颜色的贡献再加上物体本身的发光，作为结果返回。
  2. 如果深度大于 2，那么依据折射和反射光线的占比来进行概率随机，选择追踪反射光或折射光中的一支。换言之，如果反射光占比高，那么就有较大的可能性随机到反射光线，反之亦然。注意，这时由于我们放弃了另一种可能性，会导致结果的期望产生偏差，所以需要进行一次修正。

对玻璃材质，前几次反射、折射对最终的结果的贡献比较大，所以选择都进行追踪。但是这样做一下子使得某些追踪的采样次数发散为了原来的 4 倍，如果深度改为 4，可能就是 16 倍。所以代码把发散深度限制在 2 这样一个比较小的数值，也算是在效果和性能间达到一个平衡。

代码如下：

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

最后解释一下各个变量的含义：

  • a 和 b：用以计算 的临时变量。
  • R0：入射角为 0 时反射光的占比。
  • c：即 。
  • Re：当前入射角下反射光的占比。
  • Tr：折射光的占比。
  • RP：随机到反射光时，对结果的修正的因子。
  • TP：随机到折射光时，对结果的修正的因子。

这样，光追算法就全部完成了。

感想

一直觉得光线追踪是非常高端的存在，因为没有基础，想了解也不知从何看起。一个偶然的机会发现到这个 demo，觉得挺有意思，就拿过来研究了。结果发现看起来还是比较吃力，好在另外一个老外的解读 PPT 帮了很大的忙。

把学习过的东西消化了，再重新写出来是件很折磨的事情，但收获也不小。每一个公式一步一步推导过来，对一些地方「为什么这么做」也有了不一样的体会。数学还是挺重要的，而且能够给人带来很强的成就感。