具体理论学习推导参考Wikipedia,Inverse Transform Method，PBRT-P643,

下文中的CDF与PDF等概率论术语定义可在概率论与数理统计等书中找到

生成一个随机数满足一个已知解析式的分布，是一个非常常见的需求(C++11已经内部提供不少给定分布)，这里将推导如何将最常见的均匀分布随机数转换为其他已知解析式的分布。

先把结论放在这里，若有一随机变量X满足分布PDF（对应一CDF），使用一个均匀分布的随机变量Y，则Y对该CDF的反函数进行采样的结果满足PDF。

以下内容摘自Wikipedia

已知有分布$F(x)=1-e^{-\lambda x}$，其中$x\geqslant 0$，可以求得CDF的反函数，如下图所示。

$x=F^{-1}(y)=-\frac{1}{\lambda}(1-y)$

其中$y_i$是满足区间(0, 1)均匀分布的随机数(即，Y~U(0, 1))，这里可以明显看到由反函数求得的X的分布情况，在越接近0的位置分布的越多，越背离则越稀疏，满足对于结果的设想。

证明

若有CDF满足$F(x)=Pr\{X\leqslant x\}$，则起反函数$F^{-1}$满足

$F^{-1}(y)=min\{x: F(x)\geqslant y\}, y\in [0, 1]$

之所以使用下确界的原因是，$F(x)$满足右连续且单调不减

令$X=F^{-1}(U)$其中$U$满足在区间$[0,1]$上连续均匀分布，现要证

$Pr\{X\leqslant x\}=F(x), x\in \Re$

则要证$Pr\{F^{-1}(U)\leqslant x\}=F(x)$成立。此处$F$连续且单调不减，则同时对不等式套用$F$不等式成立

$\{F^{-1}(U)\leqslant x\}=\{U\leqslant F(x)\}$

则有$Pr\{U\leqslant F(x)\}$成立，下图中描述了$U\leqslant F(x)$的概率就等于$F(x)$本身

则得证
$Pr\{F^{-1}(U)\leqslant x\}=F(x)$

有了Inverse CDF这项工具，就可以很方便的将已有的满足均匀分布的随机数转化为满足一已知PDF的随机数分布。

二维随机变量中条件概率直观理解(其中关于联合以及边缘分布律/概率密度分布的概念可以在上海交大的PPT以及Wikipedia中找到)。在渲染中主要在重要性采样中会用到二维随机变量的这些概念

离散条件分布和连续条件分布的定义都可以在引用中看到，这里只关注离散型二维随机变量的分布。

首先有联合=条件*边缘的定义，当然同时要满足$\sum_{x}\sum_{y}P(X=x, Y=y)=1$的前提条件。

$P(X=x, Y=y)=P(Y=y|X=x)P(X=x)=P(X=x|Y=y)P(Y=y)$

若X,Y满足相互独立，那么也满足以下等式

$P(X=x, Y=y)=P(Y=y)P(X=x)$

在渲染中为了对图像上亮度更高的点有更高的概率采样到，那么势必要生成满足图像PDF的二维随机变量。此处易证在图像的采样上，X与Y二者之间非相互独立。

此图中$P(X=x, Y=y)=\frac{b}{b+8a}$，而$P(X=x)=P(Y=y)=\frac{b}{b+4a}$，容易看出$P(X=x, Y=y)\not = P(Y=y)P(X=x)$

那么求解上述就只剩下一种途径，就是对于给定图像求解边缘以及条件分布函数。这也是PBRT中使用边缘+条件来生成联合的原因。

以下使用OF测试得到的可视化结果，因为一维采样重建可以在二维中体现因此此处不做一维对比。具体实现办法可以在PBRT的Distribution中看到详细过程

左右为不同采样数量但可以明显发现在图像中亮度较高位置采样数量更多，采样结果明显偏向于向图像亮处。

这里图像较暗是因为没有采用Tone Mapping的关系，但因为不影响全局显示，因此没做特殊处理。此处若使用普通PNG等被截断过的图像则无法有效的对图像上高光部分进行有效采样。

考研微分方程计算中，并未包含伯努利方程的要求，但一样可以通过变量代换实现求解。下题为一个经典的变量代换求解，但自己发现却可以使用更巧妙的办法解答，张宇解析中并没有提及，因此记录在此。

求 $xdy=(1+ylnx)ydx (x>0)$

变量代换:

$\Rightarrow$ $xy'=y+y^2lnx$

$\Rightarrow$ $\frac{x}{y^2}y'=\frac{1}{y}+lnx$

令t = $\frac{1}{y}$，那么$\frac{dt}{dx}=-\frac{1}{y^2}\frac{dy}{dx}$

$\Rightarrow$ $-x\frac{dt}{dx}=t+lnx$

凑出一阶线性微分方程形式

$\Rightarrow$ $\frac{dt}{dx}+\frac{1}{x}t=-\frac{lnx}{x}$

$\Rightarrow$ $t=e^{-\int{\frac{1}{x}dx}}(\int{e^{\int{\frac{1}{x}dx}}(-\frac{lnx}{x})dx}+C)$

$\Rightarrow$ $t=\frac{1}{x}(\int{x(-\frac{lnx}{x})dx}+C)=-\frac{1}{x}(xlnx-x+C)=\frac{1}{y}$

$\Rightarrow$ $x=xy-xylnx+Cy$

凑导数形式:

$\Rightarrow$ $xy'=y+lnxy^2$

$\Rightarrow$ $\frac{y-xy'}{y^2}=-lnx$

$\Rightarrow$ $(\frac{x}{y})'=-lnx$

$\Rightarrow$ $-\frac{x}{y}=xlnx-x+C$

$\Rightarrow$ $x=xy-xylnx+Cy$

文中将会涉及到个人在内存泄露上的一些经验，主要来自于Atmos渲染器在实现可重入的过程。

一般内存溢出的问题，有较好找的，也有非常不容易发掘的。这里按照个人经验排序，从简到难记录。将会不定期更新。

0.误分配内存，这种情况最好查找，也最经常犯。Atmos在实现日志系统中，误在字符转换中，遗漏了一个char数组忘记释放，以至于全盘内存上扬.

1.只有在构造函数中分配了指定内存，经常偷懒不写析构函数的问题就在于，经常忘记释放构造函数中直接分配的内存

class A
{
public:
      A(){
            b = new B();
      }

      B* b;
};

2.Pimple所带来的弊端。Pimple能将程序实现从头文件中剥离，这算是一个好处，但也经常带来麻烦

// h
class A
{
public:
    A();
    ~A();
private:
    class B;
    B* b;
};

// cpp
A::A(){
    b = new B();
}
A:~A(){
    delete b;
}

乍一看似乎是没有问题，但这严重依赖于要在内部实现B的时候也要实现其析构函数.

// cpp
class A::B{
public:
    B(){}
    //~B(){}
};

在Atmos中，图像导入和导出都是被封装在cpp中不外漏的。因此在实现中我经常遗漏了内部实现中B析构函数的实现，或者是实现了压根没有记得真正将图像内存直接释放。

3.我称之为浅释放的问题……经常有的时候并想不起来这东西需要手动释放，例如Atmos中遇到的。实现BVH生成的二叉树，只保存了root结点

class BVH
{
public:
    ~BVH(){
        delete root;
    }
    Node* root;
};

问题在于，root只是二叉树的开始，需要手动遍历释放，或者实现Node的析构来完成释放。

void deleteNode(Node* root)
{
    if(root->left)  deleteNode(root->left);
    if(root->right)  deleteNode(root->right);
    delete root;
    root = NULL;
}

// 或者是
Node::~Node(){
    if(left) delete left;
    if(right) delete right;
}

BVH::~BVH(){
    delete root;
}

本质上说就是忘记了数据结构的定义而直接对指针本身进行了内存的释放操作

4.虚析构函数的重要性，在一般面对资源管理的类型上尤其要注意到，父类指针若想要释放实例化子类的对象，就必须要实现虚析构函数，否则子类对象的析构函数就无法被调用。
例如上例中描述的BVH，其继承primitiveSet，但因为父类未实现虚析构，导致即便其实现了析构函数也不起作用.

class father
{
public:
    //~father(){/*A*/}
    virtual ~father(){/*B*/}
}

class child:public father
{
public:
    ~child(){
        // A:not work 
        // B:it worked
    }
};

5.第三方库导致的内存未释放。这类问题非常难以察觉，以至于一开始根本不会意识到是这里的问题，除非在表象上总结出规律来。Atmos在实现图像和模型导入导出时，从未调用过库函数的释放API，这就导致了所有的内存全都未被释放。因此，在想要实现可重入之前，务必祥读文档才是真理。

6.有一些连锁的，具有内存释放前后关系的是最难以解决的。比如，对象A中包含B的指针，该B的指针由A来创建并返回给使用者。那么问题来了，在A的析构函数中可否释放B？

答案是无解，需要具体情况具体分析，有的情况做好约定或者引用计数等办法，可以在A中直接释放；而有的情况则是A仅仅只做一个指针保留，那么此时就不应该去释放。还是一句话，具体情况具体分析。

本文的三角形相交测试算法为Möller-Trumbore algorithm

这里的back-facing culling表示的是若此三角形背对着观察者，那么就直接在相交测试中跳过，返回false。这在正常逻辑思维下是正常的，直到遇到了Refraction的材质。

首先是正常参数方程定义的两个球的对比，基本上除了目前未调正确的Fresnel以外，折射部分是没有太大问题的。

Atmos的渲染结果

Mitsuba的渲染结果

一旦将其改为网格形成的球体，结果令人吃惊的不一致。

Atmos的渲染结果

Mitsuba的渲染结果

在ScatchPixel原文中有这样一段

“… Remember as well that the user might want to cull (discard) back-facing triangles. …”用户自定义背面剔除

事实上，默认的三角形理应是双面的，因为但当类似于玻璃材质模型光线在几何体内部时必然会遇到光线从背后穿过三角形的问题。一开始我并没有意识到这是三角形的背面剔除的问题，球上折射的内容让我误以为是因为法线的错误问题。

正确的不跳过背面三角形的Atmos渲染结果如下

最终，提供材质球的修复前后对比。

BUG非常隐晦，一开始的确完全没有想到是剔除的问题

错误结果

正确结果

坑爹，这BUG找了老子一整晚

Deprecated(2017.03.30)

这一篇主要记录PBRT中Path Tracing的理解笔记。一旦涉及到了代码层面的实现，会显得无从下手。

首先是路径追踪的蒙特卡洛估计量写法

$\frac{L_e(p_i \to p_{i-1})f(p_i \to p_{i-1} \to p_{i-2}|cos\theta_{i-1}|)}{p_A(p_i)}(\prod^{i-2}_{j=1}\frac{f(p_{i+1} \to p_j \to p_{j-1})|cos\theta_j|)}{p_{\omega}(p_{j+1}-p_j)})$

路径追踪加上了直接光采样，本质上就是应用重要性采样以及对于Dirac delta Distribution等的处理。

如图所示

使用一简化后的路径追踪表达式(省略了pdf,russian rullte，以及立体角到面积形式的表示，改变了累乘的上限)

$L_e(P_i \to P_{i-1})*(\prod_{j=1}^{i-1}f(P_{j+1} \to P_j \to P_{j-1})cos\theta_j)$

意味着，当前追踪到的点上求得的直接光照，是用来计算上一节点的Reflected Radiance的。令$E$表示光源，那么直接光采样得到的结果，表示当前点$P_i$位置处，向$P_i \to P_{i-1}$反射的光照，因此需要在计算直接光照时，乘以权重(仅表意，严格来讲这里是做了一次多重重要性采样)

$L_e(P_i \to P_{i-1})=L_e(E \to P_i)*\frac{f(P_i,w_{E \to i},w_{i \to i-1})cos\theta_{i-1}}{pdf}$

因此这里$P_i$点处计算得到的直接光照，被用于点$P_i$用作间接光照/反射光来传递给下一节点。由于运用了光路可逆，因此在代码实现过程中循环过程会是相反的，也就是正向追踪光路。

这一事实较难以直接从公式中get到，下方通过PBRT源代码验证。

PBRT通过保存一个pathThroughput变量，用于记录原式子右边的累乘项

$\prod^{i-2}_{j=1}\frac{f(p_{i+1} \to p_j \to p_{j-1})|cos\theta_j|)}{p_{\omega}(p_{j+1}-p_{j})}$

    for (int bounces = 0; ; ++bounces) 
    {
        // 逆向计算直接光照在上一节点的贡献/反射光照
        Vector wo = -ray.d;
        L += pathThroughput * UniformSampleOneLight(wo,...);

        // 计算当前结点BRDF以及出射方向
        Spectrum f = bsdf->Sample_f(wo, wi, pdf...);

        // 式子中的累乘项，晚于直接光照的计算
        pathThroughput *= f * AbsDot(wi, n) / pdf;

        // 新的跟踪光线
        ray = Ray(p, wi, ...);

        if (bounces == maxDepth) break;

        if (!scene->Intersect(ray,...)) 
        {
            // 无限远区域光光照
            L += pathThroughput * infiniteLight;
            break;
        }
    }
    return L;

因此在理解增量形式路径追踪的实际实现时，最关键的点为以下两点

1.当前追踪到的点上计算得到的直接光照，是作用于上一节点的反射光照上的。

2.循环中累乘项应晚于直接光照计算，因此其累乘项截止到上一节点为止。

Deprecated(2017.03.30)

对于Unidirectional Path Tracing的理解是建立在Path Tracing之上

本文用于记录当前阶段对于Path Tracing的理解。此笔记的问题来源于，对于Path Tracing使用Direct Lighting的疑惑。

很多教材上提到了渲染方程

$L_o(p,\vec{w_o}=L_e(p,\vec{w_o})+\int_{2\pi}f_r(p,\vec{w_i},\vec{w_o})L_i(p,\vec{w_i})cos\theta_i d\omega_i$

然后紧接着说到，我们将此方程改写成直接光照部分和间接光照部分，比如Advanced Global Illumination，问题出在，当真正遇到代码实现的时候，压根找不到间接光的部分，不论是BDPT还是PT都无法直观的理解到间接光的部分。或者说不清楚间接光在代码中的体现。

下图是帮助理解的核心

关键部分就是间接光照 = 反射光照，这从根本上解释了间接光照的含义就是从其他位置处反射的光照，而非从光源上。

假定所有表面都不自发光，也就是$L_e(p,\vec{w_i})=0$，那么定义X位置处上的入射辐射度为$L_x$，直接光照部分记为$D_x$，间接光照/反射光照部分记为$R_x$，BRDF项记为$f_{rx}$

$\begin{array}{ccc}{L}_x& =& f_{rX}\ast (D_X+R_X)\\ & =& f_{rX}\ast (D_X+f_{rA}\ast (D_A+R_A\left)\right)\\ & =& f_{rX}\ast (D_X+f_{rA}\ast (D_A+f_{rB}\ast (D_B+R_B)\left)\right)\\ & =& f_{rX}\ast (D_X+f_{rA}\ast (D_A+f_{rB}\ast (D_B+f_{ri}\ast (\mathrm{.}\mathrm{.}\mathrm{.}\left)\right)\left)\right)\end{array}$
\begin{array}{l}
L_x &=& f_{rX}(D_X+R_X)
\ &=& f_{rX}(D_X+f_{rA}(D_A+R_A))
\ &=& f_{rX}(D_X+f_{rA}(D_A+f_{rB}(D_B+R_B)))
\ &=& f_{rX}(D_X+f_{rA}(D_A+f_{rB}(D_B+f_{ri}(…))))
\end{array}

也就是

$L_o(p,\vec{w_o})=\int_{2\pi}f_r(p,\vec{w_i},\vec{w_o})(Le(rc(p,\vec{w_i},-\vec{w_i})+Lr(rc(p,\vec{w_i}),-\vec{w_i})))cos\theta_i d\omega_i$

这里$L_e$项表示emitted radiance，因为其他表面没有自发光项，因此就是光源部分的辐射度贡献，即直接光照direct illumination部分；$L_r$项表示reflected radiance，也就是间接光照indirect ilumination。

根据无穷积分的定义，也可以改写为

$\begin{array}{ccc}{L}_o(p,\overrightarrow{w_o})& =& {\int }_{2\pi }{f}_r(p,\overrightarrow{w_i},\overrightarrow{w_o})Le\left(rc\right(p,\overrightarrow{w_i},-\overrightarrow{w_i})cos{\theta }_{i}d{\omega }_i\\ & +& {\int }_{2\pi }{f}_r(p,\overrightarrow{w_i},\overrightarrow{w_o})Lr\left(rc\right(p,\overrightarrow{w_i}),-\overrightarrow{w_i})\frac{cos{\theta }_{i}cos{\theta }^{\mathrm{\prime }}}{\mathrm{\mid }\mathrm{\mid }{p}^{\mathrm{\prime }}-p\mathrm{\mid }{\mathrm{\mid }}^2}V(p,p^{\mathrm{\prime }})dA\end{array}$
\begin{array}{l}
L_o(p,\vec{w_o}) &=& \int_{2\pi}f_r(p,\vec{w_i},\vec{w_o})Le(rc(p,\vec{w_i},-\vec{w_i})cos\theta_i d\omega_i
\ &+& \int_{2\pi}f_r(p,\vec{w_i},\vec{w_o})Lr(rc(p,\vec{w_i}),-\vec{w_i})\frac{cos\theta_i cos\theta'}{||p'-p||^2}V(p,p')dA
\end{array}

右边的$V(p',p)$项就是shadowRay计算是否可视的结果

所谓的最简单的递归形式的Path Tracing本质上是对渲染方程的另一种解读，可以想象有以下光路

假定光路上任意点$p_i$的$BRDF(p)cos\theta_i$乘积定义为$Q_i$，若光源上出射辐射度为$L_i$

那么$L_o(p,\vec{w_o})=L_i$ Q_D Q_C Q_B Q_A

回想伪代码实现

for i : spp

    new ray()

    Li += li(ray, scene) / spp

end

等价于

$\begin{array}{ccc}{L}_o(p,\overrightarrow{w_i},\overrightarrow{w_o})& =& (Li\ast Q_D\ast Q_C\ast Q_B\ast Q_A\\ & +& Li\ast Q_E\ast Q_C\ast F_B\ast Q_G\\ & +& Li\ast Q_H\ast Q_I\mathrm{.}\mathrm{.}\mathrm{.})\\ & \ast & \frac{1}{spp}\end{array}$
\begin{array}{l}
L_o(p,\vec{w_i},\vec{w_o}) &=& (LiQ_DQ_CQ_BQ_A
\ &+& LiQ_EQ_CF_BQ_G
\ &+& LiQ_HQ_I…)
\ &*& \frac{1}{spp}
\end{array}

这一计算式本身就是$\int_{2\pi}f_r(p,\vec{w_i},\vec{w_o})L_i(p,\vec{w_i})cos\theta_idw_i$的蒙特卡洛近似实现了

结论

本身递归形式和增量循环形式的Path tracing本身都是一种渲染方程的实现方式，并无差错之分。但前者的问题在与光路击中光源的概率太低的时候，比如场景中只有点光源时那么结果会是一片黑。因此在实现效果上考虑了直接光采样的Path Tracing会更胜一筹。