加公式

学习笔记/图形学/3D高斯代码解析.md

@@ -524,13 +524,6 @@ $$\Sigma=RSS^TR^T$$
 ```cpp
 	// Compute 2D screen-space covariance matrix
 	float3 cov = computeCov2D(p_orig, focal_x, focal_y, tan_fovx, tan_fovy, cov3D, viewmatrix); //计算椭球投影成椭圆的样子[cov.x,cov.y,cov.y,cov.z]，即2D协方差矩阵[abbc]
-
-	// Invert covariance (EWA algorithm) //求2D协方差矩阵的特征值
-	float det = (cov.x * cov.z - cov.y * cov.y); //协方差矩阵[abbc]的模ac-b^2
-	if (det == 0.0f)
-		return; //协方差矩阵的模为0说明投影出来椭圆面积为0，退出
-	float det_inv = 1.f / det; //协方差矩阵[abbc]的特征值
-	float3 conic = { cov.z * det_inv, -cov.y * det_inv, cov.x * det_inv }; //协方差矩阵[abbc]的逆 det*(c,-b,a)
 ```
 
 ### 2D协方差矩阵计算函数`computeCov2D`
@@ -671,6 +664,30 @@ $$\Sigma'=JW\Sigma W^TJ^T$$
 }
 ```
 
+### 2D协方差矩阵求逆
+
+求坐标$\bm x$上的高斯函数值$G(\bm x)$公式来自3DGS原论文：
+
+$$G(\bm x)=e^{-\frac{1}{2}\bm x^T\Sigma^{-1}\bm x}$$
+
+所以为了后面渲染过程方便计算这里先算出了$\Sigma^{-1}$保存为`conic`，称为“锥体(conic)参数”（高斯球投影为实心椭圆，其边界为圆锥曲线）：
+
+$$
+\begin{aligned}
+\Sigma=\left[\begin{matrix}a&b\\b&c\end{matrix}\right]\\
+\Sigma^{-1}=\frac{\Sigma^*}{|\Sigma|}=\frac{\left[\begin{matrix}c&-b\\-b&a\end{matrix}\right]}{ac-b^2}
+\end{aligned}
+$$
+
+```cpp
+	// Invert covariance (EWA algorithm) //求2D协方差矩阵的行列式
+	float det = (cov.x * cov.z - cov.y * cov.y); //协方差矩阵[abbc]的模ac-b^2
+	if (det == 0.0f)
+		return; //协方差矩阵的模为0说明投影出来椭圆面积为0，退出
+	float det_inv = 1.f / det; //协方差矩阵[abbc]的行列式
+	float3 conic = { cov.z * det_inv, -cov.y * det_inv, cov.x * det_inv }; //协方差矩阵[abbc]的逆 = det*(c,-b,a)
+```
+
 ### 计算高斯球与成像平面上哪些tiles相交
 
 ```cpp
@@ -1010,11 +1027,15 @@ block(tile)里的每个thread(像素)都要对这个tile中的所有高斯球进
 
 ### for循环处理该像素的每个高斯球
 
-这个for循环分两个部分，首先`BLOCK_SIZE`个thread并行读取一批`BLOCK_SIZE`个高斯球进`collected_*`数组里，然后每个thread各自处理读进来的`BLOCK_SIZE`个高斯球：
 ```cpp
 	// Iterate over batches until all done or range is complete
 	for (int i = 0; i < rounds; i++, toDo -= BLOCK_SIZE)
 	{
+```
+
+这个for循环分两个部分，首先`BLOCK_SIZE`个thread并行读取一批`BLOCK_SIZE`个高斯球进`collected_*`数组里：
+
+```cpp
 		// End if entire block votes that it is done rasterizing
 		int num_done = __syncthreads_count(done); //计算一个block里done为true的数量
 		if (num_done == BLOCK_SIZE) //如果done有BLOCK_SIZE个true，则表明所有批次的高斯球全部处理完成，可以退出
@@ -1030,13 +1051,36 @@ block(tile)里的每个thread(像素)都要对这个tile中的所有高斯球进
 			collected_conic_opacity[block.thread_rank()] = conic_opacity[coll_id]; //读取待处理的高斯球2D协方差矩阵的逆和不透明度
 		}
 		block.sync(); //同步，确保读取全部完成
+```
+
+然后每个thread各自处理读进来的`BLOCK_SIZE`个高斯球：
 
+```cpp
 		// Iterate over current batch
 		for (int j = 0; !done && j < min(BLOCK_SIZE, toDo); j++) //block(tile)里的每个thread(像素)都要对这个tile中的所有高斯球进行渲染
 		{
 			// Keep track of current position in range
 			contributor++;
+```
+
+#### 求高斯函数值$G(\bm x)$中的指数部分
+
+求高斯函数值$G(\bm x)$公式来自3DGS原论文：
+
+$$G(\bm x)=e^{-\frac{1}{2}\bm x^T\Sigma^{-1}\bm x}$$
+
+其中的指数部分展开一下就是：
+
+$$\begin{aligned}
+-\frac{1}{2}\bm x^T\Sigma^{-1}\bm x&=-\frac{1}{2}\left[\begin{matrix}\Delta x&\Delta y\end{matrix}\right]\left[\begin{matrix}A&B\\B&C\end{matrix}\right]\left[\begin{matrix}\Delta x\\\Delta y\end{matrix}\right]\\
+&=A\Delta x^2 + 2B\Delta x \Delta y + C\Delta y^2
+\end{aligned}$$
+
+其中，$(\Delta x, \Delta y)^T=(x_i, y_i)^T-(x_{\text{pixel}}, y_{\text{pixel}})^T$为高斯球在成像平面上投影的高斯椭圆中心$(x_i, y_i)^T$到当前像素位置$(x_{\text{pixel}}, y_{\text{pixel}})^T$的向量。
 
+对应代码中的`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;`：
+
+```cpp
 			// Resample using conic matrix (cf. "Surface 
 			// Splatting" by Zwicker et al., 2001)
 			float2 xy = collected_xy[j]; //高斯球中心在像素平面的坐标
@@ -1045,7 +1089,21 @@ block(tile)里的每个thread(像素)都要对这个tile中的所有高斯球进
 			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; //通过相对位置、高斯球2D协方差矩阵的逆计算高斯分布的指数部分
 			if (power > 0.0f)
 				continue;
+```
+
+#### alpha-blending
+
+alpha-blending算法计算的当前像素颜色$C$由高斯球$i$的颜色$C_i$和透明度参数$\Alpha_i$按深度顺序由近及远混合而成，并最后加上背景颜色$C_{bg}$：
+
+$$\begin{aligned}
+\alpha_k&=\Alpha_kG(\bm x_k)\\
+T_j&=\prod_{k=1}^{j-1} (1 - \alpha_k)\\
+C&=T_{N+1}\cdot C_{bg}+\sum_{i=1}^{N}T_{i}\cdot\alpha_{i}\cdot C_{i}
+\end{aligned}$$
 
+对应代码中的`float alpha = min(0.99f, con_o.w * exp(power));`、`float test_T = T * (1 - alpha)`、`C[ch] += features[collected_id[j] * CHANNELS + ch] * alpha * T;`、`out_color[ch * H * W + pix_id] = C[ch] + T * bg_color[ch];`：
+
+```cpp
 			// Eq. (2) from 3D Gaussian splatting paper.
 			// Obtain alpha by multiplying with Gaussian opacity
 			// and its exponential falloff from mean.
@@ -1071,11 +1129,7 @@ block(tile)里的每个thread(像素)都要对这个tile中的所有高斯球进
 			last_contributor = contributor;
 		}
 	}
-```
 
-### 输出渲染结果
-
-```cpp
 	// All threads that treat valid pixel write out their final
 	// rendering data to the frame and auxiliary buffers.
 	if (inside)