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 training
	size_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 buffers
	int 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 sorted
	duplicateWithKeys << <(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和values
	CHECK_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个Gaussian
	if (idx >= P)
		return;

	// Generate no key/value pair for invisible Gaussians
	if (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的ID
				key <<= 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=16
	uint32_t step = msb;
	while (step > 1)
	{
		step /= 2;       //缩小2倍
		if (n >> msb)   //右移16位，相当于除以2的16次方
			msb += step;
		else
			msb -= step;
	}
	if (n >> msb)     //如果n的最高位大于0，则msb+1
		msb++;
	return msb;
}

device级别的并行基数排序

// Sort complete list of (duplicated) Gaussian indices by keys
	CHECK_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, // 排过序的keys
	uint2* 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; // 当前tile
	if (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 parallel
	const 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 space
	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 };
	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）

其中，Σ代表协方差矩阵，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_modifier
	const glm::vec4 rot, // 表示旋转的四元数
	float* cov3D) // 结果：三维协方差矩阵
{
	// Create scaling matrix
	glm::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 rotation
	glm::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 quaternion
	glm::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 Sigma
	glm::mat3 Sigma = glm::transpose(M) * M;

	// Covariance is symmetric, only store upper right
	cov3D[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, // 该线程负责第几个Gaussian
	int 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 float4
	conic_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 rasterize
	bool 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 variables
	float 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 complete
	for (int i = 0; i < rounds; i++, toDo -= BLOCK_SIZE) //代码使用 rounds 控制循环的迭代次数，每次迭代处理一批点云数据。
	{	
		// 检查是否所有线程块都已经完成渲染：
		// 通过 __syncthreads_count 统计已经完成渲染的线程数，如果整个线程块都已完成，则跳出循环。
		// End if entire block votes that it is done rasterizing
		int 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 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();

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

公式（3）

		// 迭代处理当前批次的点云数据：
		// Iterate over current batch
		for (int j = 0; !done && j < min(BLOCK_SIZE, toDo); j++) //在当前批次的循环中，每个线程处理一条点云数据。
		{
			// Keep track of current position in range
			contributor++;

			// 计算当前点的投影坐标与锥体参数的差值：
			// 计算当前点在屏幕上的坐标 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、其他参考资料待补充。