当前位置：首页 > news >正文

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

news 2025/7/9 8:41:34

一、简介

本文紧接在[图形学]smallpt代码详解（1）之后，继续详细讲解smallpt中的代码，包括自定义函数（第41到47行）和递归路径跟踪函数（第48到74行）部分。

二、smallpt代码详解

1.自定义函数（第41到47行）

smallpt源代码在第41到47行声明定义了clamp()、toInt()和intersect()三个函数：

double clamp(double x)：
```
inline double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }
```
clamp()函数将输入参数x限制在[0,1]范围内，避免浮点数计算出错，出现小于0.0或者大于1.0的计算结果。
int toInt(double x)：
```
inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); } 
```
toInt()函数将输入的在[0,1]范围内的输入参数x先进行伽马矫正（Gamma Correction），伽马校正线性颜色空间转换为 sRGB 颜色空间。伽马校正有助于提高图像在显示器上的视觉效果。，然后映射到范围[0,255]。
bool intersect(const Ray &r, double &t, int &id)：
```
inline bool intersect(const Ray &r, double &t, int &id){ 
// n 为场景中的球面物体个数
// d 用来记录场景中与光线相交的第一个物体id
// t 用来记录场景中与光线相交的第一个交点到光线原点的距离
double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20; 
// 遍历场景中的所有球面对象，然后使用 sphere[i] 对象中的函数interscet() 函数计算
// 光线是否与 球面对象 sphere[i] 相交，并找到距离光线原点最近的交点距离 t 和球面对象 i
for(int i=int(n);i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;} 
return t<inf; 
} 
```
intersect()函数用于计算光线r与整个渲染场景中所有球面对象的相交计算，如果光线r与场景存在交点，则返回true，并且将交点到光线原点的距离t，写入输入的变量double &t中，并且返回相交球面物体的id，写入输入的变量int &id中。
代码首先确定场景中的球面物体个数n：
```
double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20;
```
然后遍历场景中的所有球面对象，然后使用 sphere[i] 对象中的interscet() 函数计算光线是否与球面对象 sphere[i] 相交，对于sphere[i]对象中的intersect()函数的讲解请查看[图形学]smallpt代码详解（1）中的三、smallpt代码详解-1.自定义数据结构部分（第4到第22行）部分，并找到距离光线原点最近的交点距离 t 和球面对象 i。
```
for(int i=int(n);i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;} 
```

2.递归路径跟踪函数（第48到74行）

radiance()函数是光线跟踪（路径跟踪）中最主要的代码，该函数使用递归的方式实现路径跟踪，对场景进行渲染。
函数Vec radiance(const Ray &r, int depth, unsigned short *Xi)的返回值和参数介绍如下:

返回值 Vec：radiance() 函数计算的光线在此次传播中的辐射值；
参数 r：要处理的光线
参数 depth：此光线的传播深度
参数 Xi：随机数，作为随机采样种子值

radiance()函数整体可以分为四个部分是否结束递归（与场景无交点或者俄罗斯轮盘赌失败）（红色）、处理漫反射（绿色）、处理镜面反射（蓝色）和处理折射（橙色）部分。流程如下：

radiance()函数代码流程：
radiance()函数流程图
下面介绍各部分的详细代码实现：

2.1) `是否结束递归（与场景无交点或者俄罗斯轮盘赌失败）`：

double t;                               // distance to intersection 
int id=0;                               // id of intersected object 
// 计算当前光线 r 和场景的交点，如果没有交点，则返回 Vec()，即返回 (0,0,0)
if (!intersect(r, t, id)) return Vec(); // if miss, return black 
// 如果 r 与场景有交点，那么记录交点对应的球面对象 obj
const Sphere &obj = spheres[id];        // the hit object 
// x 为光线 r 与 obj的交点位置
// n 为交点 x 处的球面法向
// nl 为定向法向，用来描述光线r是从球外射向球内，还是从球外射向球内
//    假如 r 为从球外射向球内，则 nl 为朝向球外的法向
//    假如 r 为从球内射向球外，则 nl 为朝向球内的法向
// f 为 obj 的漫辐射颜色
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.

代码中首先判断光线r是否与场景相交：

1). 如果不相交则说明光线射出到场景之外，因此直接返回(0,0,0)。
2). 如果光线r与场景相交，则计算得到交点位置x，交点所处的球面对象obj，和交点处的法向n，与r在球面同边的法向nl。
- 2.1) 如果此时的递归深度大于5，那么使用俄罗斯轮盘赌确定是否结束递归。在俄罗斯轮盘赌策略中，smallpt使用交点处颜色f的(r,g,b)三个颜色分量中的最大值作为阈值p，然后使用erand48()函数生成一个在[0,1]范围内的随机数。此处之所以使用颜色(r,g,b)中最大值作为阈值p是考虑假如交点处材质颜色接近黑色，可以尽可能地结束递归，因为此时(r,g,b)三个分量值都很小，因此该点最终反射的radiance很小，对渲染地结果影响也很小。反之，该交点颜色接近白色，说明该交点最终反射的radiance可能很大，因此要尽可能地让光线继续传播下去。
  - 2.1.1) 假如该随机数小于p，说明赢得此次俄罗斯轮盘赌，将交点处的颜色值f变为原来的(1/p)倍数，然后继续往下运行；
  - 2.1.2) 假如该随机数大于p，说明败于此次俄罗斯轮盘赌，直接返回交点处的自发光颜色obj.e。
- 2.2) 如果此时的递归深度小于等于5，则继续往下运行。

2.2) `处理漫反射`：

if (obj.refl == DIFF){                  // Ideal DIFFUSE reflection // r1 是在范围 [0,2Pi] 范围内的随机数// r2 是在范围 [0,1] 范围内的随机数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));

在处理漫反射材质表面时，在交点处的上半球面上的cos加权采样，采样得到光线r对应的入射方向d，然后递归地计算光线下一次传播radiance(Ray(x,d),...)的结果。
半球面上的cos加权采样，即各方向的采样概率密度与 $cos(\theta)$ 成正比， $\theta$ 为采样方向与交点处法向nl的夹角。

采样步骤如下：

首先采样方位角r1和半球对应的圆平面半径长度r2s。
然后计算以交点处法向nl为(0,0,1)轴的坐标系(u,v,w)。
再之后根据采样的方位角r1和圆平面半径长度r2s在u-v平面的圆平面上采样得到采样点p，该方法即在平面u-v上的圆平面内均匀采样得到采样p。采样点p对应到半球面上的映射点P的向量即为目标方向向量d。
该方法与半球面上的cos加权采样等效的原因请查看博客[图形学]在半球面上均匀采样和cos加权采样中的3.在半球面上的cos加权采样：方法二部分。

半球面上cos加权采样示意图

采样点p在(u,v,w)坐标系下的坐标为：
$(cos (r 1) * r 2 s, s in (r 1) * r 2 s, 0)$ ，
映射点P在(u,v,w)坐标系下的坐标即为：
$cos(r1)*r2s, sin(r1)*r2s, sqrt(1-r2s^2))=(cos(r1)*r2s, sin(r1)*r2s, sqrt(1-r2))$
那么在世界坐标系(x,y,z)下，半球面上的映射点P坐标即为：
$(u * cos (r 1) * r 2 s, v * s in (r 1) * r 2 s, w * s q r t (1 - r 2))$
由于我们是使用交点处的法向nl作为坐标系w轴，暗含了以该交点为坐标系原点，因此映射点P的坐标即采样向量，即代码中的变量d。
最后返回：该交点处的自发光obj.e+该交点处的颜色值f乘以光线下一次传播计算得到的radiance值。

2.3) `处理镜面反射`：

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

处理镜面反射相对比较简单，根据光线r方向和交点处法向n可以方便地计算得到入射光线方向。
根据反射光方向和法向计算入射光方向（或者根据入射光方向和法向计算反射光方向）的算法示意图如下所示：
反射光线计算示意图
最后返回：该交点处的自发光obj.e+该交点处的颜色值f乘以光线下一次传播计算得到的radiance值。

2.4) `处理折射`：

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

处理折射时，首先判断光线是否发生全内反射，如果发生全内反射，那么就直接根据理想的反射方向继续在玻璃内传播。如果没有发生全内反射则基于斯涅尔定律计算折射方向，然后使用菲涅尔方程和Schlick近似 计算折射和反射的比例。
折射示意图

全内反射判断：
根据全内反射定律，当入射角 $\theta_{a}>\theta_{c}$ 时发生全内反射。而 $\theta_{c}=arcsin(\frac{n_b}{n_a})$ ；
即当：
$\theta_{a} > \theta_{c}=arcsin(\frac{n_b}{n_a}) \\\ \theta_{a} > arcsin(\frac{n_b}{n_a}) \\ sin(\theta_{a}) > \frac{n_b}{n_a} \\ sin^2(\theta_{a}) > (\frac{n_b}{n_a})^2 \\ 1-cos^{2}(\theta_{a})>(\frac{n_b}{n_a})^2 \\ (\frac{n_a}{n_b})^2(1-cos^{2}(\theta_{a}))>1 \\ 1 - (\frac{n_a}{n_b})^2(1-cos^{2}(\theta_{a})) < 0 \tag{1}$
时，发生全内反射。
而代码中nnt即 $\frac{n_a}{n_b}$ ，ddn等于 $-cos(\theta_{a})$ ，ddn*ddn即 $cos^{2}(\theta_{a})$ 。
因此，当cos2t=1 - nnt * nnt * (1 - ddn * ddn)<0时，发生全内反射。
计算折射方向：
根据斯涅尔定律，假如入射光方向为 $D$ ，折射光方向为 $T$ ，那么应该满足如下公式：
$n_{a}*sin(\theta_{a}) = n_{b}*sin(\theta_{b}) \tag{2}$
其中 $\theta_{a}$ 为入射角， $theta_{b}$ 为折射角， $n_a$ 为入射介质折射率， $n_b$ 为出射介质折射率。

如上图所示，我们可以得到：
$sin(\theta_a)=\sqrt{1-cos^2(\theta_a)}=\sqrt{1-(D*N)^2}\\ sin(\theta_b)=\frac{n_a}{n_b}sin(\theta_a) = \frac{n_a}{n_b}\sqrt{1-(D*N)^2} \\ cos(\theta_b) = \sqrt{1-sin^2(\theta_b)}=\sqrt{1-\frac{n_{a}^2}{n_b^2}(1-(D*N)^2)} \tag{3}$
单位向量 $B$ 等于：
$|cos(\theta_a)| N) = \frac{D - |cos(\theta_a)| N}{sin(\theta_a)}=\frac{D +(D\cdot N)N}{\sqrt{1-(D\cdot N)^2}} \tag{4}$
目标折射光向量 $T$ 等于：
$Bsin(\theta_b)-Ncos(\theta_b) = \frac{D +(D\cdot N)N}{\sqrt{1-(D\cdot N)^2}} * \frac{n_a}{n_b}\sqrt{1-(D \cdot N)^2} - N\sqrt{1-\frac{n_{a}^2}{n_b^2}(1-(D\cdot N)^2)} \\ = \frac{n_a}{n_b}*(D+D\cdot N)N-N\sqrt{1-\frac{n_{a}^2}{n_b^2}(1-(D\cdot N)^2)} \\ = \frac{n_a}{n_b}D+\frac{n_a}{n_b}N(D\cdot N)-N\sqrt{1-\frac{n_{a}^2}{n_b^2}(1-(D\cdot N)^2)} \\ = \frac{n_a}{n_b}D+N \left( \frac{n_a}{n_b}(D\cdot N)-\sqrt{1-\frac{n_{a}^2}{n_b^2}(1-(D\cdot N)^2) } \right) \tag{5}$
因为代码中nnt= $\frac{n_a}{n_b}$ ，ddn= $-cos(\theta_{a})=-(D\cdot N)$ ，cos2t=1 - nnt * nnt * (1 - ddn * ddn)= $1-\frac{n_{a}^2}{n_b^2}(1-(D\cdot N)^2)$ 。
因此，结果折射光线 $T$ 就等于代码中的：
```
	Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm(); 
```
计算折射和反射的比例：
根据菲涅尔方程和Schlick近似，可以计算得到光线发生折射和反射的比例，Schlick近似公式如下：
$n=\frac{n_a}{n_b} \\ R_{0}=\frac{(n-1)^2}{(n+1)^2}=\frac{(n_a-n_b)/n_b}{(n_a+n_b)/n_b}=\frac{n_a-n_b}{n_a+n_b} \\ R_{r}(\theta_{a}) = R_{o} + (1-R_{o})(1-cos(\theta_{a}))^{5} \tag{6}$
其中 $R_{0}$ 为光线垂直入射（即入射角 $\theta_{a}=0$ ）时的反射率，即代码中的变量R0， $R_{r}(\theta_{a})$ 为入射角为 $\theta_{a}$ 时的反射率，即代码中的变量Re。
```
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); 
```
在结束折射部分代码时，由于smallpt中规定折射材质的颜色c=(1.0,1.0,1.0)，在radiance()函数开始，使用p=max(f.r, f.g,f.b)进行俄罗斯轮盘赌时肯定无法失败、结束递归。因此在此处增加了一次俄罗斯轮盘赌，使用P=0.25+0.5*Re作为俄罗斯轮盘赌的阈值。
当递归深度大于2时使用P作为阈值进行俄罗斯轮盘赌，决定是否结束递归。如果不结束，就将Re比例的光线用于接下来的镜面反射，Tr=1-Re比例的光线用于接下来的折射。

在接下来的 [图形学]smallpt代码详解（3） 中，将继续讲解 smallpt 中的main函数部分，包括渲染场景的定义、相机的设置、光线的生成、滤波处理减少噪点和保存最终渲染结果几部分。

三、参考

[1].smallpt: Global Illumination in 99 lines of C++
[2].smallpt: Global Illumination in 99 lines of C+±Presentation slides
[3].光线跟踪smallpt详解 (一)
[4].斯涅尔定律
[5].菲涅尔方程
[6].Schlick近似

一、简介

二、smallpt代码详解

1.自定义函数（第41到47行）

2.递归路径跟踪函数（第48到74行）

2.1) 是否结束递归（与场景无交点或者俄罗斯轮盘赌失败）：

2.2) 处理漫反射：

2.3) 处理镜面反射：

2.4) 处理折射：

三、参考

相关文章：

2.1) `是否结束递归（与场景无交点或者俄罗斯轮盘赌失败）`：

2.2) `处理漫反射`：

2.3) `处理镜面反射`：

2.4) `处理折射`：