3DGS中的光栅化渲染过程（结合代码）

渲染过程主要是把3DGS投影到2D图像平面上。

3DGS中每个高斯由位置（均值）、协方差矩阵（描述形状和方向）和不透明度（alpha）表示。
投影到2D包括投影均值的位置，和投影协方差矩阵。

另外介绍一下快速光栅化，
为了提高效率，将屏幕划分为16x16像素的块（tiles）。每个3D高斯可能与多个块重叠，因此需要在这些块中实例化。
使用GPU上的Radix排序算法对所有实例化的高斯进行排序，主要依据它们的深度信息。这种排序确保了在后续的混合过程中，能够按照从远到近的顺序处理高斯，从而减少所需的混合操作数量。

具体为如下步骤：
1.投影转换：首先，每个3D高斯从世界空间或相机空间转换到屏幕空间。这涉及到计算高斯的均值（中心点）在屏幕空间的位置，以及将高斯的协方差矩阵从3D空间转换到2D屏幕空间。

2.重叠检测：对于屏幕空间中的每个像素块（tile），需要确定哪些3D高斯与该块重叠。这通常通过计算高斯的边界（例如，使用高斯的协方差矩阵来定义一个边界框）并检查这个边界框是否与像素块相交来完成。

3.实例化：如果一个3D高斯与某个像素块重叠，那么这个高斯就会在该块中被“实例化”。这意味着在渲染算法中，会为这个块创建一个或多个对应的2D高斯表示，用于后续的光栅化和混合过程。

4.分配键值：为了在后续的排序和混合过程中有效地处理这些实例化的高斯，每个实例都会被分配一个键值（key），这个键值通常包含了高斯的深度信息和它所属的像素块的索引。这样，可以对所有实例化的高斯进行全局排序，以确保在混合过程中按照正确的顺序处理它们。

5.光栅化：排序后的高斯会被用来光栅化每个像素块，这个过程涉及到计算每个高斯对块内每个像素颜色的贡献，并根据高斯的不透明度进行混合。

其中，协方差矩阵由3D投影到2D的过程如下：
3D高斯的协方差矩阵 Σ 是一个3x3的对称矩阵，描述了高斯分布在3D空间中的扩散程度和方向。它包含了高斯在三个正交轴上的方差以及这些轴之间的协方差。
世界坐标系到相机坐标系：
首先，需要将3D高斯从世界空间转换到相机空间。这通常通过乘以一个视图矩阵 V 来实现。
${\Sigma }^{{}^{\prime }}=V\Sigma V^T$，其中 $V^T$ 是视图矩阵的转置。

相机坐标系到像素坐标系：

深度剔除：首先，需要剔除高斯在相机空间中的深度信息，以获得其在屏幕空间的2D表示。这可以通过忽略协方差矩阵的第三行和第三列来实现，因为这些行和列包含了深度方向的信息。
映射：然后，应用一个相机内参将2D点从相机空间映射到屏幕空间。这个变换可以通过投影矩阵 P 实现，该矩阵将相机空间中的点转换为屏幕空间中的点。对于协方差矩阵，我们只关心其在屏幕平面上的投影，因此可以使用投影矩阵的2x2子矩阵 $P_{2D}$来实现这一转换。
最终，2D屏幕空间中的协方差矩阵
${\Sigma }^{{}^{\prime \prime }}=P_{2D}{\Sigma }^{{}^{\prime }}{P}_{2D}^T$,
这里面，${\Sigma }^{{}^{\prime }}$是映射到相机坐标系的协方差矩阵，$P_{2D}$是根据相机内参组成的一个变换矩阵。

协方差矩阵这部分的映射代码见diff-gaussian-rasterization/cuda_rasterizer/forward.cu

__device__ float3 computeCov2D(const float3& mean, float focal_x, float focal_y, float tan_fovx, float tan_fovy, const float* cov3D, const float* viewmatrix)
{
	// 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.
	//使用视图矩阵将3D均值点变换为视图空间中的点 t
	float3 t = transformPoint4x3(mean, viewmatrix);
	//根据视野范围限制变换后的 x 和 y 坐标，以防止在视口边缘出现失真。
	const float limx = 1.3f * tan_fovx;
	const float limy = 1.3f * tan_fovy;
	const float txtz = t.x / t.z;
	const float tytz = t.y / t.z;
	t.x = min(limx, max(-limx, txtz)) * t.z;
	t.y = min(limy, max(-limy, tytz)) * t.z;
	
	//计算雅可比矩阵 J，它根据焦距和变换后的坐标，将3D空间的变化关系到2D空间的变化。
	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);
	//从视图矩阵构造3x3矩阵 W，用于进一步的变换，paper中的公式(5)
	glm::mat3 W = glm::mat3(
		viewmatrix[0], viewmatrix[4], viewmatrix[8],
		viewmatrix[1], viewmatrix[5], viewmatrix[9],
		viewmatrix[2], viewmatrix[6], viewmatrix[10]);
	//计算矩阵 T，结合视图矩阵和雅可比矩阵，用于将协方差从3D变换到2D
	glm::mat3 T = W * J;
	
	//从输入的 cov3D 数组构建3D协方差矩阵 Vrk
	glm::mat3 Vrk = glm::mat3(
		cov3D[0], cov3D[1], cov3D[2],
		cov3D[1], cov3D[3], cov3D[4],
		cov3D[2], cov3D[4], cov3D[5]);
	//使用 T 变换3D协方差矩阵，得到2D协方差矩阵 cov
	glm::mat3 cov = glm::transpose(T) * glm::transpose(Vrk) * T;

	// 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;
	//返回包含2D协方差矩阵相关分量
	return { float(cov[0][0]), float(cov[0][1]), float(cov[1][1]) };
}

preprocess中完成了筛选在相机视锥体内的点，3D投影2D，
2D椭圆的计算(通过求逆2D协方差矩阵得到)，确定高斯点影响的tiles范围，颜色计算。
存储：
深度值，屏幕空间半径，2D投影坐标，2D椭圆参数和不透明度，影响的瓦片(tiles)数量

    float3 p_view;
    //检查点是否在相机视锥体内
	if (!in_frustum(idx, orig_points, viewmatrix, projmatrix, prefiltered, p_view))
		return;

	// 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 };

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

	// Compute 2D screen-space covariance matrix
	//2D协方差计算：将3D协方差矩阵转换为2D屏幕空间协方差矩阵
	float3 cov = computeCov2D(p_orig, focal_x, focal_y, tan_fovx, tan_fovy, cov3D, viewmatrix);

	// Invert covariance (EWA algorithm)
	//计算2D椭圆（conic）：通过求逆2D协方差矩阵得到
	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.
	//确定屏幕空间范围
	//计算高斯点在屏幕上的半径 
	//确定高斯点影响的图像瓦片（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;
	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;

	// If colors have been precomputed, use them, otherwise convert
	// spherical harmonics coefficients to RGB color.
	//如果颜色预计算过，直接使用
	//否则，从球谐函数（spherical harmonics）系数计算RGB颜色
	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;
	}

	// 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
	//2D椭圆参数和不透明度
	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);

全部的渲染流程在rasterizer_impl.cu的forward函数中
大致如下：

CHECK_CUDA(FORWARD::preprocess
// 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
// Sort complete list of (duplicated) Gaussian indices by keys
CHECK_CUDA(cub::DeviceRadixSort::SortPairs(
CHECK_CUDA(FORWARD::render(

渲染函数render中，涉及paper中的公式(1) ~ (3)

// 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). 
	//通过这一步限制了透明度alpha的取值在[0,1]
	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;
}