3D Gaussian Splatting前向渲染代码解读

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、其他参考资料待补充。

原文地址：https://blog.csdn.net/xiner0114/article/details/142875111

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！