本文为我自己的笔记，有错误希望指出，不胜感激。
大约十月头上，因为返校延期，我打算写一个建议的PBR光追器。
在此之前，我也就看过Games101, Ray tracing week系列的前两本，还啃过一点PBRT，所以我也只是个萌新罢了。
我的实现是基于Ray tracing in one weekend那个迷你光追器。当然，我对架构来了一场大改，断断续续大概有一个半月，直到现在才做好。整个过程磕磕绊绊，大部分时间花在了查资料，改Bug上。
如今进行复盘，发现当时查的很多都是无用资料，有些甚至起到了误导作用。我相信很多和我一样的图形学萌新估计都有这样的经历，所以我把我的学习过程和找到的资料记录在这，希望后来人能轻松一点吧。
在实现的过程里，我遇到最棘手的bug是由NaN导致的。如果你还不熟悉IEEE的浮点数规定，你可能会一头雾水；或者你接触过javaScript，会知道这是Not a Number的意思。
可以做一个简单的小实验，进行以下的计算，然后让C++输出
1/0
1.0/0
1+1.0/0
sqrt(-1.0)
acos(2)
1+sqrt(-1.0)
我直接给出结论
-1.未定义行为，输出什么都有可能
-2.+inf
-3.nan
-4.nan
-5.nan
整数除零是未定义行为，是必须避免的。
浮点数除以0等于INF，这个是IEEE754的规定，不仅有+inf还有-inf
如果进行了错误的运算，比如对-1开方，取2的arccos，就会得到nan
任何数字加上nan都等于nan，这就是nan的传递性了。当然，你也可以对inf做下实验，比如inf+nan，或者inf/inf。
C在头文件<math.h>里提供了检测inf和nan的宏：isnan()，isinf()
但不论哪种输出，都是需要极力避免的，有时候一个nan、inf就能毁了一整张渲染图，还附赠一个难以发现的bug
为了防止这些恼人因素的出现，你需要考虑到各种边界情况，甚至是浮点数的尾数误差....
我有一个简单粗暴的预防方式：断言，assert。这个函数在<assert.h>里，如果assert里面的值为1，就继续执行，如果为0，就让程序crash。
所以如果是nan/inf，就让它crash了吧，早点crash胜过事后打断点。
于是你可以这么写
float a = sqrt(-1.);
assert(!isnan(a));//Wow,it crashed here.
不止nan，任何可能发生与预期不一致行为的地方，都应该加上assert，如果你去看看pbrt的源代码，会发现里面assert满天飞。
相信我，这对你的debug大有好处
首先我们需要达成一点共识，比如渲染方程的形式如下,而且你知道这些符号的含义
$$\int_{\Omega}f(w_i\rightarrow w_o)Li(w_i)(w_n\cdot w_i)dw_i$$
众所周知，解渲染方程是通过蒙特卡洛方法来近似，蒙特卡洛的形式如下
$$\frac{1}{N}\sum_{i=1}\frac{F(x_i)}{p(x_i)}$$
我对蒙特卡洛有个极其简单的理解，那就是加权的黎曼积分，p(x)就是那个权重
先从一个最简单的积分开始$\int_{0}^{2}\frac{1}{2}xdx$，这个积分的结果是$1$，毫无疑问。但是这里让我用蒙特卡洛的方法来做一下。
首先，函数p(x)是什么？概率密度函数（后文以pdf指代）。
在这里，积分可以理解为求面积。可以把函数的积分面积进行一个划分
如果随机在三角形面积上选点，那么面积大的划分块上，点会更多，在面积小的划分块上，点会比较少。此处的pdf就可以理解为点出现在对应划分上的概率...
如果把划分的区间无限缩小，那么pdf就会变成一个连续的函数。
于是用 (1/pdf) 作为权重，相当于做了一个归一化，把过大的划分面积缩小，过小的划分面积放大，那么结果就会趋近于均值，也就是要求的积分值（此处也就是面积）
于是在积分域创建一个随机数序列x，然后带入蒙特卡洛的公式，就能求出来这个积分了
进行一个简单的c++模拟，那么应该是形如这样的代码
float monteCarlo(int N){
float res = 0;
for (int i = 0; i < N; i++) {
float x = mrand(0,2);
res += f(x)/pdf(x);
}
res = res/N;
return res;
}
此处f(x)就是0.5x，x是在积分域0，2上的随机数，那么pdf是多少？
这里暂时认为它是个均匀分布，那么均匀分布的pdf就是区间长度的倒数，也就是一个常数，0.5
运行一下这个函数，会发现当迭代的次数越多，这个积分结果就越准确，因为其本身简单，所以迭代了10次差不多也收敛了。
那换一个更复杂的函数呢？求 $y=0.5x^2$ 在0，2上的积分，身为接受过高等训练的智慧直立两脚兽，我们很容易可以算出其结果是4/3。
但电脑不这么想，在我这，迭代了500次的结果是1.32563，迭代了10000次的结果是1.36147。很明显，哪里错了。
问题出在pdf的选择上，因为pdf选的是均匀分布，而目标积分是个二次函数，两者的形状相差太大了，因而pdf不能很好的描述点在图形上的分布。
那么这里设计一个pdf，$p(x)=cx^2$
$$\int^{2}_{0}cx^2dx=1$$
$$c=\frac{3}{8}$$
$$p(x)=\frac{3}{8}x^2$$
于是我这里凑了个pdf出来：$\frac{3}{8}x^2$，把其替换掉pdf，再运行一下
在我这里，迭代10次的结果是1.33333，非常完美的结果
概率论小知识：PDF在全域上的积分为1
于是可以得到一个结论，如果pdf的分布符合被积函数的形状，那么这个蒙特卡洛就会快速收敛。
之前举的例子实际上都不难，因此不论pdf选什么收敛都蛮快的。
但是在图形学中，我们面对的是更加棘手的函数：渲染方程。在大部分的积分域上，其的值都是0（因为只有当射线击中了光源才有值）。
这时候，我们希望光线的出射方向能带来大量的贡献。因为pdf大的地方，其函数值也一定大，所以我们会想要对pdf大的方向进行采样（这也就是所谓重要性）。于是怎么让采样的方向符合pdf的形状，就是重要性采样的问题。
让出射光线直接打中光源？
详见我以前的这篇博客(在半球上均匀采样的推导)，本文主要利用逆采样进行采样。
要求蒙特卡洛，要采样，都先要得到一个pdf。
观察渲染方程，发现里面的brdf项和Li项都难以下手，先从简单的cos项开始。
$$\int_{\Omega^+}dot(w_i,w_g)dw=1$$
$$\int_{0}^{2\pi}\int_{0}^{\frac{\pi}{2}}cCos(\theta)Sin(\theta)d\theta d\phi=1$$
$$c=\frac{1}{\pi}$$
所以我们得到了$p(w)=\frac{dot(w_i,w_g)}{\pi}$.
接下来先求其边缘概率密度和联合概率密度
$$p(\theta,\phi)=\frac{cos\theta sin\theta}{\pi}$$
$$p(\theta)=2sin\theta cos\theta$$
$$p(\phi|\theta)=\frac{1}{2\pi}$$
求其分布函数
$$P(\phi|\theta)=\int_{0}^{\phi}\frac{1}{2\pi}d\phi=\frac{\phi}{2\pi}$$
$$P(\theta)=\int_{0}^{\theta}sin2\theta d\theta=-0.5cos2\theta+0.5=sin^2\theta=1-cos^2\theta$$
注：这里采用cos的原因是asin的返回值在(-0.5pi,0.5pi)之间，不方便计算
接着逆采样
$$\phi=2\pi \xi_1$$
$$\theta=cos^{-1}\sqrt{1-\xi_2}$$
其实还有一种Cos-weighted采样方法，pbrt里称为Malley's method，其思想就是先在圆盘上均匀采样，然后投影到半球上
对圆盘采样是用极坐标完成的，$r=\sqrt{\xi_1},\phi=2\pi\xi_2$，那么将其投影到半球上就是一个三维坐标，$(rcos(\phi),1-r^2,rsin(\phi))$
至于推导么....画个图看看
这个结果的正确性可以通过计算pdf来验证，在此不多计算。
微表面曾经是个光学概念，它认为在宏观世界的一块表面其实是由无数个微小的形状不同的表面组成的。
Games101举过一个例子，澳大利亚在国际空间站上看着能反光！（反光那块是陆地
最常用的微表面模型是Cook-torrance brdf
$$f(w_o)=\frac{D(w_h)G(w_o,w_i)F(w_i,w_h)}{4dot(w_g,w_i)dot(w_g,w_o)}$$
其中D是法线分布函数，G是shadowing-masking function,F是菲涅尔函数。
这个模型的推导放在后面
假设能看到的表面是宏表面，宏表面的法线也就是几何法线$w_g$，约定宏表面的面积为1；宏表面是由无数微表面构成的，每块微表面都有独立的法线$w_m$。
spatial space to statistical space
考虑从一块微表面，上面有许许多多更小更小的微微表面，我们认为这些微微表面足够光滑，并且可以视为一个点pm，那么对于每个点，其都有一个法线$w$,求此点的法线的函数就是$w(p)$，这就是空间维度。在空间维度，我们计算面积的方式很简单，可以理解为数点，每个微表面上有几个点，加起来，那么就是微表面的面积
$$microsurfacearea=\int_{\Mu}dp_m$$
但是我们不可能准确知道微表面上每个点，这时候就要将其转换到统计学维度上看，这个转换的桥梁就是NDF。
定义函数$D(w)$用来描述微表面上法线的分布（后称NDF Normal distribution function），其定义式为
$$D(w)=\int_{\Mu}\delta_w(\omega_m(p_m))dp_m$$
M代表积分域是整个微表面，其中，D(w)的单位是$\frac{m^2}{sr}$，sr来自delta函数中的参数的倒数，$\omega_m(p_m)$的意思是求微点pm在微表面上的法线，这个dleta函数的含义是：只有当$w_m(p_m)$朝向w方向时，其值为1，其余时刻为0。
那么这个积分就是对于微表面上所有 法线朝向w 的 点 进行一个统计。
假如对于所有法线可能朝向的方向进行这个统计，毫无疑问我们会得到一张直方图，横轴是法线朝向的方向，纵轴是以该方向作为法线的点的出现次数（也就是这个法线的出现次数）。把这个离散的直方图连续化，那么就是一个 类似 概率密度函数的存在（当然，只是类似，因为其积分不为1）。
于是求微表面面积就可以这么列出：
$$microsurfacearea=\int_{\Omega}D(w_i)dw_i$$
这就相当于把微表面上的点投影到了该点法线方向的球面上求和，也就等于微表面的面积。
$$microsurfacearea=\int_{\Mu}dp_m=\int_{\Omega}D(w_i)dw_i$$
于是我们完成了从空间维度到统计学维度的一个转换，而统计学维度的公式就是我们能真正使用的了。
我之前说，NDF只是一个类似概率密度的存在，这是因为，其在积分域上的积分并不为1。从直观的方面来理解，就是微表面的面积随着不同的微表面而变化，因此上式求微表面面积积分不为1，NDF不是概率密度函数
为了使其成为一个真正的概率密度函数，需要做一个归一化，这就引出了投影面积。
在这之前，我们一直是在微观的层面进行计算，这里需要引入一个更大的层面，宏表面，将其大小定义为1（$m^2$）。宏表面就是我们能直观看到的表面，它是由微表面构成，而微表面由点构成。
那么，所有微表面在宏表面上的投影面积，就是宏表面的大小，也就是1。
把微表面的面积投影到了宏表面上，于是存在这么一个关系：
$$\int_{\Omega}(\omega_m\cdot w_g)D(w_m)dw_m = 1m^2$$
$$\int_{M}(w_m(p_m)\cdot w_g)dp_m = 1m^2$$
其中，wg是宏表面的法线。
Masking function and Shadowing function
法线分布函数是描述了在单位宏表面上的法线朝向某个立体角的概率。
但是不是所有法线都能被看到，如果从某个稍微偏一些的角度去看这块微表面，就会发现某些表面被微表面上的起伏遮挡了
发生这种问题的根本原因，是NDF项只是给出了法线分布的概率，而并没有约束微表面的排列，于是如图中所示，同一个NDF可以用在不同的微表面上，但是反射光线的分布大不一样。
这种从视点看去发生的微表面自遮挡被称为Masking，于是我们需要一个Masking function来建模这种自遮挡。这个函数被写作G1。
首先先在空间维度来看，假设$G1=(w_o,p_m)$是一个二值的函数，从wo看pm点没被遮挡那么为1，反之为0；那么就能得到投影在wo方向上的可见微表面面积
$$ProjArea=\int_{M}G1(w_o,p_m)Dot_{clamp}(w_o,w_m(p_m))dp_m$$
注：此处的$Dot_{clamp}$是把点乘的结果裁剪到了(0,1)，因为此处不考虑背向于wo的微表面的面积
同样的，把这个式子转换到统计空间，那么就是
$$ProjArea=\int_{\Omega}G1(w_o,w_m)Dot_{clamp}(w_o,w_m)D(w_m)dw_m$$
此时这个投影面积是可以确定的。
首先考虑一个特殊情况，当观察方向和法线一致情况的时候，如果以上面图片中的微表面为例，那么所有的微表面都是可见的，即G1永远为1，那么这个式子的结果就是微表面的面积在法线方向上的投影，即微表面在宏表面上的投影面积1m^2。如果缓缓转动wo，那么此时在wo上的投影面积就会慢慢变小直到0，这个面积大小就是$dot(w_o,w_g)$。
$$dot(w_o,w_g)=\int_{\Omega}G1(w_o,w_m)Dot_{clamp}(w_o,w_m)D(w_m)dw_m$$
这就是对G1项的一个约束，当然想要求出这其表达式还远远不够，还需要微表面的性质（microsurface profile）进一步进行约束。
PBRT里只提到了V-Cavity模型，Smith模型。这两种遮蔽函数是符合物理的
当然最常用的还是Smith模型。
Smith G1的形式如下
$$G1=\frac{X^+(dot(w_o,w_m))}{1+\Lambda(w_o)}$$
其中的$X^+(x)$是个二值的函数，当x<0的时候值为0，当x>0的时候值为1。
$\Lambda(w_o)$是对微表面坡度的积分(an integral over the slopes of the microsurface)，这部分涉及到了坡度空间这个概念，会在对可见微表面采样这块提到。
SmithG1一般写作为
$$G1=\frac{1}{1+\Lambda(w_o)}$$
在读其他pbr文章或者论文时，经常会见到没头没脑的一句话：微表面是一个高度场。千万不要把这句话想的太复杂，他的意思只是说，微表面是一块有高低起伏的表面而已。
既然有高低起伏，那就要考虑到这些高低起伏是怎么分布的了。
Smith遮蔽函数认为这些微表面互相不相关，不连续，就像下图一样。
Masking Function只考虑了从视线往微表面看去被遮挡的情况，忽略了从光源处向微表面被遮挡，这种现象被称为Shadowing，不难发现Shadowing只是把Masing的视点位置换到了光源位置，这时候需要一个联合遮蔽-阴影函数G2同时描述这两种遮挡。
这里只记录Smith联合遮蔽-阴影函数。
Smith联合遮蔽-阴影函数分成五种，而只有前两种常用，所以我也只会介绍这两种。
分离遮蔽阴影型认为Shadowing和Masking互不相关，其形式为
$$G2(w_i,w_o)=G1(w_i)G1(w_o)$$
这个G2是最简单的，当然也是最不准确的。
高度相关遮蔽阴影型考虑了微表面高度导致的遮蔽和阴影相关性。其形式为
$$G2=\frac{X^+(dot(w_o,w_m))X^+(dot(w_i,w_m))}{1+\Lambda(w_o)+\Lambda(w_i)}$$
这个版本的G2最常用。
微表面法线函数定义了法线在微表面上的分布，自然对这个函数采样我们就能得到符合微表面的一个法线。
因前文提到以下公式。
$$\int_{\Omega}(\omega_m\cdot w_g)D(w_m)dw_m = 1m^2$$
所以可以把下式作为目标PDF并且进行采样。
$$p(w_m)=(\omega_m\cdot w_g)D(w_m)$$
Beckmann分布
各项同性版本：
$$D(w_h)=\frac{e^{tan^2\theta/\alpha^2}}{\pi\alpha^2cos^4\theta}$$
各项异性版本：
$$D(w_h)=\frac{e^{-tan^2\theta}((\frac{cos\phi}{a_u})^2+(\frac{sin\phi}{\alpha_v})^2)}{\pi\alpha_u\alpha_vcos^4\theta}$$
GGX (aka.Trowbridge-Reitz)分布
各项同性版本： $$D(w_h)=\frac{a^2}{\pi(cos^2\theta(\alpha^2 - 1)+1)^2}$$ 各项异性版本： $$D(w_h)=\frac{1}{\pi\alpha_u\alpha_vcos^4\theta(1+tan^2\theta(\frac{cos^2\phi}{a_u^2}+\frac{sin^2\phi}{a_v^2}))^2}$$
这里有必要说明一下，$\theta$和$\phi$是$w_h$在以宏表面法线为上轴的法线空间中的极坐标参数。这个可以参考PBRT Ch.8 Reflection Models 的实现。
在实现时，可以一律使用各项异性的版本（各项同性也是特殊的各项异性版本嘛）。
基于此，只要算出其各项异性版本的逆采样就行捏。
这里继续给出Beckmann和GGX各项异性的逆采样结果
Beckmann：
$$ \phi=\left\{
\begin{matrix}
arctan(\frac{\alpha_v}{\alpha_u}tan(2\pi\xi_1)) &\xi_1\in[0,0.25] \\
arctan(\frac{\alpha_v}{\alpha_u}tan(2\pi\xi_1))+\pi &\xi_1\in(0.25,0.75)\\
arctan(\frac{\alpha_v}{\alpha_u}tan(2\pi\xi_1))+2\pi &\xi_1\in[0.75,1]
\end{matrix}
\right.
$$
$$\theta=arctan(\frac{ln\xi_2}{\frac{cos^2\phi}{a_u^2}+\frac{sin^2\phi}{a_v^2}})$$
GGX：
$$ \phi=\left\{
\begin{matrix}
arctan(\frac{\alpha_v}{\alpha_u}tan(2\pi\xi_1)) &\xi_1\in[0,0.25] \\
arctan(\frac{\alpha_v}{\alpha_u}tan(2\pi\xi_1))+\pi &\xi_1\in(0.25,0.75)\\
arctan(\frac{\alpha_v}{\alpha_u}tan(2\pi\xi_1))+2\pi &\xi_1\in[0.75,1]
\end{matrix}
\right.
$$
$$\theta=arctan(\sqrt{ \frac{\xi_2}{ (1-\xi_2) (\frac{cos^2(\phi)}{\alpha_u^2}+\frac{sin^2(\phi)}{a_v^2}) } } )$$
这里的$\phi$后跟着一大串式子，这是因为arctan取值范围为$(-\frac{\pi}{2},\frac{\pi}{2})$，而$\phi$的取值范围为$(0,2\pi)$
这里就需要对值做一个重映射，把其值映射到$(0,2\pi)$。
ps:你想让他取值范围为$(-\pi,\pi)$也莫问题，毕竟sin cos一算一个样
在pbrt的代码中，这个$\phi$的映射与此处写法不同，但实际上莫有差别。
至于这个式子么...当然不是我推的，可以去看看这位大佬的博客，上面有详细的过程：Sampling Anisotropic Microfacet BRDF - A GRAPHICS GUY'S NOTE
最后，我们只要把求出的$\theta$和$\phi$转化为光线出射的方向，再算出PDF的值，就能完成对微表面BRDF的重点采样了。
那么啥不连带着菲涅尔项和遮蔽函数一起逆采样了？好问题，我不知道
坡度空间，简单来说，就是描述三维表面上梯度变化的一个空间。对于定义在高度场中的微表面模型，这种表示方式雀氏很合适。
可以这么理解坡度空间：假设存在一个平面，其法线为n，存在切向量$(\Delta x,\Delta y,\Delta z)$
把法线的起点移动到(0,0,0)，并将其向外延长，与平面z=1相交。这时候，交点的x与y分别代表了原来平面的$-\frac{\delta z}{\delta x}$和$-\frac{\delta z}{\delta x}$，也就是负的梯度。这点在梯度空间中记为$(\tilde{v}_x,\tilde{v}_y)$
图源：Slope Space in BRDF Theory - Nathan Reed 可视化坡度空间
一般使用~上标表示坡度空间中的值
这里给出几个重要结论：
关于坡度空间更详细的内容请看这个大佬的博客==> Slope Space in BRDF Theory - Nathan Reed
这部分的所有内容都是根据这俩文章来的，其中第二个是一篇的代码，抄起来很简单，PBRT里的相关部分:microfacet.cpp中的Sample_wh函数，也是论文中的实现版本
（你看现在咱有俩可以抄的代码了，不至于写不出吧？
paper - Importance Sampling Microfacet-Based BSDFs using
the Distribution of Visible Normals
supplemental1
我此处只是简单记录一下。
首先，常用的坡度空间分布（the distribution of slopes，法线分布的坡度空间形式，也就是$P^{22}$）和遮挡函数都具备有拉伸不变性，简单来说，微表面拉伸后原来会被遮挡的射线还是被遮挡，原来不被遮挡的射线还是不被遮挡。
这里以PBRT中函数TrowbridgeReitzSample(const Vector3f &wi, Float alpha_x,Float alpha_y, Float U1, Float U2)为例。
这就是文章中第一个思想：既然拉伸不变，那可以把各向异性的微表面拉伸成各项同性。
// 1. stretch wi
Vector3f wiStretched =
Normalize(Vector3f(alpha_x * wi.x, alpha_y * wi.y, wi.z));
//注：拉伸光线和拉伸微表面等价，康康上面的图
这引出了第二点：旋转：既然各项同性，那不妨在$(x',0,z')$空间上采样，得到坡度。这样做的好处是减少了一维的计算。
// 2. simulate P22_{wi}(x_slope, y_slope, 1, 1)
Float slope_x, slope_y;
TrowbridgeReitzSample11(CosTheta(wiStretched), U1, U2, &slope_x, &slope_y);
//注：对P22采样的具体实现放在了该函数的重载中，里面做了挺多优化来防止NaN出现，还是照抄吧
采样完成后再旋转相应的角度回到$(x,y,z)$上。
// 3. rotate
Float tmp = CosPhi(wiStretched) * slope_x - SinPhi(wiStretched) * slope_y;
slope_y = SinPhi(wiStretched) * slope_x + CosPhi(wiStretched) * slope_y;
slope_x = tmp;
最后，再把坡度拉伸回去，放到原来的相性上。
// 4. unstretch
slope_x = alpha_x * slope_x;
slope_y = alpha_y * slope_y;
计算，得到法线。
// 5. compute normal
return Normalize(Vector3f(-slope_x, -slope_y, 1.));
GGX的相应实现也差不多，只是P22采样的部分变了下而已。