这两天微信打算换个头像，所以打算弄个小渲染器来渲染新头像。。。以下是结果:

贴几张不同渲染滴的图片：

在传统的计算机图形学中，光栅通常用于生成像素。因为光栅化后，我们只选择我们看到的图形，视野外的几何信息被切割，导致场景的整体信息丢失。然后它会导致一些渲染 例如，水只能反映我们在屏幕上看到的像素：

此外，一些基于全球信息的渲染效果不能很好地运行，如全球光。在光栅管道下，我们通常需要通过各种类似的算法来模拟现实世界的光影现象，这是非常麻烦和有效的 “还凑合”

超越传统图形流水线的光线跟踪是一种 现代 渲染。光跟踪是一种基于物理的渲染方法，通过模拟介质中光的真实性能来输出足够接近现实的图像。光从光源开始，通过各种反射折射进入相机：

从光源开始 无数条 到达相机的路线不太现实，从灯泡开始模拟光的传播，所以曲线拯救了国家，从相机的方向开始向场景投射射线，试图找出到达光源的可能路径：

本博客代码主要实现路径跟踪技术。路径跟踪技术是光跟踪的一个分支。通过将光投射到场景中，模拟光的行径，在物体和光源之间找到可行的路径，最后返回积累的颜色。

和光栅化不同，光线追踪的第一步是向世界空间中投射光线，这一步叫做 RayCasting，对应位置的像素代表投射光遇到的实体，同时给相应位置的像素一定的颜色：

假设我们发射的光线击中一个点，我们将返回这个点的颜色作为最终像素的颜色。

那么如何描述一点颜色呢？我们介绍渲染方程 ↓

故事要从 元和二年 1986 从2000年开始，科学家们首次提出渲染方程来描述一个点受到的光的影响。渲染像素的颜色是渲染方程的过程。一般意思是:一个点的光是由 2 部分是：

1. 这是自己发出的光
2. 来自任何方向的光击中该点积累形成的光

注意这个方向，允许这个点从它那里接收 法向半球 范围内方向的入射光：

如果我们想计算法向半球任何方向的入射光的积累，我们必须计算一个 半球积分！ 以下是渲染方程的简化形式：

L ( p ) = E ( p ) ∫ L ( q ) ? c o s ( θ ) d ω i L\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) \rightarrow x \ 点的光强度 \\ E(x) \rightarrow x \ 点发出的光 \\ 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’s cosine law ，或者说 phong 光照模型中的 N dot L 这个错误我在写完博客的半年之后才发现。。。不过因为我们用的是简化版的渲染方程，并非基于物理，即使你去掉这个 cosine，也能够得到相对不错的结果，真正正确的伪代码应该是：

至于这里为啥使用错误的图，得到的图片仍然在一定程度上符合物理定律呢？我猜测是这样的：

民科警告⚠

对于 p 点的单位面积 dp，光能够发射到 q 点的比率，取决于 dp 在入射光方向上面的投影面积：

了解了渲染方程的成分，就可以进行渲染方程的求解了。在高等数学中，积分的计算需要找到被积函数的原函数，和积分变量。可是渲染方程是一个困难积分，无法精确地计算其原函数，于是需要找寻其他方法对困难积分进行计算。

那么怎样计算一个困难积分的？我们引入蒙特卡洛方法 ↓

我更愿意称之为丢豆子。考虑最直观的情况，欲求区间上一函数的积分，我们往 x 区间上面丢豆子，并且计算豆子命中的位置的 y 的值，最后把他们加起来作为积分的估计：

丢豆子的过程称之为【采样】，如果我们使用 均匀分布 的 x 进行丢豆子，就能得到上图等宽度的柱状图近似。

事实上在实际问题中，豆子的位置不会总是服从均匀分布。那么每一个豆子的贡献，除了豆子命中位置的 y 值，还取决于豆子 命中该位置的概率。蒙特卡洛方法允许我们使用 x 的任意的概率密度函数，来对积分进行估计。

假设采样 x 的概率密度函数为 P D F ( x ) PDF(x) PDF(x)，被积函数为 f ( x ) f(x) f(x)，那么 x 点采样的贡献就是：

f ( x ) P D F ( x ) \frac{f(x)}{PDF(x)} PDF(x)f(x)​

所以积分计算就很简单了。首先按照产生一堆符合概率密度函数 PDF 分布的随机变量 x，然后对每一个 x 都计算 f ( x ) P D F ( x ) \frac{f(x)}{PDF(x)} PDF(x)f(x)​ 最后求他们的 均值 即可。现在回过头来看刚刚的均匀分布的丢豆子，其中 P D F ( x ) = 1 b − a PDF(x) = \frac{1}{ b-a} PDF(x)=b−a1​，那么我们估计 x 2 x^2 x2 的积分就可以这么计算：

可以看到仅 5 次采样就可以获得还不错的结果。我们和真实的积分值十分逼近了！采样的次数越多，差异就越少，所以蒙特卡洛方法可以做到对积分结果的 无偏 估计，这是好特性。

知晓了困难积分的近似求解方式，我们开始正式求解渲染方程 ↓

渲染方程是对于被求解点 p 的法向半球的积分，那么我们的随机变量就是在法向半球内的射线的方向，假设我们取法向半球内 均匀分布的射线 方向，那么就有 P D F ( x ) = 1 2 π PDF(x) = \frac{1}{2\pi} PDF(x)=2π1​，因为半球面积就是 2 π 2\pi 2π，如图：

于是有如下的伪代码：

pathTracing(p)
{ 
        
	L = 0.0
	for(i in SAMPLE)
	{ 
        
		wi = random()	// 随机选择一个方向
		if(wi hit q)	// 射线 wi 击中 q 点
		{ 
        
			cosine = dot(nq, wi)	// q 点法向量 nq 和 wi 的夹角余弦值
			L += cosine * pathTracing(q) / PDF(wi)
		}
	}
	return L / SAMPLE	// 返回均值
}

其中 SAMPLE 就是我们的采样次数。那么问题就来了，这个递归的复杂度是指数，复杂度非常高带来的就是极大的计算资源的消耗，因为光线会有爆炸般的增长，以 SAMPLE=100 为例：

那么我们只有限制 SAMPLE=1 才能防止指数增长。而一次采样的估计肯定是不准确的，于是我们对每个像素，发射多条光线，然后平均他们的结果。每个像素的光线数目叫做 SPP，即（Sample Pre Pixel），下图演示了 SPP=3 的情况，我们找寻了 3 条到光源的路径，并且平均他们的贡献：

注意只有在第一次采样时发射若干条光线，其余的时候我们只随机选取一个方向发射光线并且递归计算。那么伪代码就改成：

// 追踪一条光线
pathTracing(p)
{ 
        
	L = 0.0
	wi = random()	// 随机选择一个方向
	if(wi hit q)	// 射线 wi 击中 q 点
	{ 
        
		cosine = dot(nq, wi)	// q 点法向量 nq 和 wi 的夹角余弦值
		L += cosine * pathTracing(q) / PDF(wi)
	}

	return L
}

// 对一个像素投射 SPP 条光线
L = 0.0
for(i in SPP)
{ 
        
	wi = random()	// 随机选择一个方向
	if(wi hit q)	// 射线 wi 击中 q 点
	{ 
        
		L += pathTracing(q) / PDF(wi)
	}
}
L /= SPP

这一步也没啥特别的，就是向每一个像素投射光线，然后求解渲染方程，没了。。。

渲染方程的伪代码有了，我们通过 c++ 实现它 ↓

着手编写一个在 Windows 10 下运行的 x64 程序，程序以图片的形式输出场景的渲染结果。我们以 Vusial Studio 2019 作为 IDE，此外我们还需要额外的帮助库。

数学运算库

首先是数学运算库，我们需要一个能够表示 三维向量，并且对向量进行加减乘除点积叉乘等操作的帮助库。你可以自己写一个简易的 class，也可以使用现成的第三方库，这里我使用的是 glm，它的网站在 这里，你也可以从它的 GitHub 上面获取。此外，也可以通过 vcpkg 包管理工具来下载，只需要运行命令：

vcpkg install glm

如果在安装时遇到任何困难，可以参考我以前的博客：传送门①，传送门②

图像输出

你可以使用任何流行的图像处理的库来进行图像输出，他们可以是 Opencv，Qt，甚至是 OpenGL，但是这里我使用非常轻量级的 svpng。svpng 不是一个 c++ 的第三方库，它仅是一个 inc文件：

你只需要把它放在你的工程目录下，然后再 #include "svpng.inc" 即可调用它。svpng 就一个非常简单的功能，就可以帮我们保存 png 图像，调用 svpng 函数即可。函数的原型长这样：

void svpng(FILE* fp, unsigned w, unsigned h, const unsigned char* img, int alpha)

其中 FILE 是文件指针，w 和 h 是图片的宽度和高度，img 是图像的像素值数组，alpha 是透明度，我们填 0 即可。通过如下的代码就可以将一个范围在 [0, 1] 之间的 double 浮点数 buffer 输出到图片上：

// 输出 SRC 数组中的数据到图像
void imshow(double* SRC)
{ 
        
    
    unsigned char* image = new unsigned char[WIDTH * HEIGHT * 3];// 图像buffer
    unsigned char* p = image;
    double* S = SRC;    // 源数据

    FILE* fp;
    fopen_s(&fp, "image.png", "wb");

    for (int i = 0; i < HEIGHT; i++)
    { 
        
        for (int j = 0; j < WIDTH; j++)
        { 
        
            *p++ = (unsigned char)clamp((*S++) * 255, 0.0, 255.0);  // R 通道
            *p++ = (unsigned char)clamp((*S++) * 255, 0.0, 255.0);  // G 通道
            *p++ = (unsigned char)clamp((*S++) * 255, 0.0, 255.0);  // B 通道
        }
    }

    svpng(fp, WIDTH, HEIGHT, image, 0);
}

clamp 是截断函数，glm 库带的，如果报错那么您可以删掉它并且换成您自己的 clamp，只是时刻注意 SRC 是 [0, 1] 范围的 double，我们习惯这么表示颜色，同时方便计算，不容易被截断。此外，svpng 默认图像的 RGB 通道是相邻的，我们直接利用指针进行遍历即可。

随便在 SRC 中写点什么，比如输出 xy 的值作为 rg 通道。如果你看到如下的图片被生成，那么很成功！

如果找不到 svpng.inc 那么检查你的 vs 工程是否配置正确，将包含 svpng 的目录添加到 vs 的 include 目录：

多线程加速

光线追踪运算量巨大，单靠简单的单线程程序无法高效执行，但是因为 每个光线的采样是相互独立的，于是我们可以利用多线程加速。Visual Studio 有自带多线程加速的 openmp 库，无需 自己手动下载，只需要引入：

#include <omp.h> // openmp多线程加速

同时在项目设置中，允许 vs 使用多线程：

然后在需要并行执行的 for 循环之前加上：

omp_set_num_threads(50); // 线程个数
#pragma omp parallel for

for()
{ 
        
	...
}

就可以享受多线程加速的福利了。此外，我墙裂建议你打开 O2 优化同时将运行模式调整到 Release，以获取最大运行速度：

一切就绪？我们准备进入光与粒的世界 ↓

光线追踪的第一步是投射光线，我们模拟相机投影与成像的规则，指定一个 [-1, 1] 范围内的投影平面和一个视点，然后根据输出图片的像素位置，计算其对应投影平面上的坐标，最后用坐标减去视点坐标，得到 视线的方向向量，如图：

值得注意的是，图片的 xy 轴原点是在图片左上方，而实际投影我们需要一个在左下方的原点（即平面几何坐标系），所以 y 要做一次 flip。此外，在世界坐标系下，我们确定相机的位置和投影平面的位置，让相机看向 z 轴负方向：

相机配置就绪，我们尝试 输出相机的射线投射方向，其中 imshow 是上面编写的显示图片的函数：

double* image = new double[WIDTH * HEIGHT * 3];
memset(image, 0.0, sizeof(double) * WIDTH * HEIGHT * 3);
double* p = image;
for (int i = 0; i < HEIGHT; i++)
{ 
        
    for (int j = 0; j < WIDTH; j++)
    { 
        
        // 像素坐标转投影平面坐标
        double x = 2.0 * double(j) / double(WIDTH) - 1.0;
        double y = 2.0 * double(HEIGHT - i) / double(HEIGHT) - 1.0;

        vec3 coord = vec3(x, y, SCREEN_Z);          // 计算投影平面坐标
        vec3 direction = normalize(coord - EYE);    // 计算光线投射方向
        
        vec3 color = direction;

        *p = color.x; p++;  // R 通道
        *p = color.y; p++;  // G 通道
        *p = color.z; p++;  // B 通道
    }
}

imshow(image);

输出如下结果说明相机有在工作：

相机投射了光线，光线和场景物体相交，我们希望描述这一过程 ↓

在计算机图形学中通常使用三角形来描述任意形状的物体，因为三角形具有很好的几何特征，并且易于进行求交，裁剪等操作。在开始之前，我们先做一些规范化的定义：

结构定义

假设我们用起点（start）和方向（direction）来描述一个射线：

// 光线
typedef struct Ray
{ 
        
    vec3 startPoint = vec3(0, 0, 0);    // 起点
    vec3 direction = vec3(0, 0, 0);     // 方向
}Ray;

在开始编写三角形类之前，我们先确定求交操作到底要返回那些信息：

1. 是否相交
2. 交点位置，用于作为我们下一次弹射的起点
3. 相交位置的表面属性：比如法向量，表面颜色，材质属性，发光度，粗糙度等

那么首先我们有表面属性的定义，使用结构体能很好的帮我们组织数据，同时易于拓展：

// 物体表面材质定义
typedef struct Material
{ 
        
    bool isEmissive = false;        // 是否发光
    vec3 normal = vec3(0, 0, 0);    // 法向量
    vec3 color = vec3(0, 0, 0);     // 颜色
}Material;

然后是光线求交结果的定义：

// 光线求交结果
typedef struct HitResult
{ 
        
    bool isHit = false;             // 是否命中
    double distance = 0.0f;         // 与交点的距离
    vec3 hitPoint = vec3(0, 0, 0);  // 光线命中点
    Material material;              // 命中点的表面材质
}HitResult;

然后是三角形的 class 的定义：

class Shape
{ 
        
public:
    Shape(){ 
        }
    virtual HitResult intersect(Ray ray) { 
         return HitResult(); }
};

// 三角形
class Triangle : public Shape
{ 
        
public:
    Triangle(){ 
        }
    Triangle(vec3 P1, vec3 P2, vec3 P3, vec3 C) 
    { 
         
        p1 = P1, p2 = P2, p3 = P3; 
        material.normal = normalize(cross(p2 - p1, p3 - p1)); material.color = C;
    }
    vec3 p1, p2, p3;    // 三顶点
    Material material;  // 材质

    // 与光线求交
    HitResult intersect(Ray ray) 
    { 
         
        HitResult res;
        // ...
        return res; 
    };
};

这里我墙裂建议使用 虚函数+指针+继承 的编程习惯，因为我们光线和任意图形求交，都有一致的返回结果，即 HitResult 结构体。我们使 c++ 的指针特性，可以通过一套代码，完成多种复杂图元的求交。此外，在添加一种新图元的时候，主代码不需要任何的改动！（虽然现在我们只有三角形。。。

求交计算

和三角形的求交分为两个步骤：首先判断光线和三角形所在平面是否有交点，然后再判断交点是否在三角形内部。思路很清晰，我们先来看光线与三角形所在平面相交的判断：

其中 t 表示了射线起点到交点的距离，如果 t 小于 0 那么表示三角形在摄像机之后！ ，然后开始判断点是否在三角形中。我们连接顶点与 P 点，然后判断连线与边的叉乘方向是否与 法向量 一致。如果三个顶点的判断都通过，说明 P 在三角形中，否则不在：

注意此处法向量 N 是通过 p 1 p 2 ⃗ × p 1 p 3 ⃗ \vec{p_1p_2} \times \vec{p_1p_3} p1​p2​