--- title: Untitled date: 2024-01-01 18:57:57 excerpt: tags: rating: ⭐ --- # 前言 - https://github.com/graphdeco-inria/gaussian-splatting/tree/main/gaussian_renderer 基于Sibr渲染器制作的3D高斯查看器。 # 项目结构 - [x] gaussian - render - sibr_gaussian - apps - SIBR_gaussianViewer_app - [x] diff-gaussian-rasterization(CUDA) # render - sibr_gaussian - picojson：JSON库 - rapidxml：XML库 - **nanoflann**：是一个c++11标准库，用于构建具有不同拓扑（R2，R3（点云），SO(2)和SO(3)（2D和3D旋转组））的KD树。 ## GaussianSurfaceRenderer >主要用于渲染椭圆体，估计是用于Debug用的。 ### GaussianData - GaussianData()：通过构造函数形参接受CPU端读取的高斯数据，再通过调用glCreateBuffers()、glNamedBufferStorage()创建GL缓存对象并且初始化，并使用GLuint进行记录（index）。 - render：给Shader绑定GL缓存，并且绘制数组实例。 ### GaussianSurfaceRenderer - GaussianSurfaceRenderer()：初始化相关变量。 - 初始化VS/Frag Shader。 - rayOrigin、MVP、alpha_limit、stage变量 - 创建idTexture、colorTexture贴图变量以及过滤器 - 创建fbo对象以及depthBuffer之后调用makeFBO()正式创建FBO - 创建清屏Shader。 - makeFBX()：创建idTexture、colorTexture、depthBuffer FBO，用于将顶点数据传递到FragShader中。 - process()：整个渲染过程逻辑处理。 1. 清屏。 2. 判断如果分辨率与FBO大小不同，则重新创建FBO。 3. 获取绘制Buffer的Index，调用glDrawBuffers() 绘制colorTexture、idTexture。 4. 开启深度测试关闭Blend模式。 5. 给Shader绑定相关`_paramMVP`、`_paramCamPos`、`_paramLimit`、`_paramStage`变量，并且调用GaussianData.render()进行一次**不透明物体**的渲染。以小方盒的形式绘制点云数据。 6. 调用glDrawBuffers() 绘制colorTexture。 7. 关闭深度测试，开启透明Blend模式。 8. GaussianData.render()进行一次**透明物体**的渲染，融合模式**additive blendnig**。以小方盒的形式绘制点云数据。 9. 开启深度测试，关闭Blend模式。 10. 将结果显示在屏幕上？ ## GaussianView 继承自sibr::ViewBase，用与调用渲染器以及显示结果。 ### GaussianView - GaussianView()： - 初始化_pointbasedrenderer渲染器 - 初始化_copyRenderer渲染器 - 载入图片并且加入debug模式（应该sibr自带的那个多视角图片debug模式） - 载入*.ply点云文件，函数为loadPly()。 - CUDA相关处理，应该是为了计算3D高斯结果所需的数据。 - 生成GaussianData指针变量gData。 - 初始化3D高斯渲染器对象_gaussianRenderer。 - 创建GL缓存对象imageBuffer。 - CUDA插值操作。 - 绑定3个geomBufferFunc、binningBufferFunc、imgBufferFunc仿函数，用来调整CUDA渲染时的缓存大小（创建或者回收内存空间） - onRenderIBR()：View的渲染函数。 - Ellipsoids（椭圆体渲染）：使用_gaussianRenderer->process() 进行渲染。(OpenGL) - Initial Points：`_pointbasedrenderer->process()`渲染点。 - Splats：使用CudaRasterizer::Rasterizer::forward()进行渲染。最后通过_copyRenderer->process()复制回imageBuffer缓存。 - onGUI()：GUI相关逻辑。 CUDA文件位于`SIBR_viewers\extlibs\CudaRasterizer\CudaRasterizer\cuda_rasterizer\rasterizer_impl.cu`以及`forward.cu`，这些为核心逻辑。 ## Shader 可以理解为将点云渲染成一个个的椭圆体，每个椭圆体的颜色与点云数据中的颜色相关。 ### VertexShader 1. 取得IndexID。 2. 使用IndexID从传入Shader的Buffer中获取的椭圆体中心、alpha、ellipsoidScale、q（四元数rotation），之后将rotation转成3x3矩阵 ellipsoidRotation。 3. 取得当前顶点Index并获得坐标。再乘以椭圆体旋转值并加上椭圆体中心坐标，取得最终的WorldPos（当前顶点的世界坐标）。 4. 使用IndexID从传入Shader的Buffer中取得**辐射照度?辐射强度?** 数据。 5. 将不符合要求的顶点堆到vec4(0,0,0,0)点。 6. 输出顶点数据到FragShader。 ### FragShader 1. 计算摄像机=>当前顶点世界坐标的方向向量dir。 2. 调用closestEllipsoidIntersection()，计算与椭圆体的相交的坐标与相交点的法线。 1. 计算椭圆体空间的localRayOrigin与localRayDirection 2. 计算椭圆与直线相交的方程。 3. 计算摄像机朝向的椭圆体的外表面。如果是内表面最终颜色值 * 0.4。 4. 将相交的世界坐标乘以MVP矩阵，得到摄像机View坐标下的的世界坐标。 5. 计算深度缓存。 6. 计算Alpha。 7. 渲染`out_color = vec4(align * colorVert, a);` 也就是colorTexture 8. 渲染`out_id = boxID;`也就是idTexture # CudaRasterizer **本人没学过CUDA，以下仅仅是对代码的猜测。** 额外需要了解Tile渲染方式（具体可以看**Tiled-Based Deferred Rendering(TBDR)**) https://zhuanlan.zhihu.com/p/547943994 - 屏幕分成`16 * 16`的tile，每个tile进行单独计算。之后对每个像素进行计算。 - 取得对应tile中Start与End的位置，对已经排序完的高斯点进行计算，求微分。 - 计算当前像素的透明度T - 2D协方差 => power => alpha。 - 每次循环都进行`float test_T = T * (1 - alpha)`，当test_T极小时（不透明）则停止循环。 - T = test_T。 - 计算当前像素的颜色，也就是计算各个方向接受的辐射照度。 - `for (int ch = 0; ch < CHANNELS; ch++)` `C[ch] += features[collected_id[j] * CHANNELS + ch] * alpha * T;` - 计算最终贡献值 - 如果当前像素在范围中则输出 - `final_T[pix_id]`最终透明度。 - `n_contrib[pix_id]`最终贡献值。 - `out_color[ch * H * W + pix_id]`最终颜色。`C[ch] + T * bg_color[ch]` 对屏幕分Tile ![[ScreenSpaceTile.jpg]] 以此减少需要遍历的点云数量。 ![[TileRange.jpg|500]] 每个点云相当于空间中当前位置空间的辐射强度分布。 ![[GS_radiation.jpg]] 一个像素的渲染会计算这个像素范围内所有的点云的辐射强度、透明度，最后求微分。下图两条横线内相当于一个像素的范围。 ![[一个像素需要计算范围内所有电源的辐射强度.png|500]] ## rasterizer_impl.cu - getHigherMsb() - checkFrustum()：判断点云是否在视锥内，返回一个bool数组。 - duplicateWithKeys() - identifyTileRanges()：确定每个Tile的工作起点与终点。 - markVisible()：标记高斯点云是否处于可视状态。 - GeometryState::fromChunk()：计算数据块的指针偏移，并且返回创建的GeometryState结构体对象。 - ImageState::fromChunk()：计算数据块的指针偏移，并且返回创建的ImageState结构体对象。 - BinningState::fromChunk()：计算数据块的指针偏移，并且返回创建的BinningState结构体对象。 - forward()：前向渲染可微分光栅化的高斯。具体见下文。 - backward()：生成优化所需的梯度数据，并传递到forward()。**该项目中目前未被调用** 相关数据结构体定义在rasterizer_impl.h中： ```c++ struct GeometryState { size_t scan_size; float* depths; char* scanning_space; bool* clamped; int* internal_radii; float2* means2D; float* cov3D; float4* conic_opacity; float* rgb; uint32_t* point_offsets; uint32_t* tiles_touched; static GeometryState fromChunk(char*& chunk, size_t P); }; struct ImageState { uint2* ranges; uint32_t* n_contrib; float* accum_alpha; static ImageState fromChunk(char*& chunk, size_t N); }; struct BinningState { size_t sorting_size; uint64_t* point_list_keys_unsorted; uint64_t* point_list_keys; uint32_t* point_list_unsorted; uint32_t* point_list; char* list_sorting_space; static BinningState fromChunk(char*& chunk, size_t P); }; ``` ### forward() 1. 创建相关变量：GeometryState、ImageState、minn、maxx。 2. FORWARD::preprocess() 3. 计算所有tile的高斯点云总量。 4. 根据需要需要渲染的高斯点云总量来调整CUDA buffer大小。 5. 创建BinningState。 6. duplicateWithKeys() 7. getHigherMsb() 8. 对高斯点运行排序。 9. cudaMemset(imgState.ranges, 0, tile_grid.x * tile_grid.y * sizeof(uint2)); 10. 调用identifyTileRanges()，确定每个Tile的工作起点与终点。 11. 取得点云颜色数组。 12. FORWARD::render() ## forward.cu ### preprocess() 在光栅化之前，对每个高斯进行初始化处理。 - 只处理在视锥中并且在盒子中的高斯。 - 使用投影矩阵对点云的点进行变换，并进行归一化，赋予给新变量p_proj。 - 计算协方差矩阵cov3D。 - 计算2D屏幕空间的协方差矩阵cov - Invert covariance - Compute extent in screen space (by finding eigenvalues of 2D covariance matrix). Use extent to compute a bounding rectangle of screen-space tiles that this Gaussian overlaps with. Quit if rectangle covers 0 tiles. - 如果没有颜色数据则从球谐函数中计算辐射照度。 - 存储当前数据。 - `depths[idx]` - `radii[idx]` - `points_xy_image[idx]` - `conic_opacity[idx]` - `tiles_touched[idx]` ```c++ // Invert covariance (EWA algorithm) float det = (cov.x * cov.z - cov.y * cov.y); if (det == 0.0f) return; float det_inv = 1.f / det; float3 conic = { cov.z * det_inv, -cov.y * det_inv, cov.x * det_inv }; // Compute extent in screen space (by finding eigenvalues of // 2D covariance matrix). Use extent to compute a bounding rectangle // of screen-space tiles that this Gaussian overlaps with. Quit if // rectangle covers 0 tiles. float mid = 0.5f * (cov.x + cov.z); float lambda1 = mid + sqrt(max(0.1f, mid * mid - det)); float lambda2 = mid - sqrt(max(0.1f, mid * mid - det)); float my_radius = ceil(3.f * sqrt(max(lambda1, lambda2))); float2 point_image = { ndc2Pix(p_proj.x, W), ndc2Pix(p_proj.y, H) }; uint2 rect_min, rect_max; if (rects == nullptr) // More conservative { getRect(point_image, my_radius, rect_min, rect_max, grid); } else // Slightly more aggressive, might need a math cleanup { const int2 my_rect = { (int)ceil(3.f * sqrt(cov.x)), (int)ceil(3.f * sqrt(cov.z)) }; rects[idx] = my_rect; getRect(point_image, my_rect, rect_min, rect_max, grid); } if ((rect_max.x - rect_min.x) * (rect_max.y - rect_min.y) == 0) return; ``` ### render() 对所有Tile进行并行计算。针对CUDA核心数量创建对应的Block以及对应数据。`int collected_id[BLOCK_SIZE]、float2 collected_xy[BLOCK_SIZE]、float4 collected_conic_opacity[BLOCK_SIZE]`。 递归所有的Block，计算透明度、Color以及贡献值（用于计算平均值）。 ```c++ // Iterate over batches until all done or range is complete for (int i = 0; i < rounds; i++, toDo -= BLOCK_SIZE) { // End if entire block votes that it is done rasterizing int num_done = __syncthreads_count(done); if (num_done == BLOCK_SIZE) break; // Collectively fetch per-Gaussian data from global to shared int progress = i * BLOCK_SIZE + block.thread_rank(); if (range.x + progress < range.y) { int coll_id = point_list[range.x + progress]; collected_id[block.thread_rank()] = coll_id; collected_xy[block.thread_rank()] = points_xy_image[coll_id]; collected_conic_opacity[block.thread_rank()] = conic_opacity[coll_id]; } block.sync(); // Iterate over current batch for (int j = 0; !done && j < min(BLOCK_SIZE, toDo); j++) { // Keep track of current position in range contributor++; // Resample using conic matrix (cf. "Surface // Splatting" by Zwicker et al., 2001) float2 xy = collected_xy[j]; float2 d = { xy.x - pixf.x, xy.y - pixf.y }; float4 con_o = collected_conic_opacity[j]; float power = -0.5f * (con_o.x * d.x * d.x + con_o.z * d.y * d.y) - con_o.y * d.x * d.y; if (power > 0.0f) continue; // Eq. (2) from 3D Gaussian splatting paper. // Obtain alpha by multiplying with Gaussian opacity // and its exponential falloff from mean. // Avoid numerical instabilities (see paper appendix).float alpha = min(0.99f, con_o.w * exp(power)); if (alpha < 1.0f / 255.0f) continue; float test_T = T * (1 - alpha); if (test_T < 0.0001f) { done = true; continue; } // Eq. (3) from 3D Gaussian splatting paper. for (int ch = 0; ch < CHANNELS; ch++) C[ch] += features[collected_id[j] * CHANNELS + ch] * alpha * T; T = test_T; // Keep track of last range entry to update this // pixel. last_contributor = contributor; }} ``` ```c++ // All threads that treat valid pixel write out their final // rendering data to the frame and auxiliary buffers. if (inside) { final_T[pix_id] = T; n_contrib[pix_id] = last_contributor; for (int ch = 0; ch < CHANNELS; ch++) out_color[ch * H * W + pix_id] = C[ch] + T * bg_color[ch]; } ``` # apps - SIBR_gaussianViewer_app 调用`gaussianviewer/renderer/GaussianView.hpp`封装的App。