从0开始学习PBRT

概况

  • 基于`Physically Based Rendering <https://www.pbr-book.org/3ed-2018/contents>`__第三版

  • 过一遍必要的知识点
    尽量记录自己的理解和推导
    记录有用的资料/资源
  • 把图形学捡起来
    做一个小小渲染器
  • 其实有点觉得从0直接看pbrt能更好学习图形学原理
    因为能从最底层的结构一步步搭建程序且不用折腾市面上的图形api
    当然一般都是从图形api开始因为更容易出结果
  • 这个书是参考书性质
    无基础的话需要颠来倒去看
    很多知识点都是一环套一环
    每个小问题都可能引出一堆问题
  • 需要学习研究大量各方面知识包括数学物理编程等等
    需要看各种相关论文和书籍
  • 要吃透所有知识点很困难
    需要反复学习思考
    是个长期过程
    有的细节要敢于暂时忽略

需要的基本知识

  • 线性代数

  • 微积分

  • 概率

  • 如果有图形学基础会比较舒服


可以先大致看看书,然后搞一下初始代码。

理论上是不需要图形api(opengl等),也不需要gui的。
但计划做成gui程序,用起来舒服点。以后可以不断加功能,有整体感。总的目标是一个可用的引擎。
gui用https://github.com/ocornut/imgui
后端用glfw(vulkan)。
glfw的作用是辅助起windos图形程序,处理基本的事件/上下文交互/各种初始化等。
vulkan只是用于glfw,并不参与pbrt。
imgui是gui库,可选glfw作为后端。

起环境

安装vs2022
安装cmake
glfw
感兴趣可以下载编译一下观察它的例子程序。
直接拿imgui代码
用imgui里的example_glfw_vulkan里的CMakeLists.txt起vs2022项目
需要事先安装vulkan sdk
实际主要编译的是glfw源码,imgui源码,imgui_impl_vulkan.cpp和main.cpp。

编译成功运行可看到imgui的demo窗口。

这个例子就是教你如何用imgui起gui程序框架。
可以在这个项目基础上搞自己的代码。
可以把CMakeLists.txt拿出来改改,把main.cpp拿出来,glfw和imgui做成git子模块。整个做成自己干净的项目。
程序入口在main.cpp的main里。
我们只需关心imgui的用法。其他vulkan相关的不用管。
当前的目标是不用图形api,只用数理逻辑学习pbrt原理。
gui用来设置参数和启动渲染。

把我们自己的界面和逻辑写在别的文件里管理好。最后在main里驱动即可。

到此gui基础框架起好了。


后话
用默认的imgui+vulkan不太容易显示图片,只能输出到文件再另外打开。
没多久就烦了,想要直接在gui里显示图片。然后我学了vulkan,对vulkan流程进行封装,再把imgui集成进去。
经过一周多的剧烈折腾知道了vulkan的基本流程,同时知道了imgui的大概原理。这样舒服多了。
也可以在任意窗口显示图形了。
以后要进一步发展的话绕不开这些东西。
所谓现代的图形api真的比较复杂,内容多功能多灵活性大 。

下面看pbrt代码。同时学习pbr理论。同时往自己项目加代码。
pbrt的代码做的比较全,各种概念/类都比较精细,是一个较大的框架。
后果就是繁琐,抄/改起代码来比较麻烦。好处就是实现下来能经手更多的代码,更好地理解原理。
计划从pbr的最小功能/必要的操作开始,用尽量少的代码跑通流程。再按需添加功能。

我的代码进程:

  • 基本几何

  • 矩阵变换

  • 透视相机

  • 球体

  • 球体与光交互

  • PointLight

  • RandomSampler

  • WhittedIntegrator

  • Lambertian bxdf

  • matte/玻璃/镜面材质

到这里可以出图了

然后

  • 错误处理

  • 多线程渲染

  • PathIntegrator 路径追踪

  • mipmap / hdr / 环境光 / infinite light

  • 三角形 / mesh / 简单glTF加载

  • BBox / BVH

  • 微小面bxdf

  • 完善gui(物体的crud / material的crud / 场景存取)

  • 简单texture

  • kd-tree

  • photon mapping 光子映射

  • 优化 / 内存池

  • 体渲染 / 次表面

初期看pbrt可能有时摸不着头脑。对比看下这个非常好,循序渐进。每个阶段都有实际输出。
抄pbrt的代码的一个问题是会陷入复杂代码里,初级阶段不容易输出。
等到能输出的时候,整个框架已经大致抄出来了。
可以同时参考这个网站出图。

我是抄到大概了解了WhittedIntegrator的流程以后开始出图。

pbrt的代码我比较懒,都没编过没跑过。就硬看书和代码,再用自己的代码一点点调起来。


1 介绍

1.2 PHOTOREALISTIC RENDERING AND THE RAY-TRACING ALGORITHM

1.2.2

计算ray和面的交点

1.2.3 光的分布

严格的点光源实际是不存在的。基于物理的光照一般是基于面光源。

我们关心的是交点的情况。实际处理上,是基于交点所在的微分面来计算。
关心交点上积存的光能。
以一点发出的光能为例
设总光能是P,在距离r上,可算出每个微小面dA上的能量dE=P/半径为r的求面积。
还可引入光线与物体表面的夹角,我们只需要垂直于表面的能量。
最后取dE在表面上的垂直分量即可。

1.2.4 阴影

shadow ray
从点向各种光源做直线,看看中间是否有遮挡。
无遮挡就计算光,到最后都没处理到的就是阴影。

1.2.5 表面散射

入射光如何散射出去

material
brdf
双向反射系数分布函数
\(brdf=f(p, 入射方向,射出方向)\)
可以看成物体的材质特性

这样有了光,方向,表面和材质。就可以算点上的光了。

radiance 辐射率

其他各种df

1.2.6 间接光

原始的光源打到物体上反射出去,又可以看作是一个光源。

光滑的物体比如镜子、金属等,为什么可以看到其他物体。
todo????因为光打到其他物体上反射出来,形成一个新的光源。而镜子材质的特性,对光的吸收较小,就保留了新光源的信息,就看到了物体。
毛糙的材质吸收大量的光

间接怎么理解。即除了本身的发光外,所有光源打过来的光,包括直接光源和其他物体反射过来的光。

一般来说,点射到眼睛的光为点本身发出的光(如果是发光体)加上所有间接光的和。
即渲染方程
间接光的和是积分的形式。对于一个微小的球面,结合brdf,用积分算出每个光源打过来的光能总和。

对于不是极为简单的场景,实际光源可以说是无数个。硬算肯定是不行的。要做简化。

最初的Whitted算法。对一个点只考虑一些重要的方向进来的光,比如镜面方向和折射方向。大大的简化。

在这基础上可以改进。比如多取一些与镜面方向相近的光,可得到glossy效果。

1.2.7 光的传播

之前考虑是真空理想状态。光不会衰减,不会被介质影响。
真实世界有各种因素影响光的传播。

两种影响。吸收和散射(发散到不同的方向)

1.3 pbrt系统概览

介绍代码结构

1.4 并行化

各个像素的计算是独立的,那么可以简单地分给多个cpu去算。速度成倍增长。
后续可研究并行计算。

多线程细节非常多,可以追究到硬件/cpu等层面。

  • 看c++的lambda相关

  • 看c++的thread相关

  • 看c++的memory模型。new/placement new等。

  • 线程里可以读公共数据,尽量不要去写。

  • 线程里不要频繁动态分配内存,不要频繁起std::vector之类。

  • 分配内存用内存池。

  • 有些数据结构需要手写,可比std的更实用高效。

  • 把场景平均分块给多线程计算。可能有的部分更复杂,造成拖慢整体进度。需要实现让完成任务的线程自动去分担没完成的任务,不让cpu空闲。


pbrt的内存池MemoryArena比较简单直接

  • 定一个block大小比如256k。每次实际分配的大小为256k的整数倍,为一个大块。

  • 分配内存要使用带对齐的api。

  • 有一个free链表存可用的大块

  • 有一个busy链表存正在使用的大块

  • 有一个独立的当前大块

  • 申请时先看当前大块剩余空间是否足够。

    • 足够直接用。

    • 不足的话当前大块进busy表,在free表中查找是否有空间足够的大块。

      • 找到就用

      • 找不到就必须新分配内存

  • 只支持单线程。在每个线程起唯一一个实例。

  • 每个线程会进行大量循环操作,每次的操作是独立线性的。

  • 每次循环之前清空block(只是逻辑清空,都进入free表,并不实际free内存)。

  • 线程结束时调用实际的free。

  • 总体只需申请1次循环中操作锁需要的内存,后续的循环大都是复用内存。


2 geometry 和 transformation

确定坐标系

定义基本的向量/点/ray之类。

空间变换。各种矩阵操作。四元数暂时不用实现。

坐标系

分左手/右手。zUp/yUp。

坐标系会影响一些变换。比如lookat。

我常用的blender是右手zup。ue是左手zup。
我先按blender来。x轴向右,z轴向上,y轴指向前方。
即右手zup。

2.5 ray

\({\Large r(t)=o+td}\)

2.5.1 RayDifferential

RayDifferential这个类在初期用不到。它在第10.1章用于Texture的相关计算。

见PerspectiveCamera::GenerateRayDifferential

目前只考虑lensRadius=0。
给原始ray的屏幕坐标x,y分别加1,生成两条辅助光。原点一样。
这样构建出距离1,后续可以方便地反向算出偏微分值。

2.7 transformation

各种矩阵变换

齐次坐标

平移
scale
分别绕xyz轴的旋转
绕任意轴的旋转
lookat
坐标是如何映射到最后的屏幕上的?

经典的mvp矩阵

  • local space 先有local坐标,即模型自己的坐标,建模时的坐标。

local坐标C

  • world space 我们创造各种程序、各种世界,每个世界都有各自的坐标系/原点。 把模型导入一个世界,就必须指定模型在这个世界中的位置、方向。 也就是把模型从原点做矩阵变换。经过平移/旋转等等放到世界里。

MODEL矩阵M

把坐标乘上M得到 M*C。world space里的坐标

  • view space 观察空间/相机空间 任何一个画面一定是基于一个观察位置和观察方向。或者比如说基于固定的相机或者眼睛。

把之前得到的世界坐标先对到相机的三个轴上。
也就是转换到相机空间。
也就是摆放相机的操作。
可以说后续的图像处理起来会比较方便。整体更统一、直观、耦合度低。
比如基于相机坐标,可以魔改出各种简单或复杂相机模式/镜头,搞出千奇百怪的照片。
View矩阵V
即lookat矩阵

此时坐标变成 V*M*C

  • clip space

到了相机空间,就要考虑怎么把各个坐标对应到最终的照片上。也就是一个投影的过程。
例如做成透视/正交图片。
Project矩阵P
本质就是相机参数的配置。
这个内容在pbr的第六章
再乘矩阵P,最后得到P*V*M*C
及常见的mvp变换。

lookat矩阵

https://learnopengl.com/Getting-started/Camera

lookat矩阵是如何推导的?找来找去没有一家给出。
都是啪一下就直接填好了矩阵。那么自己硬推试试。
本质上是一个通用的转换坐标系的公式。
给定相机在world space的坐标,dir,up,把任意world space坐标转换到相机坐标系。

坐标系是上z 右x 前y

经过相机原点的dir/right平面设为A。那么物体到A的距离,也就是物体到相机在up方向的距离,就映射为z。其他轴类似。
经过相机原点的up/dir平面设为B。那么物体到B的距离,就映射为x。
经过相机原点的up/right平面设为C。那么物体到C的距离,就映射为y。
画一下图。以up轴为例,可以发现,点到平面的距离,就是点到相机的向量在up方向上的投影长度。
投影长度正好就是向量的点积。
设C为camera,T为要转换的坐标。RUD分别为right/up/dir。
那么距离

\(L = dot(T - C, U)\)

算一下dot,得到
距离

\(L = Ux(Tx - Cx) + Uy(Ty - Cy) + Uz(Tz - Cz) = dot(U, T) - dot(U, C)\)

这个距离就是最终的z值。同理可算出最终的x,y。

\(\left\{\begin{matrix} x=dot(R, T)-dot(R, C)\\ y=dot(D, T)-dot(D, C)\\ z=dot(U, T)-dot(U, C) \end{matrix}\right.\)

这里的正负关系也是天然正确的。
T - C如果是指向up的大方向,即与up的平面夹角小于90度,即指向平面上部,是正值。
如果指向平面下部,是负值。
这和dot操作吻合,不用再处理符号了。
又发现dot操作正好和矩阵相乘非常匹配。
可以尝试往矩阵相乘上凑。如果凑出一个矩阵,能把原坐标转成最终的\(xyz\),这个矩阵就是lookat矩阵。

还是对于z,凑一个矩阵相乘。

\(\begin{bmatrix} Ux &Uy &Uz & 0 \end{bmatrix} \times \begin{bmatrix} 1& 0& 0 &-Cx \\ 0& 1& 0&-Cy \\ 0& 0& 1&-Cz \\ 0& 0& 0&1 \end{bmatrix}\)

计算结果为\(\begin{bmatrix} Ux &Uy &Uz & -dot(U, C) \end{bmatrix}\)
此时如果再乘上T即\(\begin{bmatrix} Tx \\Ty\\Tz\\1 \end{bmatrix}\)
正好能得到和上面一样的\(z=dot(U, T) - dot(U, C)\)
那么答案已经出来了。

其他两轴一样处理,最后得到lookat矩阵。

\(\begin{bmatrix} Rx& Ry& Rz& 0 \\ Dx& Dy& Dz& 0 \\ Ux& Uy& Uz& 0 \\ 0& 0& 0&1 \end{bmatrix} \times \begin{bmatrix} 1& 0& 0 &-Cx \\ 0& 1& 0&-Cy \\ 0& 0& 1&-Cz \\ 0& 0& 0&1 \end{bmatrix}\)

可以接着看看第六章Camera

2.10 相交

SurfaceInteraction的一些解释

面的normal问题

面的normal是如何确定的?blender里是可以指定每个面的normal的。
导出的gltf模型也是带normal的。我也做成模型指定normal?

3 形状

3.1 基本接口

shape基类

各种具体的形状比如三角、球、圆柱等等。
一般需要实现形状与ray的相交。

3.1.1 bbox

AABB

3.1.2 光与bbox求交

本质方法就是分别求光与bbox的三组对立面的交点。
每组得到t1/t2。如果三组t1/t2有交集,说明光穿过了bbox。
特殊情况是光的原点就在bbox内部,这样肯定是相交。

An Efficient and Robust Ray–Box Intersection Algorithm

某个垂直于x轴的平面

\(x=X_{1}\)

光打到某点的坐标

\(O+tD\)

那么光打到平面,只看x轴,有

\(O_{x} + tD_{x} = X_{1}\)

可求得

\({\Large t=\frac{X_{1}-O_{x}}{D_{x}}}\)

那么对两个对立面,按光的方向,可求得tmin和tmax。

这样子依次算三个轴,如果t有交集就是相交。

可以把\({\Large \frac{1}{D}}\)在三个轴上的值提出来先算好,再做乘法,更快。一般乘法会快一些。

可以把\({\Large \frac{1}{D}}\)提到求交函数外面,这样同一个光与大量bbox求交时,可以提速。

  • 特殊情况
    光在某个轴的方向为0。除的话为无穷大。而IEEE又有0和-0之分,除出来是正无穷和负无穷。会造成检测失败。
    不过我在VS2022里貌似没试出来,除0和-1都得到inf,并没有-inf。可能方法不对,暂忽略。
    最终应该再区分一下\(D\)的正负方向,按需交换tmin和tmax。
  • 错误处理见3.9.2

3.2 球体

球体与射线的交点

// 设s为球面上的交点
// 对于球 ||s - c||=r
// 对于ray s = o + dt
// ||o + dt - c|| = r
// c = 0, 0, 0
// (o + dt) dot (o + dt) = r*r  自己点乘自己是长度的平方
// 算出来得到 oo + 2tod + ttdd = rr
// 解t的一元二次方程
// a = dd
// b = 2od
// c = oo - rr
// att + bt + c = 0

解出t。即可算出交点。

从球坐标看

\(x = r sin(θ) cos(φ)\)
\(y = r sin(θ) sin(φ)\)
\(z = r cos(θ)\)
θ是向量与z轴夹角
φ是向量与x轴夹角

y / x = tan(φ)

可以算出交点的φ = atan(y / x)

phimax默认值是360

3.6 三角形mesh

Möller–Trumbore算法
加载mesh时就把所有三角形转到world。不同于球体等,球体是每次求交时把光转到local。

glTF模型

3.9 error处理

浮点数计算是有误差的。需要学习数值计算相关。

推荐的书
[Accuracy and Stability of Numerical Algorithms]
更正式。得看一下它前两章,相互参考。有的概念更容易懂。
误差例子
光1和球体先算出交点1,另有一个远处的光源和交点确定光2。
这个光2和球体再算相交,得到的交点2和交点1有误差,而事实上应是同一个点。
例如重复相交
先求得一个光与球体的相交,由于误差,交点可能在数值上是在球体表面以下。
这时再算反射光的话,光又会与这个球体假相交,就完蛋。
一个解决办法是加一个误差t。
比如重复相交的情况,如果算得的交点与第一个交点距离小于t,就忽略这个相交。
但如果反射方向与面越接近平行,那么这个t需要越大才行。
但t太大的话又可能把合法的相交剔除。
不可行。
误差处理很重要。因为会直接导致图像完蛋。
刚开始抄代码,做简单场景时忽略误差或者简单处理一下还能凑合。
但做到一定程度以后一定会出问题。必须严肃处理。

可以先跳过,出问题后再回来看。


3.9.1 浮点数计算

浮点数的表示

https://en.wikipedia.org/wiki/Single-precision_floating-point_format

浮点数是经过巧妙设计的。可以发展出一套方法,得到令人吃惊的小误差。

IEEE32位浮点数
1bit 存符号
8bit 存指数 -127~127
23bit 存有效数字
有效数字部分,以二进制normalized的格式来存。normalized开头一定是1。比如二进制小数1.1101011。
之前的0,没有意义,不用存。
然后因为开头一定是1,那又不用存。所以是从实际数据的第二位开始存。
所以23bit有效数字部分实际存了24bit的数据。第一位的1隐藏。

要知道IEEE浮点数和bit数据如何相互转换

比如1的最终bit数据是0x3f800000

0???
0 00000000 00000000000000000000000

1
0 01111111 (1.)00000000000000000000000  (1) * pow(2, 0)

0b0111111100000000000000000000000 = 0x3f800000

2
0 10000000 (1.)00000000000000000000000  (1) * pow(2, 1)

3
0 10000000 (1.)10000000000000000000000  (1 + 1/2) * pow(2, 1)

4
0 10000001 (1.)00000000000000000000000  (1) * pow(2, 2)

5
0 10000001 (1.)01000000000000000000000  (1 + 1/2*2) * pow(2, 2)
需要关注有效数字部分的最小粒度/最小单位。两个小数之间最小的差。
也就是二进制0.00000000000000000000001的值。
称作\(ulp\)/unit in last place/spacing。
它值为\(2^{-23}\)
约为十进制1.1920928955078125e-07
这也是cpp里std::numeric_limits<float>::epsilon()的值

https://en.wikipedia.org/wiki/Machine_epsilon

指数取255,小数全0时,表示正负无穷。
指数取255,小数非0,都表示nan。not a number
一个数值计算中如果出现了nan,结果一定是nan。
pbrt里用了大量std::isnan()防止出现nan。
需要认识到浮点数不是均匀分布的。
例如从1开始,每次加1个\(ulp\),直到把有效部分的23bit填满。再+1\(ulp\),有效数字清0,指数+1。实际值变成了2。
而从2开始,一样操作,最后再+1\(ulp\)。变成了4。
可见IEEE浮点数是在每对\([2^{x}, 2^{x+1}]\)之间划分了同样数量的间隔。整体并不是均匀分隔的。
指数越大跳跃越大。
间隔大小为\(2^{x-23}\)

算术计算
IEEE745标准保证了浮点数的加减乘除开方运算,在输入不变的情况下,一定会给出同样的结果。
并且这些结果是最接近正确值(如果有一个无限精度的算法能得出一个完美的正确值)。

就是说它保证了绝对的稳定性和极高的准确性。

以数轴来看,所有无限精度的准确值组成了整个数轴。
而有限精度的IEEE浮点数分布在数轴。
一个计算的无限精度正确结果可以落在任意一点。而IEEE保证每个计算的结果都能四舍五入到最近的一个浮点数,即rounding。

\({\large a \oplus b = round(a+b)}\)

\({\large a \ominus b = round(a -b)}\)

\({\large a \otimes b = round(a * b)}\)

\({\large a \oslash b = round(a /b)}\)

\({\large sqrt(a) = round(\sqrt a)}\)

对任意实数x进行rounding,记作\(fl(x)\)
[Accuracy and Stability of Numerical Algorithms]定理2.2给出
一定有\(fl(x)=x(1 \pm e)\),其中\(|e|<\frac{ulp}{2}\)

也就是说任意x一定可以rounding到一个相对误差在\(\frac{ulp}{2}\)内的浮点数。

证明一下
假设有x,有效数字为\(t\),指数部分为\(E\)。则\(x=tE\)。 有效数字一定从1开始。\(1\le t<2\)
-----|---------------x--------|--------------------|------------
    t1E             tE       t2E                  t3E
     |             |       |
\(tE\)的实际值在\(t1E\)\(t2E\)之间。
那么经过rounding,x一定为\(t1E\)\(t2E\)里更近的那个。

rounding的操作性质,离哪个更近就往哪边靠,决定了\(|fl(x)-x|\le\frac{t2E - t1E}{2}\)

\(t2E - t1E=ulpE\)

所以有 \(|fl(x) - x| \le \frac{ulp E} {2}\)

两边除以\(x=tE\)

\(|(fl(x) - x) / x| \le \frac{ulp}{2t}\)

因为\(t \ge1\),所以右边\(\le \frac{ulp}{2}\)

即右边的最大值是\(\frac{ulp}{2}\),而左边小于等于右边,所以左边一定\(\le\frac{ulp}{2}\)

即得到了相对误差的上限。


\(x(1 \pm e)\)是一个标准模型。

也可以用别的模型比如,\(fl(x)=\frac{x}{1 \pm e}\)也成立。

\({\large a\oplus b = round(a + b)}\)

最终会落在

\([(a + b)(1−e), (a + b)(1+e)]\)

其中\(e\)一定小于等于\(\frac{ulp}{2}\)。见上面的证明。

对于32位浮点数,\(e=2^{-24}= 5.960464477539063e-08\)
pbrt里用这个值作为MachineEpsilon
std::numeric_limits<float>::epsilon() / 2

到此可以看懂pbrt的EFloat构造函数是什么意思。
结合求一元二次方程的Quadratic函数
它把float sqrtd = std::sqrt(discriminant)转成EFloat,是这样操作的。
EFloat efloat_sqrtd(sqrtd, MachineEpsilon * sqrtd);
MachineEpsilon是相对误差上限,MachineEpsilon * sqrtd就是绝对误差上限。
最后存下了sqrtd本身和sqrtd的上下限。
EFloat的加减乘除也是直观的。就是计算时把上下限都考虑进去,取可能的最大最小值。
最后再分别NextFloatUp(加1个ulp)和NextFloatDown(减1个ulp)得到保守的上下限。

误差的传导

在IEEE标准的基础上,可以推理出一套方法来分析和限定浮点数计算中的误差。

相对误差和绝对误差

绝对误差就是round过后的结果与正确结果的差的绝对值。

相对误差就是绝对误差占正确值的百分比。

以四个浮点数相加为例

\({\large (((a ⊕ b) ⊕ c) ⊕ d)}\)

\({\large \subset (((a+b)(1 \pm e) +c)(1 \pm e) + d)(1 \pm e) }\)

\({\large = (a+b)(1±e)^{3} + c(1±e)^2 + d(1±e)}\)

为指数误差取一个更大但更好算的保守误差

\({\Large (1±e)^{n}\le 1±(n+1)e}\)

可算出最终误差为

\({\Large 4e|a+b| + 3e|c|+ 2e|d|}\)

可以看出相加的顺序对误差有影响。较小的数先加,误差更小。

\({\Large (1±e)^{n}\le 1±(n+1)e}\)

这个近似,仍不满意。

又有人发明出一个上限(见[Accuracy and Stability of Numerical Algorithms]3.1)

\({\Large 1+\theta n}\)

其中\({\Large |\theta n|\le \frac{n}{1-ne}}\)

\(n\)越大,误差累积越大。

最后取

\({\Large Gn=\frac{ne}{1-ne}}\)

之前的四数相加就变成

\({\Large \begin{aligned} & (((a ⊕ b) ⊕ c) ⊕ d) \\ & ⊂(a+b)(1±G_3) + c(1±G_2) + d(1±G_1) \\ & ⊂ a + b + c + d ± (a+b)G_3 ± cG_2 ± dG_1 \end{aligned}}\)

pbrt里的\(gamma(n)\)函数就是算\(\frac{ne}{1-ne}\)

对于乘法rounding

\({\Large \begin{aligned} & a(1±G_i) ⊗ b(1±G_j) \\ & ⊂ ab(1±G_{i+j+1}) \end{aligned}}\)

\((1 \pm G_i)(1 \pm G_j) = (1 \pm G_{i+j})\)


矩阵变换的误差

对Vector3进行tansform

\(\begin{bmatrix} m00& m01& m02& m03\\ m10& m11& m12& m13\\ m20& m21& m22& m23\\ m30& m31& m32& m44 \end{bmatrix} \times \begin{bmatrix} x\\ y\\ z\\ 1 \end{bmatrix} = \begin{bmatrix} newX\\ newY\\ newZ\\ 1 \end{bmatrix}\)

\(\begin{aligned} newX & = ((m00 ⊗ x) ⊕ (m01 ⊗ y)) ⊕ ((m02 ⊗ z) ⊕ m03) \\ & ⊂ (m00x(1±G1) ⊕ m01y(1±G1)) ⊕ (m02z(1±G1) ⊕ m03) \\ & ⊂ (m00x(1±G1) + m01y(1±G1))(1±G1) ⊕ (m02z(1±G1) + m03)(1±G1) \\ & ⊂ (m00x(1±G2) + m01y(1±G2)) ⊕ (m02z(1±G2) + m03(1±G1)) \\ & ⊂ (m00x(1±G2) + m01y(1±G2) + m02z(1±G2) + m03(1±G1))(1±G1) \\ & ⊂ m00x(1±G3) + m01y(1±G3) + m02z(1±G3) + m03(1±G2) \end{aligned}\)

可以把最后一项的G2提升为G3。把误差范围提高,不影响正确性。最后得到

\(⊂ (m00x + m01y + m02z + m03) ± G3(|m00x| + |m01y| + |m02z| + |m03|)\)

所以\(newX\)的绝对误差范围是

\(±G3(|m00x| + |m01y| + |m02z| + |m03|)\)

同理可得\(newY\)\(newZ\)的误差。


3.9.2 保守的ray-bounds相交

计算光与bounds的相交,允许假相交,但绝对不能漏掉一个。

\({\large t=(X_{1}-O_{x})\frac{1}{D_{x}}}\)

这个计算一个减一个除一个乘,误差为\(G_3\)。那么tMin和tMax都带一个\(G_3\)
可以给tMax加上两倍的\(G_3\),再做比较。
这样确保tMin和tMax之间加入了保守的误差,确保不会发生遗漏。

3.9.3 可容错的三角形相交

3.9.4 交点的误差限制

见图3.43。如果有一个计算出来的交点,以及误差上限,那么可以保证实际的交点一定在box里。

计算相交可能产生较大的误差。

例如算出某个交点,得到\(t\)值,带误差\(\delta\)

再生成光

\(P=O+tD\)

对于x轴有

\(x=O_x⊕(t±δ)⊗Dx\)

最后可到误差上限为

\(\gamma_1|O_x| + \delta_t(1 \pm \gamma_2)|D_x|+ \gamma_2|tD_x|\)

可以看出

1.这个误差不受控,我们计算出的t的误差是可控的。但\(O_x\),\(D_x\)可以是任意值,可能非常大。

2.对于\(t\)来说,通常计算相交要进行十次级别的运算,也就是能达到\(1±G10\)级别,也是很大的误差。

所以计算相交可能产生较大的误差。


需要想办法把误差控制在几个ulp之间。

有一种方法是先算一个交点a,再从交点a生成光计算第二个交点。

这样t一定非常小,趋近于0。所以3.13的误差值也会变小。

另一种方法是把交点直接投到面上。

这样不用很多计算。误差也很小。

对于不同shape,有不同的计算方法。

例如球形,球的半径r是已知的,正确的。

而算出来的点到球心的距离可能大于r或小于r,那么就把算出来的xyz以实际的球的r值,做缩放,以匹配r。

\({\Large xNew =\frac{xr}{\sqrt{ x^2 + y^2 + z^2}}}\)

\({\Large xNew = \frac{x⊗r}{\sqrt{x⊗x ⊕ y⊗y ⊕ z⊗z}}}\)

xyz本身带误差。所以有

\({\Large \begin{align} \LARGE & \subset \frac{x(1±e)⊗r}{\sqrt{x(1±e) ⊗ x(1±e) ⊕ y(1±e) ⊗ y(1±e) ⊕ z(1±e) ⊗ z(1±e)}} \\ & \subset \frac{xr(1±e)^2}{\sqrt{x^2(1±e)^2 ⊕ y^2(1±e)^2 ⊕ z^2(1±e)^2}} \\ & \subset \frac{xr(1±e)^2}{\sqrt{(x^2(1±e^2) + y^2(1±e)^2)(1±e) ⊕ z^2(1±e)^2}} \\ & \subset \frac{xr(1±e)^2}{\sqrt{(x^2(1±e)^3 + y^2(1±e)^3) ⊕ z^2(1±e)^2}} \\ & \subset \frac{xr(1±e)^2}{\sqrt{x^2(1±e)^4 + y^2(1±e)^4 + z^2(1±e)^3}} \\ \end{align}}\)

\(z^2\)\(G3\)提升为\(G4\)。有

\({\Large \subset \frac{xr(1±e)^2}{ \sqrt{ x^2 + y^2 + z^2} \sqrt{(1±e)^4} }}\)

除法也带误差。有

\({\Large \begin{align} & \subset \frac{xr(1±e)^2 } { \sqrt{ x^2 + y^2 + z^2} \sqrt{(1±e)^4}(1±e)} \\ & \subset \frac{xr(1±e)^2 } { \sqrt{ x^2 + y^2 + z^2} (1±e)^3} \\ & \subset \frac{xr(1±G5) }{ \sqrt{ x^2 + y^2 + z^2}} \end{align}}\)

最终需要带上\(1±G5\)的误差。


3.9.5 可容错的光的生成

例如要生成一个反射光线

先算出交点a,它的误差上下限可以限定出一个长方体,xyz都取大和最小作为对角。

准确的点就在这个长方体里。面是准确的,所以面一定与这个长方体相交。

如果不做处理,直接从a生成光,那么由于误差,a可能会再次与面相交。

现在的方法是先把a沿面的normal方向向外延伸,挪出可能的误差范围,让它不可能再与面相交。

再正常计算。

沿面的normal方向向外延伸的好处,挪动距离最少。而且阴影和反射的效果更好。

见OffsetRayOrigin函数

SpawnRay相关问题


加上错误处理后效果非常明显

处理之前图像有很多黑点noise,就是因为误差导致二次相交。处理完后一眼看去没有noise了。

思考

IEEE保证程序算出来的值与准确值的误差在e之内。

科学家给出一些e基础上的实用的保守近似。

我们拿来跟踪这些误差,算出累积。

在各个层面根据误差做出调整,规避问题。


4 加速

树状结构。前期可忽略,用暴力遍历。
现在主流是BVH。[STRIKEOUT:所以不看kd-tree了。]

4.3 BVH

  • 树形结构。从一个包含了所有物体的bbox为root开始。

  • 大致按质点的3d坐标划分所有物体。

  • 每个节点对应一个bbox,包住所属的物体。

  • 一般把一个节点划分为两个子节点,也就是划分为两组物体。

  • 子节点bbox一定在父节点bbox之内。所以整个子树也都在父节点bbox之内。

  • 物体最终只会存在叶子节点。

  • 一个叶子可能包含多个物体。例如质点非常近甚至重合的物体。

  • 不同节点的bbox可能相交。因为划分是按质点,无关物体实际大小。可能物体超大,就会出现相交。

  • 划分没有一个明确的最优解,一般是用概率等数学手段来得到较好的划分。 你想怎么划分都可以,只要你保证bbox那些数据计算正确,查询起来结果就正确,只是效率有差距。 不过即使是较为简单的划分也比暴力遍历强得多。

memory可控
每个基本形状只能在一个叶子里。总结点数限定在2n-1之内。
比kd-tree更容易创建,但计算光线相交相对慢一点。
数值方面更强,误差更小。
更容易更新,适合动态场景。
四种划分方法
SAH / HLBVH / Middle / EqualCounts
我们只看SAH

http://www-sop.inria.fr/members/Stefan.Popov/media/BVHPacketsOnGPU_IRT07.pdf

https://people.csail.mit.edu/amy/papers/box-jgt.pdf

https://jcgt.org/published/0007/03/04/paper-lowres.pdf

https://github.com/madmann91/bvh

https://www.scratchapixel.com/lessons/3d-basic-rendering/minimal-ray-tracer-rendering-simple-shapes/ray-box-intersection

https://github.com/jan-van-bergen/GPU-Raytracer

http://15462.courses.cs.cmu.edu/fall2016content/lectures/12_accelstructures/12_accelstructures_slides.pdf

pbrt是按前3个论文来实现的。
[Realtime Ray Tracing on GPU with BVH-based Packet Traversal]
[On fast Construction of SAH-based Bounding Volume Hierarchies]
[An Efficient and Robust Ray–Box Intersection Algorithm]

论文顺序

[Ray Tracing Deformable Scenes Using Dynamic Bounding Volume Hierarchies]
直接给出了用sah创建bvh的伪代码。查询和更新的流程也有。
0基础先看这个上手。可以直接看懂bvh是啥意思。重点看。
[Experiences with Streaming Construction of SAH KD-Trees]
对kd-tree引入binning方法。
我没看
[Realtime Ray Tracing on GPU with BVH-based Packet Traversal]
给出一个基于gpu的并行查询方法,和基于binning的bvh快速创建方法
我们要看其中的第5节binning方法。它只是讲了大概,可以看一下了解大概情况。
[On fast Construction of SAH-based Bounding Volume Hierarchies]
主要讨论了针对bvh的binning方法。
有具体流程。重点看。

pbrt默认就是用的binning方法。依据上面两个论文。

[An Efficient and Robust Ray–Box Intersection Algorithm]
包围盒求交

论文研究

[Ray Tracing Deformable Scenes Using Dynamic Bounding Volume Hierarchies]

bvh可以采用任何数量的分叉。这里采用二叉。轴对齐的AABB。

采用sah
光与某个box相交的概率=此box的表面积/root节点的box表面积
top-down顺序创建bvh
从三角形集合\(S\)开始,不断分成两组\(S_{1}\)\(S_{2}\),作为两个节点。
对每个节点再递归做处理。直到只剩最后一个节点,作为叶子。

科学家总结出一种光线与S求交的耗时表达

\({\Large T=2T_{AABB}+ \frac{A(S_{1})}{A(S)}N(S_{1})T_{tri}+ \frac{A(S_{2})}{A(S)}N(S_{2})T_{tri}}\)

\({\Large =2T_{AABB}+\frac{T_{tri}}{A(S)} (A(S_{1})N(S_{1})+ A(S_{2})N(S_{2}))}\)

\(T_{AABB}\)为光线与aabb相交的耗时
\(T_{tri}\)为光线与三角形相交的耗时
\(A(S)为S集合的包围盒的表面积\)
\(N(S)为S集合的三角形数量\)

每次分组的目的,就是要让T尽可能小。

\(S\)分为两组,有\(2^{N}\)种方法。

找到最佳的划分不现实,我们转而找一个较好的划分。

使用一系列轴对齐的平面来划分三角形。
对于一个平面,每个三角形可能在平面的左边、右边或重叠。
如果发生重叠,就不好区分左右,所以只考虑三角形的包围盒的质点。这样就能明确区分左右。
对于每个划分,都要用T来评估。
观察T,每次只需其中的变量\(A(S_{1})\)\(N(S_{1})\)\(A(S_{2})\)\(N(S_{2})\)就可以进行评估了。
## 划分S的大致伪代码。有些边界细节可能有误,暂时忽略。
## S为三角形合集,以每个三角形包围盒的质点代表三角形进行计算。
def partition(S):
    bestCost = max ## 最优开销
    bestAxis = None ## 最优的轴
    bestIndex = -1

    for axis in [x, y, z]: ## 分别对xyz轴
        sort(S, axis) ## 按当前轴排序

        ## S[i]为[三角形i/点i]的信息
        ## 设S[i].leftArea为i左边(不包含i)的总包围盒面积。
        ## 设S[i].rightArea为i右边(不包含i)的总包围盒面积。

        S1 = []
        S2 = S

        ## 遍历第一遍。算所有的leftArea
        for i = 0 to |S|:
            S[i].leftArea = Area(S1) ## 空集面积为无穷大
            S1.append(S[i]) ## 添加当前点到S1

        S1 = S
        S2 = []

        ## 遍历第二遍。算右边面积。同时评估sah。
        for i = |S| - 1 to 0:
            S[i].rightArea = Area(S2)

            ## 算SAH。只需A(S1),N(S1),A(S2),N(S2)四个值就可以进行对比。
            sah = SAH(S[i].leftArea, |S1|, S[i].rightArea, |S2|)
            if bestCost > sah: ## 如果得到更小的sah。记录下来。
                bestCost = sah
                bestAxis = axis
                bestIndex = i

            d把点i从S1挪到S2

    ## 到此已经得到了最好的划分

对于划分平面的挑选,作者尝试了平均分、用包围盒的面分、以质点分。得到的结果相近。

求交

如果只判断光与bvh是否相交。就很简单。

bool intersect(node, ray):
    if node is leaf:
        return node.triangle.intersect(ray)

    if intersect(node.left, ray) or intersect(node.right, ray):
        return True

    return False

如果要获得最近的相交,就不好搞。

首先bin的对应,是根据质点来的。
有个问题是如果某个三角形超大,它可能质点在最远处,但实际上是和光先相交。
这种怎么处理。最简单就是得检查所有相交的三角形。

pbrt貌似有优化,有点绕,一下子没看懂它怎么得到最近的,先不管。

可以想到的一种优化是。
深度优先,一旦找到一个确定的更近的相交,更新这个全局的t。
然后继续递归查询,用这个t可以把更远的bbox剔除,也就是可以提前从中间剪枝。
这样就省去了很多计算,不用找到所有相交。
还可以每次判断一下光的方向,决定先查左边还是右边。尽量顺着光的方向来查。
这样先找到正确答案的概率可能更大,剪枝可能更有效。
毕竟一般情况三角形bbox大概率还是按着bin的划分来分布的,远处的bin包含超大的三角形不是常态。
三角形发生变化时的bvh更新流程
每帧重建bvh太慢。
希望只做增量更新,更新少量数据,又能维持sah评估。但貌似还没有很好的研究成果。
那么退而求其次,三角形发生改变时,只更新bbox,不改变bvh的树形结构。
创建的bvh的时间点
一般从初始状态创建就挺好。比如游戏里人物的初始T pose。
看情况也需要找一些更好的状态来创建bvh。
讲述了一些评估bvh创建点的方法。

论文研究

[Realtime Ray Tracing on GPU with BVH-based Packet Traversal]
看第5节
基础方法需要不断在每个轴上排序所有三角形。
挺慢的。
有人针对kd-tree发明了binning方法。
现在发现binning方法在bvh上也适用。可提速、减少内存使用。
思想是对每个轴只遍历一遍三角形。把质点放进bin或者bucket里。
记录count等信息。可能延伸到多个bin。
最终的bin信息可以用来评估sah所需的值。
和原始方法类似。对bin信息也要遍历两次。
第一次从左到右。
第二次从右到左。
bin数量是可调的。量少速度快,量多更准确。
不能太多,论文的环境中是最好256个,与64kb的L1 cache匹配。
目前硬件又进化了十几年情况可能不一样。

这个论文貌似没具体讲质点如何映射到bin??只能去看它参考的论文。


论文研究

[On fast Construction of SAH-based Bounding Volume Hierarchies]
引入binning。与上一篇论文是同期的,互有参考。
bvh比kd-tree更好
节点更少
只使用三角形质心,不用管overlap。
节点限定在2n-1。内存友好。
  • 每个bin包含的数据?
    count和bbox
  • bin的数量?
    设定bin数量为\(K\),可按需设置,4和16都是可行的。pbrt定为12。
  • 选轴
    每次划分先遍历所有三角形算出整个bbox。
    根据整个bbox,选出跨度最大的轴,只在这个轴上划分。
  • 如何生成bin信息?
    遍历三角形集合。算出每个质点在选定轴上的最大跨度中,所处的位置。
    根据这个位置确定bin的index。
    比如处于跨度1/4处,有12个bin,那么就选1/4处也就是第3个bin。
    确定index后,对应bin的count+=1。
    并且把当前三角形的bbox和对应bin的bbox进行融合,
    如果bin的bbox包不住三角形bbox,就要扩大bin的bbox以包住三角形。
    所以最后每个bin的数据是质点落在对应范围内的三角形的总数和总的bbox。
  • 如何评估sah?

    ## 大致伪代码
    final_cost = max
    
    for i = 0 to |bin|: ## 遍历所有bin
    
        count_L = 0
        count_R = 0
        bbox_L
        bbox_R
    
        ## 算左边的count和总的bbox
        for j = 0 to i:
            count_L += bin[j].count
            bbox_L += bin[j].bbox
    
        ## 算右边的count和总的bbox
        for j = i to |bin|:
            count_R += bin[j].count
            bbox_R += bin[j].bbox
    
        cost = (count_L * bbox_L.area() + count_R * bbox_R.area()) / 总bbox面积
    
        if final_cost > cost: ## 得到了更好的cost
            final_cost = cost
    
    ## 到此得到了最好的cost
    
可以看到binning方法不用每次对三个轴排序。只选定跨度最大的轴。
遍历一遍三角形,按质点的跨度对应到bin,生成bin信息。
每个bin是一组三角形。
\(K^2\)复杂度,以每个bin计算sah。得到最好的划分。
总体不是找最优解,而是找一个较好的解。

4.3.1 BVH的创建

见上述论文

4.3.2 SAH

见上述论文

4.3.3 HLBVH

hierarchical linear bounding volume hierarchy


4.4 kd-tree

还是得实现kd-tree。光子映射会用到。
其实是比较容易的。还可以用暴力计算来验证正确性,比较舒服。

创建:

  • 对一堆三维空间点进行划分,实现logN复杂度查询某个点最临近的k个点。

  • 把所有点按一排数组看待,不断从中间分成两段,作为左右子树。

  • 每层按某个坐标轴划分,xyz依次来。

  • 左边的子树的坐标都小于等于当前节点,右边都大于等于当前节点。

  • 每次找位于正中间的点,用快速选择算法,见上面链接和算法导论。

  • 找中间的点,确保了树的平衡。也确定了每个节点的index。 子树index也很好算,左边就是i*2+1,右边就是i*2+2。 所有节点可存在一个数组里。

knn查询:

  • 维护一个类似priority_queue的数据,存候选点和距离,最多存k个候选。 能快速删除其中的最大距离。

  • 维护候选数据中的最大距离。

  • 从index0根节点开始。如果候选未满,插入当前节点。更新候选数据。

  • 正常情况下递归搜索左右子树。

  • 如果此时候选已满,可能有机会剪枝。 比如此时按x轴分成左右。 如果目标点在当前节点的左边,目标点与当前点在x上的距离为A。 那么根据划分的性质,当前点的右子树中所有的点与目标点的三维距离一定大于A。 如果此时候选已满,且候选最大距离比A要小,那么右子树就没必要再查了。


5 色彩和辐射度量学

spectral power distribution
光谱能量分布SPD
主要问题是如何表示spd
1.RGBSpectrum
2.SampledSpectrum

5.2

5.2.1 XYZ色彩

5.2.2 RGB色彩

5.3 RGBSpectrum

RGBSpectrum虽不完美但简单高效
用带权重的红绿蓝信息的和表示。
同样的rgb值在不同的显示器上效果可能有较大差别。也就是产出的spd不一致。

5.4 Radiometry 辐射度量学

pbr的基石。建立在物理学原理基础上。
辐射度量学+几何光学
如果深挖,能去到量子力学,好在图形学一般不用走到这个层面。

几何光学已经足够建立起pbr的光模型。

包含的性质:
Linearity
两个光系统的组合等于他们的sum
energy conservation
光发散时,能量不会凭空增加。
No polarization
不考虑光的偏振
No fluorescence or phosphorescence
不考虑荧光和磷光
Steady state
默认自然界的光已经达到一种稳定态
辐射分布不再发生变化

光的衍射和干扰很难实现

5.4.1 基本的物理量

https://en.wikipedia.org/wiki/Radiometry

energy Q

能量
单位J焦
光源发射出光子。每个光子都有波长属性,带着一定能量。
所有基本的辐射学物理量都是在测量光子。

一个光子所带的能量

flux/Radiant flux

\({\Large \Phi=\frac{dQ}{dt}}\)

辐射通量/辐射功率
单位瓦

功率按时间积分得到能量

Irradiance / Radiant exitance

辐照度/辐射出射度

单位面积上的辐射功率。一个是入,一个是出。

\({\Large E = \frac{d\Phi}{dA}}\)

A是趋近于0的,所以相当于一点上的辐射功率。

点光源

\({\Large E = \frac{\Phi}{4πr^2}}\)

这里可以看出能量按距离的平方衰减。

按面积积分得到辐射功率。

对这个积分的思考。
是对什么进行积分?是对哪个函数进行积分?
对E函数积分。E函数是什么形式?变量是什么?
E(p, a)。变量是点p和球体每块面积。
由于现在考虑的是均匀的情况,也就是E与a无关,在任何一块面积上的E值都相等。
a可以从E里拿掉。E可以看成常数。
如果是一个复杂的光源,在不同方向上的E不一样,那么E就是关于a的函数。

solid angle / intensity

立体角的单位是球面度sr
一个完整的球体是4π个sr
图5.9
一个物体c对于点p的立体角,角的值实际上等于c投影在点p为圆心的单位球体表面的面积。
想象一下如果一个物体把点p全包围,那么投影面积就是\(4πr^2\),对于单位球体就是4π。

辐射强度intensity,为单位球面度的辐射功率,和Irradiance类似。

\({\Large I=\frac{d\Phi}{dω}}\)

intensity按立体角积分得到Φ

Radiance

最后最重要的辐射率

Irradiance给出的是单位面积(一点)上的辐射功率,但是没有区分方向。

Irradiance再除一个dω,得到辐射率。
即功率在微小面积上再考虑微小立体角。

\({\Large L=\frac{dE}{dω}}\)

\({\Large L=\frac{d^2\Phi}{dAdω}}\) (这里书上写错了,漏掉了2次)

辐射率将被大量使用。
它是其他物理量的基本,从辐射率可以积分得到所有其他物理量。
且在真空是个定值。

5.4.2 入射/出射辐射率函数

在不同物体表面的连接处,辐射率函数一般不连续。
比如完全不透明的镜面,从镜面往下很微小的位置,函数就完全是两码事了。

\({\Large L_i(p, ω)}\)

5.4.3 Photometry 光度学

https://en.wikipedia.org/wiki/Photometry_(optics)

所有辐射度量学的物理量都可以转化成对应的光度学物理量。通过对spectral response curve积分。


5.5 辐射度量学的积分

渲染中最常见的操作是对辐射度量学物理量进行积分。

\({\Large L=\frac{dE}{dω}}\)

有一点p,对应一个微小面积,有一束光从某个方向ω打过来。
根据物体的材质特性,最终得到的光。就是\(L_i(p, ω)\)

\(L_i(p, ω)\)描述了点p上不同ω上的光能。

然后我们关心的是半球上所有ω的\(Li(p, ω)\)的和。
也就是方程5.4里的积分。
其中的cosθ是因为Lambert’s Law。
每束光实际产生贡献的,都只是它在normal方向的分量。
初学的话不太好理解这个积分。
到底有几束光?之前是在考虑一束光,但这里又来个积分。
只能按多束光来想。那么就是把所有光加起来。是一个和。
和就能转成积分形式,只能理解成先把所有光处理一下,得到一个统一的ω的函数?
\(Li(p, ω)\)描述的是所有的光源从任意ω方打到p点的情况,ω是变量,并不是一个定死的方向。

积分考虑的是所有光。写成积分的形式更优雅统一。并且后续可以用蒙特卡洛方法处理。

https://zhuanlan.zhihu.com/p/33464301

后面brdf也是用这种积分形式,得统一理解。

5.5.1 对立体角积分

5.5.2 球坐标积分

球坐标。用两个角度来表示三维向量
一般情况,对球坐标积分,比对立体角积分更方便。
可以用球坐标积分替代立体角积分
球坐标的微小面如何表示。
立体角的定义就是它从圆心投射在单位球体表面的面积。
所以\(dw = dA\)
\(dA\)又用微小长方形表示。
按周长的比例算
设纵向微小长度为y,有

\({\LARGE \frac{y}{2πr} = \frac{dθ}{2π}}\)

>\({\Large y=rdθ}\)

横向微小长度x。θ处的横切面半径为\({\Large rsinθ}\),有

\({\LARGE \frac{x}{2πrsinθ} = \frac{dφ}{2π}}\)

\({\Large x=rdφsinθ}\)

针对单位球体,r=1

\({\Large dw= dA = xy=r^2sin\theta d\theta d\phi = sin\theta d\theta d\phi}\)

最终的半球积分形式就是\({\Large sinθdφdθ}\)\(\theta\)从0到\(\pi/2\)积90度,\(\phi\)从0到\(2\pi\)积一圈。

\({\Large E(p, n) = \int_{0}^{2\pi } \int_{0}^{\frac{\pi}{2} } Li(p, θ, φ) cosθ sinθ dθ dφ}\)

5.5.3 转化为对面积积分

对立体角积分还可以转化为对面积积分。可用于简化计算。

例如对于方程

\({\Large E(p, n) = ∫Li(p, w)|cosθ|dw}\)

是针对一定范围的立体角积分。

如果要针对一个四边形求E。那么用立体角就不好算。

见图5.17

首先按比例算

\({\Large dw/4π = dA/4πr^2}\)

\({\Large dw = \frac{dA}{r^2}}\)

一般情况dA不是正对着点p,要加上\(cosθ_o\)

\({\Large dw = \frac{cos\theta_o dA}{r^2}}\)

替换后方程变为

\({\Large E(p, n) = ∫_{A}L |cos\theta_i| \frac{cos\theta_o dA}{r^2}}\)


5.6 表面反射

光打在表面上,表面会使其发散。部分光会返回到环境。
此模型的2个重点
1. 反射光的光谱分布
2. 方向分布
例如柠檬
柠檬会吸收大量蓝光,反射红绿光。所以显示成黄色。
这是基于光谱分布
同时对于某些角度,可看到柠檬表面的高光,即趋近于白光,并没有哪种波长的光被大量吸收。
这是基于方向分布
例如镜面
基本完全基于方向分布

透明物质比较复杂,要考虑光在物质内部(表面以下)传播,入射点和出射点可能有一定距离。

2个模型
- brdf
相对简单高效,如果不用考虑光在内部的传播,就用brdf。
  • bssrdf brdf的一般形式。 可更精确地描述透明物质。

5.6.1 BRDF (bidirectional reflectance distribution function)

先有入射光\(Li(p, ωi)\)。想要得到某个出射方向ωo的光\(Lo(p, ωo)\)

一个入射光ωi,可能产生n多出射光,能量被分散。
我们关心的是出射到眼睛的那一部分光,即ωo方向上的光。
就有f(ωo, ωi),这个是由材质的特性决定的。是一个类似百分比的东西。
不同材质都必须建模,给出一个对应的f函数,反映出材质的特性。

这里的积分又跟之前一样,考虑的是所有方向上的所有光对ωo的贡献。

第八章详细讲brdf。

重要性质
1. 函数在某点p,交换w0和w1,值相等。
2. 出射能量小于等于入射能量。即半球范围f的积分小于等于1。
>\({\Large \int_{H^2(n)}f_r(p, \omega_o, \omega')cos\theta'd\omega' \le1}\)

其中\(\theta'\)是跟随每一个\(\omega'\)的。

就是说给定入射角\({\Large \omega_o}\),半球范围内所有可能的出射方向的贡献的和小于等1。
brdf的设计必须符合这个规定。

最后形成了重要的方程

\({\Large L_o(p, \omega_o)= \int_{S^2}f(p, \omega_o, \omega_i)L_i(p, \omega_i)|cos\theta_i|d\omega_i }\)

半球范围内所有的入射光对\(\omega_o\)方向的光能,经过brdf处理,再积起来,就得到了p点向\(\omega_o\)方向的光能。


5.6.2 bssrdf


6 camera

6.2 相机投影

相机空间的坐标如何投影到最终的2d屏幕

  • screen space
    z值在0-1之间
  • ndc space
    见图6.1
    在xy平面从(0, 0)到(1, 1),在z方向的范围是(0, 1)。

ndc这个东西也是各家可以自己定的。有的xyz是从0到1,有的xyz是从-1到1。

pbr的ndc是从0到1
d3d是0到1
opengl是-1到1

我先以-1到1来计算

对于不同的ndc和相机坐标的左右手,结果矩阵是不一样的。

这个坐标系挺坑。一般资料都把相机look方向作为z轴。
我也是按这个推导的。
如果最终的场景坐标系定的是zup。那就还要转换。
可以重做矩阵。或者矩阵不变,变换完换轴。

6.2.1 正交

较为直接的mapping。没有近大远小的立体效果。

正交的输入,基本就是定义一个立方体。立方体之外的点都丢掉。

保留距离和角度关系。平行的线显示出来还是平行。

坐标在各个轴计算之于立方体范围的比例,然后把比例乘到ndc即可。

相机坐标系下。给定左右边界lr,上下边界tb,near平面n,far平面f。

\({\Large \frac{X-l}{r-l}=\frac{newX+1}{2}}\)

就是按比列投到-1和1区间

算得

\({\Large newX=2x/(r-l) - (r+l)/(r-l)}\)

同理

\({\Large newY=2y/(t-b) - (t+b)/(t-b)}\)

\({\Large newZ=2z/(f-n) - (f+n)/(f-n)}\)

那么可以凑出变换矩阵

\(\begin{bmatrix} 2/(r-l) & 0& 0& -(r+l)/(r-l)\\ 0& 2/(t-b)& 0& -(t+b)/(t-b)\\ 0& 0& 2/(f-n) & -(f+b)/(f-n)\\ 0& 0& 0& 1 \end{bmatrix}\)

这个矩阵乘 \(\begin{bmatrix} x\\ y\\ z\\ 1 \end{bmatrix}\)能得到\(\begin{bmatrix} newX\\ newY\\ newZ\\ 1 \end{bmatrix}\)

又有一般\(r/l\)\(t/b\)都取对于相机原点对称。就是取景框框是对称的。
\(r+l=0\)
\(r-l=2r\) (即立方体宽)w
\(t+b=0\)
\(t-b=2t\) (即立方体高)h

最终的正交矩阵

\(\begin{bmatrix} 2/w& 0& 0& 0\\ 0& 2/h& 0& 0\\ 0& 0& 2/(f-n)& -(f+b)/(f-n)\\ 0& 0& 0& 1 \end{bmatrix}\)

6.2.2 透视/立体

近大远小。不保留距离和角度关系。平行线也不再平行显示。

透视的输入,fov,n,f。
fov形成一个frustum椎体。
目的是把frustum里面任意一个坐标的xy投到near平面上,得到坐标Xp,Yp。
从上方视角来看,根据相似三角形关系。
>\(Xp/X = n/Z\)
>\(Xp = Xn/Z\)
同理侧面来看
>\(Yp/Y = n/Z\)
>\(Yp = Yn/Z\)
然后算Xp和Yp在near平面的比例关系。Xp和Yp代入之前的newX,newY。
>\(newX=2x/(r-l) - (r+l)/(r-l)\)
\(=2Xn/Z(r-l) - (r+l)/(r-l)\)

\(newY=2Yn/Z(t-b) - (t+b)/(t-b)\)

\(newZ\)的情况。一开始有点难理解。查了一圈有好几种说法。

如果按简单的线性映射,再凑矩阵。

\((Z-n)/(f-n) = (newZ+1)/2\)
\(newZ=2Z/(f-n) - (f+n)/(f-n)\)

发现很难像正交那样凑出来。

后来有了新的理解
首先x和y轴的映射需要有一个被广泛接受的“正确操作”,因为是直接给人看的2d画面。一般就是按比例线性映射。
而z轴的映射,并不需要存在一对一的正确值。只需确保顺序不变即可。
所以顺序不变的基础上,你能凑一个简单的矩阵就行。

\(xy\)正常凑。填上AB,再从结果反解AB。

\(\begin{bmatrix} 2n/(r-l) & 0& -(r+l)/(r-l)& 0\\ 0& 2n/(t-b)& -(t+b)/(t-b)& 0\\ 0& 0& A& B\\ 0& 0& 1& 1 \end{bmatrix}\)

注意\(newX\)\(newY\)都和Z相关,是1/Z的关系。

如果用\([X, Y, Z, 1]\)去乘,得到的结果需要/Z才是最终映射值。

换句话说得到的结果是最终映射值乘Z。

所以映射方程是\(AZ+B=(newZ)Z\)。这样映射,顺序是不会打乱的,所以正确性没问题。

对于x和y也是一样的。

就是统一构造了一个Z。

A下面的1也是这个原因,为了构造Z。

这样用\([X, Y, Z, 1]\)乘出来就是\([ZnewX, ZnewY, ZnewZ, Z]\)。最后统一除以Z,就得到了最终的(-1,1)的映射。

对于\(AZ+B=(newZ)Z\)

已知n映射到-1,f映射到1。

填上nf两个已知的映射。就能解出AB了。

\(An+B=(-1)n\)
\(Af+B=(1)f\)

解得

\(A = (f+n)/(f-n)\)
\(B = 2fn/(n-f)\)

得到Perspective矩阵

\(\begin{bmatrix} 2n/(r-l)& 0& -(r+l)/(r-l)& 0 \\ 0& 2n/(t-b)& -(t+b)/(t-b)& 0\\ 0& 0& (f+n)/(f-n)& 2fn/(n-f)\\ 0& 0& 1& 1 \end{bmatrix}\)

再有

\(r+l=0\)
\(r-l=2r\) (即n平面宽)w
\(t+b=0\)
\(t-b=2t\) (即n平面高)h

\(\begin{bmatrix} 2n/w& 0& 0& 0 \\ 0& 2n/h& 0& 0\\ 0& 0& (f+n)/(f-n)& 2fn/(n-f)\\ 0& 0& 1& 1 \end{bmatrix}\)

高和宽与fov/aspect ratio相关。fov可分水平和垂直,我们用水平fov。

\(tan(fov/2) = w / 2n\)
\(w = 2ntan(fov/2)\)
\(h = w / asp\)

得到用fov和asp表示的矩阵

\(\begin{bmatrix} 1/tan(fov/2)& 0& 0& 0 \\ 0& asp/tan(fov/2) & 0& 0\\ 0& 0& (f+n)/(f-n)& 2fn/(n-f)\\ 0& 0& 1& 1 \end{bmatrix}\)

另一种角度

http://learnwebgl.brown37.net/08_projections/projections_perspective.html

正常的坐标系,坐标可以是无限大小。
把一批跨度很大的坐标映射到(-1, 1),必然会出现误差,造成实际不同z映射出来相同。

这里他是从误差的角度去看,直接引入非线性映射。

越近的点变化幅度越大,所以误差越小。同理越远的点允许误差越大。


7 采样和重建

pbr的输出是一对离散的像素。但实际的辐射率/光在film上是连续的。
离散像素的计算直接影响到输出效果。
稍微出点错图像就会不好。反过来微小的改动也可能造成较大的改进。

采样理论。从连续函数上采集样本,再用这些样本重建相似的函数。

Sampler的实现

Filter类决定如何blend像素附近的像素,并确定最后的像素。

Film类将样本进行累积,计算像素。

7.1 采样理论

数学上我们是用连续的函数描述图像。
但最终的输出是离散的像素。
如何从连续的函数中得到离散的像素。
几乎唯一的方法就是用光线追踪对函数进行采样。

可以对每个像素进行一次精确的采样,直接作为最终输出。

但是如果想要得要较好的结果,通常要对一个像素的周边区域进行额外的采样。
进行处理后得要最后的像素值。
采样和重建必须借助各种近似,势必带来误差,最终造成锯齿/闪烁等。
根本原因就是采样的原理决定了无法完整地还原原函数。多少会有信息丢失。

7.1.1 频域和傅立叶变换

函数的频域/频率。
波动越频繁,频率越大。
大多数函数都可以分解为一系列正弦余弦函数的和。
傅立叶转换就是这一转换过程。

转换出来的信息就是原函数的频域分布。

频域可以从另一个角度看函数的性质。

傅立叶逆变换

Dirac delta distribution

7.1.2 理想化采样和重建

impulse train/shah/Ш符号
采样函数
描述采样过程。

首先有\(\delta\)函数

https://mathworld.wolfram.com/DeltaFunction.html

https://en.wikipedia.org/wiki/Dirac_delta_function

\(\delta(x) = \left\{\begin{matrix} +\infty, &x=0 \\ 0, &x \ne 0 \end{matrix}\right.\)

\(\int_{-\infty}^{\infty }\delta (x)dx=1\)

它不是一个普通意义的函数。实际上是一种分布。在0处一柱擎天。

shah函数/Dirac comb

\({\Large Ш_{T}^{}(x)=\sum_{i=-\infty}^{\infty} \delta (x-iT)}\)

https://mathworld.wolfram.com/ShahFunction.html

https://en.wikipedia.org/wiki/Dirac_comb

\(\delta\)函数这样处理就变成一个周期性的分布。就像一个梳子,所以用Ш这个符号。

\(i\)为所有整数。

\(T\)为周期,也作为采样率。

如果再乘上一个\(f(x)\),变成

\({\Large Ш_{T}^{}(x) f(x)=\sum_{i=-\infty}^{\infty} \delta (x-iT) f(iT)}\)

一般会取unit impulse即\(\delta(0) = 1\),刚好就是是对\(f(x)\)进行采样。见图7.4。

选一个重建过滤函数\(r(x)\),进行卷积操作,可以用这些采样点生成结果函数\(\tilde f\)

\((Ш_T(x)f(x))\otimes r(x)\)

convolution 卷积

\(f(x)\)\(g(x)\)

\(g(x)\)按y轴翻转,即变成\(g(-x)\)

得到两个确定的函数/图形。

这时给\(g(-x)\)一个变量\(t\)作为x轴上的位移。有\(g(t-x)\)

那么随着\(t\)的改变,\(g(-x)\)就能在x方向滑动。

对于某一个\(t\),可计算\(f(x)g(t-x)\),得到一个变量为\(x\)的函数。把这个函数对所有\(x\)积分,就得到了\(t\)对应的卷积值。

对于每个\(t\),都有一个卷积值,那么就形成了一个函数\(Juan(t)\)

就得到了卷积\(Juan(t)=f(x)\otimes g(x)\)

写成积分形式

\({\Large Juan (t)=\int_{-\infty}^{\infty}f(x)g(t-x)dx =f\ast g}\)

Nyquist frequency

https://www.shadertoy.com/view/XlcSz2

https://www.shadertoy.com/view/MdjGR1

接下来就是一些图像/信号处理的知识,涉及傅立叶变换、各种滤波等等。可以大胆地说出,看不懂。

至少学一些基本知识,有个大致概念。其他有缘再见。

可以直接转10.4.3看代码。看看怎么对图像进行resampling,然后生成mipmap。

7.1.3 锯齿


7.2 采样接口

7.3 分层采样

7.8 图形重建

7.8.1 滤波函数

https://en.wikipedia.org/wiki/Sinc_function

https://en.wikipedia.org/wiki/Lanczos_resampling

sinc滤波是一种理想的滤波,可以去除某个给定频率之上的部分,且不影响低频部分。

\(sinc(x) = \frac{sin(x)}{x}\)

实数域上的积分为\(\pi\)

进行normalize后

\(sinc(x) = \frac{sin(\pi x)}{\pi x}\)

实数域上的积分为1

\(sinc(0)=1\)

lanczos重采样

lanczos是一种低通滤波,同时可用于样本之间的平滑插值。

代码实现:

\(L(x)=\left\{\begin{matrix} \frac{sinc(\pi xa)}{\pi xa} \frac{sinc(\pi x)}{\pi x} & x \le 1\\ 0& x>1 \end{matrix}\right.\)

\(a\)取2。

Gamma Correction


8 反射模型

brdf/btdf统称bsdf

brdf对反射建模,btdf对出射建模,是对于一点的。
至于光在物体内部(表面之下)的传输,得bssrdf来。
反射数据的来源
- 测量
真实世界测得的物理量
  • 对于现象的建模
    根据现象/感觉建模。比如常见的roughness。
  • 模拟
    基于物理/几何等的模拟
  • 光波等原理。物理光学
    较为深层次的模拟。更耗资源,且通常并没有得到更好的效果。
  • 几何光学

基本术语

反射
描述一束光打到点上,出射光的情况。

主要分4类

  • diffuse
    光向所有方向发散
    如果是均匀的话是完美diffuse,实际不常见。
  • glossy specular
    光向一个较小的范围发散。高光。
    比如塑料。可显示模糊的倒影。
  • perfect specular
    光向一个极小的范围发散。 完美镜面反射。镜子
  • retro reflection
    原路返回
reflectance 反射率
发射率分布函数
isotropic/anisotropic
多数物体是各向同性的
意思是选一个点,绕normal旋转,光的分布不变。

各向异性的例子。衣物,cd表面,brushed metal。

几何设置

对于物体表面的一点,设stn三个向量,对应xyz轴。
所谓shading coordinate system
对应代码里BSDF类的ss/ts/ns。
在球坐标系。
对于任意向量W,可算出它与n之间的夹角θ。

\(cos(θ) = dot(w, n) = w.z\)

针对W可做一系列helper函数。求一些常用值。

在模型的Intersect函数里。会计算各种微分量。
例如sphere,在球坐标系。
\(x = r sin(θ) cos(φ)\)
\(y = r sin(θ) sin(φ)\)
\(z = r cos(θ)\)
θ是向量与z轴夹角
φ是向量与x轴夹角
交点处tan(phi)=y/x
float phi = std::atan2(hitPoint.y, hitPoint.x);
float theta = std::acos(Clamp(hitPoint.z / radius, -1, 1));
然后
float u = phi / phiMax;
float v = (theta - thetaMin) / (thetaMax - thetaMin);
把总量的百分比u和v当作微分量。
为什么是百分比?这个初学是不清楚的。
后续才懂了,因为就是要算百分比。uv是所谓参数化坐标,可用在texture处理中。
比如做贴图时要计算物体面上每个点对应的texture值,
也就是算每个点对应的texture上的横竖坐标上的百分比,即uv。

现在要算u变化时,xyz对应的变化率。也就是分别对xyz求u的偏导数,dp/du。

  dx/du
  = d (r sin(θ) cos(φ)) / du
  = rsin(θ) d cos(u*phiMax)/du         // float u = phi / phiMax;代入
  = -rsin(θ)sin(φ)phiMax
  = -y * phiMax

  dy/du
  = d(r sin(θ) sin(φ)) / du
  = rsin(θ) d sin(u*phiMax) / du
  = rsin(θ)cos(φ)phiMax
  = x * phiMax

  dz/du
  = d (r cos(θ)) / du
  = 0

  Vector3f dpdu(-phiMax * hitPoint.y, phiMax * hitPoint.x, 0);

同样可算出dp/dv

  Vector3f dpdv = Vector3f(hitPoint.z * cosPhi, hitPoint.z * sinPhi, -radius * std::sin(theta))
    * (thetaMax - thetaMin);

interaction的n和shading.n都初始化为Cross(dpdu, dpdv)。

BSDF的初始化中

ns(si.shading.n),
ss(Normalize(si.shading.dpdu)),
ts(Cross(ns, ss))

其实就是dpdu和dpdv和他们的Cross

一些要点
- 入射光Wi和出射光Wo都是normalize的。都从点指向外部。
Wi从点指向光源
  • n指向外部。
  • 本地坐标

  • brdf/btdf的实现不应关心Wi和Wo是否在同一半球。


8.1 基本接口

f函数。也就是brdf。

8.1.1 反射率

hemispherical-directional reflectance半球-方向反射率

\({\Large \rho_{hd}(w_o)=\int_{H^2(n)} f_r(p, w_o, w_i)|cos\theta_i|dw_i }\)

一束入射光\(w_o\)对半球范围所有方向(\(w_i\))的贡献(\(f_r(p, w_o, w_i)\)),加起来。

fr是材质的特性,给定光的入射方向wi,给出从wo方向射出的光能。wi与wo可互换。
只考虑normal方向的实际贡献,即cosθi。
见代码中的BxDF::rho()

一般都是用蒙特卡洛算积分。

不带方向的反射率为ρhh

d

8.1.2 BxDF SCALING ADAPTER


8.2 SPECULAR REFLECTION AND TRANSMISSION

镜面的反射和传输相对简单。

完美的镜面反射就是一个ωi对应唯一一个ωo。入射角和出射角相等。
\(θ_i=θ_o\)
\(φ_i=φ_o+π\)

传输就是折射

Snell’s law折射定律
index of refraction折射指数/折射率

\({\Large \eta_isin\theta_i = \eta_t sin\theta_t}\)

一般情况,折射率和光波长相关,所以一束光打进去会有不同方向的光打出来。即dispersion色散。
图形学一般忽略色散。太难算,而且忽略并不影响图像的准确。
\(\eta\)是折射指数。常见的空气为1,水1.3,玻璃1.5,钻石2.4。
根据公式,光从低折射指数的物质往高指数的物质打,方向就会靠近normal。反过来的话路线一样。

8.2.1 FRESNEL

fresnel效应
光在不同材质的交接处发生反射和折射,还会被吸收。反射和折射都会分掉能量,此消彼长。
入射光与normal的夹角对折射率有影响。
视线越靠拢normal,反射越少,折射越多。也就是越能看透材质。
反之视线与normal夹角越接近90度,反射越多,折射越少。越看不透,越趋近于镜面。

例如落地玻璃,正对着看,主要为折射。而越平着看,越体现为反光。


折射率又和物质的导电性质有关

dielectrics 绝缘体
折射率一般为1到3。
玻璃/矿物油/水/空气
conductors 导体
价电子可自由移动,从而导电。这种材质不透明,会反射大量的光。
少量的光会传进材质,被快速吸收,通常在0.1μm之内。
所以只有很薄的一部分材质能传输光。
我们忽略这个传输。
Semiconductors 半导体
硅/锗。不考虑这种材质。
一组fresnel方程可同时适用于导体和绝缘体
其中对于绝缘体可以做简化。

给定两种材质的折射率,入射角,fresnel方程可算出反射和折射的比例。

fresnel方程的原理需要光学,电磁学相关。咱忽略。

Schlick给出了一种近似算法。

8.2.2 specular reflection

8.2.3 specular transmission

镜面的btdf


8.3 LAMBERTIAN REFLECTION

Lambertian brdf是最简的一种反射模型。
将入射光完美地均匀反射到所有方向。
是一种理想化的模型。

8.4 MICROFACET MODELS 微小面模型

可以把物体表面看作一堆微小面。
面建模为高度场,用统计的方法描述。
图8.13
微小面模型中三种需要注意的问题。
  1. 点被山峰遮挡

  2. 点在阴影里

  3. 光在山谷里多次反射。

微小面模型主要组件
1. 微小面的整体分布
2. 单个微小面的brdf
单个微小面的brdf通常使用完美镜面反射。
对于一些半透明材质,会使用specular transmission。
Oren–Nayar模型对微小面使用Lambertian反射。
图8.13的问题需要重点考虑。
希望用较小的代价得到较好的近似值。

8.4.1 Oren–Nayar 漫反射

Oren和Nayar观察到完美的Lambertian反射并不真实。
Lambertian反射是把光均匀发射到所有方向。过于理想化。
而真实情况是靠近入射光方向会反射的更多。越以靠近入射光的角度去看物体,越亮。

V型微小面按高斯球形分布。参数σ,表示粗糙程度,微小面倾角的标准差。

对于V型微小面这种模型,多次反射的问题只需在相邻V型之间考虑。

同时考虑三个问题时,没法得到一个完备的表达式。
他们找到一个近似。

给定σ,θi和θo就可以算出反射率了。

8.4.2 微小面分布函数

基于微小面的模型,对于金属/塑料/毛玻璃等材质的镜面反射和传播有较好的表现。

微分面的分布函数D(ωh)。类似一个概率密度。

分布函数必须normalize,才能保证物理上正确。

图8.15。分布函数必须满足条件。给定dA,函数生成的投影面积也必须为dA。

积分形式

\({\LARGE \int_{H^2(n)}D(\omega_h)cos\theta_hd\omega_h=1}\)

\({\Large n}\)是宏观表面的normal。
\({\Large\omega_h}\)是微小面的normal。
\({\Large \theta_{h}}\)\({\Large \omega_{h}}\)\({\Large n}\)之间的夹角。
\({\Large cos\theta_h}\)是取\({\Large\omega_h}\)在宏观normal上的有效分量。
在半球上对D(ωh)d(ωh)在normal方向的分量积分,结果必须为1。
即半球上每个可能的\({\Large D(\omega_h)cos\theta_h}\)值加起来为1。
D函数说直白点,有一块表面,上面是无数杂乱无章的V形微小面。
某一个微小面的normal可以是任意值。
D函数可以看作是描述了每个可能的normal值的概率。
分布函数怎么来的,靠科学家发明。
必须积分为1。且在渲染上具有良好的性质。
两种常见的微小面分布
1. BeckmannDistribution
2. \({\LARGE D(\omega_h)=\frac{e^{-tan^2\theta_h/\alpha^2}}{\pi\alpha^2cos^4\theta_h}}\)
可以引入各项异性。之前只是竖直方向的\(\theta\)影响分布,现在让水平的转动\(\phi\)也影响分布。
可以做出更多效果。
  1. TrowbridgeReitzDistribution,即GGX。 GGX更符合实际物体的物理特性。

图8.16展示了两者的图形。
\(\alpha\)是可调的参数,用来生成不同效果的材质。
\(\alpha=0.5\)时,函数大致是钟形。
这时\(\theta\)向0集中,即光打到这个表面,得到的normal大概率接近宏观的normal。
也就是整体较为平整,到45度时概率就很小了。
如果是强烈的锯齿表面,两边的分布应该会更多。
GGX甚至超过90度还有分布。这微小面都倒扣了!
代码流程
积分器中Sample_f采样方向,走微小面的Sample_f
Sample_wh()根据\({\Large \omega_{o}}\)随机采样出一个normal\({\Large \omega_{h}}\)
每个微小面按镜面反射处理,那么可得到\({\Large \omega_{i}}\)
有了\({\Large \omega_{h}}\)可以算出D。有了D值可以得到pdf。

微小面采样见14.1.1


8.4.3 MASKING AND SHADOWING

只实现一个微小面分布D(ωh)已经可以出效果了。
但是还不够的。还有一个重点是微小面存在背对光线和被遮挡的情况。
Smith’s masking-shadowing function G1(ω, ωh)
给出了给定ωh后,某个ω方向可见面积的占比。
见图8.17
我们做微小面的一个预设是,在宏观上某些数理性质要和光滑的平面一致。
比如以\(\theta\)角度去看光滑的平面,眼睛看到的面积是\(dAcos\theta\)
那么以微小面来计算时,也得符合这个值。否则物理上就错了。
然后,我们当让可以把G1定为一个常量去迎合\(dAcos\theta\)。即从所有角度去看,被遮挡的占比相等。
但这肯定不符合现实,效果很差。
所以需要有针对角度的函数G1。所有需要有人发明出各种函数来达到真实的效果。
之前的微小面分布函数也是这个道理。先定一个前提,再凑一个模拟函数。

规则的积分形式

\({\LARGE cos\theta=\int_{H^2(n)}G_1(\omega,\omega_h)max(0,\omega \cdot \omega_h)D(\omega_h)d\omega_h}\)

给定入射方向\(\omega\)
有无数种可能的\(\omega_h\),每种\(\omega_h\)的概率为\(D(\omega_h)\)
\(\omega_h\)状态下能被看到的概率为\(G_1(\omega,\omega_h)\)
再取\(\omega_h\)方向的有效分量\(dot(\omega,\omega_h)\)
这三者乘起来再积分为\(cos\theta\)

这里有一层含义,G函数是之于D函数基础上的。是要配合的,不是独立的。

\({\large A^+(\omega)}\)为面朝入射光的面积,\({\large A^-(\omega)}\)为微小面投影到后面一个微小面上的面积。
那么被遮挡的面积就是两者相减。
那么可见的占比可以写成

\({\LARGE G_1(\omega)=\frac{A^+(\omega)-A^-(\omega)}{ A^+(\omega)}}\)

\(G_1(\omega,\omega_h)\)的区别是此时\(\omega_h\)隐含在了面积里。

再引出一个\({\LARGE \Lambda(\omega) =\frac{ A^-(\omega)}{A^+(\omega)-A^-(\omega)} }\)
看不见的面积与看得见的面积之比。

有关系\({\LARGE G_1(\omega)=\frac{1}{1+\Lambda(\omega)}}\)


`A Reflectance Graphics Model for Computer <https://graphics.pixar.com/library/ReflectanceModel/paper.pdf>`__

这个论文是82年的Cook-Torrance模型。

论文研究 `Microfacet Models for Refraction through Rough Surfaces <https://www.cs.cornell.edu/~srm/publications/EGSR07-btdf.pdf>`__

这个论文主要讲折射的微表面模型,引入了GGX。

宏观表面的BSDF记作\({\Large f_s(i, o, n)}\),描述入射i,出射o,normal n条件下得到的能量比例。

如果细分成反射和折射,那么就是BRDF(\(f_r\))和BTDF(\(f_t\))。

\(f_s\)\(f_r\)\(f_t\)的和。

可以不纠结每个微小表面,而是整体来看,用统计的方法(\(D_m\))来描述所有微表面的normal的分布。

之前学习过,发明出来的\(D_m\)必须保证物理上的正确性。

  • \({\Large 0 \le D(m)}\)
    任意m的D值大于等于0
  • \({\Large 1 \le \int D(m)d\omega_m}\)
    \(\omega_m\)为某个m对应的立体角。
    即所有方向的D(m)加起来大于等于1。
  • \({\Large (v \cdot n) = \int D(m)(v\cdot m)d\omega_m}\)
    给定任意方向v,考虑v在normal上的分量。
    左边的宏观值与右边的微观分量和要相等。
    这是一个通用的式子,对任意v有效。
    很多资料给出的是\({\Large \int D(m)(m \cdot n)d\omega_m =1}\)
    这其实是\({\Large v=n}\)的特殊情况。

Shadowing-Masking函数\({\Large G(i, o, m)}\)描述normal为m的微小表面,入射i,出射o情况下能被看到的概率。

  • \({\Large 0 \le G(i, o, m) \le 1}\)
    是一个普通的概率值
  • \({\Large G(i, o, m) = G(o, i, m)}\)
    来去对称
  • 如果\({\Large (i \cdot m) \le 0}\)\({\Large (o \cdot m) \le 0}\)\({\Large G(i, o, m) = 0}\)
    看微表面的背面,G一定为0。
    论文里带上了\({\Large (i \cdot n)}\)\({\Large (o \cdot n)}\),意思是扩展为正反面都成立。更严谨。
微表面的bsdf记作\({\Large f_s^m}\)
理论上微表面的bsdf可以用任何模型。但一般采用完美镜面模型。
论文就采用完美反射和完美折射

论文研究

`Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs <https://www.jcgt.org/published/0003/02/03/paper.pdf>`__

。。。


9 material

brdf和btdf描述了光在点上如何传播。
但渲染器需要知道点上具体包含哪些bxdf。

引入BSDF类。包含一系列bxdf。

9.1 BSDF

brdf和btdf的集合。

BSDF实例在SurfaceInteraction里

9.2 material

每种具体的material实现一个ComputeScatteringFunctions,
根据材质的性质和交点的情况,往SurfaceInteraction的BSDF里加相应的BxDF。

10 texture

pbrt里的texture比较宽泛。
本质上是一个映射。比如从一个表面的uv参数或空间的xyz,映射到其他domain比如一个颜色或一个实数。
texture可能是一个高频波动源,造成严重的锯齿。
可用第七章的方法处理。而更好的解决方法是用texture函数调整它的频谱信息,以符合采样率。
对于许多texture函数,计算一个合理的采样率,然后进行反走样,不是很难,但非常高效。

10.1 采样和抗锯齿

为了在texture函数中实现抗锯齿,需要解决两个问题。

  1. 必须计算texture空间中的采样率。
    屏幕空间的采样率是已知的,根据图像大小和像素的采样率。
    但必须确定场景中某个面上的采样率。
  2. 得到texture采样率以后,必须依据采样原理来计算texture值,
    使其不高于采样率能承载的最大频率波动。
    (除去超出奈奎斯特极限的频率)

10.1.1 找出texture的采样率

预先总结,这一节实际讲如何计算屏幕空间的xy坐标与texture的uv坐标的相对变化。

设某个面上的texture函数T(p),p为面上的点。即p点的texture值。
不考虑遮挡问题。
设f(x, y)为从屏幕空间的(x, y)到p点的转换。
可以把T(p)写成T(f(x, y))。即屏幕上(x, y)处的texture值(颜色)。
以图10.2为例。正交投影。
屏幕空间分辨率为xr和yr。
texture坐标为s,t。都在[0-1]范围。
那么屏幕坐标(x, y)对应的texture坐标(s, t)为(x/xr, y/yr)。
对于复杂的场景/相机模式/texture mapping,要找出屏幕坐标和texture坐标的精确对应是有难度的。
幸运的是,对于texture抗锯齿,我们不需要计算每个精确的f(x, y)。
只需要知道像素采样点的变化,和对应的texture坐标的变化。
方法是通过偏微分近似。
这是一个standard linear approximation
具体可参考[THOMAS' CALCULUS]13章偏微分。

\({\Large f(x, y) ≈ f(x_0, y_0) + (x - x_0)\frac{∂f}{∂x}(x_0, y_0) + (y - y_0)\frac{∂f}{∂y}(x_0, y_0)}\)

算出初始点(x0, y0)的值和对应的偏微分,就可以估算其他点(x, y)的值。

代码见
SurfaceInteraction::ComputeDifferentials()

从屏幕坐标(x,y)映射到texture(可看成一张图)的(u,v)

需要先在交点p处找出切面。
已知normal为n,已知点p。切面上的任意点(x,y,z)都应该满足,(x,y,z)到p的向量与n垂直。
所以有

\({\Large (x-p_x, y-p_y, z-p_z) \cdot n = 0}\)

\({\Large xn_x + yn_y+zn_z -n\cdot p = 0}\)

这样得到了切面方程。


PerspectiveCamera::GenerateRayDifferential
目前只考虑lensRadius=0。
给原始ray的屏幕坐标x,y分别加1,生成两条辅助光。原点一样。
这样构建出距离1,后续可以方便地反向算出偏微分值。
见SurfaceInteraction::ComputeDifferentials
两条辅助光rx和ry都打到切面,交点为px和py。
计算交点。也就是

\({\large p=o+t(dir)}\)

的t

有光与面的求交公式

\({\LARGE t=\frac{-n\cdot o - d}{n\cdot dir}}\)

其中d是之前求得的切面方程的d项也就是\(-n\cdot p\)
见代码中的tx,ty。
这里代码与书的公式对不上,是搞错了,得按公式改正代码。见https://pbrt.org/errata-3ed

tx和ty有了,px和py也有了。


p看作一个函数,从屏幕空间的xy映射到世界空间的xyz。
可用上述的线性近似

对于rx有

\({\Large p_x\approx p(x_0, y_0) + (x-x_0)\frac{∂p}{∂x}(x_0, y_0) + (y - y_0)\frac{∂p}{∂y}(x_0, y_0)}\)

对于rx,y的变化为0,那么

\({\Large p_x\approx p(x_0, y_0) + (x-x_0)\frac{∂p}{∂x}(x_0, y_0)}\)

由于辅助光的创建方式,\(x-x_0=1\)。有

\({\Large p_x\approx p(x_0, y_0) + \frac{∂p}{∂x}(x_0, y_0)}\)

\({\Large \frac{∂p}{∂x}(x_0, y_0) \approx p_x - p(x_0, y_0) }\)

ry同理。
即代码中的
dpdx = px - p;
dpdy = py - p;
按我理解这个推导跳过了一些步骤,不是表面上这么直接,它最后为什么直接就是两个点/向量相减了?
我只能按xyz三个分量来理解。
暂时只能这样,好理解一点。从后续内容看应该是对的。
p函数,从屏幕空间映射到三维世界空间,其实有三部分,分别得到xyz。
\({\large p_1(x, y) = newX}\)
\({\large p_2(x, y) = newY}\)
\({\large p_3(x, y) = newZ}\)
实际有三个计算,每个都符合上面的近似公式,最后再合起来。
看上去是在操作vector3,本质是单独计算再联写?

dpdxdpdy的本质是场景内三维坐标,对于屏幕空间二维坐标的变化率。


另一方面
dpdudpdv在之前的章节有详细的推导。在这里被用上了。
本质是场景内三维坐标,对于uv的变化率。
理论上有了p/x,p/u的关系,就可以算出u/x的关系。即uv对于屏幕xy的变化率。
记得微积分里有这个内容。

这里还是运用近似,这次xy同时移动。

\({\Large p'\approx p(x_0, y_0) + \Delta x\frac{∂p}{∂x} + \Delta y\frac{∂p}{∂y} \approx p(x_0, y_0) + \Delta u\frac{∂p}{∂u} + \Delta v\frac{∂p}{∂v} }\)

\({\LARGE \frac{\partial p}{\partial x}}\)\({\LARGE \frac{\partial p}{\partial y}}\)上面已算出。
屏幕空间\(\Delta x\)\(\Delta y\)都为1。
这时算出\(\Delta u\)\(\Delta v\),就得到了\({\LARGE \frac{\partial u}{\partial x}}\)\({\LARGE \frac{\partial v}{\partial x}}\)

解方程

\({\Large p'_x -p_x = \Delta u\frac{∂p_x}{∂u} + \Delta v\frac{∂p_x}{∂v}}\)

\({\Large p'_y -p_y = \Delta u\frac{∂p_y}{∂u} + \Delta v\frac{∂p_y}{∂v}}\)

\({\Large p'_z -p_z = \Delta u\frac{∂p_z}{∂u} + \Delta v\frac{∂p_z}{∂v}}\)

三个方程两个未知数,其中一个方程可能导致退化,解不出。
解决方法是根据n中xyz值来选取合适的方程来算。

10.1.2 滤波

可以用box filter进行AA。目前我代码抄过来有问题。


10.2 texture坐标的生成

常见的各种mapping方法

多数情况没有一个通用的完美的方法来映射texture。比如球体靠近极点的地方一定会有失真。

10.2.1 2d(uv)mapping

10.3 texture接口

10.4 image texture

是最常用的texture

10.4.1 texture存储管理

每个image都需要memory。复杂场景需要非常多的内存。
pbrt做了复用,同个image只需load一次。

10.4.3 mip maps

mipmapping是最流行的反走样方法。

频率如何理解?
目前的理解
比如有一张texture image,只考虑一边为例,长为x像素,那这个image的频率可以说是x。
然后场景里有一个面,需要贴这个texture,面的长为y。也就是需要在image上采样y次,就说采样频率为y。
走样的原因就是因为x比y大,或者说x分辨率更高。相差越大问题越大。
比如场景里很远的物体,屏幕上只占很小的面积,就100像素,却要把100万像素的图贴上去。
如果简单按比例去贴,就会出问题。
隔10000个像素抽一个像素,再拼起来,对于简单/线性的图片,貌似问题不大。
但是对于绝大多数图片,结果就是破碎的。
想个例子,把人体采样,只采10个样本,可能就会跳过眼睛之类,造成直接眼睛没了,脖子没了。。。
结果就是不成人形。
所以做texture sampling之前先要做所谓的过滤,把图像频率和采样频率统一。
必须先给我把100万像素的图处理到大致100像素,且看上去是对的,无走样。
我再去以100的频率去采样,就不会出大问题。
比如先得把100万像素的人体做成100像素的人体,而且看上去是正常的人体。
只是分辨率低一点,不能缺胳膊少腿。
那么怎么把100万像素处理到100像素而不缺胳膊少腿。
需要filtering算法
Nyquist limit/Nyquist frequency
两种filtering算法
1.trilinear interpolation
简单快速
2.elliptically weighted averaging/EWA
慢且复杂但是效果极佳

图10.10两种算法对比。EWA细节上更好。trilinear interpolation在边缘处细节丢失。

mipmap会把原始图片做成image pyramid。
如果原始图片尺寸不是2的次方,需要先resize成2的下一个次方的尺寸。

resize涉及采样原理。

因为图片尺寸只可能变大,所以不会引入走样。走样的本质起因是想把图片缩小。

详细流程见MIPMap::resampleWeights()

从一个原始图像生成mipmap的流程

    如果原始图像长宽不是2的次方。需要放大到下一级2的次方。长宽不一定要相等。

    先有重采样weight函数。
    MIPMap::resampleWeights(int oldRes, int newRes)
      //对于一串数据,从oldRes长度放大到newRes长度。

      i从0到newRes循环
        center = (i + .5f) * oldRes / newRes; //算出对应原始index

      以center为中心取4个index做Lanczos操作
            记录每个i的4个index和对应的Lanczos数据,作为weight。weight要做normalize。


    对于2d图像。横竖两个方向都要做一次重采样。

    对于每一行
        i从0到newRes循环
            用resampleWeights算出的4个index和weight,结合原始图像的值,算出新值。

    此时得到图a

    对于每一列
        i从0到newRes循环
            一样的方法,不过这次是weight结合图a的值,算出最终值。

    此时完成了放大图像到2的次方长宽。

    最后生成pyramid。下一级图像每4个值一组,就像一个个田字,每个田算平均值,作为一个点,形成了上一级的图像。
    这样不断缩小一半。就做好了mipmap。
trilinear interpolation/三线性插值

11 volumn scattering

通常假设光在真空中传播,也就意味着传输过程中radiance不变。
但是真实世界中比如烟,雾都会造成光的衰减和散射。

这一章讲光如何被participating media(空间中的大量粒子)影响。

volumn scattering模型的基础:大量的粒子形成一个统计学框架,而不是去死算每一个粒子。

用participating media可以模拟各种烟、云、激光、次表面等。


11.1 VOLUME SCATTERING PROCESSES

影响radiance分布的3个主要过程

  • 吸收:光能转换成其他能量比如热能。radiance减少。

  • 发光:发光粒子导致radiance增加。

  • 散射:radiance与粒子碰撞改变方向。

11.1.1 吸收

考虑火产生的浓烟,为什么看不到烟后面,因为其中的粒子吸收了光。
烟越浓光被吸收越多。同理烟能产生阴影。

吸收率写作\(\sigma_a\),是一个概率密度,表示光每传播一单位距离时被吸收的量。

有一些radiance到达点p

\({\Large Li(p, -\omega)}\)

经过一个微小距离\(dt\),出射能量变为

\({\Large Lo(p, \omega)}\)

吸收率的作用体现为

\({\Large Lo(p, \omega) - Li(p, -\omega) }\)

\({\Large = dLo(p, \omega) = -\sigma_a Li(p, -\omega)dt}\)

这是一个微分方程\({\Large y' = -\sigma_ay}\)

\({\Large y = y_0 e^{-\sigma t} }\)


12 光源

12.1 发光

绝对零度以上的物体有移动的原子。
根据麦克斯韦方程。带电的原子粒子的移动,会让物体发出一定波长的电磁辐射。
其中多数在室温下时为红外频率。想要得到可见频率,需要更高的温度。
把电能转化成电磁辐射有很多方法,也就是不同的光源类型。
了解一些其中物理原理有助于对光源建模。

Incandescent白炽灯 即 tungsten钨丝灯

有一根灯丝。电流通过,温度上升,发出电磁辐射。
用一个毛玻璃灯泡罩住,是为了得到合适的spd(光谱功率分布)。
白炽灯大部分电磁辐射都在红外范围。意味着大部分能量实际转为了热能而非光能。

Halogen卤素灯

也有一根灯丝。但灯泡中填充卤素气体。
随时间推移,白炽灯的灯丝会蒸发。而卤素灯的卤素气体会回到灯丝上,延长了寿命。
卤素气体也不会吸附到灯泡上,导致灯泡变黑。而白炽灯就会变黑。

Gas-discharge lamps 气体放电灯

把电流通进氢气/氖气/氩气/其他金属气体,产生光能。
通常在灯泡内测涂荧光层来产生各种颜色(波长)。

Led

electroluminescence

对上述光源,物理本质是电子和原子的碰撞,把外层的电子推向更高的能量级别。
当这些电子回到低能量时,会发出光子。
也有其他的原理,比如化学发光(荧光棒),生物发光(萤火虫)。不做研究。

12.1.1 黑体辐射源

理想的黑体会吸收所有照射到它表面的电磁辐射,无论频率和角度。完全不会反射能量。
黑体在热平衡状态会发出黑体辐射。辐射的值由普朗克定律给出,跟温度有关。
一般呈现黑色,温度很高时呈现出的其他颜色。
普朗克定律
给定温度,波长,可计算出辐射率。

12.2.1 检测光是否被阻挡


12.3 点光源

12.4 DISTANT LIGHTS/directional light

其实近似于无限远的点光源

12.5 AREA LIGHTS

绑定模型,从模型表面发光。

12.6 INFINITE AREA LIGHTS

看成一个巨大的球体包围场景,从所有方向向场景发光。
主要应用在环境光。

13 蒙特卡洛积分

我们用到的积分方程一般没有解析解。那么就得考虑数值解。

数值解
trapezoidal integration 和 Gaussian quadrature
梯形积分和高斯积分
对于低维度光滑的积分比较高效,但高维度不连续的积分收敛不理想。
蒙特卡洛积分是一种可行方法
引入随机性,使收敛速度与维度无关。

正确使用随机,在算法设计领域是革命性的。

随机算法分两大类

  1. 拉斯维加斯
    每次的输出都是一致的
    比如快速排序。每次选pivot可以是随机的,但是最终结果是精确的,只要输入相同,结果一定相同。
  2. 蒙特卡洛
    每次的输出不一致,是随机性的,但是平均结果是基本正确的。
    可理解为每次的结果基本正确,有时好一点有时坏一点。
    如果多算几次,取平均值,基本能保证较好的结果。
蒙特卡洛积分使用随机采样,只需算随机x的f(x),就能得出f(x)的积分。
代码容易实现,而且对不连续的函数同样有效。
比如算一个点上的辐射率的积分,按方程的话要考虑所有立体角。诸如此类很难算或者根本没法算。
蒙特卡洛方法直接选择一组方向进行计算。
蒙特卡洛的主要缺点
用n个样本来估计积分,收敛率,\(O(n^-0.5)\)
意思是如果要减少当前一半的错误,需要做4倍于当前的采样。
对于生成的图像来说,蒙特卡洛的错误就是noise,随机出现像素不自然地过亮或者过暗。
这一技术的研究方向基本是如何减少noise,同时尽量少增加采样数。

13.1 概率

两本推荐的概率书
[Introduction to Probability Models]
[A FIRST COURSE IN PROBABILITY]
随机变量X
给定一个概率值\(\xi\)[0,1],可以map到某个X。
cumulative distribution function(CDF)
写作\(P(x)\)

\(P(x) = Pr(X\le x)\)

即所有小于等于x的分布的和

probability density function
PDF所谓的概率密度函数,写作\(p(x)\)

13.1.1 连续随机变量

渲染领域,连续随机变量比离散随机变量更常见。

canonical uniform random variable
标准平均随机变量。拉丁字母\(\xi\)
平均分布在[0, 1)
  1. 容易生成

  2. 可先生成标准平均随机变量,再进行变换。

pdf积起来就是cdf

\({\Large P(x \in [a,b]) = \int_{a}^{b}p(x)dx}\)

cdf微分就是pdf

\({\Large p(x)=\frac {dP(x)}{dx}}\)

整个pdf下方的面积一定是1。

13.1.2 期望值和方差

函数的期望值定义

\({\Large E_p[f(x)] = \int_{D}^{} f(x)p(x)dx}\)

方差

\({\Large V[f(x)] = E[(f(x)-E[f(x)])^2]}\)


13.2 monte carlo estimator

假设要算这样的积分

\({\Large \int_{a}^{b} f(x)dx }\)

可以做一个所谓蒙特卡洛估算子。

\({\Large F_N=\frac{b-a}{N} \sum_{i=1}^{N}f(X_i)}\)

其中Xi是平均分布的随机变量。
这样把一个严格的积分转化为有限的N次计算。

看看期望值

\({\Large E[F_N] = E[\frac{b-a}{N} \sum_{i=1}^{N}f(X_i)] }\)

\({\Large = \frac{b-a}{N} \sum_{i=1}^{N}E[f(X_i)] }\)

\({\Large = \frac{b-a}{N} \sum_{i=1}^{N} \int_{a}^{b} f(x)p(x)dx }\)

因为Xi是平均分布,所以p(x)=1/(b-a)。

\({\Large = \frac{1}{N} \sum_{i=1}^{N} \int_{a}^{b} f(x)dx }\)

\({\Large = \int_{a}^{b} f(x)dx }\)

可以看到把积分问题转到概率框架下,得到的期望值是积分本身。


对于非均匀分布,可以做这样的估算子。

\({\Large F_N=\frac{1}{N} \sum_{i=1}^{N}\frac{f(X_i)}{p(X_i)}}\)

\({\Large E[F_N]=\frac{1}{N} \sum_{i=1}^{N} \int_{a}^{b}\frac{f(x)}{p(x)}p(x)dx}\)

\({\Large = \int_{a}^{b} f(x)dx }\)

同样的结果。这个就是蒙特卡洛积分的理论依据了。采样本,除pdf,做平均,期望值和积分相等。


pdf也就是p(x)到底从哪来?上面的推导其实说明了pdf可以是任意分布,随便造都可以!
只要按你造的这个pdf去采样,期望值就没问题。

可以自己手动模拟试试。给一个函数任意造pdf,再模拟多次采样,结果是一样的。

这个其实挺隐晦的。开始并没有仔细看这些理论,会认为每个场景下都会有一个”正确”的pdf。其实不对。
pdf是针对采样的,并不是某个brdf的性质或者某个函数的性质。

那么怎么选pdf。不同的pdf会造成不同的方差,效果可能非常大。

最简单的是平均采样。

如果pdf和原函数f形状相似,方差就小。见后面的重要性采样相关内容。


13.3 对随机变量采样

蒙特卡洛需要在概率分布中采样。
需要做一套1d/2d的采样方法。

13.3.1 INVERSION方法

看懂图13.1/13.2/13.3

本质是从一串离散的权重开始,例如3,2,1。
算出各自的概率,总和为1,同时算出cdf(概率累加)。
最后得到cdf=[0, 0.5, 5/6, 1]
以后就可以用一个[0, 1]的概率值来采样了。概率所在的cdf区间决定index。
代码见Distribution1D
例子是离散形式,pdf(x)其实就是(1, 3/6), (2, 2/6), (3, 1/6)这几个对子。
x轴为数组index,y轴为概率。
对于连续的分布。
\(pdf(x)\)变成连续的函数

\({\Large cdf(x)=\int_{0}^{x} pdf(x')dx'}\)

积分区域是从\(-\infty\)到指定的某个x。x对应实际随机变量的范围。例如上面取了1,2,3,那么范围就是[1, 2, 3],。
如果小于0的区域没有分布,也可以从0开始。总之是从有分布的地方开始。
积分范围总之就是所有可能的取值集合。
这时可以输入一个概率值,反过来找X。
也就是给cdf的反函数,代入概率,求得最后的X。
这样可以输入均匀分布的随机数,得到对应的X值,也就是对X均匀采样了。

\({\large cdf(X)=\xi}\)

\({\large X=cdf^{-1}(\xi)}\)

所以叫INVERSION方法。

[Introduction to Probability Models]这本书11.2.1章的内容。
对一个pdf函数均匀采样的步骤
1. 积分得到cdf 2. 计算\(cdf^{-1}\)
3. 生成均匀分布的随机数U[0, 1]
4. \(X=cdf^{-1}(U)\)

其中积分域要注意是从有分布的地方开始到x,不一定是从0开始。


pow(x, n)分布的采样

一般做成\(pdf(x) = cx^n\)

用c把pdf(x)的积分值限定在1,积分上下限可以按需取值。

算出c。那么pdf(x)就有了。再按上述方法得出cdf反函数。


指数分布采样

体渲染会用到,对指数函数在0到无穷大上采样。

\({\Large f(x) = ce^{-ax}}\)

再看下INVERSION方法。
例子中是积分范围已知,pdf已知,积分值规定为1。
此时要在0到无穷大上采样,所以积分范围就是0到无穷大。
cdf同样可以规定为1,这样采样的时候取0-1的随机数就行。
但pdf是未知的有个常数c,其实这里有一个隐含的操作,可以说是normalize。
cdf为什么可以规定为1,就是要引入系数,做normalize。让最终的pdf积分为1,以方便取随机数。
刚好问题里就有一个系数,那么可以直接算。
如果没有系数,pdf是明确的,那么就要加个系数并算出来。
否则随机数就得按cdf值来取。
(之前例子的cdf=1是造出来的特殊1,pdf是凑好的概率,所以刚好不用算系数。而一般函数采样的话,积起来的cdf肯定是五花八门,所以要做normalize。)

\({\Large \int_{0}^{\infty} ce^{-ax}dx =1}\)

\({\Large -\frac{c}{a}e^{-ax}|_{0}^{\infty} =1}\)

\({\Large 0 - (-\frac{c}{a}) =1}\)

\({\Large c = a}\)

所以normalize后的pdf函数为

\({\Large pdf(x) = ae^{-ax}}\)

积分

\({\Large \int_{0}^{x} ae^{-ax'}dx' = -e^{-ax}|_{0}^{x} = 1 - e^{-ax}}\)

求逆

\({\Large y = 1 - e^{-ax}}\)

\({\Large e^{-ax} = 1-y}\)

\({\Large y^{-1} = -\frac{ln(1-x)}{a}}\)

采样

\({\Large X = -\frac{ln(1-\xi)}{a}}\)

Piecewise-Constant 1D

就是Distribution1D的实现

13.3.2 REJECTION METHOD

有些函数不好找pdf,或者不好求cdf。
可用rejection方法来采样。这是一种类似投飞镖的模式。
见图13.5
已知f(x),但不好按上述方法采样。
这时如果有一个好求的pdf p(x),且cp(x)>f(x)。
那么可以做一些近似。

不断地对p(x)采样,得到ξ和X,如果(X, ξcp(X))在f(X)图像下方,就认为对f(x)采样成功,取(X)。


13.5 在不同分布之间变换

给定pdf函数\(p_x(x)\),随机变量为\(X_i\)
如果又给一个y函数对x进行转换,\(Y_i=y(X_i)\)。那么\(Y_i\)的分布情况是什么?
这部分一开始没理解。找了些稍详细的教学。耐心推理一下,能更好地理解。
这类问题换种说法就是:把pdf的变量转到另一种变量空间(系统/坐标系),以求计算更容易或者算出原来无法计算的式子。
有随机变量X的pdf,\(f(x)\)
有函数\(r\)把X映射到\(Y\)\(Y=r(X)\)

先体验离散的例子

所有X的分布。和为1。

\({\large \begin{matrix} f(0)=0.1 \\ f(1)=0.2 \\ f(2)=0.4 \\ f(3)=0.1 \\ f(4)=0.2 \end{matrix}}\)

r把X转换成Y

\({\large \begin{matrix} r(0)=2 \\ r(1)=1 \\ r(2)=2 \\ r(3)=1 \\ r(4)=5 \end{matrix}}\)

Y值的分布。和为1。

\({\large \begin{aligned} & g(1)=f(1) + f(3) = 0.3 \\ & g(2)=f(0) + f(2) = 0.5 \\ & g(5)=f(4)=0.2 \\ \end{aligned}}\)

把Y转回X是\({\Large r^{-1}(Y)}\)

可总结出公式

\({\Large g(y)=\sum_{x\in r^{-1}(y)}f(x)}\)

\({\large x\in r^{-1}(y)}\)意思是多个x可能映射到同个y,那么从y返回去就是一个x集合。所以是一个和的形式。

对于连续的f。就是积分了。

\({\Large g(y)=\int_{r^{-1}(y)}f(x)dx}\)

把所有y对应的x的f(x)积起来。

\({\Large g(y)}\)即是Y的pdf,所有\({\Large g(y)}\)值的和为1。

pdf是针对某个y值,如果把y值扩展为从\(-\infty\)开始,就变成cdf。

\({\Large G(y)=\int_{r^{-1}(-\infty,y)}f(x)dx}\)

之前的例子r函数是随意的。那么从y反出来的x也是一个随意的集合。
那么积分范围是乱的,可能重叠,可能杂乱分布,很难受。

这时如果规定r函数xy一一对应并且单调递增,那么r的反函数可以形成一个正常的积分范围。

且有
>\({\Large G(y)=F(r^{-1}(y))}\)
就是随便取一个Y,对应X,两个cdf相等。
比如

\({\large \begin{matrix} r(0)=2 \\ r(1)=3 \\ r(2)=4 \\ r(3)=7 \\ r(4)=18 \end{matrix}}\)

可以体会到,一一对应并且单调递增的效果就是把映射做成有序并严格隔离。这样两个cdf就能相等。

\({\Large G(y)=F(r^{-1}(y))}\)

两边微分cdf得pdf

\({\LARGE g(y)=f(r^{-1}(y))\frac{d}{dy}(r^{-1}(y))}\)

如果r是单调递减,有

\({\Large G(y)=1-F(r^{-1}(y))}\)

\({\LARGE g(y)=-f(r^{-1}(y))\frac{d}{dy}(r^{-1}(y))}\)

如果r递减,有

\({\Large \frac{d}{dy}(r^{-1}(y)) \lt 0 }\)

可以把前面的负号拿到后面。递增递减的两个\({\Large g(y)}\)可以统一写成

\({\LARGE g(y)=f(r^{-1}(y))|\frac{d}{dy}(r^{-1}(y))|}\)

\({\Large r^{-1}(y) = x}\)。最终

\({\LARGE g(y)=f(x)|\frac{dx}{dy}|}\)

这样就可以求对变量做变换以后的分布情况。
比如已知pdf函数g(y),可以算\({\large x=y^2+y}\)的pdf。
对y求导,两边一除就得到了f(x)。

13.5.1 多维的变换

涉及多重积分的换元。见[THOMAS' CALCULUS]14.8。

有的积分难算,可以尝试替换变量,实际就是换一个坐标系,变得容易算。

单变量的替换:

\({\LARGE \int_{g(a)}^{g(b)} f(x)dx=\int_{a}^{b}f(g(u))g'(u)du }\)

其实是多变量的特殊形式。

对于两个变量
\({\Large f(x, y)}\)\({\Large R}\)上积分,想要转成\({\Large f(g(u, v), h(u, v))}\)\({\Large G}\)上对\({\Large u,v}\)积分。
也需要一个类似\({\Large g'(u)}\)的因子\(J\),即Jacobian

\({\LARGE J(u, v) = \begin{vmatrix} \frac{\partial x}{\partial u} & \frac{\partial x}{\partial v} \\ \frac{\partial y}{\partial u} & \frac{\partial x}{\partial v} \end{vmatrix}}\)

最终的转换

\({\LARGE \iint_{R} f(x,y)dxdy=\iint_{G}f(g(u,v), h(u, v))J(u, v)dudv }\)

更多维度也一样适用。
理解[THOMAS' CALCULUS]里的例子。

上一节的 \({\LARGE g(y)=f(x)|\frac{dx}{dy}|}\) 也是Jacobian的一维转换特例。
多维也一样适用

\({\LARGE g(y)=f(x)|J_x(y)|}\)

比如三维数据
>\({\LARGE J_x(y) = \begin{vmatrix} \frac{\partial x_1}{\partial y_1} & \frac{\partial x_1}{\partial y_2} & \frac{\partial x_1}{\partial y_3}\\ \frac{\partial x_2}{\partial y_1} & \frac{\partial x_2}{\partial y_2} & \frac{\partial x_2}{\partial y_3}\\ \frac{\partial x_3}{\partial y_1} & \frac{\partial x_3}{\partial y_2} & \frac{\partial x_3}{\partial y_3} \end{vmatrix}}\)

13.5.2 极坐标的转换

\({\Large \left \{\begin {aligned} & x=r cos\phi \\& y=rsin\theta \end{aligned}\right.}\)

\({\Large J(r, \theta) = \begin{vmatrix} cos\theta & -rsin\theta \\ sin\theta & rcos\theta \end{vmatrix}}\)

\({\Large =r(cos^2\theta + sin^2\theta )}\)

\({\Large =r}\)

所以

\({\Large p(r, \theta)=r p(x, y)}\)

13.5.3 球面坐标的转换

\({\Large \left \{\begin {aligned} & x=rsin\theta cos\phi \\& y=rsin\theta sin\phi \\& z=rcos\theta \end{aligned}\right.}\)

\({\Large J(r, \theta, \phi) = \begin{vmatrix} sin\theta cos\phi & rcos\theta cos\phi & -rsin\theta sin\phi\\ sin\theta sin\phi & rcos\theta sin\phi & rsin\theta cos\phi\\ cos\theta & -rsin\theta & 0 \end{vmatrix}}\)

\({\Large =r^2(sin\theta cos^2\theta cos^2\phi + sin^3\theta sin^2\phi + sin\theta cos^2\theta sin^2\phi + sin^3\theta cos^2\phi)}\)

\({\Large =r^2(sin\theta cos^2\theta + sin^3\theta )}\)

\({\Large =r^2sin\theta}\)

所以笛卡尔坐标系xyz上某个区域A的积分可转成球面坐标系对应的B区域的积分

\({\Large \iiint_AF(x,y,z)dxdydz=\iiint_BG(r, \theta,\phi)(r^2sin\theta)drd\theta d\phi}\)

对于pdf有

\({\Large p(r,\theta, \phi) = p(x,y,z)r^2sin\theta}\)


另外
在 [[#5 5 2 球坐标积分]] 中推导过

\({\Large dw= sin\theta d\theta d\phi}\)

那么对于半球上的积分

\({\Large \int p(w)dw=\int p(\theta,\phi)d\theta d\phi = 1}\)

\({\Large dw}\)代入得

\({\Large \int sin\theta p(w) d\theta d\phi = \int p(\theta,\phi)d\theta d\phi}\)

所以

\({\Large sin\theta p(w)=p(\theta,\phi)}\)

这样也能从对于立体角的pdf转成对于极坐标的pdf。


13.6 多维变换的2d采样

设2d联合密度函数\({\Large p(x, y)}\),一般xy是相关的,是不能拆分的。
对于可以拆分的函数,比如\({\Large p(x, y)=p_x(x)p_y(y)}\)
那么分别采样x和y即可。

对于不可拆分的密度函数,要从\({\Large p(x, y)}\)里求marginal density function

\({\Large p(x)=\int p(x, y)dy}\)

看做只针对x的密度函数。实际这是每个x,在所有y上的平均密度。

概率上有

\({\LARGE p(y|x)=\frac{p(x, y)}{p(x)} }\)

所谓的conditional density function

这个非常厉害,能把整个采样分出先后。
先已知\(p(x, y)。\)marginal density function\(p(x)\),对\(p(x)\)采样出X。
X代入\({\Large \frac{p(x, y)}{p(x)}}\)得到只有y的函数,再采样出Y。

这样就把2d采样转成了两个分开的1d采样。

13.6.1 对半球均匀采样

半球范围对立体角采样。
对于立体角来说,pdf是常数。有\(p(w) = C\)。我们要求这个C。

那么因为pdf积分为1

\({\Large \int_{H^2} p(w)dw = 1}\)

\({\Large C\int_{H^2} dw = 1}\)

半球范围的立体角和为\(2\pi\)

\({\Large C2\pi = 1}\)

所以

\({\Large p(w)=\frac{1}{2\pi}}\)

[[#13 5 3 球面坐标的转换]] 里有

\({\Large sin\theta p(w)=p(\theta,\phi)}\)

所以

\({\Large p(\theta,\phi)=\frac{sin\theta}{2\pi}}\)

这个例子对\((\theta,\phi)\)采样是可分开的,但是用上一节的marginal density function也是可以求的。

\(\phi\)积分

\({\Large p(\theta) = \int_{0}^{2\pi} \frac{sin\theta}{2\pi} d\phi=sin\theta}\)

\({\LARGE p(\phi|\theta)=\frac{p(\theta,\phi)}{p(\theta)} = \frac{1}{2\pi}}\)

有了pdf,按[[#13 3 1 INVERSION方法]]采样。

\({\Large P(\theta)=\int_{0}^{\theta}sin\theta'd\theta' = -cos\theta|_{0}^{\theta} = 1-cos\theta}\) (积分得到cdf)

\({\Large 1-cos\theta = \xi_1}\) (给定随机值反算变量)

\({\Large \theta = acos(1- \xi_1)}\) (完成采样)

\({\Large P(\phi|\theta) = \int_{0}^{\phi} 1/2\pi d\phi' =\phi/2\pi|_{0}^{\phi}=\phi/2\pi }\) (积分得到cdf)

\({\Large \phi/2\pi = \xi_2}\) (给定随机值反算变量)

\({\Large \phi = 2\pi\xi_2}\) (完成采样)

最终转到笛卡尔坐标系即球面上的xyz

\({\Large \begin {aligned} & x=rsin\theta cos\phi =r\sqrt {1-(1-\xi_1)^2} cos(2\pi\xi_2) \\& y=rsin\theta sin\phi =r\sqrt {1-(1-\xi_1)^2} sin(2\pi\xi_2) \\& z=rcos\theta=r(1-\xi_1) \end{aligned}}\)

到这里可以看到,对一个半球均匀采样,乍看很简单,实际包含大量推导。


13.6.2 对单位disk采样

对disk采样看上去简单,对半径采样再对角度采样,实际不对。见图13.10。

针对圆形均匀采样,从面积的角度来看,采样是均匀的,采到每个单位面积的概率为

\({\Large p(x,y)=\frac{1}{\pi R^2}}\)

是一个常数,和要采样的实际圆形半径相关。
这里pbrt取R=1,得到\({\Large 1/\pi }\),是做了normalize,对单位圆采样。

我呢保留R,要注意不能写成r,容易跟后面的变量r搞混。要把R看成常数。

根据13.5.2进行极坐标的转换

\({\Large p(r, \theta)=\frac{r}{\pi R^2}}\)

marginal density,对\(\theta\)积分

\({\Large p(r)=\int_{0}^{2\pi}\frac{r}{\pi R^2}d\theta = \frac{2r}{R^2}}\)

conditional densities

\({\Large p(\theta|r) = \frac{p(r, \theta)}{p(r)} = \frac{1}{2\pi} }\) (是一个常数)

积分得cdf并采样

\({\Large P(r)=\frac{r^2}{R^2} = \xi_1}\)

\({\Large r=R\sqrt{\xi_1}}\)

\({\Large P(\theta|r)=\frac{\theta}{2\pi} = \xi_2}\)

\({\Large \theta=2\pi\xi_2}\)


虽然实现了disk的均匀采样,但是上述方法对disk有distortion。13.8.3节讲了它的坏处。

一个较好的方法是concentric mapping


13.7 RUSSIAN ROULETTE

蒙特卡洛路径追踪算法中,需要在一条路径上迭代(弹射)。每次迭代可以看成有相等的开销,但是光能是不断衰减的,对整体光能的贡献主要集中在开头的数次迭代。
这样用同样的开销去做很多不重要的计算,非常不合算。

另外最完美的迭代次数是无穷大,这也不现实。

简单的处理办法是设一个阈值,光能小到一定程度就停止计算。或者指定最大弹射次数,比如10次强行终止。
一个问题是不好决定阈值和弹射次数,可能跟具体场景等强烈相关,没有通用性。
一个问题是这样造成结果是有偏的(biased),在数学上不准确。

Russian roulette可以设置终止条件,按概率去终止一些不想要的计算,同时还能保证结果无偏(unbiased)。

这里面是有数学理论的,并不是按字面意思随便转个轮盘,1/10概率停止这样。
一开始我没有详细看,出的图都是有问题的。

首先有期望值的定义

\({\Large E[X]=\sum_{i=1}^{n}x_ip(x_i)}\)

对于一串离散随机变量,每个值乘上其概率,再加起来就是期望值。
可以理解为一个“代表值”或者“平均值”。
蒙特卡洛路径追踪里,每一次弹射我们会算beta,可以看作能量的scale。
beta会越来越小,当beta很小时,比如小于0.1。我们就开始不太愿意费劲继续去算后续的贡献值。

这时可以做一个概率q,算一个\(beta'\)作为新的beta。

\({\Large beta' =\begin{cases} \frac{1}{q} beta & (概率q) \\ 0 & (概率1-q) \end{cases} }\)

q是完全可以自定义的,比如当beta小于0.2时,我决定此时得有0.1的概率不要再继续弹射了,也就是让\(beta'\)返回0。这时q取0.9。

然后我们可以写代码。当beta小于0.2时,取一个[0, 1]随机数,如果大于0.1,继续迭代。否则终止,也就是后续的贡献取0。


为什么是\({\Large \frac{1}{q} beta}\)?简单计算一下,原始正确值是beta。我们按概率来看,进行10次操作,会有9次返回\({\Large \frac{1}{q} beta}\),1次返回0。

那么10次操作

\({\Large 平均值 = \frac{\frac{beta}{0.9} \times 9 + 0\times 1}{10}=beta}\)

平均值还是beta并没有改变。在概率框架下操作次数越多就越准。
如果不除q,也就是强行终止,平均值就偏小,能量损失了1/10。

本质上就是利用概率,随机终止操作,同时保证beta的期望值不变,即无偏。

\({\Large E[beta'] = \frac{1}{q} beta \times q + 0 \times (1-q)}\)

\({\Large = beta = beta \times 1 = E[beta]}\)


13.10 重要性采样

重要性采样是减小方差的好工具

对于这个蒙特卡洛估算子

\({\Large F_N=\frac{1}{N} \sum_{i=1}^{N}\frac{f(X_i)}{p(X_i)}}\)

如果选的p(x)与f(x)图形越相似,方差就越小,收敛速度越快。

假设

\({\Large p(x) = cf(x)}\)

那么

\({\Large \int p(x)dx = \int cf(x)dx = 1}\)

\({\Large c = \frac{1}{ \int f(x)dx}}\)

这时估算子里

\({\Large \frac{f(X_i)}{p(X_i)} = 1/c = \int f(x)dx}\)

正好已经得到了答案,方差为0。

这个是完美的假设,完全不实用。但如果能找到近似的p(x),是能减小方差的。


14 光的传播1-表面反射

这一章结合光线追踪/辐射度量学/蒙特卡洛采样,实现DirectLightingIntegrator和PathIntegrator。
为什么叫积分器。因为它本质上负责计算辐射达到稳定状态时方程的积分。

直接光的积分

\({\Large L_o(p, w_o)=\int_{S^2} f(p, w_o, w_i)L_i(p, w_i)|cosθ_i|dw_i}\)

\(L_i\)是直接光从\(w_i\)打到p上的光能。
\(f(p, w_o, w_i)\)是材质的bxdf函数,给出光从\(w_i\)\(w_o\)出时带多少比例光能。
\(cosθi\)是取normal方向分量。
最后按半球积分得到世界上所有光打到p点,并通过\(w_o\)方向射到眼睛的光能。
这描述了完美的直接光。

但这个积分一般是没法正常算的。用蒙特卡洛方法能得出近似结果。

\({\LARGE L_o(p, w_o) ≈ \frac{1}{N} \sum \frac{f(p, w_o, w_i)L_i(p, w_i)|cosθ_i|}{pdf(w_i)}}\)

比如理想状态需要计算半球上所有方向,有无数条光。但现在我只选少数N个方向算光能。
每条光的能量除以光射到这个方向的概率,得到一个近似的总能量。
类似于选一些个体来代表全体,那么选的个体越多,结果越好。

14.1 对反射函数采样

BxDF::Sample_f()
默认对

BxDF::Pdf()

14.1.1 微小面采样

Sample_wh()是从\({\Large \omega_{o}}\)随机采样出一个normal\({\Large \omega_{h}}\)

\({\Large \theta_{h}}\)\({\Large \omega_{h}}\)和Z轴之间的夹角。

经典的方法是对\(D(\omega)\)进行全局采样。
还有一种更高效的方法是只对对可见面采样。

各向同性Beckmann为例

\({\LARGE D(\omega_h)=\frac{e^{-tan^2\theta_h/\alpha^2}}{\pi\alpha^2cos^4\theta_h}}\)

\({\LARGE \xcancel{ p(\theta_h)=\frac{e^{-tan^2\theta_h/\alpha^2}}{\pi\alpha^2cos^4\theta_h}}}\)

\(p(\theta_h)\)采样,需要求积分
用wolframalpha求出来没看懂。和pbrt也不一样
有erf函数。超出了知识范围。
有缘再见

到这里必须看懂13章里的各种转换和采样。之前没有仔细看13章,所以这里搞不懂。

由[[#8 4 2 微小面分布函数]] 规定

\({\LARGE \int_{H^2}D(\omega_h)cos\theta_hd\omega_h=1 }\)

\(D(w_h)\)并不是一个标准的pdf。pdf需要计算。

\({\LARGE \int_{H^2}D(\omega_h)cos\theta_hd\omega_h=1 = \int_{H^2}p(\omega_h)d\omega_h}\)

\({\LARGE D(\omega_h)cos\theta_h= p(\omega_h)}\)

\({\LARGE p(\omega_h) = \frac{e^{-tan^2\theta_h/\alpha^2}}{\pi\alpha^2cos^3\theta_h} }\)

这个才是正确的pdf。
然后进行坐标转换。

[[#13 5 3 球面坐标的转换]] 里有

\({\Large sin\theta p(w)=p(\theta,\phi)}\)

所以

\({\LARGE p(\theta,\phi) = \frac{sin\theta e^{-tan^2\theta_h/\alpha^2}}{\pi\alpha^2cos^3\theta_h} }\)

到此就和[[#13 6 1 对半球均匀采样]]差不多的套路了。
这里用分开采样,合起来采也可以,结果一样。

\({\Large p(\theta, \phi) = p(\theta)p(\phi)}\)

其中\({\Large p(\phi)=1/2\pi}\)

所以

\({\LARGE p(\theta) = \frac{2sin\theta e^{-tan^2\theta_h/\alpha^2}}{\alpha^2cos^3\theta_h} }\)

对其积分并采样

\({\Large -e^{-tan^2\theta/\alpha^2}|_{0}^{\theta} = 1-e^{-tan^2\theta/\alpha^2} =\xi}\)

\({\Large tan^2\theta=-\alpha^2ln(1-\xi)}\)

这里得到了和pbrt一样的结果。完成了对\(\theta\)的采样,实际应用中采出\(tan^2\theta\)已经够了。
通过三角函数公式可以算出其他sin/cos值,不用算出具体角度。

\(\phi\)是独立的圆周上的均匀分布,\({\Large p(\phi)=1/2\pi}\),所以

\({\Large \phi=2\pi\xi}\)

到此完成了各向同性Beckmann采样。


各向同性GGX

\({\LARGE D(\omega_h)=\frac{\alpha^2}{\pi cos^4\theta_h(\alpha^2+tan^2\theta_h)^2}}\)

\({\LARGE p(\omega_h)=\frac{\alpha^2}{\pi cos^3\theta_h(\alpha^2+tan^2\theta_h)^2}}\)

\({\LARGE p(\theta, \phi)=\frac{sin\theta \alpha^2}{\pi cos^3\theta_h(\alpha^2+tan^2\theta_h)^2}}\)

\({\LARGE p(\theta)=\frac{2sin\theta \alpha^2}{cos^3\theta_h(\alpha^2+tan^2\theta_h)^2}}\)

\({\Large \frac{\alpha^2}{(\alpha^2 - 1)((\alpha^2-1)cos^2(t)+1)}|_{0}^{\theta} = \frac{\alpha^2}{(\alpha^2 - 1)((\alpha^2-1)cos^2(\theta)+1)}-\frac{1}{\alpha^2 - 1} = \xi}\)

最终采样

\({\Large cos^2\theta = \frac{1-\xi}{\alpha^2 \xi-\xi+1}}\)
\({\Large \phi=2\pi\xi}\)

[[#8 4 3 MASKING AND SHADOWING]]里引入了G函数,效果更真实。
之前是只针对分布函数\(D\)算pdf,如果要引入G,就得联合G一起算pdf。套路是一样的。

14.1.2 FresnelBlend

14.1.3 镜面反射和传输

14.1.4 傅立叶bsdf

14.1.5 例子:估算反射率

计算8.1.1的半球反射
ρhd(Wo)=∫fr(p, wi, wo)|cosθi|dwi
见BxDF::rho()
不过这个好像在pbrt代码里并没有调用
但是在自己实现光子映射中会用到

14.1.6 BSDF::Sample_f()

给定wo,BxDFType。
根据BxDFType选出一个bxdf,进行采样bxdf->Sample_f()。得出光能/wi/pdf等。
随机选出的bxdf决定了方向。
如果选出的bxdf不是BSDF_SPECULAR,
对于所有非镜面反射的bxdf比如LambertianReflection和MicrofacetReflection,算f(wo, wi)。

14.2 对光源采样

Light::Sample_Li()

14.2.4 infinite area light

环境map,一般情况不同方向的光差别巨大。

波动越小越好。图14.10做对比。波动越大noise越多。


14.3 直接光

14.4 光的传输方程/渲染方程

light transport equation (LTE)
描述了场景中的光达到稳定状态时的信息。

14.4.1 基本推导

LTE的基本原则是能量的守恒。微小的能量改变都需要进行追踪。
p点向\(w_o\)方向发射的光能\(L_o(p, w_o)\)
我们的终极目标就是得到这个值。

先有直接光的积分

\({\Large L_o(p, w_o)=\int_{S^2} f(p, w_o, w_i)L_i(p, w_i)|cosθ_i|dw_i}\)

p点在\(w_o\)方向射出的光能,为所有光打到p点,进行brdf计算后,在\(w_i\)方向的总和。

如果考虑p点本身也发光,在\(w_o\)方向为\(L_e(p, wo)\)
那么就有了渲染方程。就是p点本身发光加上反射光。

\({\Large L_o(p, w_o)=L_e(p, w_o) + \int_{S^2} f(p, w_o, w_i)L_i(p, w_i)|cosθ_i|dw_i}\)

再看其中的\(L_i(p, w_i)\)是从哪来的?
一定是从某点p2发出的光。那么它其实也是一个\(L_o\),与\(L_i\)反向,光能相等。

设点p2为\(t(p, wi)\),有

\({\Large L_i(p, w_i) = L_o(t(p, w_i), -w_i)}\)

替换进渲染方程得到最终形态。方程14.13

\({\Large L_o(p, w_o)=L_e(p, w_o) + \int_{S^2} f(p, w_o, w_i)L_o(t(p, w_i), -w_i)|cosθ_i|dw_i}\)

这里面包含了\(L_o\)的递归形式。
真实世界就是实时在算这个方程,非常厉害。

14.4.2 LTE的解析解

14.4.3 基于面的LTE方程

渲染方程14.13其实相对复杂,原因是场景中的几何体的关系隐藏在\(t(p, w_i)\)中。

如果把这个式子做成显式会更好。

需要把对立体角积分转为对面积分
见5.5.3章节

\({\Large E(p, n) = ∫Li(p, w)|cos\theta|dwi}\)

\({\Large dw = cos\theta_o dA/r^2}\)

\({\Large E(p, n) = ∫La |cos\theta_i| cos\theta_o dA/r^2}\)

假设3个点。光从p3打到p2再打到p。

定义p2到p的光能为\(L(p2 \top)\)

\(L(p2 \to p) = L(p2, w)\)
其中\(w = p - p2\)
\(p2-p\)\(w_o\)
\(p3-p2\)\(w_i\)

bsdf函数,\(f(p2, w_o, w_i)\),可以写成 \(f(p3 \to p2 \to p)\)

做一个G函数

\({\Large G(p \leftrightarrow p2) = V(p \leftrightarrow p2) \frac{ |cos\theta_i||cos\theta_o| }{|p-p2|^2}}\)

其中V为两点是否相互可见。如果不可见返回0。

做替换,\(dw\)转为\(dA\)

最后渲染方程14.13变为14.15

\({\Large L(p2 \to p) = Le(p2 \top) + ∫ f(p3 \to p2 \to p) L(p3 \to p2) G(p3 \leftrightarrow p2) dA(p3)}\)

积分范围是场景内所有的面(所有的p3对应的面)。

原始的反射是对所有立体角积分。现在转成了对场景内所有的微小面积分。

14.4.4 对path积分

有了方程14.15。可以推出更灵活的LTE,包含path信息,能延伸到更高维的path空间。

对于程序来说可以做成显式的积分,更好处理,而不是最初的较笨拙的递归形式。

设4个点p0/p1/p2/p3。p0为相机。有

把式子拉长、对齐,更直观。

L(p1->p0) = Le(p1->p0)
      + ∫ L(p2->p1)                                                   f(p2->p1->p0) G(p2<->p1) dA(p2)

          展开L(p2->p1)

          = Le(p1->p0)
          + ∫ (Le(p2->p1) + ∫ L(p3->p2) f(p3->p2->p1) G(p3<->p2) dA(p3))  f(p2->p1->p0) G(p2<->p1) dA(p2)

          把展开得到的Le项和∫分别乘出来

          = Le(p1->p0)
          + ∫ Le(p2->p1)                                  f(p2->p1->p0) G(p2<->p1) dA(p2)
          + ∫ ∫ L(p3->p2) f(p3->p2->p1) G(p3<->p2) dA(p3) f(p2->p1->p0) G(p2<->p1) dA(p2)

          再展开

          = Le(p1->p0)
          + ∫ Le(p2->p1)                                                                    f(p2->p1->p0) G(p2<->p1) dA(p2)
          + ∫ ∫ Le(p3->p2)                                  f(p3->p2->p1) G(p3<->p2) dA(p3) f(p2->p1->p0) G(p2<->p1) dA(p2)
          + ∫ ∫ ∫ L(p4->p3) f(p4->p3->p2) G(p4<->p3) dA(p4) f(p3->p2->p1) G(p3<->p2) dA(p3) f(p2->p1->p0) G(p2<->p1) dA(p2)

          可无限次展开,每次展开就是多一次光的弹射。
          可发现规律

          每次展开会多出一项,下一层的Le乘上层的最末项(除去本身)。
          并把最末项(除去本身)增加一重反射积分。
代码实现
可观察到\(\int f()G()dA\)是一个连乘的模式。
越内层的积分就是越下层的弹射光。用蒙特卡洛方法,除以pdf就可以轻易地得到积分。
而内层的积分可以看成常数,跟外层无关,所以每次外层的积分可以在下一层复用。代码用简单循环就行,不用递归先算内部。
可以定义一个\(beta\),每层对应一个\(beta\)值,转到下一层时更新\(beta\),乘上当前的\(\int f()G()dA\)

把L分成直接光和间接光。

∫ L(p4->p3) f(p4->p3->p2) G(p4<->p3) dA(p4)
也就是对于每一次弹射,把这个计算分成直接光和间接光。
直接光可以对所有光源采样。间接光就是对材质采样出一个path,走蒙特卡洛。

采样出来的path,在蒙特卡洛的思想中,就代表了打到场景中所有的物体的面的path,也就是近似得到了所有的间接光。

到此已经实现了全局光照!

14.5 路径追踪

一些细节解释


15 光的传播2-体渲染

我是先看[Realistic Image Synthesis Using Photon Mapping]学了光子映射。
其中包括比较完整的体渲染理论,还包含ray marching/次表面问题等等。包括了不少pbrt第11章的内容。
见文档[光子映射学习]
推荐先看这本书,先把微分方程学一下,把体渲染公式推出来。推导其实我是参考这里,非常好的资料(有些小错)。
然后这一章再从pbrt的角度学习体渲染,并实现体渲染积分器。
后续再实现光子映射的体渲染。
pbrt是从path tracing的角度求能量Li,每算一条光就用渲染方程得出一个点的能量。
光子映射是先只考虑每个光子的能量,光子每走一步都遵循体渲染方程更新能量。射完所有光子最后再算每个点的能量。

15.1 传输方程

传输方程是光在介质中传输的核心。

上一章的光的传输方程其实是传输方程的特例,省略了介质的影响。

复习体渲染方程
对于一个光子来说
能量改变 = 介质发光 + (in-scattering) - (out-scattering)
对于处于介质中的眼睛来说,我们要求的是某个方向射过来的所有能量的和。
设眼睛为p,看向某个点p0,都处于某种介质中。 这里就和path tracing不太一样了,path tracing里p和p0之间是真空,完全不用处理。
而体渲染里p和p0之间的所有点都是要考虑的。
体渲染模型本质上就是把p和p0之间的所有点,对p(眼睛)的能量贡献加起来。
重要的是可以把p和p0之间的每个点独立考虑。
每个点都有一个增加的能量,或者说一个初始能量,这个能量通过介质传到眼睛,途中按系数衰减。把所有点的能量都加起来就行了。

首先对于任何能量的衰减,都有微分方程

\({\Large L'(s) = -\sigma(s) L(s)}\)

s是移动的距离,解得

\({\LARGE L(s) = L(0) e^{-\int_{0}^{s}\sigma(t)dt}}\)

推导见[光子映射学习]
右边的exp其实就是初始值的倍数。这时令L(0)为1,就得到了距离为s的两点之间的能量衰减的程度!

\({\LARGE T_r(s) = e^{-\int_{0}^{s}\sigma(t)dt}}\)

这个在pbrt叫做beam transmittance,有的叫做optical depth。对任何能量都适用。

再回到问题,介质中某\({\Large (p, \omega)}\)上总的能量增加为介质发光(Le) + (in-scattering),对于任何一点p都可以采样计算。

\({\Large L_{gain}(p, \omega) = L_e(p, \omega) + \sigma_s(p, \omega) \int_{S^2}^{}p(p, \omega', \omega)L_i(p,\omega')d\omega'}\)

那么p和p0之间的所有点对p的能量贡献为

\({\Large L_i(p, \omega) = \int_{0}^{s} T_r(s') L_{gain}(p', \omega)ds'}\)

其中s为p0到p的距离。
如果s取无穷大,就是一个无限大小的介质。
如果s取某个正常值,更符合0距离处于某个面的情景。
比如说有一个面处于介质中,眼睛看面上某个点,眼睛到这个点距离为s。
这样的话还需要加这个点上的初始能量。
跟中途的点都是一个道理,算一次到眼睛的衰减即可。

\({\Large L_i(p, \omega) = \int_{0}^{s} T_r(s') L_{gain}(p', -\omega)ds' + T_r(s)L_o(p_0, -\omega)}\)

再回头看看体渲染方程

\({\LARGE L(s) = \int_{0}^{s} e^{-\int_{t}^{s}p(t')dt'} q(t)dt + L(0)e^{-\int_{0}^{s}p(t')dt'} }\)

惊人的相似。可以说本质上完全一致!


15.2 采样

15.2.1 均匀介质

再观察一下transmittance,距离和衰减比例的关系。
得到的值是剩余值。

\({\LARGE T_r(s) = e^{-\int_{0}^{s}\sigma(t)dt}}\)

只考虑均匀介质有

\({\Large T_r(s) = e^{-\sigma s}}\)

用图形计算器画出图像看下,是一个快速衰减的log曲线。
\(\sigma\)随便取值粗略看下效果。
取0.001,s到8000多单位时基本完全衰减。可模拟晴朗的大气。
取2时,s到5个单位时已经基本完全衰减。比较像浓雾。

s初始0时值为1,无穷大时值为0。

13.3章算过了对这种函数的采样

\({\Large X = -\frac{ln(1-\xi)}{\sigma}}\)

做简单的测试
\(\sigma\)取0.001时,采样出来的平均值为1000。
\(\sigma\)取2时,采样出来的平均值为0.5。

同时normalize后的pdf函数为

\({\Large pdf(s) = \sigma e^{-\sigma s}}\)

那么就可以得到0到s的衰减值的蒙特卡洛积分了。


问题是为什么这样采样?我也没全懂,不过有个性质是,这个pdf的期望值是\(1/\sigma\)
刚好就是光子在介质中发生一次交互前平均移动的距离。
具体使用时其实就是ray marching。每次前进时做一次采样,作为前进的距离。
这样正好符合\(\sigma\)所定义的规律。平均移动\(1/\sigma\)后做一次交互。

这也是[Realistic Image Synthesis Using Photon Mapping]10.4.1里讲的Adaptive Ray Marching

期望值:

\({\Large E[X] = \int_{0}^{\infty} x\cdot pdf(x)dx}\)

根据[托马斯微积分]第8.2的Integration by Parts方法

\({\Large \int_{a}^{b} uv'dx = uv|_{a}^{b}-\int_{a}^{b} vu'dx}\)

那么

\({\Large E[X] = \int_{0}^{\infty} s\cdot pdf(s)ds}\)

\({\Large = \int_{0}^{\infty} s(\sigma e^{-\sigma s})ds}\)

\({\Large = (-se^{-\sigma s})|_{0}^{\infty} - \int_{0}^{\infty} -e^{-\sigma s}\cdot 1 ds}\)

\({\Large = (-se^{-\sigma s})|_{0}^{\infty} - \frac{ e^{-\sigma s}}{\sigma}|_{0}^{\infty}}\)

\({\Large = \frac{1}{\sigma}}\)


采样的代码不好看,直接看看不懂。需要先推公式。
以体渲染方程为基础

\({\Large L_i(p, \omega) = \int_{0}^{s} T_r(s') L_{gain}(p', -\omega)ds' + T_r(s)L_o(p_0, -\omega)}\)

采样流程实际是做ray marching。先预先求一个不考虑介质的交点(比如打在某个物体的面上)。得到距离赋给tMax。
p0就是预计的交点。如果没有交点那tMax就是默认的无穷大。

再对transmittance采样出距离,分两种情况。tMax左边和右边。

如果距离落在tMax之外

就相当于跳过了两点之间的所有中间点。只需计算L0也就是p0的能量,带上衰减。

能量剩余值就是交点上的Tr。落在tMax之外,是一个范围,按整体看,它的总的pdf是pdf函数从tMax积分到无穷大。
13.3章推过

\({\Large pdf(s) = \sigma e^{-\sigma s}}\) (normalize后)

从0开始到某个s的积分

\({\Large \int_{0}^{s} \sigma e^{-\sigma s'}ds' = -e^{-\sigma s}|_{0}^{s} = 1 - e^{-\sigma s}}\)

那么从x到无穷大的积分就是用1减它

\({\Large pdf = e^{-\sigma s} = Tr}\)

正好与Tr相等。那么按蒙特卡洛,结果就是采样值除以pdf,等于1。
总体是对的,不过要做平均操作。具体见代码。
如果距离落在tMax之内
就是走一步后还在介质内。那么只考虑介质的计算,也就是只算方程的积分部分。
正常采样即可

采样值为Tr,pdf为\(\sigma\)Tr。同样pdf要做平均。

另外考虑in-scattering的系数

\({\Large L_i(p, \omega) = \int_{0}^{s} T_r(s') L_{gain}(p', -\omega)ds}\)

\({\Large L_{gain}(p, \omega) = L_e(p, \omega) + \sigma_s(p, \omega) \int_{S^2}^{}p(p, \omega', \omega)L_i(p,\omega')d\omega'}\)

忽略自发光Le。Sample函数返回的Beta值要再乘\({\Large \sigma_s}\)

这样就完成了各向同性介质的采样函数HomogeneousMedium::Sample


数学问题
这个采样里面我感觉数学上还是有东西。没完全懂。
整个方程是按蒙特卡洛来算,但右边是个定值,在tMax右边的整个分布上是定值。
怎么想都不自然。
平时都是按单点来算的,现在貌似有个合并的思想,把右边当成一整块来考虑,所以右边的pdf是要算积分的。

以后再研究。


15.2.3 对phase函数采样

\({\LARGE p(\theta, \phi) = \frac{1-g^2}{4\pi(1+g^2-2gcos\theta)^{1.5}}}\)

首先对于立体角积分,它本质是有两个变量,因为和\(\phi\)无关所以可以简写。

\({\LARGE p(\theta) = \frac{1-g^2}{4\pi(1+g^2-2gcos\theta)^{1.5}}}\)

这个函数如何采样。得搞懂它的原理。
并且我觉得实际有不少细节,并不直观。而书里面都是直接出结果,这不好。
它给定了积分值是1。是对立体角积分,得按第5.5.2章的积分方法。
第5.5.2章是半球积分,而这里是所有方向所以\(\theta\)积分域是0到\(\pi\)

\({\Large \int_{0}^{2\pi } \int_{0}^{\pi } p(\theta) sin\theta d\theta d\phi = 1}\)

\({\Large \int_{0}^{\pi } p(\theta) sin\theta d\theta = 1 / 2\pi}\) (\(\phi\)无关,可以分开算出来)

\({\Large \int_{0}^{\pi } \frac{1-g^2}{4\pi(1+g^2-2gcos\theta)^{1.5}} sin\theta d\theta = 1 / 2\pi}\)

\({\Large u = cos\theta}\),有 \({\Large du/d\theta = -sin\theta}\),可以对积分进行替换,积分域变为[1, -1]。
刚好能把积分里的\(sin\theta\)消掉,我感觉是这个函数的发明者故意这样凑的?
这样变换成了更容易采样的\({\Large p_u}\)

\({\Large \int_{1}^{-1 } - \frac{1-g^2}{4\pi(1+g^2-2gu)^{1.5}} du = 1/2\pi}\)

处理负号

\({\Large \int_{-1}^{1 } \frac{1-g^2}{4\pi(1+g^2-2gu)^{1.5}} du = 1/2\pi}\)

(到这里其实可以算积分,验证这个积分值是否正确,也就是验证这个phase函数的cdf是否为1。我是在后面算出积分时验证了一下,是对的。)

这其实已经到了inversion方法的求积分步骤

\({\Large \int_{-1}^{1 } p(u) du = 1/2\pi}\)

本质仍然是对原函数采样,不过变量换成了\(u=cos\theta\),那么采样出来的值就是\(cos\theta\)。积分域也发生变化。

仍然需要做normalize。 因为右边是确定的值,所以可以直接除成1。

\({\Large \int_{-1}^{1 } \frac{1-g^2}{2(1+g^2-2gu)^{1.5}} du = 1}\)

这样做好了准备工作,按inversion方法求导。

\({\Large cdf(u)=\int_{-1}^{u} pdf(u')du'}\)

要注意是从-1开始积分!

\({\Large cdf(u)=\int_{-1}^{u} \frac{1-g^2}{2(1+g^2-2gu')^{1.5}}du'}\)

\({\Large cdf(u)= \frac{1-g^2}{2g\sqrt{1+g^2-2gu'}}|_{-1}^{u}}\)

\({\Large cdf(u)= \frac{1-g^2}{2g\sqrt{1+g^2-2gu}} - \frac{1-g}{2g}}\) (这里可以验证一下,u取1时,确实积分为1)

\({\Large cdf(u) + \frac{1-g}{2g} = \frac{1-g^2}{2g\sqrt{1+g^2-2gu}}}\)

\({\Large \sqrt{1+g^2-2gu} = \frac{1-g^2}{2g(cdf(u) + \frac{1-g}{2g})}}\)

\({\Large 1+g^2-2gu = (\frac{1-g^2}{2g(cdf(u)) + 1-g})^2}\)

得到反函数

\({\Large u = \frac{1}{2g}( 1+g^2- (\frac{1-g^2}{2g(cdf(u)) + 1-g})^2)}\)

采样

\({\Large cos\theta = \frac{1}{2g}( 1+g^2- (\frac{1-g^2}{2g\xi_1 + 1-g})^2)}\)

\({\Large \phi = 2\pi\xi_2}\)

其他参考

https://www.astro.umd.edu/~jph/HG_note.pdf


15.3 体光传输


16 光的传播3-双向方法

16.1

16.2 STOCHASTIC PROGRESSIVE PHOTON MAPPING

pbrt直接讲sppm,从理论开始讲,看不懂。先看别的专门讲光子映射的材料可能更容易。

光子映射是一步步进化的,从pm到ppm到sppm。内容很多。
我从最基础的论文开始。
资料:
[Realistic Image Synthesis Using Photon Mapping]
这是光子映射发明者Henrik Wann Jensen的书,专讲Photon Mapping。
这是Jensen的一个课件,基本是书的精简版。
[Photorealistic Rendering using the PhotonMap]
Wald

光子映射属于粒子追踪算法。从光源开始。创建path,最终连到camera。

16.2.1 partical-tracing的理论基础

大致概念是,从光源射出一堆能量粒子,有的附着在物体表面,有的发生弹射。
最后计算这些例子对应的能量,得到最终的图形。

直觉上比较简单,但问题也很多。比如能量如何衰减?cos项如何处理?

需要采用Veach的一个框架。它能给出一个坚实的理论基础。
的4.A。

粒子追踪算法在场景中的点\(p_j\)生成N个样本的光照。

\({\Large (p_j, \omega_j, \beta_j)}\)

每个样本记录入射角度\(\omega_j\),吞吐量权重\(\beta_j\)
吞吐量权重包含一个吞吐量函数T,和对应的pdf。

给定importance function \(W_e(p, \omega)\)

。。。



17 总结与展望


资源/参考

google.com
wikipedia.org
[Physically Based Rendering]
!!!
[THOMAS' CALCULUS]
微积分大部头。有完整的习题答案。
[Real-Time Rendering]
!!!
[Fundamentals of Computer Graphics]
!!!
[Mathematics for Computer Graphics]
[3D Math Primer for Graphics and Game Development]
[Mathematics for 3D Game Programming and Computer Graphics]
3d数学参考
[Introduction to Probability Models]
[A FIRST COURSE IN PROBABILITY]
概率大部头
基本图形学知识
非常好。每讲一些概念都有对应的实际输出。
GAMES101-现代计算机图形学入门-闫令琪
非常好的视频教学
99行代码的全局光照渲染器
Learn Computer Graphics From Scratch
从头学图形学。包含很多文章/知识点。

[Realistic Image Synthesis Using Photon Mapping]

[Advanced global illumination]

[Accuracy and Stability of Numerical Algorithms]

[Ray Tracing Deformable Scenes Using Dynamic Bounding Volume Hierarchies]

ROBUST MONTE CARLO METHODS FOR LIGHT TRANSPORT SIMULATION

[A Practical Guide to Global Illumination using Photon Maps]

是一本技术整理。不少内容是照搬/翻译论文等资料。

Vulkan学习

渲染器也需要交互,是朝着游戏引擎的形态走。图形api是绕不开的。
这里初期甚至pbr出图都没法在窗口显示,太差劲。
必须先把基本的流程和交互搞清楚。
正好学习一下新一代的图形api。
先看[Learning Vulkan]第一章
再看https://vulkan-tutorial.com。
同时看imgui的vulkan代码,和https://github.com/SaschaWillems/Vulkan这个代码。
同时自己新起一个代码,把流程一点点抄过来。实现画一个三角形。
概念比较多,只能硬看。初期就是各种query/fetch/create。

如果api报错,得把validation layer整上。有些不会导致崩溃的错误也能显示出来。

常见错误是某个配置里数量填了非0,但指针填了nullptr。

各家教程和代码都不一样,抄歪了的话,很可能造成配置对不上。
shader代码与配置不符也会出问题。

画三角形大概要一千行代码。包含基本的vulkan流程和概念。

再嵌入imgui。同时把三角形的vulkan流程和数据做抽象和封装。非常有助于了解vulkan。
可参考https://github.com/SaschaWillems/Vulkan 的结构。

[Learning Vulkan]

physical device
每个物理设备可以有多个queue。
queue可属于多种family。每个family都有各自的功能。
比如图形相关,计算相关,数据转移,内存管理。

每个family提供。。???

显式地管理memory。暴露device上不同种类的heap(属于不同memory region)。
memory两大类。host/device。
app进行command recording,到command buffer。
command buffer提交到queue。再被device处理。
command类型
action/set state/Synchronization

queue被驱动按提交顺序读取。

command buffer的创建expensive。
可以被cache,可重复提交。
多个command buffer可并行创建。
传统api是单线的。

app的职责

  1. 为command的执行做准备
    比如准备各种resource,编译shader。
    设置render状态。创建pipeline。draw call。
  2. 管理memory
  3. 同步
  4. 异常处理

queue

app要确保把cb提交到正确的queue。

queue类型
1.single queue
单线。保证顺序。cb顺序执行。
2.multiple queue
多线。不保证顺序。app必须负责同步。
同步机制
1. semaphore
  1. event

  2. fence

  3. pipeline barrier

object model

app层面。各种实例,比如device,queue,cb等等都称为object。

api层面,object以handle来识别。

  1. dispatchable handle

  2. non-dispatchable handle

object生命周期

app必须负责创建和销毁。

create/destroy
allocate/free

vulkan app的结构

hardw/softw init
presentation init <- wsi(window system integration)
resource setup
pipeline setup
descriptor/pool
shader with spir-v <- spir(glsl/hlsl/llvm) 提供预编译的shader
pipeline management
command recording
command submission

硬件初始化

loader进行工作
找driver
。。。
初始化完后,app负责:
1. create instance
2. 向device查询可用的queue
3. 查询extension,并存好。
4. 使能injectable layer。用作error checking/debugging/validating

Window presentation surfaces

以wsi管理window。
提供swapchain。first image/second image。
显示的同时准备下一步的数据。

wsi可看作display和app之间的接口。

resource setup

setup意味着存数据。
vulkan提供最底层的access control。可对memory进行精细的操控。
memory heap分两大类。
host local
host可见
device local
host可见
host不可见

分配内存expensive。所以要先分配内存再suballocate objects。

app负责:
1. 创建resource obj
2. 查询对应的memory instance。创建memory obj。
3. Get the memory requirements for the allocation.
4. 分配内存。存数据。
5. bind内存与resource obj。

pipeline setup

Descriptor sets and descriptor pools

Descriptor set是resource和shader之间的接口。
把shader和shader要使用的memory进行绑定。
比如shader与image,shader与buffer。

Shaders with SPIR-V

可用vulkan的自带工具glslc.exe编译glsl等shader代码到spv格式。
再用vkCreateShaderModule生成

Pipeline management

Recording commands

Queue submission

code

取extension相关数据
vkEnumerateInstanceExtensionProperties

vkEnumerateInstanceLayerProperties // 可选 if (settings.validation)

vkCreateInstance

vkEnumeratePhysicalDevices

对每个device
vkGetPhysicalDeviceProperties
选VK_PHYSICAL_DEVICE_TYPE_DISCRETE_GPU作为gpu
vkGetPhysicalDeviceFeatures // 可选
vkGetPhysicalDeviceMemoryProperties // 可选
获取GRAPHICS queue
vkGetPhysicalDeviceQueueFamilyProperties
选VK_QUEUE_GRAPHICS_BIT
用gpu和queue index创建logical device
vkCreateDevice
vkCreateCommandPool
vkCreateDescriptorPool
创建surface
glfwCreateWindowSurface

vkGetPhysicalDeviceSurfaceSupportKHR

获取surface的cap/format/mode
vkGetPhysicalDeviceSurfaceCapabilitiesKHR
vkGetPhysicalDeviceSurfaceFormatsKHR
vkGetPhysicalDeviceSurfacePresentModesKHR

vkCreateSwapchainKHR

render pass
vkCreateRenderPass

vkCreateImageView

vkCreateFramebuffer

vkCreateCommandPool

vkAllocateCommandBuffers

vkCreateFence
vkCreateSemaphore

ImGui::CreateContext();