diff --git "a/\345\255\246\344\271\240\347\254\224\350\256\260/\345\233\276\345\275\242\345\255\246/3D\351\253\230\346\226\257\344\273\243\347\240\201\350\247\243\346\236\220.md" "b/\345\255\246\344\271\240\347\254\224\350\256\260/\345\233\276\345\275\242\345\255\246/3D\351\253\230\346\226\257\344\273\243\347\240\201\350\247\243\346\236\220.md" index 453c06b..9cd2ec9 100644 --- "a/\345\255\246\344\271\240\347\254\224\350\256\260/\345\233\276\345\275\242\345\255\246/3D\351\253\230\346\226\257\344\273\243\347\240\201\350\247\243\346\236\220.md" +++ "b/\345\255\246\344\271\240\347\254\224\350\256\260/\345\233\276\345\275\242\345\255\246/3D\351\253\230\346\226\257\344\273\243\347\240\201\350\247\243\346\236\220.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)