GAMES101 学习笔记 Lecture 13~15

最新推荐文章于 2025-05-18 21:34:06 发布

hijackedbycsdn

最新推荐文章于 2025-05-18 21:34:06 发布

阅读量511

点赞数

分类专栏： Games笔记文章标签：学习笔记人工智能

本文链接：https://blog.csdn.net/hijackedbycsdn/article/details/130904967

版权

Games笔记专栏收录该内容

23 篇文章

订阅专栏

往期作业汇总帖

https://games-cn.org/forums/topic/allhw/

别人的总结

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

有很多是抄这里的，毕竟懒得手打

作业 5 框架代码详解

https://blog.csdn.net/qq_41835314/article/details/124969379

Lecture 13: Ray Tracing 1(Whitted-Style Ray Tracing)

Why Ray Tracing

为什么要光线追踪

做软阴影 soft shadows
光泽反射 glossy relection
间接光照 Indirect illumination

光栅化一般就是用 Phong 模型只计算一次光照，所以要实现多次光照就很难，但是也不是不行，只是比较麻烦，也很难保证物理上的正确性

一般认为光栅化是一种很快的近似方法

光栅化和光线追踪的对比

（1）光栅化不能很好的处理全局光照：阴影、光线多次弹射

（2）光栅化处理速度快，但是画面质量相对较差

（3）光线追踪结果精确，但是速度太慢

（4）光栅化更适合实时渲染，光线追踪更适合离线渲染

全局光照和环境光照的区别

（1）全局光照是指光线的多次反射的渲染

（2）环境光照是指环境中任一点均可以当成自发光源，照亮其他点

（3）两者有一定的交叉

光线追踪和路径追踪的简单对比

引用自：https://zhuanlan.zhihu.com/p/29418992

Ray Tracing：这其实是个框架，而不是个方法。符合这个框架的都叫raytracing。这个框架就是根据光路的可逆性，从视点发射ray，与物体相交就根据规则反射、折射或吸收。遇到光源或者走太远就停住。一般来说运算量不小。

Ray Casting：其实这个和volumetric可以脱钩。它就是ray tracing的第一步，发射光线，与物体相交。这个可以做的很快，在Doom 1里用它来做遮挡。

Path Tracing：是ray tracing + 蒙特卡洛法。在相交后会选一个随机方向继续跟踪，并根据BRDF计算颜色。运算量也不小。还有一些小分类，比如Bidirectional path tracing。

Ray Casting

首先要从人眼经过屏幕空间的一个像素点打出一条射线，然后取这个射线与所有物体的第一个交点

在这里插入图片描述
这里就完成了相当于光栅化中的深度测试的功能，因为我们确实是取第一个交点

获得了第一个交点之后，就考虑这个点会不会被光源照亮

因此将该交点与光源之间连线

在这里插入图片描述
如果中间没有东西挡着，这一点就会被照亮，否则就是阴影

然后就进入着色器，写回像素的值

Recursive Ray Tracing 递归光线追踪

光路会发生投射、反射、折射

在这里插入图片描述

二次反射之后的光线与物体的交点都与光源连线（也就是连一条 shadow ray），判断是否能被照亮，将所有能被照亮的交点的颜色值累加起来

折射反射的时候能量会减弱，所以不会发生那种折射多很多光路时，累加起来过曝的情况

光线与几何表面求交

光线的定义

在这里插入图片描述
r(t) 表示一个射线

光线与球相交

联立球和射线的方程

在这里插入图片描述
现在要解 t，就是解一个一元二次方程

在这里插入图片描述

光线与一般的隐式表面求交

在这里插入图片描述

实数，要有物理意义
t>0，表示传播时间

光线与显式表面求交

在一个多边形所定义的平面内随便点一个点，从这个点出发向任意方向打一条射线

如果该射线与多边形的交点为奇数个，说明射线的起点在多边形内

如果该射线与多边形的交点为偶数个，说明射线的起点在多边形外

感觉是因为射线的起点在多边形内的话，射线就是只出不进，所以是奇数……？

相切的话应该不算相交吧

这个规律推广到 3d 依然成立，要求物体不能有洞

具体怎么计算呢

简单的方法就是

对每一个三角形与给定的光线求交点，只可能有零个或者一个

可以有一些方法来加速

光线与三角形求交

前面提到了对每一个三角形与给定的光线求交点，那么具体来说怎么光线与三角形求交？

分为两步：

光线与三角形所定义的平面求交
判断光线与三角形所定义的平面的交点是否在三角形内

不考虑光线与三角形所定义的平面完全平行的情况

平面由平面上的一点 P’ 与平面的法向量所定义

那么平面上任意一点 P 的定义式就是，(P - P’)·N = 0

设 P(x,y,z) P’(x0,y0,z0)，最终定义式可以写为 ax+by+cz+d=0

在这里插入图片描述

不分步的光线与三角形求交 Moller Trumbore Algorithm

Moller Trumbore Algorithm

直接将射线方程与重心坐标联立

在这里插入图片描述推导：

$O - P_0 = -tD + b_1(P_1 - P_0) + b_2(P_2 - P_0) \\ S = -tD + b_1E_1 + b_2E_2$

使用克拉默法则计算系数

$\dfrac{[S, E_1, E_2]}{[-D,E_1,E_2]} \\ b_1 = \dfrac{[-D, S,E_2]}{[-D,E_1,E_2]} \\ b_2 = \dfrac{[-D, E_1,S]}{[-D,E_1,E_2]} \\$

三个向量组成的矩阵的行列式等于他们的混合积

混合积前面加一个负号，相当于原混合积的叉乘的两个向量的位置对换

其中

$E_1, E_2] = (S \times E_1) \cdot E_2 = S_2 \cdot E_2$

$S,E_2] = [E_2, -D, S] = (E_2 \times -D) \cdot S = (D \times E_2) \cdot S = S_1 \cdot S$

$E_1,S] = -D \cdot (E_1 \times S) = D \cdot (S \times E_1) = D \cdot S_2$

$[-D,E_1,E_2] = [E_2, -D, E_1] = (E_2 \times -D) \cdot E_1 = (D \times E_2) \cdot E_1 = S_1 \cdot E_1$

代入就得到公式

不知道为啥写不出来他那个形式的话，可以试试轮换混合积

要判断点是否在三角形，需要 t,b1,b2 都是非负的

Bounding Volumes 包围盒

先让光线与包围盒求交，如果光线与包围盒不相交，则不会与包围盒内任意面片相交，则不需要求光线与包围盒内面片的交点。

AABB包围盒（Axis-Aligned Box Bounding）：可以快速求取得到光线与平面交点

在这里插入图片描述
（1）t_enter = max(t_x_min, t_y_min, t_z_min)；进包围盒的点，为所有入点的最大值

（2）t_exist = min(t_x_max, t_y_max, t_z_max)；出包围盒的点，为所有出点的最小值

（3）相交条件：t_enter < t_exist，且 t_exist > 0.

对于三维情况，就是在三个方向上有三组 t_min t_max

如果 t_exist < 0 说明光线在出发点背后与物体相交

如果 t_enter <= 0, t_exist > 0 说明光线出发点在包围盒内

为什么包围盒的三组对面是水平或者竖直的？因为好求

在这里插入图片描述

Lecture 14: Ray Tracing 2(Acceleration & Radiometry 辐射度量学)

GTC

DLSS 2.0 Deep Learning Super Sampling

使用深度学习做超采样

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

RTXGI

加速结构

Uniform grids 均匀网格法

（1）寻找最外层AABB包围盒

（2）建立均匀分割的网格

（3）将三角面片保存到所在网格（跨网格的需要存储多份）

只需要保存物体的表面，不需要保存物体的空心的内部

在这里插入图片描述
这个图中的右上角的好像少涂了一个网格阴影？

（4）光线按光路顺序找网格，若当前网格没有存储三角面片，则继续下一个网格

光线与网格求交，相比于光线与网格内所有三角面片求交，更加简单

假设光线已经打到某个网格，怎么知道下一次与哪个格子求交？

最简单的想法就是，根据光线的方向来判断，假设这是在二维平面中，光线斜率为正，那么下一次打到的格子一定是位于当前格子的上面一格或者右边一格

对两个格子分别求交，看看实际与哪个格子有交点就好了

这个方法也可以拓展到三维

实际上可能不是这么做的，但是这个简单的想法是基本的

（5）若当前网格存储有三角面片，则与该网格内所有三角面片求交，若无交点，则继续下一个网格

（6）若光线与当前网格的三角面片有交点，则取第一个交点，结束遍历。

在这里插入图片描述

什么时候使用均匀网格法

网格密度过小，则加速效果不明显；网格密度过大，则与网格求交也会消耗大量的资源，过犹不及。

如果场景中出现了大面积的空白，那么也会有很多时间浪费在与格子求交上

空间分区方式

因为均匀网格法可能遇到场景中出现大面积空白的情况，我明明一直知道了这一大片是空白的，但是还是在这里有很多光线与网格求交

所以希望有一个方法将空间切成不均匀的网格

Oct-Tree 八叉树：每个划分都需要继续切成8块，直至块内的对象个数较少或者切割次数到达一定数值。

八叉树在维度比较大的时候（比如 3），划分出来的网格也比较多，是 8 的指数

KD-Tree 二叉树：可以先x再y后z，每个划分都需要被化

二维中，交替地水平、竖直、水平、竖直地划分，使空间尽可能均匀

三维中就是 x y z 这样轮换

BSP-Tree

1、首先，选择一个分区超平面。为了方便讨论，我们将使用一个二维世界，而我们的根节点将是一条线。

2、使用初始分区超平面对世界上的所有多边形进行分区，并将它们存储在正面或背面多边形列表中。(Partition all polygons in the world with the initial partition hyperplane, storing them in either the front or back polygon list.)

3、在前后多边形列表中进行递归或迭代，创建一个新的树节点，并将其附加到父节点的左侧或右侧叶子。(recurse or iterate through the front and back polygon list, creating a new tree node and attaching it the the left or right leaf of the parent node.)
————————————————
版权声明：本文为CSDN博主「packdge_black」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/packdge_black/article/details/114681992

计算超平面还是有点复杂的，所以这里先介绍更简单的使用维度来切分的 KD Tree

BVH

KD Tree

KD-Tree 结构

非叶子结点存储：

（1）分割的坐标轴x或y或z

（2）分割的具体坐标位置，包括xmin、xmax、ymin、ymax、zmin、zmax即切割的坐标位置p

（3）子节点指针

（4）不存储三角面片

叶子结点存储：

（1）内部三角面片集合

在这里插入图片描述

KD-Tree构建

切割方向，第1刀x，第2刀y，第3刀z。

实际上，可以不用正中间切割，可以根据需要切割。

KD-Tree遍历

伪代码如下

function traversing(node, light) 
{
    计算光线与node结点包围盒的交点：tmin、tmax;
    if 与包围盒不存在交点:
        return null;
    if 孩子结点是叶节点:
        while 遍历tmin所在孩子结点的三角面片:
            if 有交点:
                记录多个交点的最小的t_result_min所在点信息
        if 存在交点t_result_min:
            找到交点
            return t_result_min所在点信息
        while 遍历另一个孩子结点的三角面片:
            if 有交点:
                记录多个交点的最小的t_result_min所在点信息
        if 存在交点t_result_min:
            找到交点
            return t_result_min所在点信息
        与此包围盒内三角面片不存在交点
        return null

    result = traversing(node的tmin所在的子节点)
    if result != null:
        return result
    result = traversing(node的另一个子节点)
    if result != null:
        return result
    与此包围盒内子包围盒也不存在交点
    return null
}

最终结果result = traversing(root, light)

emmm 感觉这个遍历孩子节点的面片这里是有点怪的

因为如果某一个节点 A 的孩子节点是叶节点的话，这个节点 A 的另一个孩子节点并不一定也是叶节点吧……

课件的例子：

给定 KD Tree 的根节点

那么光线先与根节点求交

在这里插入图片描述
判断到光线确实与这个节点相交了，那么就与这个节点的孩子节点求交

在这里插入图片描述

如果与光线相交的孩子节点为叶节点，那么光线就与这个叶节点中的所有三角形面片求交

在这里插入图片描述
如果与光线相交的孩子节点为非叶节点，那么就在这个非叶节点递归判断

缺点

对于有一个三角形构成的场景中，很难判定一个三角形和包围盒有交集

一个包围盒被划分成了很多个小格子后，一个物体可能会在多个格子里，那么在这些格子都要保存这个物体，最后判断这个物体和光线是否相交也要执行很多次。

层次包围盒 Bounding Volume Hierarchy（BVH）

之前的方法都是按照空间来划分的，这个是按照物体划分的。

这个方法现在广泛应用在光线追踪中

结构

非叶子结点存储：

（1）包围盒

（2）子节点列表

叶子结点存储：

（1）包围盒

（2）三角形面片列表

优点：一个物体只能出现在一个包围盒

缺点：两个不同的包围盒可能相交

构建

构建流程：

（1）寻找包围盒

（2）将三角形面片集合按一定规则划分为两个子集合

（3）计算两个子集合的包围盒

（4）递归对子集合进行划分

（5）直至子集合满足一定要求

切分子集合方式：

（1）选择一个方向去切割

（2）总是选择最长的边进行切割

（3）选择这个被切分的包围盒的“中间”的那个物体切分

我看到别人的笔记写的是“选择中间的三角面片切割，保证两侧的三角形个数一致”

虽然我的目的确实是“保证两侧的三角形个数一致”

但是“选择中间的三角面片切割”这种描述听上去就像是划分空间来建树的，这是 KD Tree 的做法

但是 BVH 是以物体作为包围盒中的基础单元来建树的

所以感觉这个说法会须会有一些混淆？

这其中涉及到了选择中位数的问题

更进一步，这就是一个 topK 问题，典型的解决思路是快速选择算法，快选就是对快排的修改，在快排的排序过程中，只排序一边

遍历

Intersect(Ray ray, BVH node) {
	if (ray misses node.bbox) return;		//光线和包围盒不相交
	if (node is a leaf node)				//该节点时叶子节点
		test intersection with all objs;		//测试光线和该节点里面的所有物体是否相交
		return closest intersection;			//返回最近的交点
	hit1 = Intersect(ray, node.child1);	//递归左子树
	hit2 = Intersect(ray, node.child2);	//递归右子树
	return the closer of hit1, hit2;
}

BVH 中的物体移动了该怎么办

只能重新计算了……？

KD Tree 和 BVH 的对比

在这里插入图片描述
KD Tree 划分空间一个物体可能被包含在多个叶子节点中

BVH 划分物体的集合一个物体只存在于一个集合中

Radiometry 辐射度量学

物理正确的光照计算。

Q：光的能量 energy

辐射能是电磁辐射的能量，以焦耳 J 为单位

Φ：Radiant Flux 辐射通量光通量 power

辐射通量是单位时间内或发射，或反射，或传送，或接受的能量，以瓦特 W 为单位

$\mathrm \Phi \equiv \dfrac{\mathrm d Q}{\mathrm d t}$

光学中有另一个单位 lm lumen 流明

Radiant Intensity 辐射强度 Light Emitted From A Source

Irradiance 辐[射]照度 Light Falling On A Surface

Radiance 辐[射]亮度 Light Traveling ALong A Ray

在这里插入图片描述

Radiant Intensity

单位立体角所发射的能量

$\equiv \dfrac{\mathrm d \Phi}{\mathrm d w}$

单位为 $\dfrac{W}{sr}$ 其中 sr 是立体角的单位

或者 $\dfrac{lm}{sr} = cd = candela$

角度相关的定义

对于二维

$\theta = \dfrac{l}{r}$

角度 = 弧长/半径

那么一个圆的角度为 $2\pi$

对于三维

$KaTeX parse error: Undefined control sequence: \Omgea at position 1: \̲O̲m̲g̲e̲a̲ ̲= \dfrac{A}{r^2…$

立体角 = 面积/半径的平方

那么一个球的角度为 $4\pi$

单位面积在球坐标系中的表达式为：

在这里插入图片描述
进而得到单位立体角在球坐标系中的表达式

可见，面积不仅仅是 $\mathrm d \theta \mathrm d \phi$ ，还有一个 $\sin\theta$ 的系数，这就说明面积在球面上不是一个均匀的划分，越靠近极点越小

可以验证 $\Omega = \int_{S^2}\mathrm d w = \int_{0}^{2\pi}\int_{0}^{\pi}\sin\theta\mathrm d \theta \mathrm d \phi = 4\pi$

Radiant Intensity 与 Radiant Flux

已知辐射强度的定义为 $\equiv \dfrac{\mathrm d \Phi}{\mathrm d w}$

并且我们也知道单位立体角在球面上的积分为 $\Omega = \int_{S^2}\mathrm d w = \int_{0}^{2\pi}\int_{0}^{\pi}\sin\theta\mathrm d \theta \mathrm d \phi = 4\pi$

那么假设一个球面上的辐射强度是常量，那么就可以计算辐射通量 $\Phi = \int_{S^2}I\mathrm d w = 4\pi I$

那么辐射通量与辐射强度之间的关系就是 $\dfrac{\Phi}{4\pi}$

作业 5

cmake 构建报错：无效的数值参数

报错类似：

cl : 命令行 error D8021: 无效的数值参数“/Wextra” [D:\Work\VisualStudioProject\Games101Homework\Games101Homework\05\build\RayTraci
ng.vcxproj]

cl : 命令行 error D8021: 无效的数值参数“/Wshadow” [D:\Work\VisualStudioProject\Games101Homework\Games101Homework\05\build\RayTrac
ing.vcxproj]

cl : 命令行 error D8021: 无效的数值参数“/Wreturn-type” [D:\Work\VisualStudioProject\Games101Homework\Games101Homework\05\build\Ra
yTracing.vcxproj]

感觉应该是因为 Linux 和 Windows 上有所不同吧

在 CMakelist 中删掉这三个参数就行了

Render()

之前我以为这个距离是不知道的

然后我还以为这个 scene 的宽高就是整个屏幕空间的宽高了

所以我就用 scene 的宽和 fov 来算摄像机距离屏幕的距离了

            // generate primary ray direction
            float x;
            float y;
            // TODO: Find the x and y positions of the current pixel to get the direction
            // vector that passes through it.
            // Also, don't forget to multiply both of them with the variable *scale*, and
            // x (horizontal) variable with the *imageAspectRatio*            

            // 一般认为 FOV 是水平视场角
            float distance = scene.height / 2.0f / scale;
            Vector3f dir = Vector3f(i - scene.width/2.0f, j - scene.height/2.0f, -distance);
            dir = normalize(dir);
            
            framebuffer[m++] = castRay(eye_pos, dir, scene, 0);

之后看到别人写的，才发现原来这里是默认了镜头距离屏幕的距离是 1

看上去应该是因为他在 dir 那里写的 z 值是 -1

那么根据这个 1 就能写出屏幕的实际宽高

那么就能得到计算方法

            // generate primary ray direction
            float x;
            float y;
            // TODO: Find the x and y positions of the current pixel to get the direction
            // vector that passes through it.
            // Also, don't forget to multiply both of them with the variable *scale*, and
            // x (horizontal) variable with the *imageAspectRatio*            

            // 一般认为 FOV 是水平视场角
            float screen_height = 1.0f * scale * 2.0f;
            float screen_width = screen_height * imageAspectRatio;

            // (i,j) 表示的像素从左到右，从上到下
            // 首先加 0.5 是要求像素中心
            // 然后减 scene.width / 2.0f 或减 scene.height / 2.0f 是要转换到世界空间中
            // 值 / scene.width * screen_width 或 / scene.height * screen_height 是将目标转换到屏幕空间中
            // y 加一个负号是因为 j 表示的像素是从上到下的
            x = (i + 0.5 - scene.width / 2.0f) / scene.width * screen_width;
            y = -(j + 0.5 - scene.height / 2.0f) / scene.height * screen_height;
    
            Vector3f dir = Vector3f(x, y, -1); // Don't forget to normalize this direction!
            dir = normalize(dir);
            
            framebuffer[m++] = castRay(eye_pos, dir, scene, 0);

地板渲染错误

原因可能是在 bool rayTriangleIntersect(const Vector3f& v0, const Vector3f& v1, const Vector3f& v2, const Vector3f& orig, const Vector3f& dir, float& tnear, float& u, float& v) 中没有注意到要返回三个系数 tnear, u, v

在这里插入图片描述

地板上的蓝色小点

https://games-cn.org/forums/topic/zuoye5dibanduijiaoxianshangdexiangsuweibeizhengquexuanranwenti/

在 rayTriangleIntersect() 中判断四个系数的正负的时候，对系数加一个小量，避免将边界上的点判断为 true

作业 5 中的框架结构

从 castRay 开始

在作业 6 中写 Render() in Renderer.cpp 的时候，按照作业 5 我写成

    for (uint32_t j = 0; j < scene.height; ++j) {
        for (uint32_t i = 0; i < scene.width; ++i) {
            // generate primary ray direction
            float x = (2 * (i + 0.5) / (float)scene.width - 1) *
                      imageAspectRatio * scale;
            float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;
            // TODO: Find the x and y positions of the current pixel to get the
            // direction
            //  vector that passes through it.
            // Also, don't forget to multiply both of them with the variable
            // *scale*, and x (horizontal) variable with the *imageAspectRatio*

            // Don't forget to normalize this direction!

            Vector3f dir = Vector3f(x, y, -1);
            dir = normalize(dir);

            framebuffer[m++] = castRay(eye_pos, dir, scene, 0);
        }
        UpdateProgress(j / (float)scene.height);
    }

但是却发现作业 6 中没有 castRay 了

于是这里就需要看作业 5 中的 castRay 是怎么样的，在作业 6 中怎么对应实现它

获取 castRay 的 payload

作业 5 中的 Vector3f castRay() in Renderer.cpp

Vector3f castRay(
        const Vector3f &orig, const Vector3f &dir, const Scene& scene,
        int depth)
{
    if (depth > scene.maxDepth) {
        return Vector3f(0.0,0.0,0.0);
    }

    Vector3f hitColor = scene.backgroundColor;
    if (auto payload = trace(orig, dir, scene.get_objects()); payload)
    {

开头考虑深度是否大于最大深度，大于则返回全黑颜色

然后初始化 hitColor

然后获得 payload

if (auto payload = trace(orig, dir, scene.get_objects()); payload) 这里使用了 C++17 的新特性，允许在 if/switch 语句中初始化一个在 if/switch 作用域内的局部变量，离开作用域就析构

payload 的具体意义

payload 这个命名其实是比较没有意义的……它的本意是“有效载荷”，但是这个意义可以指很多事的载荷……就像写 unity 的时候给脚本命名一个 controller manager 之类的hhh

所以还是要具体看 trace()

std::optional<hit_payload> trace(
        const Vector3f &orig, const Vector3f &dir,
        const std::vector<std::unique_ptr<Object> > &objects)
{
    float tNear = kInfinity;
    std::optional<hit_payload> payload;
    for (const auto & object : objects)
    {
        float tNearK = kInfinity;
        uint32_t indexK;
        Vector2f uvK;
        if (object->intersect(orig, dir, tNearK, indexK, uvK) && tNearK < tNear)
        {
            payload.emplace();
            payload->hit_obj = object.get();
            payload->tNear = tNearK;
            payload->index = indexK;
            payload->uv = uvK;
            tNear = tNearK;
        }
    }

    return payload;
}

写代码时，有时候想要返回某种类型的对象，也就是说，我们可以有某个类型的值，也可以没有任何值。需要有一种方法来模拟类似指针的定义。c++17 给出了 std::optional<T> 来定义这样的对象

将 std::optional<T> 转化为 bool 时，包含值则为 true，不包含值则为 false，因此我们可以很方便地使用

if(auto opt_var = get(); opt_var)

来获得一个可能不存在的东西，并且方便地限定它的作用域，方便地写出它的判空条件语句

接下来的部分就是判断光线是否与三角面相交，不相交则返回背景颜色，相交则根据物体的具体材质有不同的表现

intersect 判断光线与不同类型的物体之间的交点

看上去他应该是使用 trace() 来判断光线是否与三角面相交

其中 object->intersect 是判断相交的函数，其他部分是获得返回值的逻辑

根据他这个循环可以看出，他是对场景中的每一个物体都判断一遍，这个物体是否会和光线相交

intersect 是 Object 类的虚函数，对于不同子类有不同的实现

作业 5 中实现了两个子类，球和三角面物体

其中，球的相交函数就是解那个一元二次方程

三角面物体的相交函数就是遍历这个物体的每一个三角面，对每一个三角面都使用 Moller Trumbore Algorithm 算一遍是否相交

    bool intersect(const Vector3f& orig, const Vector3f& dir, float& tnear, uint32_t& index,
                   Vector2f& uv) const override
    {
        bool intersect = false;
        for (uint32_t k = 0; k < numTriangles; ++k)
        {
            const Vector3f& v0 = vertices[vertexIndex[k * 3]];
            const Vector3f& v1 = vertices[vertexIndex[k * 3 + 1]];
            const Vector3f& v2 = vertices[vertexIndex[k * 3 + 2]];
            float t, u, v;
            if (rayTriangleIntersect(v0, v1, v2, orig, dir, t, u, v) && t < tnear)
            {
                tnear = t;
                uv.x = u;
                uv.y = v;
                index = k;
                intersect |= true;
            }
        }

        return intersect;
    }

看上去这里没有提前停止之类的……但是有一个取最近的 tnear 的判断，这样保证取到的是光线与物体的第一个交点

这里也可以看出返回的 uv 是三角形的重心坐标

getSurfaceProperties 获得光线击中点的表面属性

然后回到 Vector3f castRay() in Renderer.cpp

    if (auto payload = trace(orig, dir, scene.get_objects()); payload)
    {
        Vector3f hitPoint = orig + dir * payload->tNear;
        Vector3f N; // normal
        Vector2f st; // st coordinates
        payload->hit_obj->getSurfaceProperties(hitPoint, dir, payload->index, payload->uv, N, st);

这里是，如果相交，那么就能得到交点位置，并且得到这个物体在这个交点处的性质，比如物体在这个点的法线和纹理坐标

一般我们确实是使用 UV 表示纹理坐标

但是他这个是，他已经用 uv 表示光线与物体的交点在击中的这个三角形上的重心坐标了，所以就使用 st 来表示纹理左边，免得弄混，我猜

可以看他是怎么获取纹理坐标的，就是用重心坐标来插值的

    void getSurfaceProperties(const Vector3f&, const Vector3f&, const uint32_t& index, const Vector2f& uv, Vector3f& N,
                              Vector2f& st) const override
    {
        const Vector3f& v0 = vertices[vertexIndex[index * 3]];
        const Vector3f& v1 = vertices[vertexIndex[index * 3 + 1]];
        const Vector3f& v2 = vertices[vertexIndex[index * 3 + 2]];
        Vector3f e0 = normalize(v1 - v0);
        Vector3f e1 = normalize(v2 - v1);
        N = normalize(crossProduct(e0, e1));
        const Vector2f& st0 = stCoordinates[vertexIndex[index * 3]];
        const Vector2f& st1 = stCoordinates[vertexIndex[index * 3 + 1]];
        const Vector2f& st2 = stCoordinates[vertexIndex[index * 3 + 2]];
        st = st0 * (1 - uv.x - uv.y) + st1 * uv.x + st2 * uv.y;
    }

REFLECTION_AND_REFRACTION

回到 Vector3f castRay() in Renderer.cpp，接下来就是根据表面类型，光线有不同的行为

switch (payload->hit_obj->materialType) {

REFLECTION_AND_REFRACTION 即发生反射也发生折射

REFLECTION 只有反射

default 默认情况使用 Phong

            case REFLECTION_AND_REFRACTION:
            {
                Vector3f reflectionDirection = normalize(reflect(dir, N));
                Vector3f refractionDirection = normalize(refract(dir, N, payload->hit_obj->ior));

反射的计算

这个反射的计算很神奇……

其中 I 是入射光，R 是反射光，A，B 类似入射光的分量，N 是着色点的法向量的单位向量

那么计算就简单了

Vector3f reflect(const Vector3f &I, const Vector3f &N)
{
    return I - 2 * dotProduct(I, N) * N;
}

折射的计算

大部分是看 https://www.cnblogs.com/starfallen/archive/2012/11/05/2754992.html

关于全反射角的理解

一开始我还以为这个博文里面有一点说错了……

他说完全可以让 $\dfrac{1}{\eta^2}(1-\cos^2\theta_1)$ 大于 1，使得 $\dfrac{1}{\eta^2}(1-\cos^2\theta_1)$ 小于 0……

我一开始感觉这是不可能的……因为 $|\vec{t_1}| = \sin\theta_2 = \sqrt{\dfrac{1}{\eta^2}(1-\cos^2\theta_1)}$

已知 $\leq \sin\theta_2 \leq 1$ ，即 $\leq \sqrt{\dfrac{1}{\eta^2}(1-\cos^2\theta_1)} \leq 1$

即 $\leq \dfrac{1}{\eta^2}(1-\cos^2\theta_1) \leq 1$

所以推出 $\leq 1- \dfrac{1}{\eta^2}(1-\cos^2\theta_1) \leq 1$

又或者说，这个式子本身就是 $\leq \sin\theta_2 \leq 1$ 啊

之后我才发现，原来所谓的全反射指的是不遵循折射定律了，或者也不是不遵循吧……就是入射角满足的条件使得折射定律不能成立了

就是看 $\dfrac{1}{\eta^2}(1-\cos^2\theta_1)$ 这个来说， $\leq 1-\cos^2\theta_1 \leq 1$ ，那么乘一个 $\dfrac{1}{\eta^2}$ 的话，范围就是 $\leq \dfrac{1}{\eta^2}(1-\cos^2\theta_1) \leq \dfrac{1}{\eta^2}$

本来按照我的想法， $\dfrac{1}{\eta^2}(1-\cos^2\theta_1)$ 其实就是 $\sin\theta_2$ ，那么也应该就有 $\leq \sin\theta_2 \leq 1$ ，那么 $\dfrac{1}{\eta^2}(1-\cos^2\theta_1)$ 的取值范围应该是 $[0, 1]$ 和 $[0,\dfrac{1}{\eta^2}]$ 之间的交集才对

而这个交集我之前默认是得到 $[0, 1]$

其实这么想是不对的……这里的相对折射率 $\eta = \dfrac{\eta_2}{\eta_1}$ 是两个材质之间的比较，他的取值是可能大于 1 或者小于 1 的，一般我们说的一定大于 1 的折射率是单独某一个材质的绝对折射率

所以说 $\dfrac{1}{\eta^2}$ 是有可能大于 1 的，那么也就是说，在入射角 $\theta_1$ 大于某一值时， $\dfrac{1}{\eta^2}(1-\cos^2\theta_1)$ 有可能会大于 1，而这时不是要与 $[0, 1]$ 取值范围相交，而是意味着不能存在满足条件的折射角 $\theta_2$ ，也就是说不能发生折射，那就只剩下反射了

折射的计算代码

// [comment]
// Compute refraction direction using Snell's law
//
// We need to handle with care the two possible situations:
//
//    - When the ray is inside the object
//
//    - When the ray is outside.
//
// If the ray is outside, you need to make cosi positive cosi = -N.I
//
// If the ray is inside, you need to invert the refractive indices and negate the normal N
// [/comment]
Vector3f refract(const Vector3f &I, const Vector3f &N, const float &ior)
{
    float cosi = clamp(-1, 1, dotProduct(I, N));
    float etai = 1, etat = ior; //1 近似为空气折射率，ior 为介质的折射率
    Vector3f n = N;
    //如果入射光从介质 1(etai)->介质 2(etat)，则夹角 > 90°；
    //如果入射光从介质 2->介质 1，则夹角 < 90°,用于提供方向向量的 N 也要反过来，折射率之比也需要换一下
    if (cosi < 0) { cosi = -cosi; } else { std::swap(etai, etat); n= -N; }
    float eta = etai / etat;
    float k = 1 - eta * eta * (1 - cosi * cosi);
    return k < 0 ? 0 : eta * I + (eta * cosi - sqrtf(k)) * n;
}

这里他是计算了入射角和法向量的点积

那么根据两个向量之间的角度，还需要判断到底是从物体折射到空气中还是从空气折射到物体中

光线击中点的计算误差的处理

回到 Vector3f castRay() in Renderer.cpp，接下来是对 hitPoint

                Vector3f reflectionRayOrig = (dotProduct(reflectionDirection, N) < 0) ?
                                             hitPoint - N * scene.epsilon :
                                             hitPoint + N * scene.epsilon;
                Vector3f refractionRayOrig = (dotProduct(refractionDirection, N) < 0) ?
                                             hitPoint - N * scene.epsilon :
                                             hitPoint + N * scene.epsilon;

由于计算问题会出现误差，击中点如果在物体的另外一面，就会导致光线与物体表面多相交一次

这个所谓的“击中点如果在物体的另外一面”，可以用入射方向或者反射方向来看，比如这里，我从左上角向这个圆入射，那么我的击中点应该在圆的上方；那我要是从左下角打，击中点应该在圆的下方

这样观察，就可以发现击中点期望的位置的规律，那么我们只需要将计算出来的击中点位置向着我们期望的位置移动就好了

这个代码里使用了反射方向或者折射方向与法向量的点积，也就是 $\cos(\theta_r)$ ，也都是可以的

菲涅尔

获得了光线反射折射的方向，也获得了反射或折射的光线的起点

那么就可以把光线传播的过程递归下去，直到得到反射或折射的光线最终返回的颜色

                Vector3f reflectionColor = castRay(reflectionRayOrig, reflectionDirection, scene, depth + 1);
                Vector3f refractionColor = castRay(refractionRayOrig, refractionDirection, scene, depth + 1);
                float kr = fresnel(dir, N, payload->hit_obj->ior);
                hitColor = reflectionColor * kr + refractionColor * (1 - kr);
                break;

使用菲涅尔定律混合反射和折射的颜色，就可以得到这个反射或折射的光线混合得到的颜色，就可以返回本光线的颜色了

或者说着色点的颜色？

可以看到，离我们越近的水面，投射现象越强烈，离我们越远的水面反射现象越强烈。反射光与折射光的数量可以使用菲涅耳方程来计算。光由两个垂直波组成，我们称其为平行偏振光和垂直偏振光。我们需要使用两个不同的方程式计算这两个波的反射光比率，并对结果求平均值。两个菲涅耳方程为：

$F_{R\parallel} = (\dfrac{\eta_2\cos\theta_1-\eta_1\cos\theta_2}{\eta_2\cos\theta_1+\eta_1\cos\theta_2})^2$

$F_{R\perp} = (\dfrac{\eta_2\cos\theta_1-\eta_1\cos\theta_2}{\eta_2\cos\theta_1+\eta_1\cos\theta_2})^2$

通过求平均值我们可以得到反射的比率：

$F_R = \dfrac{1}{2}(F_{R\parallel} + F_{R\perp})$

项η1，η2是两种介质的折射率。项cosθ1和cosθ2分别是入射角和折射角。如前所述，由于能量守恒，折射光的比率可以简单地计算为：

$F_T = 1 - F_R$

计算代码：

// [comment]
// Compute Fresnel equation
//
// \param I is the incident view direction
//
// \param N is the normal at the intersection point
//
// \param ior is the material refractive index
// [/comment]
float fresnel(const Vector3f &I, const Vector3f &N, const float &ior)
{
    float cosi = clamp(-1, 1, dotProduct(I, N));
    float etai = 1, etat = ior;
    if (cosi > 0) {  std::swap(etai, etat); }
    // Compute sini using Snell's law
    float sint = etai / etat * sqrtf(std::max(0.f, 1 - cosi * cosi));
    // Total internal reflection
    if (sint >= 1) {
        return 1;
    }
    else {
        float cost = sqrtf(std::max(0.f, 1 - sint * sint));
        cosi = fabsf(cosi);
        float Rs = ((etat * cosi) - (etai * cost)) / ((etat * cosi) + (etai * cost));
        float Rp = ((etai * cosi) - (etat * cost)) / ((etai * cosi) + (etat * cost));
        return (Rs * Rs + Rp * Rp) / 2;
    }
    // As a consequence of the conservation of energy, transmittance is given by:
    // kt = 1 - kr;
}

这里也是先判断是否是全反射，是的话 kr 直接 = 1，不是的话就按照菲涅尔公式计算

REFLECTION

只有反射的情况就相当于在 REFLECTION_AND_REFRACTION 中去掉了折射这一项

default

默认情况的材质是使用 phong 着色模型

与之前的作业 1~3 有所区别的是，这里还加入了对阴影的计算

                    auto shadow_res = trace(shadowPointOrig, lightDir, scene.get_objects());
                    bool inShadow = shadow_res && (shadow_res->tNear * shadow_res->tNear < lightDistance2);

                    lightAmt += inShadow ? 0 : light->intensity * LdotN;

如果从着色点到光源的光线没有击中物体，那么 shadow_res == false，自然 inShadow == false

如果从着色点到光源的光线击中了物体，那么 shadow_res == true，再判断从着色点到被击中的点的距离是否小于从着色点到光源的距离，如果确实小于，那么光源就确实被一个物体挡住了，那么 inShadow == true，否则，从着色点到光源的光线击中的物体在光源的身后，光源没有被挡住，所以没有阴影

作业 6

Renderer::Render()

之前本来以为 castRay 没了，之后发现是搬到 Scene 中了

新的 castRay 也差不多，所以很简单

            // generate primary ray direction
            float x = (2 * (i + 0.5) / (float)scene.width - 1) *
                      imageAspectRatio * scale;
            float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;
            // TODO: Find the x and y positions of the current pixel to get the
            // direction
            //  vector that passes through it.
            // Also, don't forget to multiply both of them with the variable
            // *scale*, and x (horizontal) variable with the *imageAspectRatio*

            // Don't forget to normalize this direction!

            Vector3f dir = Vector3f(x, y, -1);
            dir = normalize(dir);

            framebuffer[m++] = scene.castRay(Ray(eye_pos, dir), 0);

Triangle::getIntersection()

与之前的一对比，别人已经把三个系数都算好了，我们只需要返回一个值，包装一下相交结果就行了，只是要理解一下他为什么这么写

比如他一开始判断光线是否是从三角面的背面打过来的，如果是背面则忽略，认为不相交

还有，判断 $S_1 \cdot E_1$ 的值是否接近 0，因为他要作为分母，如果接近 0 的话就会很大误差，所以当他的值小于某一个小量时就认为不相交

inline Intersection Triangle::getIntersection(Ray ray)
{
    Intersection inter;

    if (dotProduct(ray.direction, normal) > 0)  // if ray is cast from bottom of face, ingore it
        return inter;

    double u, v, t_tmp = 0;

    Vector3f pvec = crossProduct(ray.direction, e2);  // s1
    double det = dotProduct(e1, pvec); // det = s1 dot e1
    if (fabs(det) < EPSILON)  // avoid denominator is too small
        return inter;

    double det_inv = 1. / det;

    Vector3f tvec = ray.origin - v0;  // s
    u = dotProduct(tvec, pvec) * det_inv;  // u = s1 dot s / (s1 dot e1)
    if (u < 0 || u > 1)
        return inter;

    Vector3f qvec = crossProduct(tvec, e1);  // s2 = s cross e1
    v = dotProduct(ray.direction, qvec) * det_inv;  // v = s2 dot dir / (s1 dot e1)
    if (v < 0 || u + v > 1)
        return inter;

    t_tmp = dotProduct(e2, qvec) * det_inv;  // tnear = s2 dot e2 / (s1 dot e1)

    // TODO find ray triangle intersection

    inter.happened = true;
    inter.coords = ray(t_tmp);
    inter.normal = normal;
    inter.distance = t_tmp;
    inter.obj = this;
    inter.m = this->m;

    return inter;
}

剩下的就是填充返回值……可以看球与光线相交的函数是怎么做的

Bounds3::IntersectP()

粗略地看了一下这个包围盒的头文件，都是一些简单的求交函数

不知道怎么用这个 dirIsNeg？

inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,
                                const std::array<int, 3>& dirIsNeg) const
{
    // invDir: ray direction(x,y,z), invDir=(1.0/x,1.0/y,1.0/z), use this because Multiply is faster that Division
    // dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x>0),int(y>0),int(z>0)], use this to simplify your logic
    // TODO test if ray bound intersects
    
    Vector3f t_min = (pMin - ray.origin) * invDir;
    Vector3f t_max = (pMax - ray.origin) * invDir;

    float t_enter = std::max(std::max(t_min.x, t_min.y), t_min.z);
    float t_exist = std::min(std::min(t_max.x, t_max.y), t_max.z);

    return (t_enter < t_exist) && (t_exist > 0);
}

无法解析的外部符号

严重性代码说明项目文件行禁止显示状态
错误 LNK2001 无法解析的外部符号 “public: double & __cdecl Vector3f::operator” (??AVector3f@@QEAAAEANH@Z) Games101Homework D:\Work\VisualStudioProject\Games101Homework\Games101Homework\BVH.obj 1

我去搜到了别人的讲解这种错误的文章

https://blog.csdn.net/pongba/article/details/19130

他讲的是模板定义了之后，如果某个 h 头文件或 cpp 源文件没有对这个模板实例化，编译出的目标文件 obj 文件中就不会有这个模板示例的二进制代码，也就不会有这个模板示例的符号，连接器在符号表中查不到这个模板示例的符号，就会报错

实际上我是发现我对这个 Vector3f 使用了 []……实际上直接 .x .y .z 就行了

Debug

后来我发现我这个写错了，求不出正确的结果

但是别人的我一开始觉得复杂就没有看……

inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,
                                const std::array<int, 3>& dirIsNeg) const
{
    // invDir: ray direction(x,y,z), invDir=(1.0/x,1.0/y,1.0/z), use this because Multiply is faster that Division
    // dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x>0),int(y>0),int(z>0)], use this to simplify your logic
    // TODO test if ray bound intersects
    //注：dirIsNeg我觉得没必要用，我感觉我这种逻辑更容易理解，更小的t肯定是min
    float t1,t2,t_min_x, t_max_x, t_min_y, t_max_y, t_min_z, t_max_z, t_enter, t_exit;
    t1 = (pMin.x - ray.origin.x) * invDir.x; t2 = (pMax.x - ray.origin.x) * invDir.x;
    t_min_x = fminf(t1, t2); t_max_x = fmaxf(t1, t2);
    t1 = (pMin.y - ray.origin.y) * invDir.y; t2 = (pMax.y - ray.origin.y) * invDir.y;
    t_min_y = fminf(t1, t2); t_max_y = fmaxf(t1, t2);
    t1 = (pMin.z - ray.origin.z) * invDir.z; t2 = (pMax.z - ray.origin.z) * invDir.z;
    t_min_z = fminf(t1, t2); t_max_z = fmaxf(t1, t2);

    t_exit = fminf(fminf(t_max_x, t_max_y), t_max_z);
    if (t_exit < 0) return false;
    t_enter = fmaxf(fmaxf(t_min_x, t_min_y), t_min_z);
    if (t_enter >= t_exit)return false;
    return true;
}

后面才发现，确实，他需要在对 x y z 单独找 min 和 max

BVHAccel::getIntersection()

仔细看他这个 BVH，它是以 Object 划分的，而 Object 又是球和三角形的基类

所以其实一个三角面就是一个物体

所以其实一个三角面组成的物体的 BVH 是以三角面为基本单位的，而不是我之前理解的，BVH 是以每一个三角面组成的物体构成的

BVH 的构建

BVHBuildNode* BVHAccel::recursiveBuild(std::vector<Object*> objects)
{
    BVHBuildNode* node = new BVHBuildNode();

    // Compute bounds of all primitives in BVH node
    Bounds3 bounds;
    for (int i = 0; i < objects.size(); ++i)
        bounds = Union(bounds, objects[i]->getBounds());
    if (objects.size() == 1) {
        // Create leaf _BVHBuildNode_
        node->bounds = objects[0]->getBounds();
        node->object = objects[0];
        node->left = nullptr;
        node->right = nullptr;
        return node;
    }
    else if (objects.size() == 2) {
        node->left = recursiveBuild(std::vector{objects[0]});
        node->right = recursiveBuild(std::vector{objects[1]});

        node->bounds = Union(node->left->bounds, node->right->bounds);
        return node;
    }
    else {
        Bounds3 centroidBounds;
        for (int i = 0; i < objects.size(); ++i)
            centroidBounds =
                Union(centroidBounds, objects[i]->getBounds().Centroid());
        int dim = centroidBounds.maxExtent();
        // 按照 x 或 y 或 z 来排序，排完之后就方便找到中间的那个
        switch (dim) {
        case 0:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().x <
                       f2->getBounds().Centroid().x;
            });
            break;
        case 1:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().y <
                       f2->getBounds().Centroid().y;
            });
            break;
        case 2:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().z <
                       f2->getBounds().Centroid().z;
            });
            break;
        }

        auto beginning = objects.begin();
        auto middling = objects.begin() + (objects.size() / 2);
        auto ending = objects.end();

        auto leftshapes = std::vector<Object*>(beginning, middling);
        auto rightshapes = std::vector<Object*>(middling, ending);

        assert(objects.size() == (leftshapes.size() + rightshapes.size()));

        node->left = recursiveBuild(leftshapes);
        node->right = recursiveBuild(rightshapes);

        node->bounds = Union(node->left->bounds, node->right->bounds);
    }

    return node;
}

他这里使用递归建树，在输入的物体列表小于等于两个的时候创建叶子节点，否则根据物体列表在空间中的分布来判断怎么划分物体列表

这里获得到物体列表在空间中的分布

        int dim = centroidBounds.maxExtent();
        // 按照 x 或 y 或 z 来排序，排完之后就方便找到中间的那个
        switch (dim) {

这里具体判断空间分布

    int maxExtent() const
    {
        Vector3f d = Diagonal();
        // 在 x 方向上更长
        if (d.x > d.y && d.x > d.z)
            return 0;
        // 在 y 方向上更长
        else if (d.y > d.z)
            return 1;
        // 在 z 方向上更长
        else
            return 2;
    }

这种方式很聪明啊……可以发现，如果不满足第一个条件，那么第二条件自然就代表了 d.y >= d.x && d.y > d.z，同理第三个条件也是

他在递归的最后对 bound 属性赋值，可以看出这个才是最关键的部分，就是说，一个 node 本身不论是否保存物体，其实 BVH 是使用这个 bound 属性来判断求交的

然后抄了别人的

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    // TODO Traverse the BVH to find intersection
    Intersection isect;

    std::array<int, 3> dirIsNeg = { (ray.direction.x > 0), (ray.direction.y > 0), (ray.direction.z > 0) };
    isect.happened = node->bounds.IntersectP(ray, ray.direction_inv, dirIsNeg);

    if (!isect.happened) return isect;

    if (node->left == nullptr && node->right == nullptr)
        return node->object->getIntersection(ray);

    Intersection isect1, isect2;
    isect1 = getIntersection(node->left, ray);
    isect2 = getIntersection(node->right, ray);

    if (isect1.distance < isect2.distance)
        return isect1;
    else
        return isect2;
}

如果当前节点是叶子节点，就执行光线对物体求交

否则选择出距离最近的那一侧的孩子节点递归

SAH

参考了：

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

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

https://www.zhihu.com/question/364497131

BVH 一般的划分方法

之前我们使用的划分方法是找到包围盒在哪个维度上更长，在这个最长的维度上对物体排序，使用中间的物体划分

也就是先找到总包围盒最长的维度，然后在这个维度上等量划分

但是等量划分不一定是最优的

还有一种方法是，确定了划分的维度之后，在这个维度上使用中点划分，也就是求出各个物体在这个维度上的 min max 然后取中点，按照中点的值划分

但是中点划分也不一定是最优的

左：中点划分或等量划分；右：比较理想的划分

启发式的划分方法

一种更为常用且效果更好的方法是启发式的方法，这种方法通过对求交代价和遍历代价进行评估，给出了每一种划分的代价（Cost），而我们的目的便是去寻找代价最小的划分

for(第 i 种划分情况 : n 种划分情况){
	cur_cost = 第 i 种划分情况的 cost;
	if(min_cost > cur_cost){
		min_cost = min_cost;
		state = i;
	}
}
按照第 state 种划分情况来划分物体，得到 node->left, node->right

具体到这个情景中，我们还要确定在哪个轴上划分，而对于启发式的方法，其实也就是猜的方法，他也不知道在哪个轴上划分，所以都要尝试一遍

for(第 i 个轴 : 3 个坐标轴){
	for(第 j 种划分情况 : n 种划分情况){
		cur_cost = 第 j 种划分情况的 cost;
		if(min_cost > cur_cost){
			min_cost = min_cost;
			state = i;
		}
	}
}
按照第 state 种划分情况，在第 i 个轴上划分物体，得到 node->left, node->right;

基于表面积的启发式划分方法 SAH

还剩下的细节就是，具体怎么计算这个 cost

因此有人提出了基于表面积的启发式评估划分方法（Surface Area Heuristic，SAH），这里使用包围盒的表面积计算 cost

它的基本思想是包围盒的表面积越大，越容易与光线相交，所以希望包围盒的表面积尽可能小，所以使用包围盒的表面积计算 cost，相当于惩罚表面积大的项，自然就能得到表面积小的项

使用包围盒的表面积计算 cost

具体来说，BVH 树分为两种结点，一种是叶子节点，一种是非叶节点

叶子节点的求交代价： $\sum C_{intersect}(i)$ 其中 $C_{intersect}(i)$ 是每一个物体的求交代价， $C$ 是总物体的求交代价

假设 $C_{intersect}(i) = C_{intersect}$ 是一个常数的话，那么 $\sum C_{intersect}(i) = NC_{intersect}$ ，其中 $N$ 是叶子节点中物体的数量

虽然其实我不知道为什么要考虑叶子节点的求交代价……感觉问题主要是怎么对非叶节点中的物体划分

非叶节点的求交代价： $C = C_{trav} + p_AC_A + p_BC_B$ ，其中 $C_{trav}$ 是对这个非叶节点本身求交的代价；物体划分到 A, B 桶中， $C_A,C_B$ 是对 A, B 桶中每一个物体求交的代价之和， $p_A,p_B$ 是光线与 A, B 桶相交的概率

假设我们仍然认为每一个物体的求交代价都是常数，那么“包围盒的表面积计算 cost”的具体表现就是计算这个 $p_A,p_B$ ，因为在遍历每一种划分情况的时候，一旦划分情况确定，也就是物体划分到 A, B 桶的数量确定了，也就是 $C_A,C_B$ 确定了，所以自然剩下的问题就是怎么计算 $p_A,p_B$ 了

这里我们使用 $p_A = \dfrac{S_A}{S_N}, p_B = \dfrac{S_B}{S_N}$ ，其中 $S_A, S_B$ 为划分到 A, B 桶的物体的包围盒的表面积， $S_N$ 为所有物体的总包围盒的表面积

更进一步，在比较的时候，我们都是比较非叶节点之间的代价，所以可以不考虑 $C_{trav}$

还可以认为每一个物体的求交代价都是常数 1，那么 $C_A,C_B$ 就等于 A, B 桶中物体的数量

伪代码：

得到所有物体的包围盒 bounds_total;

for(第 i 个轴 : 3 个坐标轴){
	for(第 j 种划分情况 : n 种划分情况){
		确定 A, B 桶中的物体;
		得到 A, B 桶中的物体数量 C_A, C_B, A, B 桶中的物体组成的包围盒 bounds_A, bounds_B;
		p_A = bounds_A.SurfaceArea()/bounds_total.SurfaceArea();
		p_B = bounds_B.SurfaceArea()/bounds_total.SurfaceArea();
		cur_cost = p_A * C_A + p_B * C_B;
		if(min_cost > cur_cost){
			min_cost = min_cost;
			state = i;
		}
	}
}
按照第 state 种划分情况，在第 i 个轴上划分物体，得到 node->left, node->right;

SAH 的划分方法

这里还需要知道 SAH 是怎么划分 A, B 桶的

PPT 里面说是，推荐在空间上划分

那么设图元质心在划分坐标轴上位置坐标的分量为 t，则该图元所处桶的索引号为：

$\mathrm{clamp}(\mathrm{int}(\dfrac{t - t_{min}}{t_{max} - t_{min}}n_{bucket}), 0, n_{bucket}-1)$

其中桶的数量为 $n_{bucket}$ ， $t_{max}$ 为各个图元在 t 坐标轴方向上的最大值， $t_{min}$ 为各个图元在 t 坐标轴方向上的最小值

看了别人的实现，感觉别人没有似乎没有按照空间划分 A, B 桶，而是直接按照物体在 objects 数组中的下标来划分了

BVHBuildNode* BVHAccel::recursiveBuildSAH(std::vector<Object*> objects)
{
    int Bucket_num = 10;
    if (objects.size() < Bucket_num)
        return recursiveBuild(objects);

    int item_num_per_bucket = objects.size() / Bucket_num;

    BVHBuildNode* node = new BVHBuildNode();
    double minCost = std::numeric_limits<double>::infinity();
    int mindim, split_index;

    // get total bound
    Bounds3 bounds;
    for (int i = 0; i < objects.size(); ++i)
        bounds = Union(bounds, objects[i]->getBounds());

    // try 3 dim, find mindim which cost is min 
    for (int dim = 0; dim < 3; dim++)
    {
        // sort according to x or y or z axis
        switch (dim) {
        case 0:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().x <
                    f2->getBounds().Centroid().x;
                });
            break;
        case 1:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().y <
                    f2->getBounds().Centroid().y;
                });
            break;
        case 2:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().z <
                    f2->getBounds().Centroid().z;
                });
            break;
        }

        Bounds3 bucketA;
        int item_num_in_bucketA, item_num_in_bucketB;
        
        item_num_in_bucketA = 0;
        item_num_in_bucketB = 0;
        // bucketA is define outside, so it would keep its value after each loop body
        // so define a bucket_num_have_been_in_A to indicate how many buckets have been in bucketA
        for (int bucket_num_have_been_in_A = 0; bucket_num_have_been_in_A < Bucket_num - 1; bucket_num_have_been_in_A++)
        {
            int i;
            for (i = 0; i < item_num_per_bucket; i++)
            {
                bucketA = Union(bucketA, objects[bucket_num_have_been_in_A * item_num_per_bucket + i]->getBounds());
                item_num_in_bucketA++;
            }
            Bounds3 bucketB;
            for (int j = bucket_num_have_been_in_A * item_num_per_bucket + i; j < objects.size(); j++)
            {
                bucketB = Union(bucketB, objects[j]->getBounds());
            }
            item_num_in_bucketB = objects.size() - item_num_in_bucketA;
            double cost = (bucketA.SurfaceArea() * item_num_in_bucketA + bucketB.SurfaceArea() * item_num_in_bucketB) / bounds.SurfaceArea();
            if (minCost > cost)
            {
                minCost = cost;
                mindim = dim;
                split_index = item_num_in_bucketA;
            }
        }
    }
    
    switch (mindim) {
    case 0:
        std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
            return f1->getBounds().Centroid().x <
                f2->getBounds().Centroid().x;
            });
        break;
    case 1:
        std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
            return f1->getBounds().Centroid().y <
                f2->getBounds().Centroid().y;
            });
        break;
    case 2:
        std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
            return f1->getBounds().Centroid().z <
                f2->getBounds().Centroid().z;
            });
        break;
    }
    auto beginning = objects.begin();
    auto middling = objects.begin() + (split_index);
    auto ending = objects.end();

    auto leftshapes = std::vector<Object*>(beginning, middling);
    auto rightshapes = std::vector<Object*>(middling, ending);

    assert(objects.size() == (leftshapes.size() + rightshapes.size()));

    node->left = recursiveBuild(leftshapes);
    node->right = recursiveBuild(rightshapes);

    node->bounds = Union(node->left->bounds, node->right->bounds);
    return node;
}

又看了别人的一些代码，他们在划分上的部分代码如下

            auto l=(float)objects.size();
            float nums[]={1.0/6,2.0/6,3.0/6,4.0/6,5.0/6};
            for(int i=0;i<5;i++)
                nums[i]*=l;
            for(int i=0;i< 5;i++){
                auto middling = objects.begin() +(int)nums[i];
                auto leftshapes = std::vector<Object*>(beginning, middling);
                auto rightshapes = std::vector<Object*>(middling, ending);

        for(int i=1; i<B_size; ++i) {// search minCost and minCostIdx
            auto middling = objects.begin() + std::max(1, (int)(objects.size() * i / B_size));

            auto leftshapes = std::vector<Object*>(beginning, middling);
            auto rightshapes = std::vector<Object*>(middling, ending);

感觉他们也是按下标比例划分的……？

于是我自己写的是

BVHBuildNode* BVHAccel::recursiveBuildSAH(std::vector<Object*> objects)
{
    int Bucket_num = 10;
    if (objects.size() < Bucket_num)
        return recursiveBuild(objects);

    int item_num_per_bucket = objects.size() / Bucket_num;

    BVHBuildNode* node = new BVHBuildNode();
    double minCost = std::numeric_limits<double>::infinity();
    int mindim, split_index;

    // get total bound
    Bounds3 bounds;
    for (int i = 0; i < objects.size(); ++i)
        bounds = Union(bounds, objects[i]->getBounds());

    // try 3 dim, find mindim which cost is min 
    for (int dim = 0; dim < 3; dim++)
    {
        float t_min = 0, t_max = 0, t_len = 0;

        // sort according to x or y or z axis
        switch (dim) {
        case 0:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().x <
                    f2->getBounds().Centroid().x;
                });
            t_min = (*objects.begin())->getBounds().Centroid().x;
            t_max = (*(objects.end() - 1))->getBounds().Centroid().x;
            break;
        case 1:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().y <
                    f2->getBounds().Centroid().y;
                });
            t_min = (*objects.begin())->getBounds().Centroid().y;
            t_max = (*(objects.end() - 1))->getBounds().Centroid().y;
            break;
        case 2:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().z <
                    f2->getBounds().Centroid().z;
                });
            t_min = (*objects.begin())->getBounds().Centroid().z;
            t_max = (*(objects.end() - 1))->getBounds().Centroid().z;
            break;
        }
        t_len = t_max - t_min;

        std::vector<std::vector<Object*>> buckets(Bucket_num, std::vector<Object*>());
        int idx = 0;
        for (auto obj : objects) {
            switch (dim) {
            case 0:
                idx = clamp(0, Bucket_num - 1, (int)((obj->getBounds().Centroid().x - t_min) / t_len * Bucket_num));
                break;
            case 1:
                idx = clamp(0, Bucket_num - 1, (int)((obj->getBounds().Centroid().y - t_min) / t_len * Bucket_num));
                break;
            case 2:
                idx = clamp(0, Bucket_num - 1, (int)((obj->getBounds().Centroid().z - t_min) / t_len * Bucket_num));
                break;
            }
            buckets[idx].push_back(obj);
        }

        Bounds3 bucketA;
        int item_num_in_bucketA, item_num_in_bucketB;
        
        item_num_in_bucketA = 0;
        item_num_in_bucketB = 0;

        // bucketA is define outside, so it would keep its value after each loop body
        for (int new_bucket_idx_A = 0; new_bucket_idx_A < Bucket_num - 1; new_bucket_idx_A++)
        {
            for (int i = 0; i < buckets[new_bucket_idx_A].size(); i++)
            {
                bucketA = Union(bucketA, buckets[new_bucket_idx_A][i]->getBounds());
            }
            item_num_in_bucketA += buckets[new_bucket_idx_A].size();

            Bounds3 bucketB;
            for (int bucket_idx_B = new_bucket_idx_A + 1; bucket_idx_B < Bucket_num; bucket_idx_B++)
            {
                for (int i = 0; i < buckets[bucket_idx_B].size(); i++)
                {
                    bucketB = Union(bucketB, buckets[bucket_idx_B][i]->getBounds());
                }
            }
            item_num_in_bucketB = objects.size() - item_num_in_bucketA;

            double cost = (bucketA.SurfaceArea() * item_num_in_bucketA + bucketB.SurfaceArea() * item_num_in_bucketB) / bounds.SurfaceArea();
            if (minCost > cost)
            {
                minCost = cost;
                mindim = dim;
                split_index = item_num_in_bucketA;
            }
        }
    }

    switch (mindim) {
    case 0:
        std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
            return f1->getBounds().Centroid().x <
                f2->getBounds().Centroid().x;
            });
        break;
    case 1:
        std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
            return f1->getBounds().Centroid().y <
                f2->getBounds().Centroid().y;
            });
        break;
    case 2:
        std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
            return f1->getBounds().Centroid().z <
                f2->getBounds().Centroid().z;
            });
        break;
    }

    auto beginning = objects.begin();
    auto middling = objects.begin() + (split_index);
    auto ending = objects.end();

    auto leftshapes = std::vector<Object*>(beginning, middling);
    auto rightshapes = std::vector<Object*>(middling, ending);

    assert(objects.size() == (leftshapes.size() + rightshapes.size()));

    node->left = recursiveBuild(leftshapes);
    node->right = recursiveBuild(rightshapes);

    node->bounds = Union(node->left->bounds, node->right->bounds);
    return node;
}