于是我打算自己写一篇博客。以【把大象装进冰箱】进行比喻：如果说网络上的资料是【打开门，送大象进去，关门】的话，那本篇博客就是讨论【怎么装】的。比如冰箱门把手在哪里，朝哪个方向拉才能打开冰箱，要用香蕉还是鞭子把大象赶进去，怎么关上冰箱的门…
不多说批话，直接开始文本的内容！
在传统的计算机图形学中，光栅通常用于生成像素。因为光栅化后，我们只选择我们看到的图形，视野外的几何信息被切割，导致场景的整体信息丢失。然后它会导致一些渲染例如，水只能反映我们在屏幕上看到的像素：
此外，一些基于全球信息的渲染效果不能很好地运行，如全球光。在光栅管道下，我们通常需要通过各种类似的算法来模拟现实世界的光影现象，这是非常麻烦和有效的“还凑合”
超越传统图形流水线的光线跟踪是一种现代渲染。光跟踪是一种基于物理的渲染方法，通过模拟介质中光的真实性能来输出足够接近现实的图像。光从光源开始，通过各种反射折射进入相机：
从光源开始无数条到达相机的路线不太现实，从灯泡开始模拟光的传播，所以曲线拯救了国家，从相机的方向开始向场景投射射线，试图找出到达光源的可能路径：
本博客代码主要实现路径跟踪技术。路径跟踪技术是光跟踪的一个分支。通过将光投射到场景中，模拟光的行径，在物体和光源之间找到可行的路径，最后返回积累的颜色。
和光栅化不同，光线追踪的第一步是向世界空间中投射光线，这一步叫做RayCasting，对应位置的像素代表投射光遇到的实体，同时给相应位置的像素一定的颜色：
假设我们发射的光线击中一个点，我们将返回这个点的颜色作为最终像素的颜色。
那么如何描述一点颜色呢？我们介绍渲染方程↓
故事要从元和二年1986从2000年开始，科学家们首次提出渲染方程来描述一个点受到的光的影响。渲染像素的颜色是渲染方程的过程。一般意思是:一个点的光是由2部分是：
注意这个方向，允许这个点从它那里接收法向半球范围内方向的入射光：
如果我们想计算法向半球任何方向的入射光的积累，我们必须计算一个半球积分！以下是渲染方程的简化形式：
L(p)=E(p)∫L(q)cos(θ)dωiL\left(p\right)=E\left(p\right)\int_{}L\left(q\right)*cos(\theta)\\mathrm{d}\omega_{i}L(p)=E(p)+∫L(q)cos(θ)dωi
其中符号的解释如下：
L(x)→x点的光强度E(x)→x点发出的光q→从p点出发,方向为ωi的光命中了q点的物体θ→ωi与q点法向量的夹角\begin{array}{l}L(x)\rightarrowx\点的光强度\\E(x)\rightarrowx\点发出的光\\q\rightarrow从\p\点出发,方向为\\omega_{i}\的光命中了\q\点的物体\\\theta\rightarrow\omega_{i}\与\q\点法向量的夹角\end{array}L(x)→x点的光强度E(x)→x点发出的光q→从p点出发,方向为ωi的光命中了q点的物体θ→ωi与q点法向量的夹角
此外，∫\int_{}∫就是对法向半球入（出）射光线方向的积分。注意到渲染方程的递归形式，要想求解p点的光照我们必须先递归计算q点的光照值，如图：
递归式告诉我们：q点也反射了光到p点，反射光的强度等于q点的光强度乘以反射到p点的百分比cos(θ）
注：这里有个错误，我写的是θ是出射光wi和命中点q的法线的夹角，实际上渲染方程中，θ是出射光wi和出发点p点的法线的夹角，也就是Lambert’scosinelaw，或者说phong光照模型中的NdotL这个错误我在写完博客的半年之后才发现。。。不过因为我们用的是简化版的渲染方程，并非基于物理，即使你去掉这个cosine，也能够得到相对不错的结果，真正正确的伪代码应该是：
民科警告
对于p点的单位面积dp，光能够发射到q点的比率，取决于dp在入射光方向上面的投影面积：
了解了渲染方程的成分，就可以进行渲染方程的求解了。在高等数学中，积分的计算需要找到被积函数的原函数，和积分变量。可是渲染方程是一个困难积分，无法精确地计算其原函数，于是需要找寻其他方法对困难积分进行计算。
那么怎样计算一个困难积分的？我们引入蒙特卡洛方法↓
我更愿意称之为丢豆子。考虑最直观的情况，欲求区间上一函数的积分，我们往x区间上面丢豆子，并且计算豆子命中的位置的y的值，最后把他们加起来作为积分的估计：
丢豆子的过程称之为【采样】，如果我们使用均匀分布的x进行丢豆子，就能得到上图等宽度的柱状图近似。
事实上在实际问题中，豆子的位置不会总是服从均匀分布。那么每一个豆子的贡献，除了豆子命中位置的y值，还取决于豆子命中该位置的概率。蒙特卡洛方法允许我们使用x的任意的概率密度函数，来对积分进行估计。
假设采样x的概率密度函数为PDF(x)PDF(x)PDF(x)，被积函数为f(x)f(x)f(x)，那么x点采样的贡献就是：
f(x)PDF(x)\frac{f(x)}{PDF(x)}PDF(x)f(x)
所以积分计算就很简单了。首先按照产生一堆符合概率密度函数PDF分布的随机变量x，然后对每一个x都计算f(x)PDF(x)\frac{f(x)}{PDF(x)}PDF(x)f(x)最后求他们的均值即可。现在回过头来看刚刚的均匀分布的丢豆子，其中PDF(x)=1baPDF(x)=\frac{1}{b-a}PDF(x)=ba1，那么我们估计x2x^2x2的积分就可以这么计算：
可以看到仅5次采样就可以获得还不错的结果。我们和真实的积分值十分逼近了！采样的次数越多，差异就越少，所以蒙特卡洛方法可以做到对积分结果的无偏估计，这是好特性。
知晓了困难积分的近似求解方式，我们开始正式求解渲染方程↓
渲染方程是对于被求解点p的法向半球的积分，那么我们的随机变量就是在法向半球内的射线的方向，假设我们取法向半球内均匀分布的射线方向，那么就有PDF(x)=12πPDF(x)=\frac{1}{2\pi}PDF(x)=2π1，因为半球面积就是2π2\pi2π，如图：
于是有如下的伪代码：
pathTracing(p){ L=0.0 for(iinSAMPLE) { wi=random() //随机选择一个方向 if(wihitq) //射线wi击中q点 { cosine=dot(nq,wi) //q点法向量nq和wi的夹角余弦值 L+=cosine*pathTracing(q)/PDF(wi) } } returnL/SAMPLE //返回均值}其中SAMPLE就是我们的采样次数。那么问题就来了，这个递归的复杂度是指数，复杂度非常高带来的就是极大的计算资源的消耗，因为光线会有爆炸般的增长，以SAMPLE=100为例：
那么我们只有限制SAMPLE=1才能防止指数增长。而一次采样的估计肯定是不准确的，于是我们对每个像素，发射多条光线，然后平均他们的结果。每个像素的光线数目叫做SPP，即（SamplePrePixel），下图演示了SPP=3的情况，我们找寻了3条到光源的路径，并且平均他们的贡献：
注意只有在第一次采样时发射若干条光线，其余的时候我们只随机选取一个方向发射光线并且递归计算。那么伪代码就改成：
//追踪一条光线pathTracing(p){ L=0.0 wi=random() //随机选择一个方向 if(wihitq) //射线wi击中q点 { cosine=dot(nq,wi) //q点法向量nq和wi的夹角余弦值 L+=cosine*pathTracing(q)/PDF(wi) } returnL}//对一个像素投射SPP条光线L=0.0for(iinSPP){ wi=random() //随机选择一个方向 if(wihitq) //射线wi击中q点 { L+=pathTracing(q)/PDF(wi) }}L/=SPP这一步也没啥特别的，就是向每一个像素投射光线，然后求解渲染方程，没了。。。
渲染方程的伪代码有了，我们通过c++实现它↓
着手编写一个在Windows10下运行的x64程序，程序以图片的形式输出场景的渲染结果。我们以VusialStudio2019作为IDE，此外我们还需要额外的帮助库。
首先是数学运算库，我们需要一个能够表示三维向量，并且对向量进行加减乘除点积叉乘等操作的帮助库。你可以自己写一个简易的class，也可以使用现成的第三方库，这里我使用的是glm，它的网站在这里，你也可以从它的GitHub上面获取。此外，也可以通过vcpkg包管理工具来下载，只需要运行命令：
vcpkginstallglm如果在安装时遇到任何困难，可以参考我以前的博客：传送门①，传送门②
你可以使用任何流行的图像处理的库来进行图像输出，他们可以是Opencv，Qt，甚至是OpenGL，但是这里我使用非常轻量级的svpng。svpng不是一个c++的第三方库，它仅是一个inc文件：
你只需要把它放在你的工程目录下，然后再#include"svpng.inc"即可调用它。svpng就一个非常简单的功能，就可以帮我们保存png图像，调用svpng函数即可。函数的原型长这样：
voidsvpng(FILE*fp,unsignedw,unsignedh,constunsignedchar*img,intalpha)其中FILE是文件指针，w和h是图片的宽度和高度，img是图像的像素值数组，alpha是透明度，我们填0即可。通过如下的代码就可以将一个范围在[0,1]之间的double浮点数buffer输出到图片上：
//输出SRC数组中的数据到图像voidimshow(double*SRC){unsignedchar*image=newunsignedchar[WIDTH*HEIGHT*3];//图像bufferunsignedchar*p=image;double*S=SRC;//源数据FILE*fp;fopen_s(&fp,"image.png","wb");for(inti=0;i 随便在SRC中写点什么，比如输出xy的值作为rg通道。如果你看到如下的图片被生成，那么很成功！ 如果找不到svpng.inc那么检查你的vs工程是否配置正确，将包含svpng的目录添加到vs的include目录： 光线追踪运算量巨大，单靠简单的单线程程序无法高效执行，但是因为每个光线的采样是相互独立的，于是我们可以利用多线程加速。VisualStudio有自带多线程加速的openmp库，无需自己手动下载，只需要引入： #include//openmp多线程加速同时在项目设置中，允许vs使用多线程： 然后在需要并行执行的for循环之前加上： omp_set_num_threads(50);//线程个数#pragmaompparallelforfor(){ ...}就可以享受多线程加速的福利了。此外，我墙裂建议你打开O2优化同时将运行模式调整到Release，以获取最大运行速度： 一切就绪？我们准备进入光与粒的世界↓ 光线追踪的第一步是投射光线，我们模拟相机投影与成像的规则，指定一个[-1,1]范围内的投影平面和一个视点，然后根据输出图片的像素位置，计算其对应投影平面上的坐标，最后用坐标减去视点坐标，得到视线的方向向量，如图： 值得注意的是，图片的xy轴原点是在图片左上方，而实际投影我们需要一个在左下方的原点（即平面几何坐标系），所以y要做一次flip。此外，在世界坐标系下，我们确定相机的位置和投影平面的位置，让相机看向z轴负方向： 相机配置就绪，我们尝试输出相机的射线投射方向，其中imshow是上面编写的显示图片的函数： double*image=newdouble[WIDTH*HEIGHT*3];memset(image,0.0,sizeof(double)*WIDTH*HEIGHT*3);double*p=image;for(inti=0;i 相机投射了光线，光线和场景物体相交，我们希望描述这一过程↓ 在计算机图形学中通常使用三角形来描述任意形状的物体，因为三角形具有很好的几何特征，并且易于进行求交，裁剪等操作。在开始之前，我们先做一些规范化的定义： 假设我们用起点（start）和方向（direction）来描述一个射线： //光线typedefstructRay{vec3startPoint=vec3(0,0,0);//起点vec3direction=vec3(0,0,0);//方向}Ray;在开始编写三角形类之前，我们先确定求交操作到底要返回那些信息： 那么首先我们有表面属性的定义，使用结构体能很好的帮我们组织数据，同时易于拓展： //物体表面材质定义typedefstructMaterial{boolisEmissive=false;//是否发光vec3normal=vec3(0,0,0);//法向量vec3color=vec3(0,0,0);//颜色}Material;然后是光线求交结果的定义： //光线求交结果typedefstructHitResult{boolisHit=false;//是否命中doubledistance=0.0f;//与交点的距离vec3hitPoint=vec3(0,0,0);//光线命中点Materialmaterial;//命中点的表面材质}HitResult;然后是三角形的class的定义： classShape{public:Shape(){}virtualHitResultintersect(Rayray){returnHitResult();}};//三角形classTriangle:publicShape{public:Triangle(){}Triangle(vec3P1,vec3P2,vec3P3,vec3C){p1=P1,p2=P2,p3=P3;material.normal=normalize(cross(p2-p1,p3-p1));material.color=C;}vec3p1,p2,p3;//三顶点Materialmaterial;//材质//与光线求交HitResultintersect(Rayray){HitResultres;//...returnres;};};这里我墙裂建议使用虚函数+指针+继承的编程习惯，因为我们光线和任意图形求交，都有一致的返回结果，即HitResult结构体。我们使c++的指针特性，可以通过一套代码，完成多种复杂图元的求交。此外，在添加一种新图元的时候，主代码不需要任何的改动！（虽然现在我们只有三角形。。。 和三角形的求交分为两个步骤：首先判断光线和三角形所在平面是否有交点，然后再判断交点是否在三角形内部。思路很清晰，我们先来看光线与三角形所在平面相交的判断： 其中t表示了射线起点到交点的距离，如果t小于0那么表示三角形在摄像机之后！，然后开始判断点是否在三角形中。我们连接顶点与P点，然后判断连线与边的叉乘方向是否与法向量一致。如果三个顶点的判断都通过，说明P在三角形中，否则不在：