背景
最近因为找到关于光线追踪相关不错的教程，所以边学习边做记录并希望将相关资料进行分享。
光线追踪作为计算机图形学中一种可以获得良好的效果的渲染算法，有着非常广泛的应用。历史背景相关的介绍可参考百度百科或者维基百科。本文中的参考资料来自于Peter Shirley分享的资料1，这个系列资料主要有三册，本文中主要介绍第一份资料相关的内容并且补充一些资料中缺少的细节内容。这份资料当中对于从基础知识开始到学习光线追踪算法的介绍都非常详细，本文中只做简单叙述。
光线追踪算法
光线追踪，从算法名称上来考虑简单来说就是从视点开始发射一条射线通过视平面上的每一个像素点并不断进行光线与物体交叉判定，同时考虑反射折射等光学现象来渲染三维场景。算法的终止条件主要由是否碰撞到物体，衰减，追踪的深度决定的。因此单纯从算法思想上来讲并不复杂，但是细节上其实仍然有许多点值得留意学习。
本文先从渲染最简单的图像开始介绍：
渲染图像
从一张200＊100的图像来讲，假设这张图像有3条通道（RGB）那么从数据量来讲就是200＊100＊3，再将图像简单抽象一下可以看成是一个200＊100大小的矩阵，3个通道就可看为3个200＊100大小的矩阵。要是渲染这样一张图像只要确定这个矩阵之中每一个对应下标的RGB数值便可以生成一张图像。我此次学习过程中尝试使用由c++开发的stb2类库生成图像。下面是一个使用stb来渲染图像的简单案例(注：并不列出所有代码所以想尝试的朋友可以通过我最后分享的github地址来参考)：
#include <iostream>#include "vec3.h"#define STB_IMAGE_WRITE_IMPLEMENTATION#include "stb/stb_image_write.h"#define STB_IMAGE_IMPLEMENTATION#include "stb/stb_image.h"using namespace std;int main(){int nx = 200;int ny = 100;int n = 4;unsigned char *data = new unsigned char[nx * ny * n];for (int j = ny - 1; j >= 0; j--){for (int i = 0; i < nx; i++){vec3 col(float(i) / float(nx), float(ny - 1 - j) / float(ny), 0.2f);data[j * nx * n + i * n + 0] = int(255.99 * col.r());data[j * nx * n + i * n + 1] = int(255.99 * col.g());data[j * nx * n + i * n + 2] = int(255.99 * col.b());data[j * nx * n + i * n + 3] = 255;}}cout << "write png to file!" << endl;stbi_write_png("cpt1_1.png", nx, ny, n, data, nx * 4);stbi_image_free(data);return 0;}
得到以下结果：
从光线到渲染图像
既然是光线追踪算法，那么最重要的就是这条光线的定义。我们可以将定义光线 p p p如下：
p ( t ) = A + t B p(t) = A + tB p(t)=A+tB
A A A是该条光线的起点， B B B是该条光线的方向，由 t t t来控制这条光线。
有了光线定义我们可以简单尝试一下利用光线来渲染图像，假设在一个三维空间之中设置一个点当作我们的眼睛，即视点。从这个视点出发我们可以看到眼前的图像信息，那么计算机中我们模拟这个情况需要确定一下图像的解析度，就是开始提到的图像的长宽，再加上左上角顶点的位置信息便可以确定眼前渲染的这个图像。如果说要通过光线来渲染一张图像的话，那么可以通过视点和视平面(图像)上所对应一点的坐标便可以构筑成一条光线，而根据这条光线所携带不同的信息以及一定的规则便可以渲染出一张由光线渲染出的图像。举一个简单例子可以参考如下代码：
#include <iostream>#include "ray.h"#define STB_IMAGE_WRITE_IMPLEMENTATION#include "stb/stb_image_write.h"#define STB_IMAGE_IMPLEMENTATION#include "stb/stb_image.h"using namespace std;vec3 color(const ray&r){vec3 unit_direction = unit_vector(r.direction());float t = 0.5*(unit_direction.y() +1.0);return (1.0-t)*vec3(1.0,1.0,1.0) + t*vec3(0.5,0.7,1.0);}int main(){int nx = 800;int ny =400;int n = 4;vec3 lower_left_corner(-2.0,-1.0,-1.0);vec3 horizontal(4.0,0.0,0.0);vec3 vertical(0.0,2.0,0.0);vec3 origin(0.0,0.0,0.0);unsigned char *data = new unsigned char[nx * ny * n];for (int j = ny - 1; j >= 0; j--){for (int i = 0; i < nx; i++){float u = float(i) / float(nx);float v = float(ny - 1 - j) / float(ny);ray r(origin, lower_left_corner + u*horizontal + v*vertical);vec3 col = color(r);data[j * nx * n + i * n + 0] = int(255.99 * col.r());data[j * nx * n + i * n + 1] = int(255.99 * col.g());data[j * nx * n + i * n + 2] = int(255.99 * col.b());data[j * nx * n + i * n + 3] = 255;}}cout << "write png to file!" << endl;stbi_write_png("cpt3.png", nx, ny, n, data, nx * 4);stbi_image_free(data);return 0;}
结果如下：
渲染球体到图像
三维球体经常被用于渲染一些简单的场景，主要还是得益于球体简易又经典的方程：
( x − c x ) 2 + ( y − c y ) 2 + ( z − c z ) 2 = R 2 (x-cx)^2+(y-cy)^2+(z-cz)^2 = R^2 (x−cx)2+(y−cy)2+(z−cz)2=R2
此处球体的中心为 C = ( c x , c y , c z ) C=(cx,cy,cz) C=(cx,cy,cz),半径为 R R R。
现在假设球面上有一点 p p p,那么我们可以来改写上述方程为:
( p − C ) ⋅ ( p − C ) = R 2 (p-C)\cdot(p-C) = R^2 (p−C)⋅(p−C)=R2
接下来我们就可以考虑如果在三维场景中放置一个球体，要通过发射光线方式来渲染出这个球体的话，首要考虑的是确定球体的位置。换句话来说，如果场景某处有球体的话那么发射出的光线应该会和球体有交叉。则接下来的问题就是如何判定光线与场景中某球体是否有交叉。通过下述公式可以推导出判定公式：
( p ( t ) − C ) ⋅ ( p ( t ) − C ) = R 2 (p(t)-C)\cdot(p(t)-C) = R^2 (p(t)−C)⋅(p(t)−C)=R2
( A + t B − C ) ⋅ ( A + t B − C ) = R 2 (A+tB-C)\cdot(A+tB-C) = R^2 (A+tB−C)⋅(A+tB−C)=R2
t 2 ( B ⋅ B ) + 2 t ( B ⋅ ( A − C ) ) + ( ( A − C ) ⋅ ( A − C ) ) − R 2 = 0 t^2(B\cdot B)+2t(B\cdot (A-C))+((A-C)\cdot (A-C)) -R^2= 0 t2(B⋅B)+2t(B⋅(A−C))+((A−C)⋅(A−C))−R2=0
此处 p ( t ) p(t) p(t)即为上文已定义的光线，接下来通过判定式便可知光线是否与球体有交叉，即：
bool hitSphere(const vec3 ¢er, float radius, const ray &r){vec3 oc = r.origin() - center;float a = dot(r.direction(), r.direction());float b = 2 * dot(r.direction(), oc);float c = dot(oc, oc) - radius * radius;float discriminant = b * b - 4 * a * c;return (discriminant > 0);}
通过如此方式可以得到如下结果：
如果尝试将球体表面的法线向量进行可视化，有如下结果：
抗锯齿
当单纯渲染图像时物体边缘普遍锯齿感较为强烈，主要是由于图像绘制算法而导致的问题。因而需要采取一些抗锯齿的算法来减轻物体边缘的锯齿感。此处采取一个简单的抗锯齿方法，在渲染某一点像素时通过随机采样方式，取得渲染点附近多个点的像素值并求得均值即可。
效果对比如下（上图未开启，下图开启）：
漫反射材质(Diffuse Material)
漫反射现象直白地描述一下，当多条平行光线射到某一平面之后光线会朝着各不相同的方向反射出去。需要模拟这样一个现象的话，我们首先考虑到的是各不相同的方向也就是不好确定一个具体方向，在计算机中也就是可以通过随机方式来选定一个方向进行反射。同时也要考虑到光线并不是完全反射出去，考虑到衰减的情况即其中有一部分的光线会被吸收掉。将这两者的特性同时考虑并进行编码的话：
vec3 random_not_in_unit_sphere(){vec3 p;do{p=2.0*vec3(_drand48(),_drand48(),_drand48()) - vec3(1,1,1);}while(p.squared_length()>=1.0);return p;}vec3 color(const ray&r, hitable *world){hit_record rec;if(world->hit(r,0.0, FLT_MAX,rec)){//p+normal = center, center+not_in = random directionvec3 target = rec.p + rec.normal + random_not_in_unit_sphere();return 0.5*color(ray(rec.p,target-rec.p),world);}else{vec3 unit_direction = unit_vector(r.direction());float t = 0.5*(unit_direction.y() +1.0);return (1.0-t)*vec3(1.0,1.0,1.0) + t*vec3(0.5,0.7,1.0);}}
可得到如下结果：
对图像进行伽马校正之后可得到：
光线追踪中的反射(Reflection)与折射(Refraction)
本章中一部分参考资料来自于Bram de Greve3, 主要是为了补充Peter Shirley资料中省略去的一部分推导。
到这一步，我们已经可以通过基本的光线来渲染简单的图像。但是由于没有考虑到渲染场景中很多物体的各种属性，因而无法渲染出较为精致的图像。这一章节我们先从反射和折射这两个物理现象入手：
这张图源自Bram de Greve的资料，通过这张图我认为会更容易理解反射和折射现象，也便于我们理解相关数学公式。
反射向量推导
要想模拟漫反射现象，那么首先需要通过入射向量来得到反射向量，此处我们首先需要对一个向量进行分解。对于一个方向向量来说我们可以分解其为一个水平方向向量和一个垂直方向向量：
v = v ⊥ + v ∥ v = v_{\perp}+v_{\parallel} v=v⊥+v∥
接下来我们尝试计算垂直方向的向量，如上图所示的法线向量是垂直与平面的，那么垂直方向的分量可以看作是原向量向法线向量的投影，可如下表示：
v ⊥ = v ⋅ n ∣ n ∣ 2 n = ( v ⋅ n ) n v_{\perp} = \frac{v\cdot n}{|n|^2}n = (v\cdot n)n v⊥=∣n∣2v⋅nn=(v⋅n)n
同时，考虑向量与法线夹角时，我们可以得到以下的结论：
cos θ v = ∣ v ⊥ ∣ ∣ v ∣ = ∣ v ⊥ ∣ = ± v ⋅ n \cos\theta_v = \frac{|v_{\perp}|}{|v|}=|v_{\perp}|=\pm v\cdot n cosθv=∣v∣∣v⊥∣=∣v⊥∣=±v⋅n
sin θ v = ∣ v ∥ ∣ ∣ v ∣ = ∣ v ∥ ∣ \sin\theta_v = \frac{|v_{\parallel}|}{|v|}=|v_{\parallel}| sinθv=∣v∣∣v∥∣=∣v∥∣
得到了以上的推导结果我们可以开始推导反射向量：
由已知条件：
θ i = θ r \theta_i = \theta_r θi=θr
r ⊥ = − i ⊥ , r ∥ = i ∥ r_{\perp}=-i_{\perp}, r_{\parallel}=i_{\parallel} r⊥=−i⊥,r∥=i∥
可得：
r = r ⊥ + r ∥ = i ∥ − i ⊥ = ( i − i ⊥ ) − i ⊥ = [ i − ( i ⋅ n ) n ] − ( i ⋅ n ) n = i − 2 ( i ⋅ n ) n \begin{aligned} r &= r_{\perp}+r_{\parallel} \\ &=i_{\parallel}-i_{\perp} \\ &=(i-i_{\perp})-i_{\perp} \\ &=[i-(i\cdot n)n]-(i\cdot n)n \\ &=i-2(i\cdot n)n \end{aligned} r=r⊥+r∥=i∥−i⊥=(i−i⊥)−i⊥=[i−(i⋅n)n]−(i⋅n)n=i−2(i⋅n)n
得到了反射向量 r r r之后我们就可以实现类似金属的材质。
vec3 reflect(const vec3 &v, const vec3 &n){return v - 2 * dot(v, n) * n;}
当物体表面反射出光线时再加上衰减的特性，最终可以渲染出下图：
折射向量推导
在实现了类似金属渲染之后，要模拟水，玻璃，钻石类似的材质时我们需要去计算入射向量的折射向量，这是由于光线需要穿过不同材质的物体而导致的。
这里我们最常用的是斯涅耳定律（Snell’s law），即：
η i sin θ i = η t sin θ t \eta_i\sin\theta_i=\eta_t\sin\theta_t ηisinθi=ηtsinθt
sin θ t = η i η t sin θ i \sin\theta_t = \frac{\eta_i}{\eta_t}\sin\theta_i sinθt=ηtηisinθi
此处注意我们需要约束：
sin θ i ⩽ η t η i \sin\theta_i \leqslant \frac{\eta_t}{\eta_i} sinθi⩽ηiηt
否则就会变成全内反射TIR(total internal reflection)现象。
通过上面公式我们可以继续变形：
∣ t ∥ ∣ = η i η t ∣ i ∥ ∣ |t_{\parallel}| = \frac{\eta_i}{\eta_t}|i_{\parallel}| ∣t∥∣=ηtηi∣i∥∣
由于 t ∥ t_{\parallel} t∥与 i ∥ i_{\parallel} i∥在同一方向上：
t ∥ = η i η t i ∥ = η i η t ( i − i ⊥ ) = η i η t ( i − i ⊥ ) t_{\parallel} = \frac{\eta_i}{\eta_t}i_{\parallel}=\frac{\eta_i}{\eta_t}(i-i_{\perp})=\frac{\eta_i}{\eta_t}(i-i_{\perp}) t∥=ηtηii∥=ηtηi(i−i⊥)=ηtηi(i−i⊥)
又由于存在：
∣ t ∣ 2 = ∣ t ⊥ ∣ 2 + ∣ t ∥ ∣ 2 |t|^2 = |t_{\perp}|^2+|t_{\parallel}|^2 ∣t∣2=∣t⊥∣2+∣t∥∣2
则我们还需要求 t ⊥ t_{\perp} t⊥：
∣ t ⊥ ∣ = 1 − ∣ t ∥ ∣ 2 |t_{\perp}|=\sqrt{1-|t_{\parallel}|^2} ∣t⊥∣=1−∣t∥∣2
由此数值可以推出向量 t ⊥ t_{\perp} t⊥：
t ⊥ = − 1 − ∣ t ∥ ∣ 2 n t_{\perp}=-\sqrt{1-|t_{\parallel}|^2}n t⊥=−1−∣t∥∣2 n
将 t t t水平和垂直的两个分量相加得到：
t ⊥ + t ∥ = η i η t ( i − i ⊥ ) − 1 − ∣ t ∥ ∣ 2 n = η i η t i − ( η i η t i ⊥ + 1 − ∣ t ∥ ∣ 2 ) n \begin{aligned} t_{\perp} + t_{\parallel}&= \frac{\eta_i}{\eta_t}(i-i_{\perp})-\sqrt{1-|t_{\parallel}|^2}n \\ &=\frac{\eta_i}{\eta_t}i-(\frac{\eta_i}{\eta_t}i_{\perp}+\sqrt{1-|t_{\parallel}|^2})n \end{aligned} t⊥+t∥=ηtηi(i−i⊥)−1−∣t∥∣2 n=ηtηii−(ηtηii⊥+1−∣t∥∣2 )n
此时可以将 ∣ t ∥ ∣ 2 |t_{\parallel}|^2 ∣t∥∣2转化为：
∣ t ∥ ∣ 2 = sin θ t 2 = ( η i η t ) 2 sin θ i 2 |t_{\parallel}|^2 = \sin\theta_t^2= (\frac{\eta_i}{\eta_t})^2\sin\theta_i^2 ∣t∥∣2=sinθt2=(ηtηi)2sinθi2
将此式代入上式可得：
t ⊥ + t ∥ = η i η t i − η i η t ( i ⋅ n ) n − 1 − ( η i η t ) 2 sin θ i 2 n = η i η t ( i − ( i ⋅ n ) n ) − 1 − ( η i η t ) 2 ( 1 − cos θ i 2 ) n = η i η t ( i − ( i ⋅ n ) n ) − 1 − ( η i η t ) 2 ( 1 − ( i ⋅ n ) 2 ) n \begin{aligned} t_{\perp} + t_{\parallel}&=\frac{\eta_i}{\eta_t}i-\frac{\eta_i}{\eta_t}(i\cdot n)n-\sqrt{1-(\frac{\eta_i}{\eta_t})^2\sin\theta_i^2}n \\ &=\frac{\eta_i}{\eta_t}(i-(i\cdot n)n)-\sqrt{1-(\frac{\eta_i}{\eta_t})^2(1-\cos\theta_i^2)}n \\ &=\frac{\eta_i}{\eta_t}(i-(i\cdot n)n)-\sqrt{1-(\frac{\eta_i}{\eta_t})^2(1-(i\cdot n)^2)}n \end{aligned} t⊥+t∥=ηtηii−ηtηi(i⋅n)n−1−(ηtηi)2sinθi2 n=ηtηi(i−(i⋅n)n)−1−(ηtηi)2(1−cosθi2) n=ηtηi(i−(i⋅n)n)−1−(ηtηi)2(1−(i⋅n)2) n
得到最终化简的结果后我们可以通过入射光线得到最终的折射光线，编码如下：
bool refract(const vec3 &v, const vec3 &n, float ni_over_nt, vec3 &refracted){vec3 uv = unit_vector(v);float dt = dot(uv, n);float discriminant = 1.0 - ni_over_nt * ni_over_nt * (1 - dt * dt);if (discriminant > 0){refracted = ni_over_nt * (uv - n * dt) - n * sqrt(discriminant);return true;}elsereturn false;}
此处还需要补充的是，当光线射到该类材质表面时何时会发生折射现象以及何时会发生反射现象。换句话来说，就是我们需要考虑到"物体表面接收到不同角度入射的光线时该点的反射率会不相同"这样的特性。因而此处我们将使用Christophe Schlick所提出的多项式近似方法来计算反射率4。 计算公式如下： R 0 = ( n 1 − n 2 n 1 + n 2 ) 2 R_0=(\frac{n_1-n_2}{n_1+n_2})^2 R0=(n1+n2n1−n2)2 R ( θ ) = R 0 + ( 1 − R 0 ) ( 1 − cos θ ) 5 R(\theta)=R_0+(1-R_0)(1-\cos\theta)^5 R(θ)=R0+(1−R0)(1−cosθ)5 参考代码如下： 经过渲染可以得到下图： 最终我们根据Peter Shirley提供的资料在场景中随机生成各种材质的球体并进行渲染，下图是最终渲染结果： 由于没有并行处理以及使用一些加速算法最终渲染相当消耗时间，可以作为改进点。 本文针对实现光线追踪算法中一些细节上数学推导进行了补充，详细关于各个细节实现部分请参考Peter Shirley的书籍。本人在学习过程中也自己动手实现了一遍基础光线追踪算法，下面是GitHub地址，大家也可以参考其中部分代码。 Ray Tracing 另外一提此系列一共有三本书，我感觉都很值得学习，如果有机会看完后面两本，我还会继续分享相关资料的总结。 如果有任何错误或者问题请指出，谢谢！ /petershirley/raytracinginoneweekend ↩︎ /nothings/stb ↩︎ De Greve, B. (). Reflections and refractions in ray tracing. URL http://users. skynet. be/bdegreve/writings/reflection_transmission. pdf (accessed -05-30). ↩︎ Schlick, C. (1994). "An Inexpensive BRDF Model for Physically-based Rendering" . Computer Graphics Forum. 13 (3): 233–246. CiteSeerX 10.1.1.12.5173. doi:10.1111/1467-8659.1330233 ↩︎float schlick(float cosine, float ref_idx){float r0 = (1 - ref_idx) / (1 + ref_idx);r0 = r0 * r0;return r0 + (1 - r0) * pow((1 - cosine), 5);}结语
后续
