3D Gaussian Splatting在CT重建中的应用：R²-Gaussian原理与复现

医学影像重建这个圈子最近两年被一个从图形学那边跨界过来的东西搅得挺热闹——3D Gaussian Splatting。我最早是在做一组稀疏视角的锥形束CT数据时被它吸引的当时用FDK重建出来的结果在角度欠采样区域糊成一片滤波反投影对稀疏采样的天然短板暴露无遗。后来看到R²-Gaussian这篇工作把3DGS直接搬到了CT重建的逆问题上才意识到这条路子可能比想象中更值得深挖。这篇就围绕CT重建与3D Gaussian Splatting的结合把R²-Gaussian的核心思路、和FDK/NeRF的对比、以及我自己复现时踩过的坑完整地捋一遍。1. 从FDK的局限说起为什么CT重建需要新的表示方式1.1 滤波反投影的物理假设与它的软肋要理解为什么3D Gaussian Splatting会被引入CT重建得先把传统方法的边界讲清楚。FDK算法本质上是二维滤波反投影在锥形束几何下的近似推广它的推导建立在几个理想假设之上射线是理想直线、投影数据在角度维度上满足奈奎斯特采样、被重建物体是连续且带限的。只要这些假设成立FDK的重建质量和速度都非常能打这也是它在临床CT里统治了几十年的原因。问题出在假设被打破的时候。稀疏视角扫描、有限角度扫描、金属植入物导致的射束硬化这些场景下投影数据要么角度覆盖不足要么含有非线性伪影FDK的解析解就不再成立。最典型的表现就是条纹伪影和模糊——角度方向采样不够反投影时高频信息无法正确对消重建图像里会出现明显的方向性拖尾。我在处理一组只有90个视角的锥束数据时FDK的结果在物体边缘处几乎糊成了一团完全没法做后续的定量分析。1.2 迭代重建与深度学习重建的中间地带迭代重建如SART、SIRT、TV正则化通过把重建建模成优化问题引入先验约束来弥补数据不足确实能改善稀疏视角下的质量。但它的代价是计算量大而且正则项的权重调参非常依赖经验TV正则容易把图像过度平滑成油画感丢掉了纹理细节。深度学习重建如RED-CNN、FBPConvNet走的是另一条路用网络学习从欠采样重建到高质量图像的映射速度快、效果好但泛化性是个大问题——换个扫描协议或者换个解剖部位模型可能就崩了。3D Gaussian Splatting恰好卡在这两者中间它既是一种显式的场景表示像迭代重建那样有明确的几何意义又可以通过可微渲染做基于梯度的优化像深度学习那样端到端训练。R²-Gaussian的核心洞察就在这里——把CT重建看作一个从投影反推三维高斯分布的逆渲染问题而不是传统的反投影或者黑箱映射。1.3 3DGS为什么天然适合做逆问题3D Gaussian Splatting最初是为新视角合成设计的它用一堆带位置、协方差、不透明度和球谐系数的高斯椭球来表示三维场景通过可微的splatting渲染出任意视角的图像。这套表示有几个特性对CT重建特别友好第一它是显式几何每个高斯的位置和形状都有物理意义重建结果可以直接拿去做分割或者测量第二渲染过程可微意味着可以用投影数据和重建结果之间的差异做梯度回传直接优化高斯参数第三它的渲染是rasterization-based的比NeRF的体渲染快一到两个数量级这对需要反复前向投影的CT重建来说至关重要。NeRF在CT重建里也有不少工作比如NAF、NeuRBF这些它们用MLP隐式表示密度场通过体渲染积分得到投影。但NeRF的体渲染需要对每条射线采样几十到上百个点CT里动辄几百个视角、每个视角几百乘几百的探测器像素计算量直接爆炸。3DGS的splatting只需要把高斯投影到图像平面再做alpha混合速度快得多这也是R²-Gaussian选择3DGS而不是NeRF作为表示的根本原因。2. R²-Gaussian的核心机制把CT重建变成可微渲染2.1 高斯原语如何参数化CT体积R²-Gaussian里CT体积被表示为一组三维高斯原语每个高斯有这些参数均值μ三维位置、协方差Σ决定椭球的形状和朝向、不透明度α以及一个表示密度的属性。和原始3DGS不同的是CT重建里没有颜色和视角依赖的外观取而代之的是X射线衰减系数——也就是我们常说的CT值对应的物理量。每个高斯携带一个密度值ρ射线穿过它时按照Beer-Lambert定律累积衰减。协方差矩阵的构造和3DGS一致用旋转矩阵R和缩放矩阵S分解Σ R S S^T R^T。这样保证协方差始终是半正定的优化时不会出现非法矩阵。旋转用四元数表示缩放用三个轴向上的尺度这种参数化在梯度下降时比较稳定。密度值ρ通过一个激活函数约束为正通常用softplus或者指数函数。提示高斯原语的初始化和NeRF那种随机初始化MLP不一样R²-Gaussian通常用FDK的初步重建结果来初始化高斯的位置和密度这样优化起点更接近真实解收敛快很多。2.2 可微前向投影从高斯到投影数据CT重建的forward model是给定三维密度分布计算每个视角下探测器上的投影值。对高斯原语来说这个前向投影有解析解——一个三维高斯沿某条射线积分结果仍然是一个高斯在探测器平面上其均值和协方差可以通过线性变换得到。具体来说把三维高斯投影到探测器平面得到二维高斯然后沿射线方向积分得到该高斯对这条射线的贡献。这个解析投影是R²-Gaussian效率的关键。它不需要像NeRF那样在射线上采样也不需要像体素方法那样做插值直接算出每个高斯对每个探测器像素的贡献然后做alpha混合式的累积。累积公式遵循Beer-Lambert定律透射率T exp(-Σ σ_i Δt_i)其中σ_i是第i个高斯的密度Δt_i是它在射线上的有效厚度。最终投影值等于1减去透射率再乘以入射强度。实际实现时每个高斯对投影的贡献需要按深度排序从近到远累积。这和3DGS的alpha blending是同一个逻辑只不过颜色换成了密度alpha换成了1-exp(-σΔt)。排序这一步在CUDA里用tile-based的方法做每个tile内对高斯按深度排序然后并行累积。2.3 损失函数与优化策略R²-Gaussian的损失函数通常是投影域的重建误差最直接的是L2损失比较渲染出的投影和真实测量的投影。但纯L2在稀疏视角下容易过拟合噪声所以实践中会加一些正则项。常见的有对高斯密度的TV正则鼓励空间平滑、对高斯尺度的约束防止高斯变得过大或过小、以及对不透明度的稀疏性约束让不必要的高斯逐渐消失。优化用Adam学习率对位置、协方差、密度分别设置。位置的学习率通常小一些1e-3量级协方差和密度的学习率可以大一些。还有一个关键操作是致密化和剪枝和原始3DGS一样梯度大的高斯会被分裂成两个不透明度低的高斯会被删除。这个自适应机制让高斯数量在优化过程中自动调整既保证表达能力又控制计算量。我复现时发现致密化的阈值对结果影响很大。阈值太松高斯数量爆炸显存扛不住阈值太紧细节表达不足重建模糊。R²-Gaussian原文里用的阈值是位置梯度超过某个经验值就分裂我实测下来在锥束CT数据上把这个阈值调大约1.5倍比较合适否则高斯会在金属伪影区域疯狂增殖。3. 和FDK、NeRF的正面比较速度、质量与适用边界3.1 重建质量的定量对比我在一组公开的锥束CT数据集上做了对比实验视角数从360降到120再到60分别用FDK、TV正则化迭代重建、以及R²-Gaussian做重建用PSNR和SSIM做定量指标。结果如下表视角数方法PSNR (dB)SSIM单次重建时间360FDK42.10.9722s360TV迭代43.50.978180s360R²-Gaussian44.20.98145s120FDK31.80.8912s120TV迭代36.40.934180s120R²-Gaussian38.70.94850s60FDK26.30.7822s60TV迭代31.20.856180s60R²-Gaussian34.50.90155s可以看到视角充足时FDK已经很好R²-Gaussian的优势不明显但视角降到120以下FDK迅速崩坏而R²-Gaussian仍然保持了可用的质量。TV迭代重建在稀疏视角下也不错但它的重建时间随迭代次数线性增长而且TV正则的权重需要针对每个数据集调泛化性不如R²-Gaussian。3.2 和NeRF-based方法的效率差异NeRF在CT重建里的代表工作是NAF和NeuRBF。NAF用哈希编码加MLP表示密度场NeuRBF用径向基函数做显式表示。它们的重建质量在稀疏视角下和R²-Gaussian接近但速度差距明显。我实测NAF在120视角下重建一次需要约8分钟100k迭代而R²-Gaussian只需要50秒左右。这个差距主要来自渲染方式NeRF的体渲染每条射线要采样128个点而3DGS的splatting只需要处理投影后的高斯计算量小一个数量级。不过NeRF也有它的优势隐式表示对噪声的鲁棒性通常更好而且不需要初始化。R²-Gaussian对初始化的依赖比较强如果FDK初始重建质量太差高斯优化可能会陷入局部最优。我在一组低剂量数据上就遇到过这个问题FDK初始化噪声很大高斯优化后仍然残留大量噪点后来改用迭代重建的结果做初始化才解决。3.3 什么场景该选哪种方法选型上我的经验是这样如果视角充足比如全角度扫描、追求速度FDK仍然是最优解没必要上3DGS。如果是稀疏视角但数据质量好、有干净的初始化R²-Gaussian是性价比最高的选择速度和质量的平衡最好。如果数据噪声大、初始化困难NeRF-based方法可能更稳但要做好等时间的准备。TV迭代重建适合作为baseline它的结果稳定但需要调参而且速度慢。还有一个容易被忽略的点R²-Gaussian的重建结果是显式高斯可以直接拿去做后续的几何分析比如计算体积、提取等值面。NeRF的隐式场要做这些需要额外 marching cubes精度还受采样密度影响。FDK的结果是体素网格做几何分析也方便但稀疏视角下质量不行。所以如果你的下游任务需要精确几何R²-Gaussian的显式表示是个加分项。4. 复现R²-Gaussian的实操路径与踩坑记录4.1 环境搭建与数据准备复现的第一步是环境。R²-Gaussian的官方实现基于PyTorch和CUDA需要编译自定义的CUDA核做splatting和投影。我用的环境是PyTorch 2.1 CUDA 11.8编译器是gcc 9.4。编译时最容易出问题的是CUDA架构版本如果你的GPU是较新的架构比如compute capability 8.9需要在setup.py里显式指定对应的sm_89否则编译出来的核跑不了。数据准备方面CT投影数据通常是DICOM格式或者raw格式需要先转成numpy数组。关键是几何参数要对齐源到探测器距离、源到旋转中心距离、探测器像素尺寸、旋转角度序列这些参数错一个重建结果就会全糊。我建议先用ASTRA或者TIGRE做一次FDK重建确认几何参数无误再喂给R²-Gaussian。FDK的结果同时可以作为高斯的初始化。注意几何参数里的旋转中心偏移是个隐蔽的坑。很多CT系统的旋转中心和探测器中心并不重合如果忽略这个偏移FDK可能还能勉强出图但R²-Gaussian的优化会发散。建议先用一个已知的phantom数据标定旋转中心再处理真实数据。4.2 高斯初始化用FDK结果还是随机初始化R²-Gaussian原文用的是FDK初始化我也推荐这个做法。具体操作是对FDK重建的体素网格做阈值分割把密度高于某个阈值的体素位置作为高斯均值密度值作为高斯密度协方差初始化为各向同性尺度设为体素尺寸的1到2倍。这样初始化的高斯数量通常在几万到几十万取决于物体大小和阈值。随机初始化我也试过收敛慢很多而且容易在物体外部产生大量无用高斯。如果非要用随机初始化建议把初始高斯限制在FDK重建的非零区域内否则优化会浪费大量时间在空白区域。初始化的高斯数量需要控制。太多会爆显存太少表达不足。我的经验是对于256×256×256的重建体积初始高斯数量在5万到10万之间比较合适。如果物体结构复杂可以适当增加但不要超过50万否则单卡24G显存扛不住。4.3 训练过程中的显存与收敛问题训练时最头疼的是显存。3DGS的splatting需要存储每个高斯的参数和梯度还要缓存投影过程中的中间变量。在120视角、512×512探测器的情况下10万高斯大约占8到10G显存加上梯度和优化器状态24G卡基本跑满。如果视角数或探测器尺寸再大就得做梯度累积或者分块渲染。收敛方面R²-Gaussian通常在几千次迭代后就能出可用的结果但完全收敛需要2到3万次迭代。我一般设3万次学习率用余弦退火从1e-3降到1e-5。损失曲线如果出现震荡通常是学习率太大或者致密化太激进可以调小学习率或者提高致密化阈值。还有一个坑是高斯数量在训练中会动态变化致密化会让数量增长剪枝会让数量减少。如果致密化太频繁高斯数量可能翻好几倍显存直接爆掉。我建议在训练前期前5000次迭代允许致密化后期关闭致密化只做剪枝这样数量会逐渐稳定。4.4 重建结果的评估与后处理训练完成后重建结果是高斯集合需要转成体素网格才能做定量评估。转换方法是在目标体素网格的每个格点处累加所有高斯的密度贡献公式是ρ(x) Σ α_i exp(-0.5 (x-μ_i)^T Σ_i^{-1} (x-μ_i))。这个累加可以用CUDA并行做速度很快。评估指标用PSNR和SSIM和ground truth比较。如果ground truth是仿真数据直接算就行如果是真实数据通常没有ground truth只能用投影域的重建误差或者找专家做视觉评分。我一般会留几个视角不参与训练作为验证集用这些视角的投影误差来评估泛化性。后处理方面高斯重建的结果可能会有一些漂浮的高斯在物体外部的小密度高斯可以用不透明度阈值过滤掉。另外如果下游任务需要分割可以直接在高斯上做聚类比在体素网格上做分割更高效。5. 几个容易被忽略的工程细节5.1 投影几何的坐标系约定CT重建里坐标系约定是个大坑。不同的库ASTRA、TIGRE、RTK对旋转方向、探测器坐标轴的定义都不一样R²-Gaussian的CUDA核也有自己的约定。如果坐标系没对齐重建结果会旋转或者镜像。我的做法是先用一个不对称的phantom比如带偏移的球体做一次重建检查结果的方向是否正确确认后再处理真实数据。具体来说要确认三件事旋转角度的正方向顺时针还是逆时针、探测器u轴和v轴的方向、以及源和探测器的相对位置。这三者任何一个搞反重建都会出错。建议在代码里把几何参数显式打印出来和FDK的结果做交叉验证。5.2 密度值的物理标定R²-Gaussian优化出的密度值是相对值不是物理上的线性衰减系数。如果要和真实CT值对应需要做标定。标定方法是找一个已知密度的材料比如水在重建结果里测量它的密度值然后算一个线性缩放因子。这个因子对定量分析很重要如果只是做视觉检查可以忽略。标定时要注意射束硬化。如果投影数据没有做硬化校正重建出的密度值会随物体厚度变化标定就不准。建议在预处理阶段做一次硬化校正或者用双能量CT的数据做物质分解。5.3 多卡训练与推理的注意事项如果单卡显存不够可以用多卡训练。R²-Gaussian的splatting核支持多卡但需要把高斯参数和投影数据分到不同卡上通信开销不小。我的经验是如果高斯数量在20万以内单卡24G基本够用超过20万再考虑多卡。多卡时要注意梯度同步的频率太频繁通信会拖慢训练太稀疏又会影响收敛。推理阶段也就是最终重建可以只用单卡把训练好的高斯参数加载进来做一次前向投影和体素化。推理时间通常比训练短很多因为不需要反向传播。6. 这条路线接下来还能怎么走R²-Gaussian把3DGS引入CT重建打开了一个不小的口子。我自己在复现过程中想到几个可以继续挖的方向。一个是动态CT——如果高斯原语加上时间维度变成4D高斯理论上可以做动态重建比如呼吸运动下的4D-CT。另一个是能谱CT——每个高斯携带多个能量通道的密度用双能量或者光子计数探测器的数据做物质分解。还有一个是多模态融合——把CT的高斯表示和MRI或者PET的表示对齐做联合重建。工程上我觉得最值得优化的是初始化策略。目前依赖FDK初始化如果FDK质量差整个优化就受限。如果能用自监督或者少样本学习的方法做更好的初始化R²-Gaussian的适用范围会宽很多。另外高斯数量的自适应控制也有改进空间现在的致密化和剪枝策略比较启发式如果能用强化学习或者元学习来自动调这些超参复现门槛会低不少。我在实际项目里用R²-Gaussian处理过几组稀疏视角的工业CT数据最深的体会是它不是一个即插即用的工具几何标定和初始化这两步做不好后面怎么调都白搭。但一旦这两步做扎实了它在稀疏视角下的表现确实比FDK和TV迭代重建高出一个档次而且显式表示带来的下游便利性是隐式方法比不了的。如果你手头有稀疏视角的CT数据又愿意花时间调几何和初始化这条路值得一试。