配置

i7-11800h
笔记本版的 RTX3070
CUDA 11.6 vs2019

Baseline

基线为多线程版，我参考的这两个博客：
GAMES101作业7-路径追踪实现过程&代码框架超全解读
GAMES101作业7及课程总结（重点实现多线程加速，微表面模型材质）
在768 $\times$ 768 $\times$ 16(spp)的设置下，8线程release模式渲染时间为16s

CUDA化

CUDA基础就不介绍了，我也只是了解的程度，迁移这个作业代码到CUDA，注意这几个要点就行：

• 所有渲染相关函数，包括类成员函数用__device__修饰，其他初始化的函数，如BVH构建过程在主机端跑就行（__host__修饰）
• 类里面动态分配的内存需要单独用cudaMemcpy拷贝，不能直接拷贝类的对象
• Triangle，MeshTriangle这两个类的初始化要在设备端进行，不能直接cudaMemcpy拷贝到设备端，这两个类重写了虚函数，而设备端和主机端编译出来的虚表地址是不一样的。

__global__ void cpy2device(Object* d_object_ptr_)
{new (d_object_ptr_) MeshTriangle(*(MeshTriangle*)d_object_ptr_);
}
__global__ void create_d_triangles(Triangle* h_triangle_data, Object** d_triangles, int numObjects)
{int index = blockIdx.x * blockDim.x + threadIdx.x;if (index < numObjects) {d_triangles[index] = new Triangle(h_triangle_data[index]);}
}

• 由于设备端不支持使用STL容器，MeshTriangle这个类用于存储Triangle的vector需要动态初始化，即成员变量保留一个vector指针，然后在构造函数中new，否则在设备端调用MeshTriangle构造函数时会调用不支持的std::vector的构造函数

 std::vector<Triangle> *triangles = nullptr;std::vector<Object*> *d_triangles_h_ptr = nullptr;Object** d_triangles;

• 由于我BVH的构建函数recursiveBuild保留了作业本身的版本，也就是在主机端进行的，BVH构建后，BVH叶节点指向的Object对象也是在主机端对象，因此需要一个额外的更新过程来保证BVH加速结构在设备端可用，为实现这一点，还需要将原来的指针型二叉树改为线性存储的二叉树

struct BVHBuildNode {Bounds3 bounds;int leftChildIndex;int rightChildIndex;Object* object;float area;__host__ __device__ BVHBuildNode():bounds(Bounds3()),leftChildIndex(-1), rightChildIndex(-1),object(nullptr), area(0){}
};
__host__ void BVHAccel::updateDeviceNodePointer(BVHBuildNode& node, std::vector<Object*>& d_primitives) {// 如果当前节点是叶节点，则更新 object 指针if (node.leftChildIndex == -1 && node.rightChildIndex == -1) {size_t index = std::distance(primitives.begin(), std::find(primitives.begin(), primitives.end(), node.object));// 更新指针为设备内存中的对象node.object = d_primitives[index];}// 递归更新子节点的指针if (node.leftChildIndex != -1)this->updateDeviceNodePointer(nodes[node.leftChildIndex], d_primitives);if (node.rightChildIndex != -1)this->updateDeviceNodePointer(nodes[node.rightChildIndex], d_primitives);
};

优化

递归->迭代

如果只是像上面那样，保留源代码框架的算法逻辑放到GPU上跑，效率是很低的，我第一次改完代码渲染时间在30s（同样是768 $\times$ 768 $\times$ 16spp的设置）左右，接近CPU端的两倍了。关于这点，如果在性能探查器看原算法的CPU占用，就会发现我们的光追函数castRay中的光线求交和非直接光照的递归占了很多：



CUDA核函数虽然现在支持递归调用，但函数调用需要保存很多临时变量，返回地址，寄存器状态什么的，所以递归对GPU上非常宝贵的寄存器资源很不友好，如果寄存器被迅速耗尽，我们在算法执行中就需要频繁对全局内存进行读写，这通常需要几十到几百个时钟周期，相比只需要几个时钟周期的寄存器读写是非常费时间的

所以第一个优化方向就是递归转迭代，原算法的递归主要有这两个函数Scene::castRay和BVHAccel::getIntersection，先看光追部分Scene::castRay：

Vector3f Scene::castRay(const Ray& ray, int depth) const
{Vector3f accumulatedColor(0.0f);Vector3f throughput(1.0f);Ray currentRay = ray;while (true) {Intersection inter = intersect(currentRay);if (inter.happened) {//if (inter.m->hasEmission()) {//    accumulatedColor += throughput * inter.m->getEmission();//    break;//}accumulatedColor += throughput * inter.m->getEmission();// Direct lightingIntersection lightInter;float pdfLight = 0;sampleLight(lightInter, pdfLight);Vector3f wo = currentRay.direction, x = lightInter.coords, p = inter.coords, N = inter.normal, NN = lightInter.normal;Vector3f emit = lightInter.emit;Vector3f ws = normalize(x - p);float length_sq = length_square(x - p);Ray lightRay(p, ws);Intersection lightTest = intersect(lightRay);if (lightTest.happened && length_square(x - lightTest.coords) <= EPSILON) {accumulatedColor += throughput * emit * inter.m->eval(wo, ws, N) * dotProduct(ws, N) * dotProduct(-ws, NN) / length_sq / pdfLight;}// Indirect lightingif (get_random_float() <= RussianRoulette) {Vector3f wi = normalize(inter.m->sample(wo, N));float pdf = inter.m->pdf(wo, wi, N);throughput *= inter.m->eval(wo, wi, N) * dotProduct(wi, N) / pdf / RussianRoulette;Ray reflectRay(p, wi);Intersection reflectInter = intersect(reflectRay);if (reflectInter.happened && !reflectInter.m->hasEmission()) {currentRay = reflectRay;}else {break;}}else {break;}}else {break;}}return accumulatedColor;
}

原递归的终止条件就是我们迭代的终止条件，然后跟踪三个值：当前参与渲染方程计算的入射光线currentRay，作为返回值的累计颜色量accumulatedColor，以及在非直接光照计算中的当前光强的衰减throughput。由于原算法实现中，只在半球内采样一次而且是尾递归，所以这里我们可以不用栈来额外模拟递归过程。

第二个部分是BVH求交查询BVHAccel::getIntersection，刷过题的可以看出来，这就是实现一个二叉树的非递归先序遍历，注意这里我们保持了之前线性二叉树的设计，左右节点用索引值保存：

    __device__ void BVHAccel::getIntersection(const Ray& ray,BVHBuildNode* &d_nodes_on_chip, Intersection& inter)const{if(!root)return;Intersection curr_inter;Vector3f tmin, tmax, subtractor = ray.origin * ray.direction_inv;float texit, tenter;BVHBuildNode* currentNode, *stack[16];int l_idx, r_idx;Vector3f dirIsNeg(ray.direction.x < 0, ray.direction.y < 0, ray.direction.z < 0);int stackTop = 0;stack[stackTop++] = &d_nodes_on_chip[root_idx];while (stackTop > 0) {currentNode = stack[--stackTop];l_idx = currentNode->leftChildIndex;r_idx = currentNode->rightChildIndex;if (l_idx == -1 && r_idx == -1) {currentNode->object->getIntersection(ray, curr_inter);if (curr_inter.happened && curr_inter.distance < inter.distance) {inter = curr_inter;}}else {if (l_idx != -1) {currentNode = &d_nodes_on_chip[l_idx];tmin = currentNode->bounds.pMin * ray.direction_inv - subtractor;tmax = currentNode->bounds.pMax * ray.direction_inv - subtractor;if (dirIsNeg.x) device_swap(tmin.x, tmax.x);if (dirIsNeg.y) device_swap(tmin.y, tmax.y);if (dirIsNeg.z) device_swap(tmin.z, tmax.z);texit = fmin(tmax.x, fmin(tmax.y, tmax.z));tenter = fmax(tmin.x, fmax(tmin.y, tmin.z));if (tenter <= texit && texit >= 0 && (tenter < 0 || tenter < inter.distance)) {stack[stackTop++] = currentNode;}}if (r_idx != -1) {currentNode = &d_nodes_on_chip[r_idx];tmin = currentNode->bounds.pMin * ray.direction_inv - subtractor;tmax = currentNode->bounds.pMax * ray.direction_inv - subtractor;if (dirIsNeg.x) device_swap(tmin.x, tmax.x);if (dirIsNeg.y) device_swap(tmin.y, tmax.y);if (dirIsNeg.z) device_swap(tmin.z, tmax.z);texit = fmin(tmax.x, fmin(tmax.y, tmax.z));tenter = fmax(tmin.x, fmax(tmin.y, tmin.z));if (tenter <= texit && texit >= 0 && (tenter < 0 || tenter < inter.distance)) {stack[stackTop++] = currentNode;}}}}}

类似之前的，这里需要跟踪的变量就是当前节点curr_inter，然后有几个细节说明下：

• 这里的栈大小我设置成了16，其实从BVH的构建过程可以看出来，它已经非常接近平衡二叉树了，16的深度意味着最多支持30k左右的面片数量或者场景内mesh数量，完全够用
• 我是在节点入栈前判断包围盒是否相交，能减少些迭代次数
• 所有循环体内的变量，能移到循环体外就尽量拿出去，特别是Vector3f，Intersection 这种自定义类，构造函数相比赋值运算还是更昂贵些
• 在判断当前节点的子节点是否入栈的一个剪枝条件tenter < 0 || tenter < inter.distance，理论上应该是有用的，但我们的测试数据比较玩具，都是个位数的面片数量和场景mesh数量，所以我加了之后发现并没有什么性能提升
• 设备端代码中不能用std::array，dirIsNeg用Vector3f保存也一样能工作

这里我们已经将原算法中CPU占用最高的两个函数改为了非递归，算法运行时间可以从30s减少到2.1s，还是有非常大的提升的

线程组织形式

我们是在像素级做并行化，在最初的版本中，我是按一维方式来设置核函数的块大小和网格大小的，这导致一个block或者一个wrap内的线程其实是在图像的水平方向上铺开，然后分别计算，它虽然有一定的空间局部性，但还不充分，因为图像是二维上有局部性的。所以需要像下面这样组织两个大小

dim3 threadsPerBlock(BLOCK_SIZE, BLOCK_SIZE);
dim3 numBlocks(scene.width / threadsPerBlock.x + 1, scene.height / threadsPerBlock.y + 1);
// in kernel...
int xIndex = threadIdx.x + blockIdx.x * blockDim.x;
int yIndex = threadIdx.y + blockIdx.y * blockDim.y;

这样我们一个wrap内的线程在访存时能尽量地收束，在各种条件分支的判断上也不会那么发散，但它的提升相对比较有限，从2.1s降到了1.8s左右

性能分析

到这里其实已经比较完整了，很难想到能有什么优化的地方，只能借性能分析工具来看，对于我们的代码，Nsight Compute给出的报告如下，这里我们只关心执行渲染的核函数RenderPixel：


也是意料之中，算法中执行的逻辑运算其实没多少，反而是一些赋值，函数调用占了大部分，因此报告显示内存的存取比计算执行得更多，也是目前的性能瓶颈。更进一步地，我们看看Memory Workload Analysis这部分：

到这里我突然意识到光追优化起来真的好难。CUDA本身面向过程，对oop支持不那么友好。
然后和其他的需要用GPU的场景对比下，像神经网络的前向传播这种，矩阵的计算都是有非常好的访存局部性的，访存操作在一个wrap内可以合并，大部分关于CUDA的优化教程，尤其是中文社区的资源，也都是以矩阵运算为例子。但是对于光追这种，算是内存密集型，而不是计算密集型的任务，每个线程的内存存取只能在运行时确定，还有半球内随机采样这种会明显导致线程发散的操作，如果保留原有的算法结构的话，访存的优化空间好像很小。

上图可以明显看到两级Cache和全局内存之间进行了大量的内存传输，这里平摊到每个线程，大概渲染一个像素需要在L1和L2之间传输1MB的数据。L1命中率只有60%左右说明了访存确实是现在的性能瓶颈

在wrap的调度数据方面，报告提示GPU一个周期最大调度12个wrap，然而我们的算法目前每个周期平均活跃的数量为3.3，而真正发射出去执行指令的wrap，每个周期只有0.25。那阻塞的wrap在干什么呢

stall long scoreboard是主要因素，就是在等全局内存的数据准备就绪，这说明寄存器是严重不足的，比如在访问最频繁的三角形与光线求交部分，活跃寄存器数量一直处于饱和状态：

到这里就比较清晰了，下一个优化目标就是减少频繁调用的函数的局部变量，之前我关于代码的说明也提到了这点

场景节点放入共享内存

从之前的求交算法可以看到，每次求交需要对场景BVH树的节点进行大量内存访问，而RenderPixel核函数是以场景为单位进行渲染的，因此可以将一个场景的所有BVHBuildNode放入共享内存，我们BVHBuildNode占70个字节左右，共享内存大小一般为32KB，可以放几百个节点，在更复杂的场景下也够用了，因为场景的BVH树是按mesh来组织的。在报告中可以看到共享内存和L1是平级的，应该可以节省一些宝贵的L1缓存：

dim3 threadsPerBlock(BLOCK_SIZE, BLOCK_SIZE);
dim3 numBlocks(scene.width / threadsPerBlock.x + 1, scene.height / threadsPerBlock.y + 1);
int dynamicSharedMemorySize = scene.bvh->node_num * sizeof(BVHBuildNode);
RenderPixel << <numBlocks, threadsPerBlock, dynamicSharedMemorySize >> > (d_framebuffer, scene.d_scene, scene.d_state, eye_pos);
__global__ void RenderPixel(Vector3f* framebuffer, Scene* scene, curandState* states, Vector3f eye_pos)
{int xIndex = threadIdx.x + blockIdx.x * blockDim.x;int yIndex = threadIdx.y + blockIdx.y * blockDim.y;extern  __shared__ BVHBuildNode scene_bvh_nodes[];if (xIndex < s_width && yIndex < s_height){int tid = xIndex + yIndex * s_width;int local_tid = threadIdx.x + blockDim.x * threadIdx.y;if(local_tid < scene->d_bvh->node_num){scene_bvh_nodes[local_tid] = scene->d_bvh->d_nodes[local_tid];}__syncthreads();float x = (2 * (xIndex + 0.5) / (float)s_width - 1) *imageAspectRatio * scale;float y = (1 - 2 * (yIndex + 0.5) / (float)s_height) * scale;for (int k = 0; k < spp; k++) {framebuffer[tid] += scene->castRay(Ray(eye_pos, normalize(Vector3f(-x, y, 1))), 0, states[tid], scene_bvh_nodes) / spp;}}
}

然而很不幸，这一操作并没带来性能提升，只能说明场景BVH的节点访问量，相比于每个mesh的面片BVH节点访问量还是太少。

返回值改为引用传递

这时候回头看原算法的代码框架：

Intersection BVHAccel::Intersect(const Ray& ray) const
{Intersection isect;if (!root)return isect;isect = BVHAccel::getIntersection(root, ray);return isect;
}Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{Intersection res;std::array<int, 3>dirIsNeg = { int(ray.direction.x < 0), int(ray.direction.y < 0), int(ray.direction.z < 0) };if (!node->bounds.IntersectP(ray, ray.direction_inv, dirIsNeg)) {return res;}if (node->left == nullptr && node->right == nullptr) {res = node->object->getIntersection(ray);return res;}Intersection left, right;left = getIntersection(node->left, ray);right = getIntersection(node->right, ray);return left.distance < right.distance ? left : right;
}

这两个需要频繁调用的函数都在开始时初始化了一个Intersection对象，然后作为返回值，而这个对象其实是在函数调用链的末端，也就是Triangle::getIntersection这里才被真正赋值，前面的BVHAccel::getIntersection，BVHAccel::Intersect，MeshTriangle::getIntersection，Scene::intersect都是在进行无效的临时对象构造或者返回值传递。

如果熟悉c++特性的，应该知道，c++17开始强制要求编译器实现返回值优化(return value optimization，RVO)，就是函数内生成的局部变量对象或者结构体作为返回值时，通常情况下会在函数结束时拷贝到一个临时对象，赋值给调用该函数语句处的一个左值，然后销毁该局部的临时对象，而实现RVO时编译器会合并这两个操作。

由于CUDA对c++标准支持是比较滞后的，目前也最多支持到c++17的特性，还需要在编译选项中声明，所以我猜nvcc可能没有做这个事情，google也查不到对nvcc rvo的相关讨论。要知道我们的Intersection对象长这样：

class Intersection
{
public:__device__ __host__ Intersection():happened(false),coords(Vector3f(0.0f)),normal(Vector3f(0.0f)),distance(kInfinity),obj(nullptr),m(nullptr){}bool happened;Vector3f coords;Vector3f tcoords;Vector3f normal;Vector3f emit;double distance;Object* obj;Material* m;
};

不算内存对齐，一个对象也有72个字节，而实际上编译器还会对Vector3f进行内存对齐，虽然它包含3个单精浮点，但却占了24字节，也就是说一个Intersection对象超过100字节，每构造一个这个对象可能就要用掉十几个寄存器~在前面的性能分析报告就能看到，用150多个寄存器的情况下就已经接近饱和状态

所以这里我们的优化就是将光线场景求交的这个调用链上的临时对象+返回值的形式改为引用传递，之前我贴出的非递归BVHAccel::getIntersection是修改后的代码，可以看到其中相比原来的函数，在形参上多了个Intersection&，其他几个的修改类似

    __device__ void Scene::intersect(const Ray& ray, BVHBuildNode* &d_nodes, Intersection &inter) const{inter.happened = false;inter.distance = kInfinity;this->d_bvh->Intersect(ray,//this->d_bvh->d_nodesd_nodes, inter);}

这里Scene::intersect需要交点对象的distance重新初始化，因为在BVHAccel::getIntersection中我们用到了这个值判断当前节点是否更新。

经过这次修改后，再来看访存的情况

直接把L1cache命中率拉高了10个百分点，到了70%，然后L1，L2和全局内存的传输量也几乎减少为原来的一半。这里也可以看到，之前将场景BVH节点放入共享内存后只分担了500MB左右的访存量，相比上面L1cache的10GB还是杯水车薪，没有性能提升也是正常的。

最终版的算法执行时间如下：

总结

这个东西，改出第一版代码拖了两个月，包括踩坑，学习那些，然后debug一周才正常跑通，优化其实只花了3天。一个经验是，专业性或者时效性较强的问题，不要无脑chatgpt，想起改代码那段时间，我问gpt怎么在设备端调用虚函数时，居然说CUDA不支持虚函数，早期的CUDA版本，好像5.0之前确实不支持虚函数，后面就可以了，它的信息在这个问题上严重滞后，然后给我指了三个替代虚函数的方案：模板特化，函数指针和静态成员函数，条件判断。第三个就不说了，前两个在调试的过程中遇到各种各样的问题，都放弃了，直到openai账号被封，回到搜索引擎时，才知道CUDA早就可以运行时多态了。

运行时间来到1s以内的话，终于能展望下实时光追的可能了，相比8线程CPU版的16s加速了20倍，还是比较低于预期的，访存方面应该还有很多优化空间

另外，实验的设置还比较玩具，场景mesh数量，mesh面片数量，spp，分辨率，相对正常值都低很多，纹理也还没做，但增加这些值也会带来更多的优化可能。

比如在场景复杂度这方面，有BVH和合理的剪枝策略，增加面片数量影响就不会太大，我将场景中的两个只有6个面片的盒子，换成了给的4000面片的斯坦福兔子，仅增加了不到0.2s的渲染时间：

纹理方面，GPU有硬件级的纹理缓存和插值加速，额外的纹理、材质计算量也不会成为性能瓶颈；spp可以从优化光追算法本身入手之类的

最后希望所有学图形学的人，都能得到不会辜负已有努力的结果