简介InDuDoNet模型的Python复现源码包面向需要复现或深入研究该模型的科研人员、算法工程师及医学图像处理方向学习者聚焦CT图像器官与病变分割任务。项目覆盖数据加载、网络搭建、训练、推理与结果评估全流程并集成可视化与评价指标模块方便观察分割效果按Deeplesion、RatFemur等数据集分别提供训练与推理入口可据此重建实验环境、对照论文结果。压缩包内共61个文件、918KB44个Python脚本构成核心实现9个YAML配置管理多任务参数文本与Markdown说明、CSV记录实验细节MAT文件提供测试数据.gitignore维持仓库整洁目录按网络、训练器、工具等模块划分便于定位与扩展。已有385人学习浏览。资源不仅给出InDuDoNet、InDuDoNet_Plus、WNet等网络定义还配套requirements、开发日志和推理入口方便复现论文、开展消融实验或二次开发是理解与改进该模型时一个完整的代码基线。1. InDuDoNet 复现这件事到底值不值得花时间低剂量 CT 在临床里越来越常见可 sinogram投影域里一旦混进泊松噪声重建出来的 CT 图就会冒出条纹伪影和 CT 值漂移。InDuDoNet 走的是双域路线把投影域和图像域联合起来做数据一致性约束比单域 U-Net 后处理更能保住细小血管和骨边缘。很多做医学影像重建的工程师会把基于论文复现 InDuDoNet 的 Python 实现源码当作研究基线或预实验模型——它的网络不算大数据却要求得比较刁真正动手时会发现坑基本都在 Radon 变换的精度和两个子网络的协调上。这篇文章就按我自己的复现路径讲一遍适合想搭反向问题基线、准备论文对比实验或者想把手头 DICOM 数据跑通去噪管线的读者。2. InDuDoNet 的模型设计拆解双域网络结构与 Radon 变换的取舍2.1 为什么单域后处理不够投影域与图像域各自能修什么低剂量 CT 的退化过程很直接CT 扫描得到的是 X 射线穿过人体后的衰减数据也就是 sinogram投影域的数据经过滤波反投影FBP才能得到图像域。低剂量条件意味着每个角度上的射线计数都偏低投影域噪声近似服从泊松分布在重建图像里会表现为放射状条纹和整体低对比度模糊。单域网络通常只在图像域做后处理。你可以把这类方法理解成“拿一张已经被毒化的图像去修”U-Net 学到的大多是抑制条纹和拉升对比度可投影域里那些低频噪声早就被 FBP 放大成了伪影图像域很难无中生有地复原信噪比。反过来如果只在投影域做去噪sinogram 被修得太平滑重建后图像里的组织边缘会跟着损失。InDuDoNet 的核心思想就是让网络同时操作这两个域并保证两个域的中间结果互相约束。常见做法是把两个子网络都设计成 U-Net 变体投影域子网络先处理 sinogram输出修好的投影再走一层滤波反投影得到图图像域子网络在这张图上继续修细节最后把修过的图重新投影回 sinogram 域与原始测量值做一致性比对。我复现的时候会把“一致性”理解成一个软约束——它不要求网络输出机械地等于原始测量而是让修正量残差被限制在可信范围内。2.2 从 DuDoNet 到 InDuDoNet数据一致性层到底加了什么DuDoNet 系列最早把“双域”这个概念带进了金属伪影消除任务InDuDoNet 继承并强化了数据一致性层Dual-Domain Data Consistency简称 DDC。这个层的位置通常在图像域子网络的输出端。计算过程可以简化为投影域子网络输出修正后的 sinogramy_hat对 y_hat 做滤波反投影得到 x_init图像域子网络输出 x_out对 x_out 做前向 Radon 变换得到 y_proj一致性输出为y_dc y_hat alpha *y_orig - y_proj其中 y_orig 是原始低剂量 sinogram最后再把 y_dc 重建为图像送入后续网络块迭代。这里的 alpha 是残差权重。我一开始把 alpha 设为固定值 0.5结果训练曲线一直震荡后来改成可学习参数并用 0.2 初始化收敛稳定多了。工程上不建议把 alpha 设得太大否则网络会过度依赖原始测量双域网络容易退化成“一层残差校正 一层普通去噪 U-Net”。2.3 复现时的网络结构选型U-Net 变体与通道数、深度对照表InDuDoNet 的两个子网络不必完全一样。投影域处理的是角度方向的 sinogram分辨率通常也是 512x512对应探测器像素数和角度数但感受野需求更大所以通道数可以给足图像域子网络则更关注高频细节深度不能太深否则小病灶会被整体抹平。我复现时用的配置如下模块输入尺寸滤波起始通道深度归一化方式下采样方式投影域 U-Net1x180x512324InstanceNorm2x2 卷积 stride2图像域 U-Net1x512x512643InstanceNorm2x2 卷积 stride2DDC 层———无参数层—中间激活ReLU————投影域网络起始通道选 32是考虑到 sinogram 数据量本身就比图像域少显存压力相对可控图像域起始通道选 64是因为 512x512 图像里的纹理复杂度更高网络需要更多表达容量。实例归一化几乎可以说是这类重建任务的默认选择——用 BatchNorm 在小 batch 训练时会出现投影噪声的缩放偏差损失函数下得去重建图像却发灰。还有一个常被忽略的选型细节两个子网络的卷积核大小都用 3x3但在 DDC 层前面的最后一层卷积我故意用 1x1目的是让输出不要引入额外的空间相位偏移保证残差对齐。3. 复现前先把数据准备好低剂量 CT 数据集的配对切片与 sinogram 生成3.1 找数据来源与目录组织AAPM、LoDoPaB-CT 与 DICOM 序列的处理训练 InDuDoNet 需要成对数据同一扫描对象的正常剂量重建图和低剂量 sinogram。公开数据里AAPM 低剂量 CT 挑战赛的数据集梅奥诊所提供最常用但需要申请LoDoPaB-CT 数据集可以公开访问里面的投影数据和重建图都是现成的 HDF5 格式适合直接用来做模型基线对比。如果你手头只有临床 DICOM 序列就得先自己合成低剂量 sinogram——把清晰的 CT 图像做 Radon 变换得到干净投影再往投影上加泊松噪声然后重建出低剂量图。这样虽然引入一定的近似但在基线实验里完全够用。我一般把数据按下面这个结构组织data/ ct_clean/ # 512x512 归一化后的正常剂量图像 ct_low/ # 重建后的低剂量图像 sino_clean/ # 由 ct_clean 生成的干净 sinogram sino_low/ # 加噪后的 sinogram train_list.txt # 训练样本文件名列表训练时只读 train_list 指向的样本不维护内存里的全量数据集。这样切分数据、换数据集、做交叉验证都比较灵活。3.2 用 Python 把 DICOM 转成 HDF5 并生成配对 sinogram 的脚本我自己的实现里把 DICOM 转成 HDF5 是一个独立的预处理脚本不在训练循环里做。整个过程分四步读序列、转 HU、裁窗、生成 sinogram。import pydicom import numpy as np import h5py from glob import glob from scipy.ndimage import rotate def load_dicom_series(dicom_dir): 读取一个 DICOM 序列按切片位置排序并返回 3D 数组单位HU。 files glob(f{dicom_dir}/*.dcm) slices [pydicom.dcmread(f) for f in files] slices.sort(keylambda s: float(s.ImagePositionPatient[2])) array np.stack([s.pixel_array for s in slices]) rescale np.vectorize(lambda s: float(s.RescaleSlope)) intercept np.vectorize(lambda s: float(s.RescaleIntercept)) array array * rescale(slices) intercept(slices) return array def normalize_ct(image, low-1000, high400): 将 HU 值裁剪并归一化到 [0,1]低剂量去噪常用窗。 image np.clip(image, low, high) image (image - low) / (high - low) return image def radon_2d(image, n_angles180): 用旋转积分生成 sinogram角度范围 [0, 180)。 h, w image.shape sinogram np.zeros((n_angles, h), dtypenp.float32) for i, angle in enumerate(np.linspace(0, 180, n_angles, endpointFalse)): rotated rotate(image, angle, reshapeFalse, order1) sinogram[i] rotated.sum(axis0) return sinogram with h5py.File(ct_pairs.h5, w) as f: ct_3d load_dicom_series(./dicom_patient_01) ct_2d normalize_ct(ct_3d[80]) # 取中间切片演示 ct_2d ct_2d[::2, ::2] # 降采样到 256x256 以节省演示时间 sino radon_2d(ct_2d, n_angles180) f.create_dataset(ct_clean, datact_2d) f.create_dataset(sino_clean, datasino)逻辑说明这里 sinogram 的 shape 是 (180, 256)180 是角度数256 是探测器方向采样数。角度方向我固定为 0 度对应竖直方向投影代码里用rotate(image, angle)把图像旋转到该角度后沿 axis0 求和所以坐标原点在图像中心。siogram 的行的顺序对应投影角度和后面的 Radon 模块保持一致。参数说明low-1000, high400是骨窗/软组织窗常用的范围如果任务注重骨结构high 可以调到 1000 以上但会把软组织对比度压小。这里rescale/ intercept是直接从 DICOM 头读的实际临床数据里这两项必须读不能默认存在任何一个固定文件。order1表示线性插值Radon 变换里用 order1 已经足够没必要用三次样条反而会引入额外平滑。3.3 训练/验证/测试切分与数据增强随机裁剪、翻转、角度遮挡数据切分有个容易被忽略的原则同一个患者的所有切片必须落在同一份集合里否则模型会在验证集上得到非常乐观的指标。我按患者目录做分组切分而不是按切片打散。增强方面InDuDoNet 训练时的主流做法是随机裁剪 patch。全图 512x512 训练时 GPU 显存开销大且每个 batch 包含的解剖结构太单一不利于网络学习噪声分布。裁剪到 256x256 后batch size 可以开到 4 而显存压力可控。翻转增强只做水平和垂直翻转因为医学图像做旋转增强会改变投影角度语义除非你同步旋转 sinogram 和 CT 图否则容易把网络的几何先验搞乱。我也是只做翻转不做强度变换避免破坏低剂量噪声的统计特性。4. 用 PyTorch 实现 InDuDoNet 核心模块从 DDC 到联合损失4.1 可微的 Radon/iRadon 模块两种实现路线与代码InDuDoNet 的反向传播必须穿过 Radon/iRadon 模块所以这两个变换必须是可微的。skimage 自带 radon/iradon 但梯度传不过去。常见的做法是自己在 PyTorch 里用空间变换实现一个近似可微版本或者封装 ASTRA Toolbox 的 GPU 算子。封装 ASTRA 的性能确实好但部署环境配置非常折腾而且论文复现阶段往往需要快速迭代模块设计。我用的是 PyTorch 原生实现基于grid_sample的旋转加积分import torch import torch.nn.functional as F import math class RadonProjection(torch.nn.Module): 前向投影把图像转成 sinogram平行束近似。 def __init__(self, n_angles180): super().__init__() self.n_angles n_angles def forward(self, images): # images: (B, 1, H, W)值为 0~1 B, C, H, W images.shape angles torch.linspace(0, 180, self.n_angles, deviceimages.device) sino [] for theta in angles: rad torch.deg2rad(theta) cos_t torch.cos(rad) sin_t torch.sin(rad) # 旋转矩阵用于 grid_sample 的 affine_grid rot_mat torch.tensor( [[cos_t, -sin_t, 0], [sin_t, cos_t, 0]], dtypeimages.dtype, deviceimages.device ) rot_mat_batched rot_mat.unsqueeze(0).repeat(B, 1, 1) grid F.affine_grid(rot_mat_batched, images.shape, align_cornersTrue) rotated F.grid_sample( images, grid, modebilinear, padding_modezeros, align_cornersTrue ) # 沿宽度方向求和得到该角度下的投影向量 sino.append(rotated.sum(dim-1)) return torch.stack(sino, dim-1) # (B, 1, H, n_angles) class FilteredBackProjection(torch.nn.Module): 滤波反投影对 sinogram 做 ramp 滤波后反投影累加。 def __init__(self, n_angles180): super().__init__() self.n_angles n_angles def forward(self, sino): # sino: (B, 1, H, A)H 是探测器维度A 是角度数 B, C, H, A sino.shape filtered torch.fft.rfft(sino, dim2) freq torch.fft.rfftfreq(H, devicesino.device) ramp torch.abs(freq) * H # ramp filter filtered filtered * ramp proj_filtered torch.fft.irfft(filtered, nH, dim2).real recon torch.zeros_like(sino[:, :, :, 0]) # (B,1,H) rec_size H angles torch.linspace(0, 180, A, devicesino.device) for i, theta in enumerate(angles): rad torch.deg2rad(theta) cos_t torch.cos(rad) sin_t torch.sin(rad) rot_mat torch.tensor( [[cos_t, -sin_t, 0], [sin_t, cos_t, 0]], dtypesino.dtype, devicesino.device ) rot_mat_batched rot_mat.unsqueeze(0).repeat(B, 1, 1) ones torch.ones((B, 1, H, rec_size), devicesino.device) grid F.affine_grid(rot_mat_batched, ones.shape, align_cornersTrue) # 将一维投影向量平铺成二维平面旋转到对应角度 proj_vec proj_filtered[:, :, :, i] # (B,1,H) proj_2d proj_vec.unsqueeze(-1).expand(-1, -1, -1, rec_size) rotated F.grid_sample( proj_2d, grid, modebilinear, padding_modezeros, align_cornersTrue ) recon rotated[:, 0] # 角度密度修正乘以 pi / A recon recon * (math.pi / A) return recon.unsqueeze(1)逻辑说明前向投影的复杂度是 A 次旋转 积分A180 时单张 512x512 图像的投影在 GPU 上约几十毫秒训练可以接受。滤波反投影里的 ramp filter 我做了幅度修正abs(freq) * H的系数不能省否则重建图整体会偏暗或者偏亮。滤波后每个角度一维投影平铺成二维平面旋转累加得到重建图。最后的pi / A是平行束 FBP 的归一化系数。参数说明n_angles180是最常用配置对应医学 CT 扫描的常见角度抽样。如果数据里 sinogram 是 360 度范围记得把角度序列改成 0~360。自定义视角时n_angles越大条纹伪影越少但训练显存和时间线性增加。4.2 数据一致性层 DDC 的 PyTorch 实现DDC 层是 InDuDoNet 最有工程价值的部分。网上很多源码实现把它放在最后做个简单的加法但没有把 Radon 变换接上等于让两个域在中间断开。我这里给出一个带残差的完整版本class DataConsistencyLayer(torch.nn.Module): 双域数据一致性层图像域输出回到投影域与原始投影做残差约束。 def __init__(self, alpha0.2, n_angles180): super().__init__() self.alpha torch.nn.Parameter(torch.tensor(float(alpha))) self.radon RadonProjection(n_anglesn_angles) self.fbp FilteredBackProjection(n_anglesn_angles) def forward(self, sino_input, sino_enhanced, img_out): # sino_input: 原始低剂量 sinogram # sino_enhanced: 投影域子网络输出 # img_out: 图像域子网络输出 x_init self.fbp(sino_enhanced) # 投影域结果重建图 x_refine x_init img_out # 残差式图像细化 sino_refine self.radon(x_refine) # 重新投影回 sinogram # 一致性残差修正量 alpha * (原始 - 重投影) y_dc sino_enhanced self.alpha * (sino_input - sino_refine) x_dc self.fbp(y_dc) return x_dc, y_dc逻辑说明这里的img_out来自图像域子网络的输出但实际落地时我会把图像域 U-Net 的输出设计成残差模式即 U-Net 学的是修复增量而不是直接输出完整图像。x_init img_out合成最终图像后再做 Radon计算原始投影与重投影的差异。这个差异乘上可学习的 alpha 后回到 sinogram 域保证投影域信息没有被丢弃。alpha 初始 0.2 的经验是低剂量噪声越强alpha 越要往小调。alpha 太大时网络会过度相信原始投影中的噪声样本输出图像会有高频抖动alpha 太小时 DDC 层退化成恒等映射双域约束消失。训练过程中可以把 alpha 的输出值打出来观察它在 0.1~0.3 之间浮动都是正常的。4.3 联合损失函数与配比感知损失之外还要盯住投影域InDuDoNet 的训练损失通常由几部分构成图像域重建损失、投影域重建损失、以及可选的感知损失。我只用 L1 作为图像域重建损失因为 L2 会让输出图像过于平滑低剂量去噪任务里 CT 值偏置更容易出现在 L2 优化的结果里。损失项监督对象权重备注L1_img网络最终图像 vs 正常剂量图像1.0主损失SSIM_img网络最终图像 vs 正常剂量图像0.1提升结构相似度L1_sinoy_dc vs 干净 sinogram0.05保证投影域不去过头L1_img_slicex_dc vs 正常剂量图像0.1DDC 输出再补一次监督投影域损失权重不能超过 0.1否则网络会把注意力放在把 sinogram 修漂亮而 sinogram 里的微小误差经 FBP 放大后反而毁掉图像边缘。我在训练时每 5 个 epoch 会把各损失项分别打印出来如果 L1_sino 下降很快但 SSIM 停滞说明图像域子网络没有吃到有效梯度优先检查 DDC 层里 alpha 是否退化成了 0或者 Radon/iRadon 的角度方向是否颠倒。图像域 L1 和 SSIM 的配比按上面表格基本够用。感知损失不是必须的加不加对最终 PSNR 影响不大但会让纹理看起来更自然——如果做主观评价建议保留一个低权重的 VGG 感知损失。5. InDuDoNet 训练常见问题与排查显存失控、结果发糊与投影域噪声5.1 显存爆掉为什么 sinogram 的 batch size 不能照搬普通 2D 网络现象把 batch size 设成 8训练刚开始就 OOM报错提示是 CUDA out of memory。原因InDuDoNet 虽然看起来只有两个子网络但 DDC 层里的 RadonProjection 和 FilteredBackProjection 会在计算图中保存每个角度的中间旋转图180 个角度的变量叠加起来比网络本身还占显存。尤其grid_sample的梯度回传需要保存输入和网格这块开销经常被忽略。解决batch size 降到 2~4patch 从 512 降到 256。还不够的话用梯度累积模拟大 batch# 以 batch2、累积 4 步模拟 batch8 的梯度 python train_indunonet.py --batch-size 2 --grad-accum-steps 4 --patch-size 256 --amp同时在训练脚本里加两行配置torch.backends.cudnn.benchmark True把 AMP 打开。AMP 对 InDuDoNet 尤其有效因为 Radon 变换里的矩阵乘法占了大量算力半精度能把显存占用砍掉将近一半而且重建精度损失肉眼不可见。如果开了 AMP 仍然爆显存把投影域 U-Net 的起始通道从 32 降到 24换来的是投影域的表达能力略降一般可接受。5.2 loss 降了但图像发糊雾感从哪里来怎么定位现象训练 loss 曲线很漂亮PSNR 也不低但输出图像像蒙了一层薄雾组织边缘对比度明显下降。原因这是低剂量去噪复现里最典型的翻车现场。第一个嫌疑是 DDC 层的 alpha 过大压制了网络学习高频细节的能力第二个嫌疑是滤波反投影里的 ramp filter 幅度不对重建图像本身就被压低第三个原因是投影域 L1 权重过高网络学成了“把 sinogram 做均值模糊”就能降低损失。排查时我把网络中间输出可视化保存下来def dump_debug_images(model, sino_low, tag): 保存 DDC 层前后的中间结果用于判断发糊出自哪一步。 model.eval() with torch.no_grad(): sino_enh model.proj_unet(sino_low) x_init model.fbp(sino_enh) x_refine model.img_unet(x_init) y_dc, x_dc model.ddc(sino_low, sino_enh, x_refine) torch.save({ sino_enh: sino_enh, x_init: x_init, x_refine: x_refine, y_dc: y_dc, x_dc: x_dc, }, fdebug_{tag}.pt)然后在本地写个小脚本把每个张量做 min-max 归一化并存成 PNG。如果x_init已经发糊问题在投影域子网络和 FBP如果x_init清晰但x_final发糊问题在图像域子网络。我踩过的坑是 x_init 清晰、最终输出发糊最后定位到图像域 U-Net 深度太深3x3 卷积堆叠让感受野覆盖了太多无关联组织改成 3 层之后雾感明显缓解。5.3 测试时 Filtered Back Projection 与训练时不一致伪影从哪冒出来现象训练集指标很好测试集上一张真实 DICOM 图输出出现放射状发散伪影甚至图像方向都反了。原因训练时用的自定义 Radon 模块和测试时用的 skimage 或 ASTRA 重建模块角度约定不一致。最常见的是角度起点不一致有的库 0 度对应水平投影有的库 0 度对应竖直投影还有的 sinogram 第一行是 0 度有的最后一行才是 0 度。InDuDoNet 在训练时学了自定义模块的角度约定测试用另一个库就等于输入分布错位。解决训练和测试必须共用一个 Radon/iRadon 实现。我建议在项目里固定写一个geometry.py里面只保留自己实现的模块不要在推理脚本里用skimage.transform.iradon。每改一次角度参数就用 Shepp-Logan 模体做一次 round-trip 校验def validate_geometry(): from skimage.data import shepp_logan_phantom phantom shepp_logan_phantom().astype(np.float32) phantom torch.from_numpy(phantom).unsqueeze(0).unsqueeze(0) radon RadonProjection(n_angles180) fbp FilteredBackProjection(n_angles180) sino radon(phantom) recon fbp(sino) error (recon - phantom).abs().mean().item() print(fround-trip MAE: {error:.4f}) # 视觉检查recon 纹理应与 phantom 一致无方向偏转参数说明round-trip MAE对 512x512 的 phantom 一般应在 0.05 以下超过 0.1 就要检查 ramp filter 的系数或者旋转方向是否反了。我在一次换 PyTorch 版本后遇到过align_cornersTrue/False改变导致的角度偏移 0.5 像素当时就是靠这个校验函数发现的。6. 把复现的 InDuDoNet 变成可交付的推理管线torchscript 导出与验证复现的价值在能稳定地跑新数据。训练完模型后我用 torch.jit.trace 把整个推理链路固定下来避免训练脚本和推理脚本因为代码版本不同产生微妙差异。trace 时要注意输入尺寸固定InDuDoNet 的 DDC 层包含角度循环动态长度输入会触发重新 trace有时会错误地把已有权重覆盖掉所以我在 trace 前会先把批处理尺寸固定成 1x1x512x512。model.eval() dummy_input torch.randn(1, 1, 512, 512) traced torch.jit.trace(model, dummy_input) traced.save(indunonet_ct.pt)若要再提速可以用半精度推理和 ONNX Runtime 替代原生 PyTorch。半精度推理前必须验证重投影再重建的误差有些模型在 fp16 下 DDC 残差 loss 下降不稳定最好做一个 PSNR 对比实验再决定是否默认开启。推理管线里我每次还会带一个后置校验输出图的 CT 值分布范围是否和输入一致如果均值漂移超过 20 HU基本可以断定推理路径用错了重建模块。这些细节在复现 InDuDoNet 的源码工程里很容易忽略。我自己第一次交付时就是忘了训练和推理用不同 Radon 实现导致同一份数据线上线下结果对不上此后固定了一套 geometry 模块任何实验都先跑 phantom 校验再跑临床数据。希望这些踩坑记录能帮你少走一段弯路。复现 InDuDoNet 并最终跑通一份 Python 实现源码的过程本质上就是在磨你对双域数据流的掌控力。把这套模块沉淀成你自己的工具库以后再做其他 CT 重建网络的对比实验会顺畅很多。本文还有配套的精品资源点击获取 
长河君 · 2026/10/5 2:09:22 · 阅读 1.2k InDuDoNet复现指南：双域网络实现低剂量CT去噪与Radon变换
简介InDuDoNet模型的Python复现源码包面向需要复现或深入研究该模型的科研人员、算法工程师及医学图像处理方向学习者聚焦CT图像器官与病变分割任务。项目覆盖数据加载、网络搭建、训练、推理与结果评估全流程并集成可视化与评价指标模块方便观察分割效果按Deeplesion、RatFemur等数据集分别提供训练与推理入口可据此重建实验环境、对照论文结果。压缩包内共61个文件、918KB44个Python脚本构成核心实现9个YAML配置管理多任务参数文本与Markdown说明、CSV记录实验细节MAT文件提供测试数据.gitignore维持仓库整洁目录按网络、训练器、工具等模块划分便于定位与扩展。已有385人学习浏览。资源不仅给出InDuDoNet、InDuDoNet_Plus、WNet等网络定义还配套requirements、开发日志和推理入口方便复现论文、开展消融实验或二次开发是理解与改进该模型时一个完整的代码基线。1. InDuDoNet 复现这件事到底值不值得花时间低剂量 CT 在临床里越来越常见可 sinogram投影域里一旦混进泊松噪声重建出来的 CT 图就会冒出条纹伪影和 CT 值漂移。InDuDoNet 走的是双域路线把投影域和图像域联合起来做数据一致性约束比单域 U-Net 后处理更能保住细小血管和骨边缘。很多做医学影像重建的工程师会把基于论文复现 InDuDoNet 的 Python 实现源码当作研究基线或预实验模型——它的网络不算大数据却要求得比较刁真正动手时会发现坑基本都在 Radon 变换的精度和两个子网络的协调上。这篇文章就按我自己的复现路径讲一遍适合想搭反向问题基线、准备论文对比实验或者想把手头 DICOM 数据跑通去噪管线的读者。2. InDuDoNet 的模型设计拆解双域网络结构与 Radon 变换的取舍2.1 为什么单域后处理不够投影域与图像域各自能修什么低剂量 CT 的退化过程很直接CT 扫描得到的是 X 射线穿过人体后的衰减数据也就是 sinogram投影域的数据经过滤波反投影FBP才能得到图像域。低剂量条件意味着每个角度上的射线计数都偏低投影域噪声近似服从泊松分布在重建图像里会表现为放射状条纹和整体低对比度模糊。单域网络通常只在图像域做后处理。你可以把这类方法理解成“拿一张已经被毒化的图像去修”U-Net 学到的大多是抑制条纹和拉升对比度可投影域里那些低频噪声早就被 FBP 放大成了伪影图像域很难无中生有地复原信噪比。反过来如果只在投影域做去噪sinogram 被修得太平滑重建后图像里的组织边缘会跟着损失。InDuDoNet 的核心思想就是让网络同时操作这两个域并保证两个域的中间结果互相约束。常见做法是把两个子网络都设计成 U-Net 变体投影域子网络先处理 sinogram输出修好的投影再走一层滤波反投影得到图图像域子网络在这张图上继续修细节最后把修过的图重新投影回 sinogram 域与原始测量值做一致性比对。我复现的时候会把“一致性”理解成一个软约束——它不要求网络输出机械地等于原始测量而是让修正量残差被限制在可信范围内。2.2 从 DuDoNet 到 InDuDoNet数据一致性层到底加了什么DuDoNet 系列最早把“双域”这个概念带进了金属伪影消除任务InDuDoNet 继承并强化了数据一致性层Dual-Domain Data Consistency简称 DDC。这个层的位置通常在图像域子网络的输出端。计算过程可以简化为投影域子网络输出修正后的 sinogramy_hat对 y_hat 做滤波反投影得到 x_init图像域子网络输出 x_out对 x_out 做前向 Radon 变换得到 y_proj一致性输出为y_dc y_hat alpha *y_orig - y_proj其中 y_orig 是原始低剂量 sinogram最后再把 y_dc 重建为图像送入后续网络块迭代。这里的 alpha 是残差权重。我一开始把 alpha 设为固定值 0.5结果训练曲线一直震荡后来改成可学习参数并用 0.2 初始化收敛稳定多了。工程上不建议把 alpha 设得太大否则网络会过度依赖原始测量双域网络容易退化成“一层残差校正 一层普通去噪 U-Net”。2.3 复现时的网络结构选型U-Net 变体与通道数、深度对照表InDuDoNet 的两个子网络不必完全一样。投影域处理的是角度方向的 sinogram分辨率通常也是 512x512对应探测器像素数和角度数但感受野需求更大所以通道数可以给足图像域子网络则更关注高频细节深度不能太深否则小病灶会被整体抹平。我复现时用的配置如下模块输入尺寸滤波起始通道深度归一化方式下采样方式投影域 U-Net1x180x512324InstanceNorm2x2 卷积 stride2图像域 U-Net1x512x512643InstanceNorm2x2 卷积 stride2DDC 层———无参数层—中间激活ReLU————投影域网络起始通道选 32是考虑到 sinogram 数据量本身就比图像域少显存压力相对可控图像域起始通道选 64是因为 512x512 图像里的纹理复杂度更高网络需要更多表达容量。实例归一化几乎可以说是这类重建任务的默认选择——用 BatchNorm 在小 batch 训练时会出现投影噪声的缩放偏差损失函数下得去重建图像却发灰。还有一个常被忽略的选型细节两个子网络的卷积核大小都用 3x3但在 DDC 层前面的最后一层卷积我故意用 1x1目的是让输出不要引入额外的空间相位偏移保证残差对齐。3. 复现前先把数据准备好低剂量 CT 数据集的配对切片与 sinogram 生成3.1 找数据来源与目录组织AAPM、LoDoPaB-CT 与 DICOM 序列的处理训练 InDuDoNet 需要成对数据同一扫描对象的正常剂量重建图和低剂量 sinogram。公开数据里AAPM 低剂量 CT 挑战赛的数据集梅奥诊所提供最常用但需要申请LoDoPaB-CT 数据集可以公开访问里面的投影数据和重建图都是现成的 HDF5 格式适合直接用来做模型基线对比。如果你手头只有临床 DICOM 序列就得先自己合成低剂量 sinogram——把清晰的 CT 图像做 Radon 变换得到干净投影再往投影上加泊松噪声然后重建出低剂量图。这样虽然引入一定的近似但在基线实验里完全够用。我一般把数据按下面这个结构组织data/ ct_clean/ # 512x512 归一化后的正常剂量图像 ct_low/ # 重建后的低剂量图像 sino_clean/ # 由 ct_clean 生成的干净 sinogram sino_low/ # 加噪后的 sinogram train_list.txt # 训练样本文件名列表训练时只读 train_list 指向的样本不维护内存里的全量数据集。这样切分数据、换数据集、做交叉验证都比较灵活。3.2 用 Python 把 DICOM 转成 HDF5 并生成配对 sinogram 的脚本我自己的实现里把 DICOM 转成 HDF5 是一个独立的预处理脚本不在训练循环里做。整个过程分四步读序列、转 HU、裁窗、生成 sinogram。import pydicom import numpy as np import h5py from glob import glob from scipy.ndimage import rotate def load_dicom_series(dicom_dir): 读取一个 DICOM 序列按切片位置排序并返回 3D 数组单位HU。 files glob(f{dicom_dir}/*.dcm) slices [pydicom.dcmread(f) for f in files] slices.sort(keylambda s: float(s.ImagePositionPatient[2])) array np.stack([s.pixel_array for s in slices]) rescale np.vectorize(lambda s: float(s.RescaleSlope)) intercept np.vectorize(lambda s: float(s.RescaleIntercept)) array array * rescale(slices) intercept(slices) return array def normalize_ct(image, low-1000, high400): 将 HU 值裁剪并归一化到 [0,1]低剂量去噪常用窗。 image np.clip(image, low, high) image (image - low) / (high - low) return image def radon_2d(image, n_angles180): 用旋转积分生成 sinogram角度范围 [0, 180)。 h, w image.shape sinogram np.zeros((n_angles, h), dtypenp.float32) for i, angle in enumerate(np.linspace(0, 180, n_angles, endpointFalse)): rotated rotate(image, angle, reshapeFalse, order1) sinogram[i] rotated.sum(axis0) return sinogram with h5py.File(ct_pairs.h5, w) as f: ct_3d load_dicom_series(./dicom_patient_01) ct_2d normalize_ct(ct_3d[80]) # 取中间切片演示 ct_2d ct_2d[::2, ::2] # 降采样到 256x256 以节省演示时间 sino radon_2d(ct_2d, n_angles180) f.create_dataset(ct_clean, datact_2d) f.create_dataset(sino_clean, datasino)逻辑说明这里 sinogram 的 shape 是 (180, 256)180 是角度数256 是探测器方向采样数。角度方向我固定为 0 度对应竖直方向投影代码里用rotate(image, angle)把图像旋转到该角度后沿 axis0 求和所以坐标原点在图像中心。siogram 的行的顺序对应投影角度和后面的 Radon 模块保持一致。参数说明low-1000, high400是骨窗/软组织窗常用的范围如果任务注重骨结构high 可以调到 1000 以上但会把软组织对比度压小。这里rescale/ intercept是直接从 DICOM 头读的实际临床数据里这两项必须读不能默认存在任何一个固定文件。order1表示线性插值Radon 变换里用 order1 已经足够没必要用三次样条反而会引入额外平滑。3.3 训练/验证/测试切分与数据增强随机裁剪、翻转、角度遮挡数据切分有个容易被忽略的原则同一个患者的所有切片必须落在同一份集合里否则模型会在验证集上得到非常乐观的指标。我按患者目录做分组切分而不是按切片打散。增强方面InDuDoNet 训练时的主流做法是随机裁剪 patch。全图 512x512 训练时 GPU 显存开销大且每个 batch 包含的解剖结构太单一不利于网络学习噪声分布。裁剪到 256x256 后batch size 可以开到 4 而显存压力可控。翻转增强只做水平和垂直翻转因为医学图像做旋转增强会改变投影角度语义除非你同步旋转 sinogram 和 CT 图否则容易把网络的几何先验搞乱。我也是只做翻转不做强度变换避免破坏低剂量噪声的统计特性。4. 用 PyTorch 实现 InDuDoNet 核心模块从 DDC 到联合损失4.1 可微的 Radon/iRadon 模块两种实现路线与代码InDuDoNet 的反向传播必须穿过 Radon/iRadon 模块所以这两个变换必须是可微的。skimage 自带 radon/iradon 但梯度传不过去。常见的做法是自己在 PyTorch 里用空间变换实现一个近似可微版本或者封装 ASTRA Toolbox 的 GPU 算子。封装 ASTRA 的性能确实好但部署环境配置非常折腾而且论文复现阶段往往需要快速迭代模块设计。我用的是 PyTorch 原生实现基于grid_sample的旋转加积分import torch import torch.nn.functional as F import math class RadonProjection(torch.nn.Module): 前向投影把图像转成 sinogram平行束近似。 def __init__(self, n_angles180): super().__init__() self.n_angles n_angles def forward(self, images): # images: (B, 1, H, W)值为 0~1 B, C, H, W images.shape angles torch.linspace(0, 180, self.n_angles, deviceimages.device) sino [] for theta in angles: rad torch.deg2rad(theta) cos_t torch.cos(rad) sin_t torch.sin(rad) # 旋转矩阵用于 grid_sample 的 affine_grid rot_mat torch.tensor( [[cos_t, -sin_t, 0], [sin_t, cos_t, 0]], dtypeimages.dtype, deviceimages.device ) rot_mat_batched rot_mat.unsqueeze(0).repeat(B, 1, 1) grid F.affine_grid(rot_mat_batched, images.shape, align_cornersTrue) rotated F.grid_sample( images, grid, modebilinear, padding_modezeros, align_cornersTrue ) # 沿宽度方向求和得到该角度下的投影向量 sino.append(rotated.sum(dim-1)) return torch.stack(sino, dim-1) # (B, 1, H, n_angles) class FilteredBackProjection(torch.nn.Module): 滤波反投影对 sinogram 做 ramp 滤波后反投影累加。 def __init__(self, n_angles180): super().__init__() self.n_angles n_angles def forward(self, sino): # sino: (B, 1, H, A)H 是探测器维度A 是角度数 B, C, H, A sino.shape filtered torch.fft.rfft(sino, dim2) freq torch.fft.rfftfreq(H, devicesino.device) ramp torch.abs(freq) * H # ramp filter filtered filtered * ramp proj_filtered torch.fft.irfft(filtered, nH, dim2).real recon torch.zeros_like(sino[:, :, :, 0]) # (B,1,H) rec_size H angles torch.linspace(0, 180, A, devicesino.device) for i, theta in enumerate(angles): rad torch.deg2rad(theta) cos_t torch.cos(rad) sin_t torch.sin(rad) rot_mat torch.tensor( [[cos_t, -sin_t, 0], [sin_t, cos_t, 0]], dtypesino.dtype, devicesino.device ) rot_mat_batched rot_mat.unsqueeze(0).repeat(B, 1, 1) ones torch.ones((B, 1, H, rec_size), devicesino.device) grid F.affine_grid(rot_mat_batched, ones.shape, align_cornersTrue) # 将一维投影向量平铺成二维平面旋转到对应角度 proj_vec proj_filtered[:, :, :, i] # (B,1,H) proj_2d proj_vec.unsqueeze(-1).expand(-1, -1, -1, rec_size) rotated F.grid_sample( proj_2d, grid, modebilinear, padding_modezeros, align_cornersTrue ) recon rotated[:, 0] # 角度密度修正乘以 pi / A recon recon * (math.pi / A) return recon.unsqueeze(1)逻辑说明前向投影的复杂度是 A 次旋转 积分A180 时单张 512x512 图像的投影在 GPU 上约几十毫秒训练可以接受。滤波反投影里的 ramp filter 我做了幅度修正abs(freq) * H的系数不能省否则重建图整体会偏暗或者偏亮。滤波后每个角度一维投影平铺成二维平面旋转累加得到重建图。最后的pi / A是平行束 FBP 的归一化系数。参数说明n_angles180是最常用配置对应医学 CT 扫描的常见角度抽样。如果数据里 sinogram 是 360 度范围记得把角度序列改成 0~360。自定义视角时n_angles越大条纹伪影越少但训练显存和时间线性增加。4.2 数据一致性层 DDC 的 PyTorch 实现DDC 层是 InDuDoNet 最有工程价值的部分。网上很多源码实现把它放在最后做个简单的加法但没有把 Radon 变换接上等于让两个域在中间断开。我这里给出一个带残差的完整版本class DataConsistencyLayer(torch.nn.Module): 双域数据一致性层图像域输出回到投影域与原始投影做残差约束。 def __init__(self, alpha0.2, n_angles180): super().__init__() self.alpha torch.nn.Parameter(torch.tensor(float(alpha))) self.radon RadonProjection(n_anglesn_angles) self.fbp FilteredBackProjection(n_anglesn_angles) def forward(self, sino_input, sino_enhanced, img_out): # sino_input: 原始低剂量 sinogram # sino_enhanced: 投影域子网络输出 # img_out: 图像域子网络输出 x_init self.fbp(sino_enhanced) # 投影域结果重建图 x_refine x_init img_out # 残差式图像细化 sino_refine self.radon(x_refine) # 重新投影回 sinogram # 一致性残差修正量 alpha * (原始 - 重投影) y_dc sino_enhanced self.alpha * (sino_input - sino_refine) x_dc self.fbp(y_dc) return x_dc, y_dc逻辑说明这里的img_out来自图像域子网络的输出但实际落地时我会把图像域 U-Net 的输出设计成残差模式即 U-Net 学的是修复增量而不是直接输出完整图像。x_init img_out合成最终图像后再做 Radon计算原始投影与重投影的差异。这个差异乘上可学习的 alpha 后回到 sinogram 域保证投影域信息没有被丢弃。alpha 初始 0.2 的经验是低剂量噪声越强alpha 越要往小调。alpha 太大时网络会过度相信原始投影中的噪声样本输出图像会有高频抖动alpha 太小时 DDC 层退化成恒等映射双域约束消失。训练过程中可以把 alpha 的输出值打出来观察它在 0.1~0.3 之间浮动都是正常的。4.3 联合损失函数与配比感知损失之外还要盯住投影域InDuDoNet 的训练损失通常由几部分构成图像域重建损失、投影域重建损失、以及可选的感知损失。我只用 L1 作为图像域重建损失因为 L2 会让输出图像过于平滑低剂量去噪任务里 CT 值偏置更容易出现在 L2 优化的结果里。损失项监督对象权重备注L1_img网络最终图像 vs 正常剂量图像1.0主损失SSIM_img网络最终图像 vs 正常剂量图像0.1提升结构相似度L1_sinoy_dc vs 干净 sinogram0.05保证投影域不去过头L1_img_slicex_dc vs 正常剂量图像0.1DDC 输出再补一次监督投影域损失权重不能超过 0.1否则网络会把注意力放在把 sinogram 修漂亮而 sinogram 里的微小误差经 FBP 放大后反而毁掉图像边缘。我在训练时每 5 个 epoch 会把各损失项分别打印出来如果 L1_sino 下降很快但 SSIM 停滞说明图像域子网络没有吃到有效梯度优先检查 DDC 层里 alpha 是否退化成了 0或者 Radon/iRadon 的角度方向是否颠倒。图像域 L1 和 SSIM 的配比按上面表格基本够用。感知损失不是必须的加不加对最终 PSNR 影响不大但会让纹理看起来更自然——如果做主观评价建议保留一个低权重的 VGG 感知损失。5. InDuDoNet 训练常见问题与排查显存失控、结果发糊与投影域噪声5.1 显存爆掉为什么 sinogram 的 batch size 不能照搬普通 2D 网络现象把 batch size 设成 8训练刚开始就 OOM报错提示是 CUDA out of memory。原因InDuDoNet 虽然看起来只有两个子网络但 DDC 层里的 RadonProjection 和 FilteredBackProjection 会在计算图中保存每个角度的中间旋转图180 个角度的变量叠加起来比网络本身还占显存。尤其grid_sample的梯度回传需要保存输入和网格这块开销经常被忽略。解决batch size 降到 2~4patch 从 512 降到 256。还不够的话用梯度累积模拟大 batch# 以 batch2、累积 4 步模拟 batch8 的梯度 python train_indunonet.py --batch-size 2 --grad-accum-steps 4 --patch-size 256 --amp同时在训练脚本里加两行配置torch.backends.cudnn.benchmark True把 AMP 打开。AMP 对 InDuDoNet 尤其有效因为 Radon 变换里的矩阵乘法占了大量算力半精度能把显存占用砍掉将近一半而且重建精度损失肉眼不可见。如果开了 AMP 仍然爆显存把投影域 U-Net 的起始通道从 32 降到 24换来的是投影域的表达能力略降一般可接受。5.2 loss 降了但图像发糊雾感从哪里来怎么定位现象训练 loss 曲线很漂亮PSNR 也不低但输出图像像蒙了一层薄雾组织边缘对比度明显下降。原因这是低剂量去噪复现里最典型的翻车现场。第一个嫌疑是 DDC 层的 alpha 过大压制了网络学习高频细节的能力第二个嫌疑是滤波反投影里的 ramp filter 幅度不对重建图像本身就被压低第三个原因是投影域 L1 权重过高网络学成了“把 sinogram 做均值模糊”就能降低损失。排查时我把网络中间输出可视化保存下来def dump_debug_images(model, sino_low, tag): 保存 DDC 层前后的中间结果用于判断发糊出自哪一步。 model.eval() with torch.no_grad(): sino_enh model.proj_unet(sino_low) x_init model.fbp(sino_enh) x_refine model.img_unet(x_init) y_dc, x_dc model.ddc(sino_low, sino_enh, x_refine) torch.save({ sino_enh: sino_enh, x_init: x_init, x_refine: x_refine, y_dc: y_dc, x_dc: x_dc, }, fdebug_{tag}.pt)然后在本地写个小脚本把每个张量做 min-max 归一化并存成 PNG。如果x_init已经发糊问题在投影域子网络和 FBP如果x_init清晰但x_final发糊问题在图像域子网络。我踩过的坑是 x_init 清晰、最终输出发糊最后定位到图像域 U-Net 深度太深3x3 卷积堆叠让感受野覆盖了太多无关联组织改成 3 层之后雾感明显缓解。5.3 测试时 Filtered Back Projection 与训练时不一致伪影从哪冒出来现象训练集指标很好测试集上一张真实 DICOM 图输出出现放射状发散伪影甚至图像方向都反了。原因训练时用的自定义 Radon 模块和测试时用的 skimage 或 ASTRA 重建模块角度约定不一致。最常见的是角度起点不一致有的库 0 度对应水平投影有的库 0 度对应竖直投影还有的 sinogram 第一行是 0 度有的最后一行才是 0 度。InDuDoNet 在训练时学了自定义模块的角度约定测试用另一个库就等于输入分布错位。解决训练和测试必须共用一个 Radon/iRadon 实现。我建议在项目里固定写一个geometry.py里面只保留自己实现的模块不要在推理脚本里用skimage.transform.iradon。每改一次角度参数就用 Shepp-Logan 模体做一次 round-trip 校验def validate_geometry(): from skimage.data import shepp_logan_phantom phantom shepp_logan_phantom().astype(np.float32) phantom torch.from_numpy(phantom).unsqueeze(0).unsqueeze(0) radon RadonProjection(n_angles180) fbp FilteredBackProjection(n_angles180) sino radon(phantom) recon fbp(sino) error (recon - phantom).abs().mean().item() print(fround-trip MAE: {error:.4f}) # 视觉检查recon 纹理应与 phantom 一致无方向偏转参数说明round-trip MAE对 512x512 的 phantom 一般应在 0.05 以下超过 0.1 就要检查 ramp filter 的系数或者旋转方向是否反了。我在一次换 PyTorch 版本后遇到过align_cornersTrue/False改变导致的角度偏移 0.5 像素当时就是靠这个校验函数发现的。6. 把复现的 InDuDoNet 变成可交付的推理管线torchscript 导出与验证复现的价值在能稳定地跑新数据。训练完模型后我用 torch.jit.trace 把整个推理链路固定下来避免训练脚本和推理脚本因为代码版本不同产生微妙差异。trace 时要注意输入尺寸固定InDuDoNet 的 DDC 层包含角度循环动态长度输入会触发重新 trace有时会错误地把已有权重覆盖掉所以我在 trace 前会先把批处理尺寸固定成 1x1x512x512。model.eval() dummy_input torch.randn(1, 1, 512, 512) traced torch.jit.trace(model, dummy_input) traced.save(indunonet_ct.pt)若要再提速可以用半精度推理和 ONNX Runtime 替代原生 PyTorch。半精度推理前必须验证重投影再重建的误差有些模型在 fp16 下 DDC 残差 loss 下降不稳定最好做一个 PSNR 对比实验再决定是否默认开启。推理管线里我每次还会带一个后置校验输出图的 CT 值分布范围是否和输入一致如果均值漂移超过 20 HU基本可以断定推理路径用错了重建模块。这些细节在复现 InDuDoNet 的源码工程里很容易忽略。我自己第一次交付时就是忘了训练和推理用不同 Radon 实现导致同一份数据线上线下结果对不上此后固定了一套 geometry 模块任何实验都先跑 phantom 校验再跑临床数据。希望这些踩坑记录能帮你少走一段弯路。复现 InDuDoNet 并最终跑通一份 Python 实现源码的过程本质上就是在磨你对双域数据流的掌控力。把这套模块沉淀成你自己的工具库以后再做其他 CT 重建网络的对比实验会顺畅很多。本文还有配套的精品资源点击获取