渲染过程主要是把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 来实现。
Σ
′
=
V
Σ
V
T
\Sigma ^{'} = V \Sigma V^{T}
Σ′=VΣVT,其中
V
T
V^{T}
VT 是视图矩阵的转置。
相机坐标系到像素坐标系:
深度剔除:首先,需要剔除高斯在相机空间中的深度信息,以获得其在屏幕空间的2D表示。这可以通过忽略协方差矩阵的第三行和第三列来实现,因为这些行和列包含了深度方向的信息。
映射:然后,应用一个相机内参将2D点从相机空间映射到屏幕空间。这个变换可以通过投影矩阵 P 实现,该矩阵将相机空间中的点转换为屏幕空间中的点。对于协方差矩阵,我们只关心其在屏幕平面上的投影,因此可以使用投影矩阵的2x2子矩阵
P
2
D
P_{2D}
P2D来实现这一转换。
最终,2D屏幕空间中的协方差矩阵
Σ
′
′
=
P
2
D
Σ
′
P
2
D
T
\Sigma^{''} = P_{2D} \Sigma^{'} P_{2D}^{T}
Σ′′=P2DΣ′P2DT,
这里面,
Σ
′
\Sigma ^{'}
Σ′是映射到相机坐标系的协方差矩阵,
P
2
D
P_{2D}
P2D是根据相机内参组成的一个变换矩阵。
协方差矩阵这部分的映射代码见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;
}