3D Gaussian Splatting前向渲染简介

3DGS前向渲染流程

3DGS前向渲染流程介绍：
(a)泼溅步骤将 3D 高斯投射到图像空间。
(b)3D高斯将图像划分为多个不重叠的块（tiles）。
(c)3D GS复制覆盖多个块的高斯，为每个副本分配一个标识符 ID。
(d) 通过渲染有序高斯，我们可以获得所有像素的值。渲染过程相互独立。

3DGS前向渲染特点：

• 视锥剔除
• 泼溅（splatting）
• 以像素为单位进行渲染
• 瓦片（图像块）
• 并行化渲染
  3D GS将空间中的3D高斯投影到基于像素的图像平面上，这个过程被称为泼溅（splatting）。随后，3D GS对这些高斯进行排序并计算每个像素的值。

伪代码

伪代码解释：
1、将屏幕划分为16x16的tiles（对于Gaussian点来说就是bins）；
2、计算每个Gaussian点所处的tiles和相对视平面的深度；
3、根据Gaussian点相交的tiles和深度对所有Gaussian点进行排序；
排序方法：GPU Radix sort，每个bins里按Gaussian点深度进行排序；
排序完成后，每个tile都有一个list(bins of Gaussian点)，和这个tile相交的所有Gaussian点在这个list里面从近到远依次存放；
4、给每个tile在GPU里开一个thread block，将tile对应的list加载进block的shared memory里；
5、thread block中的每个thread对应tile中的一个像素，执行α-Blending；
6、计算list里下一个高斯点在当前像素投影出的颜色和α值（很显然这样无法处理两个高斯点相交的情况，所以作者强调了这个α-Blending是approximate的）；
7、将颜色与frame buffer中的颜色混合；
8、将α与透明度buffer中的透明度值相加；
9、如果透明度值大于阈值则退出计算，否则回到步骤1。

代码解读

主要是结合论文和代码进行解读。
论文名称：《3D Gaussian Splatting for Real-Time Radiance Field Rendering》
论文地址：https://repo-sam.inria.fr/fungraph/3d-gaussian-splatting/3d_gaussian_splatting_high.pdf
代码地址：https://github.com/graphdeco-inria/diff-gaussian-rasterization/tree/59f5f77e3ddbac3ed9db93ec2cfe99ed6c5d121d/cuda_rasterizer

栅格化主流程

大致流程：
1、初始化常量和变量(例如焦距、内存、tile size、image相关的变量)；
2、预处理：将3D gaussian点集投影到图像； 
3、生成Idx:为每个projected 2d tile生成idx；
4、生成key-value：为每个projected 2D tile生成key-value;
5、排序：对上一步生成key-value排序(先对tile排序，相同tile id的再按depth排序)；
6、分配range:为每一个tile分配一个range(因为一个tile会对应多个projected 2D tile，那么需要知道2D projected tile id的起始以及终止)。
7、渲染。
备注：projected 2d tile，即是3D gaussian投影到tile 网格上的坐标，多个3Dgaussian可能投影到一个2D tile上，这里projected 2d tile只是指一个3D gaussian的投影tile。

入口：
CudaRasterizer::Rasterizer::forward
代码文件目录：submodules\diff-gaussian-rasterization\cuda_rasterizer\rasterizer_impl.cu

// Forward rendering procedure for differentiable rasterization
// of Gaussians.
int CudaRasterizer::Rasterizer::forward(std::function<char* (size_t)> geometryBuffer,std::function<char* (size_t)> binningBuffer,std::function<char* (size_t)> imageBuffer,//上面的三个参数是用于分配缓冲区的函数。在submodules/diff-gaussian-rasterization/rasterize_points.cu中定义const int P, // Gaussian的数量int D, // 对应于GaussianModel.active_sh_degree，是球谐度数int M, // RGB三通道的球谐傅里叶系数个数，应等于3 × (D + 1)²const float* background,const int width, int height, // 图片宽高const float* means3D, // Gaussians的中心坐标const float* shs, // 球谐系数const float* colors_precomp, // 预先计算的RGB颜色const float* opacities, // 不透明度const float* scales, // 缩放const float scale_modifier, // 缩放的修正项const float* rotations, // 旋转const float* cov3D_precomp, // 预先计算的3维协方差矩阵const float* viewmatrix, // W2C矩阵const float* projmatrix, // 投影矩阵const float* cam_pos, // 相机坐标const float tan_fovx, float tan_fovy, // 视场角一半的正切值const bool prefiltered,float* out_color, // 输出的颜色int* radii, // 各Gaussian在像平面上用3σ原则截取后的半径bool debug)

初始化常量和变量

1、计算焦距；
2、根据3D高斯个数初始化几何相关变量内存；
3、根据固定block size，计算tile size；
4、根据H，W以及tile size初始化image 相关变量；

	/*gaussian_renderer/__init__.py的render函数，定义了tanfovx = math.tan(viewpoint_camera.FoVx * 0.5)tanfovy = math.tan(viewpoint_camera.FoVy * 0.5)*/const float focal_y = height / (2.0f * tan_fovy); // y方向的焦距const float focal_x = width / (2.0f * tan_fovx); // x方向的焦距/*注意tan_fov = tan(fov / 2) 。而tan(fov / 2)就是图像宽/高的一半与焦距之比。以x方向为例，tan(fovx / 2) = width / 2 / focal_x，故focal_x = width / (2 * tan(fovx / 2)) = width / (2 * tan_fovx)。*/// 下面初始化一些缓冲区size_t chunk_size = required<GeometryState>(P); // GeometryState占据空间的大小char* chunkptr = geometryBuffer(chunk_size);GeometryState geomState = GeometryState::fromChunk(chunkptr, P);if (radii == nullptr){radii = geomState.internal_radii;}dim3 tile_grid((width + BLOCK_X - 1) / BLOCK_X, (height + BLOCK_Y - 1) / BLOCK_Y, 1);// BLOCK_X = BLOCK_Y = 16，准备分解成16×16的tiles。// 之所以不能分解成更大的tiles，是因为对于同一张图片的离得较远的像素点而言// Gaussian按深度排序的结果可能是不同的。// （想象一下两个Gaussians离像平面很近，一个靠近图像左边缘，一个靠近右边缘）// dim3是CUDA定义的含义x,y,z三个成员的三维unsigned int向量类。// tile_grid就是x和y方向上tile的个数。dim3 block(BLOCK_X, BLOCK_Y, 1);// Dynamically resize image-based auxiliary buffers during trainingsize_t img_chunk_size = required<ImageState>(width * height);char* img_chunkptr = imageBuffer(img_chunk_size);ImageState imgState = ImageState::fromChunk(img_chunkptr, width * height);if (NUM_CHANNELS != 3 && colors_precomp == nullptr){throw std::runtime_error("For non-RGB, provide precomputed Gaussian colors!");}

预处理

preprocess

// Run preprocessing per-Gaussian (transformation, bounding, conversion of SHs to RGB)CHECK_CUDA(FORWARD::preprocess(P, D, M,means3D,(glm::vec3*)scales,scale_modifier,(glm::vec4*)rotations,opacities,shs,geomState.clamped,cov3D_precomp,colors_precomp,viewmatrix, projmatrix,(glm::vec3*)cam_pos,width, height,focal_x, focal_y,tan_fovx, tan_fovy,radii,geomState.means2D, // Gaussian投影到像平面上的中心坐标geomState.depths, // Gaussian的深度geomState.cov3D, // 三维协方差矩阵geomState.rgb, // 颜色geomState.conic_opacity, // 椭圆二次型的矩阵和不透明度的打包向量tile_grid, // geomState.tiles_touched,prefiltered), debug) // 预处理，主要涉及把3D的Gaussian投影到2D

对于FORWARD::preprocess，详细看预处理章节。

生成Idx

这步是为duplicateWithKeys做准备，计算出每个Gaussian对应的keys和values在数组中存储的起始位置Idx。
使用cub的InclusiveSum实现。

// Compute prefix sum over full list of touched tile counts by Gaussians
// E.g., [2, 3, 0, 2, 1] -> [2, 5, 5, 7, 8]
CHECK_CUDA(cub::DeviceScan::InclusiveSum(geomState.scanning_space, geomState.scan_size, geomState.tiles_touched, geomState.point_offsets, P), debug)

InclusiveSum
计算数组前缀和。所谓"Inclusive"就是第i个数被计入第i个和中.

template<typename InputIteratorT, typename OutputIteratorT>
static inline cudaError_t InclusiveSum(void *d_temp_storage, // 额外需要的临时显存空间size_t &temp_storage_bytes, // 临时显存空间的大小InputIteratorT d_in, // 输入指针OutputIteratorT d_out, // 输出指针int num_items, // 元素个数cudaStream_t stream = 0)

为排序做准备

生成key-value,其中key是 [ tile | depth ]，key是一个uint64_t，前32位表示tile id，后32位表示投影深度；value是3D gaussian的id。

// Retrieve total number of Gaussian instances to launch and resize aux buffersint num_rendered;CHECK_CUDA(cudaMemcpy(&num_rendered, geomState.point_offsets + P - 1, sizeof(int), cudaMemcpyDeviceToHost), debug); // 东西塞到GPU里面去size_t binning_chunk_size = required<BinningState>(num_rendered);char* binning_chunkptr = binningBuffer(binning_chunk_size);BinningState binningState = BinningState::fromChunk(binning_chunkptr, num_rendered);

// For each instance to be rendered, produce adequate [ tile | depth ] key // and corresponding dublicated Gaussian indices to be sortedduplicateWithKeys << <(P + 255) / 256, 256 >> > (P,geomState.means2D,geomState.depths,geomState.point_offsets,binningState.point_list_keys_unsorted,binningState.point_list_unsorted,radii,tile_grid) // 生成排序所用的keys和valuesCHECK_CUDA(, debug)

duplicateWithKeys
计算2d高斯椭圆中心点points_xy在2d像素平面上占据的tile的tileID，并将tileID|depth组合成64位的key值，value值为高斯球的编号。

// Generates one key/value pair for all Gaussian / tile overlaps. 
// Run once per Gaussian (1:N mapping).
__global__ void duplicateWithKeys(int P,const float2* points_xy,const float* depths,const uint32_t* offsets, 			//累计的tiles数量的数组uint64_t* gaussian_keys_unsorted,	 //未排序的key（tileID|depth）uint32_t* gaussian_values_unsorted, 	//未排序的valu（depth）int* radii,  	//高斯球的半径dim3 grid)		//block编号的xy两个极大值
{auto idx = cg::this_grid().thread_rank(); // 线程索引，该显线程处理第idx个Gaussianif (idx >= P)return;// Generate no key/value pair for invisible Gaussiansif (radii[idx] > 0){// Find this Gaussian's offset in buffer for writing keys/values.uint32_t off = (idx == 0) ? 0 : offsets[idx - 1];uint2 rect_min, rect_max;getRect(points_xy[idx], radii[idx], rect_min, rect_max, grid);// 因为要给Gaussian覆盖的每个tile生成一个(key, value)对，// 所以先获取它占了哪些tile// For each tile that the bounding rect overlaps, emit a // key/value pair. The key is |  tile ID  |      depth      |,// and the value is the ID of the Gaussian. Sorting the values // with this key yields Gaussian IDs in a list, such that they// are first sorted by tile and then by depth. for (int y = rect_min.y; y < rect_max.y; y++){for (int x = rect_min.x; x < rect_max.x; x++){uint64_t key = y * grid.x + x; // tile的IDkey <<= 32; // 放在高位key |= *((uint32_t*)&depths[idx]); // 低位是深度gaussian_keys_unsorted[off] = key;gaussian_values_unsorted[off] = idx;off++; // 数组中的偏移量}}}
}

查找最高有效位

int bit = getHigherMsb(tile_grid.x * tile_grid.y);

getHigherMsb
查找最高有效位(most significant bit），输入变量n表示tile编号最大值x、y的乘积。

// Helper function to find the next-highest bit of the MSB
// on the CPU.
uint32_t getHigherMsb(uint32_t n)
{uint32_t msb = sizeof(n) * 4; //4*4=16uint32_t step = msb;while (step > 1){step /= 2;       //缩小2倍if (n >> msb)   //右移16位，相当于除以2的16次方msb += step;elsemsb -= step;}if (n >> msb)     //如果n的最高位大于0，则msb+1msb++;return msb;
}

device级别的并行基数排序

// Sort complete list of (duplicated) Gaussian indices by keysCHECK_CUDA(cub::DeviceRadixSort::SortPairs(binningState.list_sorting_space,binningState.sorting_size,binningState.point_list_keys_unsorted, binningState.point_list_keys,binningState.point_list_unsorted, binningState.point_list,num_rendered, 0, 32 + bit), debug)// 进行排序，按keys排序：每个tile对应的Gaussians按深度放在一起；value是Gaussian的ID

SortPairs
该函数根据key将(key, value)对进行升序排序。这是一种稳定排序。

template<typename KeyT, typename ValueT, typename NumItemsT>
static inline cudaError_t SortPairs(void *d_temp_storage, // 排序时用到的临时显存空间size_t &temp_storage_bytes, // 临时显存空间的大小const KeyT *d_keys_in,         KeyT *d_keys_out, // key的输入和输出指针const ValueT *d_values_in,     ValueT *d_values_out, // value的输入和输出指针NumItemsT num_items, // 对多少个条目进行排序int begin_bit = 0, // 低位int end_bit = sizeof(KeyT) * 8, // 高位cudaStream_t stream = 0)// 按照[begin_bit, end_bit)内的位进行排序

排序后处理

CHECK_CUDA(cudaMemset(imgState.ranges, 0, tile_grid.x * tile_grid.y * sizeof(uint2)), debug);

// Identify start and end of per-tile workloads in sorted list
if (num_rendered > 0)identifyTileRanges << <(num_rendered + 255) / 256, 256 >> > (num_rendered,binningState.point_list_keys,imgState.ranges); // 计算每个tile对应排序过的数组中的哪一部分
CHECK_CUDA(, debug)

identifyTileRanges
一个thread处理一个point_list_keys中的tile，总共L个tile；point_list_keys：已经排序过的key列表，tileID从小到大排列（优先），depth从小到大排列；
ranges：每一项存储对应tile的的id范围[0，L-1]，这个id表示的是在point_list_keys中的索引，通过binningState.point_list找到对应高斯球编号。
例如：point_list_keys值如下：tileID：0 0 0 0 1 1 1 2 2 3 4 4…
depth: 1 2 3 4 1 4 5 3 4 2 3 5… ，那么point_list_keys[0]中的tileID即为0，ranges[0].x = 0

// Check keys to see if it is at the start/end of one tile's range in 
// the full sorted list. If yes, write start/end of this tile. 
// Run once per instanced (duplicated) Gaussian ID.
__global__ void identifyTileRanges(int L, // 排序列表中的元素个数uint64_t* point_list_keys, // 排过序的keysuint2* ranges)// ranges[tile_id].x和y表示第tile_id个tile在排过序的列表中的起始和终止地址
{auto idx = cg::this_grid().thread_rank();if (idx >= L)return;// Read tile ID from key. Update start/end of tile range if at limit.uint64_t key = point_list_keys[idx];uint32_t currtile = key >> 32; // 当前tileif (idx == 0)ranges[currtile].x = 0; // 边界条件：tile 0的起始位置else{uint32_t prevtile = point_list_keys[idx - 1] >> 32;if (currtile != prevtile)// 上一个元素和我处于不同的tile，// 那我是上一个tile的终止位置和我所在tile的起始位置{ranges[prevtile].y = idx;ranges[currtile].x = idx;}}if (idx == L - 1)ranges[currtile].y = L; // 边界条件：最后一个tile的终止位置
}

渲染

// Let each tile blend its range of Gaussians independently in parallelconst float* feature_ptr = colors_precomp != nullptr ? colors_precomp : geomState.rgb;CHECK_CUDA(FORWARD::render(tile_grid, block, // block: 每个tile的大小imgState.ranges,binningState.point_list,width, height,geomState.means2D,feature_ptr,geomState.conic_opacity,imgState.accum_alpha,imgState.n_contrib,background,out_color), debug) // 最后，进行渲染return num_rendered;

FORWARD::render详细看渲染章节。

预处理

// Perform initial steps for each Gaussian prior to rasterization.
template<int C>
__global__ void preprocessCUDA(int P,  //高斯分布的点的数量。int D,  //高斯分布的维度。int M, //点云数量。const float* orig_points, //三维坐标。const glm::vec3* scales, //缩放。const float scale_modifier, //缩放调整因子。const glm::vec4* rotations, //旋转。const float* opacities, //透明度。const float* shs, //球谐函数（SH）特征。bool* clamped, //用于记录是否被裁剪。const float* cov3D_precomp, //预计算的三维协方差。const float* colors_precomp, //预计算的颜色。const float* viewmatrix, //视图矩阵。const float* projmatrix, //投影矩阵const glm::vec3* cam_pos, //相机位置。const int W, int H, //输出图像的宽度和高度。const float tan_fovx, float tan_fovy, //水平和垂直方向的焦距切线。const float focal_x, float focal_y, //焦距。int* radii, //输出的半径。float2* points_xy_image, //输出的二维坐标。float* depths, //输出的深度。float* cov3Ds, //输出的三维协方差。float* rgb, // 输出的颜色。float4* conic_opacity, //锥形透明度。const dim3 grid, //CUDA 网格的大小。uint32_t* tiles_touched,bool prefiltered) //是否预过滤。

获取3D高斯点的id，变量初始化

auto idx = cg::this_grid().thread_rank();if (idx >= P)return;// Initialize radius and touched tiles to 0. If this isn't changed,// this Gaussian will not be processed further.// 首先，初始化了一些变量，包括半径（radii）和触及到的瓦片数量（tiles_touched）。radii[idx] = 0;tiles_touched[idx] = 0;

检查3D高斯点是否在视锥体范围内

	// Perform near culling, quit if outside.// 使用 in_frustum 函数进行近裁剪，如果点在视锥体之外，则退出。float3 p_view;if (!in_frustum(idx, orig_points, viewmatrix, projmatrix, prefiltered, p_view))return;

in_frustum
具体实现在auxiliary.h文件中。
代码路径：submodules\diff-gaussian-rasterization\cuda_rasterizer\auxiliary.h

__forceinline__ __device__ bool in_frustum(int idx,const float* orig_points,const float* viewmatrix,const float* projmatrix,bool prefiltered,float3& p_view)
{float3 p_orig = { orig_points[3 * idx], orig_points[3 * idx + 1], orig_points[3 * idx + 2] };// Bring points to screen spacefloat4 p_hom = transformPoint4x4(p_orig, projmatrix);float p_w = 1.0f / (p_hom.w + 0.0000001f);float3 p_proj = { p_hom.x * p_w, p_hom.y * p_w, p_hom.z * p_w };p_view = transformPoint4x3(p_orig, viewmatrix);if (p_view.z <= 0.2f)// || ((p_proj.x < -1.3 || p_proj.x > 1.3 || p_proj.y < -1.3 || p_proj.y > 1.3))){if (prefiltered){printf("Point is filtered although prefiltered is set. This shouldn't happen!");__trap();}return false;}return true;
}

计算高斯中心点的2D投影

// Transform point by projecting// 对原始点进行投影变换，计算其在屏幕上的坐标。float3 p_orig = { orig_points[3 * idx], orig_points[3 * idx + 1], orig_points[3 * idx + 2] };float4 p_hom = transformPoint4x4(p_orig, projmatrix);float p_w = 1.0f / (p_hom.w + 0.0000001f);float3 p_proj = { p_hom.x * p_w, p_hom.y * p_w, p_hom.z * p_w };

transformPoint4x4
具体实现在auxiliary.h文件中

__forceinline__ __device__ float4 transformPoint4x4(const float3& p, const float* matrix)
{float4 transformed = {matrix[0] * p.x + matrix[4] * p.y + matrix[8] * p.z + matrix[12],matrix[1] * p.x + matrix[5] * p.y + matrix[9] * p.z + matrix[13],matrix[2] * p.x + matrix[6] * p.y + matrix[10] * p.z + matrix[14],matrix[3] * p.x + matrix[7] * p.y + matrix[11] * p.z + matrix[15]};return transformed;
}

计算3D协方差

// If 3D covariance matrix is precomputed, use it, otherwise compute// from scaling and rotation parameters. // 根据输入的缩放和旋转参数，计算或使用预计算的3D协方差矩阵。const float* cov3D;if (cov3D_precomp != nullptr){cov3D = cov3D_precomp + idx * 6;}else{computeCov3D(scales[idx], scale_modifier, rotations[idx], cov3Ds + idx * 6);cov3D = cov3Ds + idx * 6;}

computeCov3D

3D协方差的计算公式，对应论文中的公式（6）
3d060e0c6e4d0abfbb330163843cca.png#pic_center" alt="在这里插入图片描述" />
其中，Σ代表协方差矩阵，R为旋转矩阵，S为缩放矩阵，上标T表示转置矩阵。

// Forward method for converting scale and rotation properties of each
// Gaussian to a 3D covariance matrix in world space. Also takes care
// of quaternion normalization.
__device__ void computeCov3D(const glm::vec3 scale, // 表示缩放的三维向量float mod, // 对应gaussian_renderer/__init__.py中的scaling_modifierconst glm::vec4 rot, // 表示旋转的四元数float* cov3D) // 结果：三维协方差矩阵
{// Create scaling matrixglm::mat3 S = glm::mat3(1.0f);S[0][0] = mod * scale.x;S[1][1] = mod * scale.y;S[2][2] = mod * scale.z;// Normalize quaternion to get valid rotationglm::vec4 q = rot;// / glm::length(rot);float r = q.x;float x = q.y;float y = q.z;float z = q.w;// Compute rotation matrix from quaternionglm::mat3 R = glm::mat3(1.f - 2.f * (y * y + z * z), 2.f * (x * y - r * z), 2.f * (x * z + r * y),2.f * (x * y + r * z), 1.f - 2.f * (x * x + z * z), 2.f * (y * z - r * x),2.f * (x * z - r * y), 2.f * (y * z + r * x), 1.f - 2.f * (x * x + y * y));glm::mat3 M = S * R;// Compute 3D world covariance matrix Sigmaglm::mat3 Sigma = glm::transpose(M) * M;// Covariance is symmetric, only store upper rightcov3D[0] = Sigma[0][0];cov3D[1] = Sigma[0][1];cov3D[2] = Sigma[0][2];cov3D[3] = Sigma[1][1];cov3D[4] = Sigma[1][2];cov3D[5] = Sigma[2][2];
}

计算2D协方差（3D协方差在2D的投影）

	// Compute 2D screen-space covariance matrix// 根据3D协方差矩阵、焦距和视锥体矩阵，计算2D屏幕空间的协方差矩阵。float3 cov = computeCov2D(p_orig, focal_x, focal_y, tan_fovx, tan_fovy, cov3D, viewmatrix);

computeCov2D
相机视角下的协方差矩阵，计算公式对应论文中的公式（5）

其中，J为雅可比矩阵，W为视点变换矩阵，Σ代表3D协方差矩阵。

// Forward version of 2D covariance matrix computation
__device__ float3 computeCov2D(const float3& mean, // Gaussian中心坐标float focal_x, // x方向焦距float focal_y, // y方向焦距float tan_fovx,float tan_fovy,const float* cov3D, // 已经算出来的三维协方差矩阵const float* viewmatrix) // W2C矩阵
{// The following models the steps outlined by equations 29// and 31 in "EWA Splatting" (Zwicker et al., 2002). // Additionally considers aspect / scaling of viewport.// Transposes used to account for row-/column-major conventions.float3 t = transformPoint4x3(mean, viewmatrix);// W2C矩阵乘Gaussian中心坐标得其在相机坐标系下的坐标const float limx = 1.3f * tan_fovx;const float limy = 1.3f * tan_fovy;const float txtz = t.x / t.z; // Gaussian中心在像平面上的x坐标const float tytz = t.y / t.z; // Gaussian中心在像平面上的y坐标t.x = min(limx, max(-limx, txtz)) * t.z;t.y = min(limy, max(-limy, tytz)) * t.z;glm::mat3 J = glm::mat3(focal_x / t.z, 0.0f, -(focal_x * t.x) / (t.z * t.z),0.0f, focal_y / t.z, -(focal_y * t.y) / (t.z * t.z),0, 0, 0); // 雅可比矩阵（用泰勒展开近似）glm::mat3 W = glm::mat3( // W2C矩阵viewmatrix[0], viewmatrix[4], viewmatrix[8],viewmatrix[1], viewmatrix[5], viewmatrix[9],viewmatrix[2], viewmatrix[6], viewmatrix[10]);glm::mat3 T = W * J;glm::mat3 Vrk = glm::mat3( // 3D协方差矩阵，是对称阵cov3D[0], cov3D[1], cov3D[2],cov3D[1], cov3D[3], cov3D[4],cov3D[2], cov3D[4], cov3D[5]);glm::mat3 cov = glm::transpose(T) * glm::transpose(Vrk) * T;// transpose(J) @ transpose(W) @ Vrk @ W @ J// Apply low-pass filter: every Gaussian should be at least// one pixel wide/high. Discard 3rd row and column.cov[0][0] += 0.3f;cov[1][1] += 0.3f;return { float(cov[0][0]), float(cov[0][1]), float(cov[1][1]) };// 协方差矩阵是对称的，只用存储上三角，故只返回三个数
}

计算2D协方差的逆(EWA algorithm)

相关公式待补充。

// Invert covariance (EWA algorithm)// 对协方差矩阵进行求逆操作，用于EWA（Elliptical Weighted Average）算法。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 };

计算2D协方差矩阵的特征值（转换到像素坐标系，计算投影半径）

计算2D协方差矩阵的特征值，用于计算屏幕空间的范围，以确定与之相交的瓦片。
高斯投影半径的计算公式，待补充。

// 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. // 计算2D协方差矩阵的特征值，用于计算屏幕空间的范围，以确定与之相交的瓦片。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;getRect(point_image, my_radius, rect_min, rect_max, grid);if ((rect_max.x - rect_min.x) * (rect_max.y - rect_min.y) == 0)return;

getRect
具体实现在auxiliary.h文件中。

__forceinline__ __device__ void getRect(const float2 p, int max_radius, uint2& rect_min, uint2& rect_max, dim3 grid)
{rect_min = {min(grid.x, max((int)0, (int)((p.x - max_radius) / BLOCK_X))),min(grid.y, max((int)0, (int)((p.y - max_radius) / BLOCK_Y)))};rect_max = {min(grid.x, max((int)0, (int)((p.x + max_radius + BLOCK_X - 1) / BLOCK_X))),min(grid.y, max((int)0, (int)((p.y + max_radius + BLOCK_Y - 1) / BLOCK_Y)))};
}

ndc2Pix

__forceinline__ __device__ float ndc2Pix(float v, int S)
{return ((v + 1.0) * S - 1.0) * 0.5;
}

根据高斯球谐系数计算RGB

// If colors have been precomputed, use them, otherwise convert// spherical harmonics coefficients to RGB color.// 如果预计算颜色未提供，则使用球谐函数（SH）系数计算颜色。if (colors_precomp == nullptr){glm::vec3 result = computeColorFromSH(idx, D, M, (glm::vec3*)orig_points, *cam_pos, shs, clamped);rgb[idx * C + 0] = result.x;rgb[idx * C + 1] = result.y;rgb[idx * C + 2] = result.z;}

computeColorFromSH
该函数从球谐系数相机观察每个Gaussian的RGB颜色。

// Forward method for converting the input spherical harmonics
// coefficients of each Gaussian to a simple RGB color.
__device__ glm::vec3 computeColorFromSH(int idx, // 该线程负责第几个Gaussianint deg, // 球谐的度数int max_coeffs, // 一个Gaussian最多有几个傅里叶系数const glm::vec3* means, // Gaussian中心位置glm::vec3 campos, // 相机位置const float* shs, // 球谐系数bool* clamped) // 表示每个值是否被截断了（RGB只能为正数），这个在反向传播的时候用
{// The implementation is loosely based on code for // "Differentiable Point-Based Radiance Fields for // Efficient View Synthesis" by Zhang et al. (2022)glm::vec3 pos = means[idx];glm::vec3 dir = pos - campos;dir = dir / glm::length(dir);	// dir = direction，即观察方向glm::vec3* sh = ((glm::vec3*)shs) + idx * max_coeffs;glm::vec3 result = SH_C0 * sh[0];if (deg > 0){float x = dir.x;float y = dir.y;float z = dir.z;result = result - SH_C1 * y * sh[1] + SH_C1 * z * sh[2] - SH_C1 * x * sh[3];if (deg > 1){float xx = x * x, yy = y * y, zz = z * z;float xy = x * y, yz = y * z, xz = x * z;result = result +SH_C2[0] * xy * sh[4] +SH_C2[1] * yz * sh[5] +SH_C2[2] * (2.0f * zz - xx - yy) * sh[6] +SH_C2[3] * xz * sh[7] +SH_C2[4] * (xx - yy) * sh[8];if (deg > 2){result = result +SH_C3[0] * y * (3.0f * xx - yy) * sh[9] +SH_C3[1] * xy * z * sh[10] +SH_C3[2] * y * (4.0f * zz - xx - yy) * sh[11] +SH_C3[3] * z * (2.0f * zz - 3.0f * xx - 3.0f * yy) * sh[12] +SH_C3[4] * x * (4.0f * zz - xx - yy) * sh[13] +SH_C3[5] * z * (xx - yy) * sh[14] +SH_C3[6] * x * (xx - 3.0f * yy) * sh[15];}}}result += 0.5f;// RGB colors are clamped to positive values. If values are// clamped, we need to keep track of this for the backward pass.clamped[3 * idx + 0] = (result.x < 0);clamped[3 * idx + 1] = (result.y < 0);clamped[3 * idx + 2] = (result.z < 0);return glm::max(result, 0.0f);
}

保存信息

    // 存储计算得到的深度、半径、屏幕坐标等结果，用于下一步继续处理。// 为每个高斯分布进行预处理，为后续的高斯光栅化做好准备。// Store some useful helper data for the next steps.depths[idx] = p_view.z;radii[idx] = my_radius;points_xy_image[idx] = point_image;// Inverse 2D covariance and opacity neatly pack into one float4conic_opacity[idx] = { conic.x, conic.y, conic.z, opacities[idx] };tiles_touched[idx] = (rect_max.y - rect_min.y) * (rect_max.x - rect_min.x);

渲染

renderCUDA的核心逻辑如下：
1、通过计算当前线程所属的 tile 的范围，确定当前线程要处理的像素区域。
2、判断当前线程是否在有效像素范围内，如果不在，则将 done 设置为 true，表示该线程不执行渲染操作。
3、使用 __syncthreads_count 函数，统计当前块内 done 变量为 true 的线程数，如果全部线程都完成，跳出循环。
4、在每个迭代中，从全局内存中收集每个线程块对应的范围内的数据，包括点的索引、2D 坐标和锥体参数透明度。
5、对当前线程块内的每个点，进行基于锥体参数的渲染，计算贡献并更新颜色。
6、所有线程处理完毕后，将渲染结果写入 final_T、n_contrib 和 out_color。

// Main rasterization method. Collaboratively works on one tile per
// block, each thread treats one pixel. Alternates between fetching 
// and rasterizing data.
template <uint32_t CHANNELS>
__global__ void __launch_bounds__(BLOCK_X * BLOCK_Y)// 这是 CUDA 启动核函数时使用的线程格和线程块的数量。
renderCUDA(const uint2* __restrict__ ranges, //包含了每个范围的起始和结束索引的数组。const uint32_t* __restrict__ point_list, //包含了点的索引的数组。int W, int H, //图像的宽度和高度。const float2* __restrict__ points_xy_image, //包含每个点在屏幕上的坐标的数组。const float* __restrict__ features, //包含每个点的颜色信息的数组。const float4* __restrict__ conic_opacity, //包含每个点的锥体参数和透明度信息的数组。float* __restrict__ final_T, //用于存储每个像素的最终颜色的数组。uint32_t* __restrict__ n_contrib, //用于存储每个像素的贡献计数的数组。const float* __restrict__ bg_color, //如果提供了背景颜色，将其作为背景。float* __restrict__ out_color) //存储最终渲染结果的数组。

确定当前像素范围

// 这部分代码用于确定当前线程块要处理的像素范围，包括 pix_min 和 pix_max，并计算当前线程对应的像素坐标 pix。// Identify current tile and associated min/max pixel range.auto block = cg::this_thread_block();uint32_t horizontal_blocks = (W + BLOCK_X - 1) / BLOCK_X;uint2 pix_min = { block.group_index().x * BLOCK_X, block.group_index().y * BLOCK_Y };uint2 pix_max = { min(pix_min.x + BLOCK_X, W), min(pix_min.y + BLOCK_Y , H) };uint2 pix = { pix_min.x + block.thread_index().x, pix_min.y + block.thread_index().y };uint32_t pix_id = W * pix.y + pix.x;float2 pixf = { (float)pix.x, (float)pix.y };

判断当前线程是否在有效像素范围内

// 根据像素坐标判断当前线程是否在有效的图像范围内，如果不在，则将 done 设置为 true，表示该线程无需执行渲染操作。// Check if this thread is associated with a valid pixel or outside.bool inside = pix.x < W&& pix.y < H;// Done threads can help with fetching, but don't rasterizebool done = !inside;

加载点云数据处理范围

// 这部分代码加载当前线程块要处理的点云数据的范围，即 ranges 数组中对应的范围，并计算点云数据的迭代批次 rounds 和总共要处理的点数 toDo。// Load start/end range of IDs to process in bit sorted list.uint2 range = ranges[block.group_index().y * horizontal_blocks + block.group_index().x];const int rounds = ((range.y - range.x + BLOCK_SIZE - 1) / BLOCK_SIZE);int toDo = range.y - range.x;

初始化共享内存

// 分别定义三个共享内存数组，用于在每个线程块内共享数据。
// Allocate storage for batches of collectively fetched data.
__shared__ int collected_id[BLOCK_SIZE];
__shared__ float2 collected_xy[BLOCK_SIZE];
__shared__ float4 collected_conic_opacity[BLOCK_SIZE];

初始化渲染相关变量

// 初始化渲染所需的一些变量，包括当前像素颜色 C、贡献者数量等。// Initialize helper variablesfloat T = 1.0f;uint32_t contributor = 0;uint32_t last_contributor = 0;float C[CHANNELS] = { 0 };

迭代处理点云数据

在每个迭代中，处理一批点云数据。内部循环迭代每个点，进行基于锥体参数的渲染计算，并更新颜色信息。

// Iterate over batches until all done or range is completefor (int i = 0; i < rounds; i++, toDo -= BLOCK_SIZE) //代码使用 rounds 控制循环的迭代次数，每次迭代处理一批点云数据。{	// 检查是否所有线程块都已经完成渲染：// 通过 __syncthreads_count 统计已经完成渲染的线程数，如果整个线程块都已完成，则跳出循环。// End if entire block votes that it is done rasterizingint num_done = __syncthreads_count(done);if (num_done == BLOCK_SIZE)break;// 共享内存中获取点云数据：// 每个线程通过索引 progress 计算要加载的点云数据的索引 coll_id，然后从全局内存中加载到共享内存 collected_id、collected_xy 和 collected_conic_opacity 中。block.sync() 确保所有线程都加载完成。// Collectively fetch per-Gaussian data from global to sharedint 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();

以下内容涉及论文中的公式（2）和（3）
公式（2）


3dfa72ce07695eba655.png#pic_center" alt="在这里插入图片描述" />

公式（3）

		// 迭代处理当前批次的点云数据：// Iterate over current batchfor (int j = 0; !done && j < min(BLOCK_SIZE, toDo); j++) //在当前批次的循环中，每个线程处理一条点云数据。{// Keep track of current position in rangecontributor++;// 计算当前点的投影坐标与锥体参数的差值：// 计算当前点在屏幕上的坐标 xy 与当前像素坐标 pixf 的差值，并使用锥体参数计算 power。// 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;// 计算论文中公式2的 alpha：// 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;}// 使用高斯分布进行渲染计算：更新颜色信息 C。// 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;}}

写入最终渲染结果

// 如果当前线程在有效像素范围内，则将最终的渲染结果写入相应的缓冲区，包括 final_T、n_contrib 和 out_color。// 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];}

参考资料：
1、https://iks-ran.me/2023/10/17/3d_gaussian_splatting/
2、https://github.com/graphdeco-inria/gaussian-splatting?tab=readme-ov-file
3、https://github.com/graphdeco-inria/diff-gaussian-rasterization/tree/main/cuda_rasterizer
4、其他参考资料待补充。