光线追踪

光线追踪

模拟光路成像

光线的假设

1. 光线是沿直线传播
2. 光线之间不会发生碰撞
3. 认为光线可以模拟为从摄像机射出到物体 -- 原本是从物体到摄像机在真实物理情况上

光线投射 Ray Casting

摄像机和物体之间存在于一屏幕，从摄像机开始向物体每个点连线，然后对应在屏幕上可以进行光栅化。
然后从光源出发向对应的每个点连线，可以计算出阴影

Ray Casting -- 生成从摄像机发出的光线

1. 假设摄像机就是一个点
2. 物体接受到光线会进行完美的反射和折射
3. 从摄像机到屏幕上的一个像素块形成一条发出光线，只考虑该发出光线 eye ray和物体的最近交点 -- 深度
4. 当产生交点，需要考虑这个点是否能够被照亮，将该点与光源连接形成第二条光路shadow ray，考虑这条光路上是否有其他点。

Whitted-Style Rat Tracing

eye ray遇见的第一个点具有折射和反射

计算着色时，需要计算每一个交点与光源的连线是否能够被照亮。
对于这一个像素，将是所有能够被照亮交点的着色结果

Ray-Surface Intersection -- 求eye ray 与物体表面的交点

对于一个点o，传播方向d，那么ray上的任意一点$r(t)=o+td(t>=0)$

ray和球的交点

球的数学表达式：$p:(p-c{)}^2-R^2=0$
相当于求p和r(t)的公共点，有$(o+td-c{)}^2-R^2=0$ 然后求解t。
根据t大于0的值的数量 可以判断交点情况

ray 和所有表面$p:f(p)=0$的交点

$$f(o+td)=0$$

求解上式，t>=0的值，选取有效值中最小值即可

ray和三角形求交点

tip:
判断一个点是否在几何体内部还是外部
以该点为起点，做一条向几何体边界方向的射线，如果交点数为奇数说明在几何体内部，如果为偶数则是外部

简化：ray和三角形所在平面的交点+判断交点是否在三角形内
平面的定义：一个平面的法线n + 平面所过的一个点$p^{{}^{\prime }}$，那么该平面可以被写作

$$p:(p-p^{{}^{\prime }})\cdot n=0$$

而p又作为ray上的一点所以有$(o+td-p^{{}^{\prime }})\cdot n=0$，可以解得t

$$t=\frac{(p^{{}^{\prime }}-o)\cdot n}{d\cdot n}$$

也可以直接求解，可以利用重心坐标来表示三角形内部一点，所以有

$$o+td=(1-b_1-b_2)P_0+b_1{P}_1+b_2{P}_2$$

求解$t,b_1,b_2$即可，然后判断合理性

ray和多个三角形组成的表面求交点

1. 基本方法：ray和路径上每个三角形求交，然后选择最近一点 -- 消耗大
2. 优化：包围盒 -- 与包围盒没有交点就不需要考虑
   3D空间中，一般使用长方体来充当包围盒，可以将长方体理解为3组对面形成的交集(AABB)
   包围盒的每组对面一般为平行于坐标平面 -- 主要是可以简化法向量
   

ray 和AABB求交点

3D空间：

ray进入了所有3组对面 才能 说明ray进入了包围盒
ray只要离开了任意1组对面 就说明ray离开了包围盒
对每组对面计算$t_{min},t_{max}$，那么ray进入和离开包围盒的时间为$t_{enter}=max(t_{min}),t_{exit}=min(t_{max})$

如果$t_{enter}<t_{exit}$说明光线会进入这个包围盒
如果ray和AABB有交点当且仅当$t_{enter}<t_{exit}$ && $t_{exit}>=0$

AABB如何加速交点的求取

基本思路

1. 构建包围盒
2. 将包围盒划分为多个格子
3. 对格子进行分类：包含空间中物体边界的格子A和没有包含的格子B
4. 检测ray与格子的相交情况，如果遇到了A类格子，那么再进行判断ray是否与物体相交

如何划分格子

1. 以空间划分（KD tree）
2. 以物体划分 (BVH)
   
   KD-Tree 以2D为例 每次划分以x,y方向交替进行，将前一步所得的格子划分。
   前一步得到的格子，可以根据判断条件决定是否继续进行划分，这样一来就相当于形成了一个二叉树。
   对于该二叉树的叶子节点，存储的是待判断是否有交点的物体的格子
   对于该二叉树的非叶子节点，存储的是自己这个格子和划分的子格子。
   然后进行ray的相交判断。对于每个节点都需要根据包围盒求得的$t_{min},t_{max}$来判断是否与该节点所表示的本身格子相交。如果不相交那么对于其下面的子树都不用判断了。如果相交那么继续判断这个子树。
   对于叶子节点的判断，当ray与叶子节点所代表的格子相交时，要继续进行判断ray是否与物体有交点

BVH划分方式

以物体为基础来划分，对物体进行分类，一类的物体用一个包围盒约束在一起（每次划分都需要重新计算包围盒），这些包围盒可能被重叠。
但对于每个叶子节点，同一类物体只可能出现在一个叶子节点中。 -- 降低重叠区域

SAH加速BVH划分

划分技巧

1. 每次划分选择当前最长边进行分割
2. 分割节点选择当前最长边上节点的中间节点
3. 当某个节点中的三角形较少时，可以停止使得该节点成为叶子节点

AABB求交伪代码

//递归
Intersect(Ray ray, BVH node)
{
  //如果光线没有和当前节点所表示的包围盒相交直接跳过
  if(ray misses node.bbox) return;
  //如果当前节点是叶子节点，那么需要判断ray和该bbox中所有物体的交点，选择最小值
  if(node is a leaf node)
  {
    findhit (ray,node.bbox.objs);
    return closest hit;
  }
  //如果不是，递归求解两个子树
  hit1 = Intersect(ray,node.leftchild);
  hit2 = Intersect(ray,node.rightchild);

  return closest(hit1,hit2);
}

辐射度量学

Radiant Energy and Flux(Power)

Radiant Energy: 能量 $Q[J=Joule]$
Radiant Flux(Power): 单位时间上的Radiant Energy$\mathrm{\Phi }=\frac{dQ}{dt}[W=Watt][lm=lumen]$
Radiant Intensity: 单位立体角所具有的Flux$I(\omega )=\frac{d\mathrm{\Phi }}{d\omega }[\frac{lm}{sr}=cd]$，具体体现就是光源在某个方向上的Flux，并不会衰减
立体角:

1. 2D平面上一个角$\theta =\frac{l(\mathrm{弧}\mathrm{长})}{r(\mathrm{半}\mathrm{径})}$
2. 3D平面上的立体角$\mathrm{\Omega }=\frac{A(\mathrm{球}\mathrm{面}\mathrm{上}\mathrm{的}\mathrm{面}\mathrm{积})}{r^2}$
3. 微分立体角$d\omega$：微分面积(dA)除以半径平方

图片
如果点光源是均匀的向四周发出Flux 那么$I(\omega )=\frac{\mathrm{\Phi }}{4\pi }$

Irradiance: 单位物体表面接受到的Flux(Power)$E(x)=\frac{d\mathrm{\Phi }(x)}{dA}$，Flux只计算和当前表面垂直的分量，随着距离衰减

Radiance: 光在传播过程中的能量，Power(Flux)在单位立体角上的单位面积的大小(per unit solid angle and per projected unit area)$L(p,\omega )=\frac{d^2\mathrm{\Phi }(p,\omega )}{d\omega dA\cos \theta }$，可以这样理解，一个单位面积向某一个方向辐射的Power

1. Irradiance per unit solid angle
2. Intensity per unit area

BRDF 双向分布反射函数

反射对于一点

某一点接受光线，再从这个点发射出去 -- 反射的另一种理解
对于某一点接收到一个方向的能量 相当于就是当前方向的 radiance 再乘以当前方向的立体角，得到的Irradiance。
反射到摄像机相当于就是当前点的Irradiance，对于这个方向辐射的能量。

即变换为这样一个函数，对于一个点的面积$dA$接收到来自一个立体角$d\omega$的Irradiance会如何分配到各个不同的立体角上

BRDF

BRDF就是定义的这样一个函数，它直接给出的是每个反射出去方向的radiance占某一个立体角射入的Irradiance的比例

$$f_r({\omega }_i->{\omega }_r)=\frac{d{L}_r({\omega }_r)}{d{E}_i({\omega }_i)}=\frac{d{L}_r({\omega }_r)}{L_i({\omega }_i)\cos {\theta }_{i}d{\omega }_i}$$

反射方程

BRDF给出了是从一个立体角射入，反射到所有立体角的关系，
实际应用求反射相当于是空间中每个光源所在的立体角入射然后反射到同一个立体角的能量

$$L_r(p,{\omega }_r)={\int }_{H^2}{f}_r(p,{\omega }_i->{\omega }_r)L_i(p,{\omega }_i)\cos {\theta }_{i}d{\omega }_{i}d{\omega }_i$$

$L_i(p,{\omega }_i)$的困境

由于光不会只反射1次所以出射的radiance可能作为其他点入射的radiance -- 导致递归

基础渲染方程

对于一个点的渲染，可以由自己发出和其他物体包括光源发出的所有radiance之和

$$L_r(p,{\omega }_o)=L_e(p,{\omega }_o){\int }_{H^2}{f}_r(p,{\omega }_i->{\omega }_r)L_i(p,{\omega }_i)(n\cdot {\omega }_i)d{\omega }_{i}d{\omega }_i$$

通过简化 -- 可以得到全局光照的表达式

$$L=E(\mathrm{光}\mathrm{照})+KE(\mathrm{直}\mathrm{接}\mathrm{光}\mathrm{照})+K^{2}E(\mathrm{间}\mathrm{接}\mathrm{光}\mathrm{照})+K^{3}E+...$$

可以发现光栅化的着色就是光照+直接光照的结果

蒙特卡罗积分-路径追踪

蒙特卡罗积分

求解定积分的一种方法，利用其来求解BRDF的那个积分式。对于定积分${\int }_a^{b}f(x)dx$，有概率密度函数$X_{i}p(x)$，那么根据蒙特卡洛积分有

$${\int }_a^{b}f(x)dx=\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{p(X_i)}$$

路径追踪

Whi-RayTracing 对于所有物体都是沿着一个方向进行折射或者反射，即对于漫反射表面存在错误的结果。
所以采用渲染方程来求解着色对于渲染方程中的积分可以使用蒙特卡洛方法求解，对于待着色点o，和观测方向${\omega }_o$

$$L_o(o,{\omega }_o)=\frac{1}{N}\sum _{i=1}^N\frac{L_i(o,w_i)\cdot f_r(o,w_i,w_o)\cdot (n\cdot {\omega }_i)}{p({\omega }_i)}$$

那么对于直接光照有如下伪代码

//对点o着色
shade(o,wo):
  //初始化着色结果
  Lo = 0.0f;
  //随机选择N个点WI 和定义一个pdf函数
  for(auto& wi)
  {
    //从o点到WI构建wi的光线
    Trace a ray r(o,wi)
    //如果光线能够到达光源说明是这条光线属于直接光照
    if(ray r hit the light)
    {
      Lo += (1/N)*Li*fr*cos(n,wi)/pdf(wi);
    }
  }
  return Lo;

对于间接光照，即相当于光先照射到q点，然后q点反射出来照射到p点的光线对于p点的着色，而这部分光的能量相当于从o点观测到q点的着色

  //q为构建的光线wi打到的其他非光源物体的点
  //shade(q,-wi) 即为从o点观测能够接受到q点的反射能量值，相当于从q点对o点的光照
  Lo += (1/N)*shade(q,-wi)*fr*cos(n,wi)/pdf(wi);

所以全局光照的路径追踪算法

shade(o,wo):
  //初始化着色结果
  Lo = 0.0f;
  //随机选择N个点WI 和定义一个pdf函数
  for(auto& wi)
  {
    //从o点到WI构建wi的光线
    Trace a ray r(o,wi)
    //如果光线能够到达光源说明是这条光线属于直接光照
    if(ray r hit the light)
    {
      Lo += (1/N)*Li*fr*cos(n,wi)/pdf(wi);
    }
    //间接光照
    else if(ray r hit the obj-point q)
    {
      Lo += (1/N)*shade(q,-wi)*fr*cos(n,wi)/pdf(wi);
    }
  }
  return Lo;

存在的问题

1. N过大-->弹射光线数量指数级增加
   只有N=1才避免这种情况 -- 路径追踪 -- 噪声过大
   对于一个像素使用多条Path，即对于一个像素内部再选取点

//对于一个像素的路径追踪 N=1
ray_generation(camPos,pixel)
  //从pixel中均匀选择N个点
  //像素的着色
  pixel_color = 0.0f;
  //对于像素中的N个采样点
  for(each sample in the pixel)
  {
    //从相机到采样点构建一条ray
    ray r(camPos,sample);
    //如果打到了物体
    if r hit obj-point p
    {
      //计算方向相反
      pixel_color += (1/N)*shade(p,-r);
    }
  }
  return pixel_color;

2. 无法终止
   概率小，即通过弹射到达光线停止
   通过概率控制停止 -- 俄罗斯轮盘赌
   1. 定义概率P，通过P来决定是否发出这条光线
   2. 发出光线的返回值为$\frac{L_o}{P}$
   3. 未发出光线的返回值是0
   4. 最终返回的期望仍然是Lo，只不过具有噪声

最终路径追踪代码

shade(o,wo):
  //初始化着色结果
  Lo = 0.0f;
  //概率P_RR控制光线是否发出
  float P_RR;
  //每次随机选择一个0~1的值
  random(ksi);
  if(ksi > P_RR) return 0.0f;

  //随机选择1个点WI 和定义一个pdf函数
  for(auto& wi)
  {
    //从o点到WI构建wi的光线
    Trace a ray r(o,wi)
    //如果光线能够到达光源说明是这条光线属于直接光照
    if(ray r hit the light)
    {
      Lo += Li*fr*cos(n,wi)/pdf(wi);
    }
    //间接光照
    else if(ray r hit the obj-point q)
    {
      Lo += shade(q,-wi)*fr*cos(n,wi)/pdf(wi);
    }
  }
  return Lo/P_RR;

修改采样方法

之前是通过在着色点上随机选择N条光线，会导致可能由于光源过小出现无法打入光源的情况。
现在从光源上采样，那么对应要将渲染方程写为在光源上的微表面积分，即考虑单位立体角$d\omega$和光源上微表面$dA$的关系
$dA$对应投影于以着色点为球心的单位球的面积就是其所对应的立体角

$$d\omega =\frac{dAcos{\theta }^{{}^{\prime }}}{||x^{{}^{\prime }}-x|{|}^2}$$

${\theta }^{{}^{\prime }}$为光源的法线和着色点与微表面连线的夹角

修改算法

1. 直接光照
   对光源采样，不再需要P_RR控制
   
2. 间接光照
   不变，需要P_RR控制
   
3. 直接光照判断光线是否会被遮挡