简介这份 MATLAB 实现包聚焦 ADMM 与 Total Variation 正则化在 CT 图像重建中的应用面向医学影像算法学习者和科研人员可改善传统滤波反投影FBP在低剂量或稀疏投影下的伪影问题同时借助 TV 正则化保持边缘、抑制噪声。算法将原问题拆分为数据拟合、TV 正则项与乘子更新三个子问题通过交替迭代逼近最优解既降低噪声又能保留细节。资源共 10 个文件以 .m 脚本为主体涵盖一阶/二阶 TV、去模糊、椒盐噪声、一维压缩感知等示例程序另附 README.md 与 Git 配置文件便于查阅说明与版本管理压缩包仅 10KB代码紧凑适合直接运行和修改。目前已有 511 人学习。通过阅读完整实现可理解 ADMM 迭代求解的交替优化逻辑并参照演示代码调整迭代次数、正则化参数等关键设置为后续 CT 重建研究提供扩展性较强的算法范例。1. 从滤波反投影到迭代重建为什么ADMM-TV能同时压噪声与保边缘CT在低剂量扫描或者稀疏角度采集时滤波反投影FBP重建出来的图像会出现明显的星芒状伪影这些细条纹理和真实解剖结构混在一起医生很难判断是病灶还是噪声。传统做法是加后处理滤波但滤完噪声边缘也糊了。ADMM-TV的思路完全不同它把“重建”本身当作一个最优化问题在目标函数里同时放数据拟合项和全变分正则项然后通过分裂变量交替迭代来求解。这套代码给我的感觉是单轮迭代的计算量只比FBP多几次卷积和阈值操作却能把条纹伪影压低一个量级边缘对比度还能保持住。它内部的MATLAB实现覆盖了去模糊、椒盐噪声恢复、二阶TV和一维压缩感知适合正在做图像重建算法对比或者想把传统解析重建方法迁移到优化框架下的工程师和研究组。2. ADMM-TV的算子分裂拆解与MATLAB文件包定位2.1 目标函数与变量分裂的关键CT图像重建要解的是线性逆问题已知投影数据 y希望恢复断层图像 x噪声和欠采样让 y Ax e 成为病态方程。FBP直接做逆投影在数据完备且信噪比高时够用一旦投影角度稀疏逆矩阵的奇异值接近零高频噪声会被放大。TV正则化的做法是在最小二乘项后面加一项梯度幅度的L1范数让重建图像尽量平滑但又允许强边缘跳变。TV项本身是凸但非光滑的直接对目标函数做梯度下降要么震荡要么收敛极慢。ADMM把原问题拆成两个能高效求解的子问题中间用线性约束把它们耦合起来这是这套资源的核心逻辑。具体到实现通常会引入辅助变量 z约束 x z然后写出增广拉格朗日函数数据拟合项 0.5 ||Ax - y||² 负责让重建结果符合测量值TV 正则项 λ ||∇z||₁ 负责抑制噪声与伪影约束 x z 通过乘子 u 和惩罚参数 rho 拉紧。迭代里交替更新 x、z、ux 子步骤处理带 L2 项的二次最小化z 子步骤处理 TV 去噪u 子步骤只管累积误差。这样 AV 项和 TV 项各自用自己最擅长的方式被求解而不是强行放在同一个代价函数里做梯度下降。2.2 压缩包内的文件功能对照拿到ADMM-Total-Variation-master.rar解压后先把 README 过一遍再按文件用途分组会比较快。文件解决的问题关键参数ADMM_DeblurTV.m线性模糊 高斯噪声下的图像去模糊lambda、rho、maxIterADMM_DeblurTV_Demo.m去模糊演示脚本包含焦散噪声添加和PSNR计算blur kernel、噪声强度ADMM_SnP.m椒盐噪声下的L1-TV恢复rho、lambda、软阈值系数DemoSnP_ADMM.m椒盐噪声演示跑完直接出对比图噪声密度、正则系数Demo_Generalized.m通用线性算子反问题演示适合CT重建A算子句柄、初值ADMM_DeblurTV_SecondOrderTV.m二阶TV去模糊减少阶梯伪影order、lambdaADMM_1D_CAPL1.m一维复值信号的L1稀疏恢复测量矩阵、稀疏度这份文件映射表在实际调试中很有用如果目标只是CT重建真正要反复改的是Demo_Generalized.m和ADMM_DeblurTV.m其他文件更像是边界情况下的对比实现。2.3 核心循环的MATLAB骨架把ADMM_DeblurTV.m的主循环剥出来结构大致是下面这段代码的样子这里省略了角标细节保留了框架% 核心ADMM-TV迭代骨架 % y: 观测图像A: 正投影算子AH: A的共轭算子 % lambda: TV正则化权重rho: ADMM惩罚参数 x AH(y); % 用反投影/共轭投影做初值 z x; % 辅助变量 u zeros(size(x)); % 拉格朗日乘子 for k 1:maxIter % x 子问题求解 (AA rho*I) x A(y) rho*(z - u) rhs AH(y) rho * (z - u); x pcg((v) AH(A(v)) rho * v, rhs, 1e-5, 30); % z 子问题对 x u 做一次TV去噪阈值由 lambda/rho 控制 z tv_denoise(x u, lambda / rho); % 乘子更新把约束误差累计到 u 上 u u x - z; end上面这段代码里pcg是预条件共轭梯度法用来解 x 子问题。当 A 是卷积算子时这个子问题可以换成 FFT 一步求解速度会快很多当 A 是 CT 系统矩阵时则每次迭代都要跑一次正反投影这也是整个算法最耗时的部分。lambda / rho是 z 子问题里 TV 去噪的阈值折算比值越大图像越平滑但边缘保留会变弱比值越小去噪能力越弱。rho还负责调节迭代稳定性和收敛速度一般需要实验确定。2.4 参数初始化的经验区间对于 256×256 的 CT 重构图像我一般先把rho设成1e-2lambda从1e-4到1e-2之间做对数扫描。rho太大时 x 子问题会变成完全依赖 TV 项的平滑滤波导致图像过磨rho太小时乘子更新又跟不上约束误差收敛非常慢。maxIter在去模糊 demo 里设为 80 足够但迁移到稀疏角度 CT 时 120 轮以上才稳定。3. 用Demo_Generalized把去模糊循环改写成CT重建3.1 椒盐噪声模型为什么要换L1保真项DemoSnP_ADMM.m演示的是椒盐噪声恢复这和 CT 重建里偶尔遇到的强脉冲伪影很像。椒盐噪声的像素值要么接近 0要么接近最大灰度如果继续用最小二乘保真项0.5||Ax-y||²离群点会在梯度里被平方放大TV 项再怎么压都压不干净。ADMM_SnP.m把保真项换成了 L1 范数目标函数变成||Ax-y||₁ λTV(x)。L1 保真项的 ADMM 会在原始变量之外再引入一个残差辅助变量残差子步骤用软阈值操作完成% L1保真项的残差更新rho1是残差惩罚参数 r max(abs(Ax - y u_r) - 1/rho1, 0) .* sign(Ax - y u_r); u_r u_r Ax - y - r;这里sign保证方向不变max(abs - 1/rho1, 0)把幅度小于阈值的残差直接清零。对这种非对称脉冲噪声L1 保真项能防止单个异常像素把重建带偏代价是收敛速度比 L2 慢通常要多跑 30% 的迭代次数。3.2 把去模糊算符替换成CT正反投影Demo_Generalized.m的设计目的就是承接这种替换。在去模糊场景里A 是模糊核卷积在 CT 里A 是 Radon 变换AH 是反投影算子。只要把 ADMM 内部的A和AH从矩阵/卷积函数句柄换成 CT 正反投影主循环代码一行都不用改。下面是一个可运行的 CT 重建示例我用 MATLAB 内置的radon和iradon作为算子把 Demo 里的模糊核换掉% CT重建用ADMM-TV替代FBP P phantom(256); % Shepp-Logan体模 theta 0:2:178; % 稀疏角度隔2度采样 y radon(P, theta); % 投影数据正弦图 A (x) radon(x, theta); % 正投影 AH (p) iradon(p, theta, Linear, None, 1, size(P,1)); % 无滤波反投影近似共轭 x0 AH(y); % FBP无滤波初值 opts.lambda 0.05; % TV权重需调节 opts.rho 0.02; % ADMM惩罚参数 opts.maxIter 100; % 这里把算子句柄和初值传给通用求解器 x_hat Demo_Generalized(A, AH, y, x0, opts);这里面几个参数容易踩坑iradon的None滤波器不保证严格满足共轭关系但作为 AH 的近似在 ADMM 里仍然可用因为迭代会逐步补偿系统误差。lambda要随角度稀疏度增大而增大角度从 0:2:178 变成 0:5:178 时我一般把lambda乘以 2 到 3 倍。rho则要保持在lambda的 1/2 到 1 倍左右否则 x 子问题和 z 子问题的尺度差太多会振荡。3.3 触发条件与迭代终止判断直接拿 FBP 做初值时前 10 轮迭代的残差下降非常明显到 40 轮后曲线趋于平缓。一个简单可靠的停止条件是在每轮计算||x - z||和||rho*(z-z_prev)||两个残差两个都小于初始值的1e-3就停这样可以避免多余的反投影计算。注意不要单独看||x - z||因为小rho会导致这个残差永久偏高此时对偶残差已经很小继续迭代只是在空转。4. 二阶TV与一维复稀疏不是只有一层差分才能保边缘4.1 一阶TV的阶梯伪影问题标准 TV 正则化里的差分是一阶差分它倾向于把重建图像中的平滑渐变区域压成分段常数块。在 CT 软组织切片上这种分段常数经常表现为“阶梯”状尤其在血管和脂肪组织过渡区最明显。ADMM_DeblurTV_SecondOrderTV.m出现的原因就是要在保持边缘的同时允许真实存在的灰度斜坡继续存在。二阶 TV 使用二阶差分算子惩罚的不是灰度跳变本身而是跳变的变化率。数学上看TV 从sum(|D x|)变成sum(|D₂ x|)其中D₂是二阶差分矩阵。这样一来边缘仍然能保持锐利但是边缘之间的过渡允许线性渐变不会被削成平顶。代价是 z 子问题不再是经典 TV 去噪而是一个高阶泛函最小化每次迭代要解更复杂的线性系统内存占用和耗时明显上升。如果只是验证效果可以在去模糊脚本里把ADMM_DeblurTV.m的主循环替换成二阶版本然后比较同一张图的横向灰度剖面% 比较一阶TV和二阶TV恢复后的灰度剖面 % x1: 一阶TV结果x2: 二阶TV结果 figure; plot(1:size(x1,2), x1(128,:), LineWidth, 1.5); hold on; plot(1:size(x2,2), x2(128,:), LineWidth, 1.5); legend(一阶TV, 二阶TV); xlabel(像素); ylabel(灰度值);运行后会看到一个明显区别一阶 TV 的剖面上会有平台和突变二阶 TV 的剖面多出斜坡过渡但真实边缘的锐度并不会下降太多。代价是二阶版本的参数调节更敏感lambda需要比一阶小一半左右否则迭代过程会出现低频振荡。4.2 一维复值信号中的L1恢复ADMM_1D_CAPL1.m处理的是另一类场景一维复值信号比如雷达回波或磁共振自由衰减信号。这类信号在稀疏字典下的表示系数既有实部又有虚部直接对模值加 L1 惩罚会丢失相位信息。这个文件里用的策略是将实部和虚部分开做软阈值或者对复数幅值做收缩具体取决于测量矩阵的结构。复值软阈值和实值软阈值有一点关键区别实值软阈值直接对每个数作用而复值软阈值通常写成max(|c| - tau, 0) * c / |c|这个操作同时保持复数符号方向并在模小于阈值时整体置零。下面这段代码展示了这种收缩操作function s soft_threshold_complex(c, tau) % c: 复数向量tau: 阈值 mag abs(c); s max(mag - tau, 0) .* (c ./ (mag eps)); end加eps是为了防止除零这是复值 L1 恢复里最常见的数值崩溃点。实际运行时如果看到NaN优先检查这里是不是没有加小常数保护。4.3 几种变体之间的梯度差异从ADMM_DeblurTV.m到ADMM_DeblurTV_SecondOrderTV.m再到ADMM_1D_CAPL1.m差分算子和保真项的范数都在变但 ADMM 的外层交替结构几乎没有变。调试时可以保留统一的打印函数把每轮的原始残差和对偶残差都记下来比较不同变体下残差下降的形状。二阶 TV 的残差曲线往往比一阶更平缓这是正常的不是发散一维复值问题则经常出现残差先快速下降再缓慢回升这是乘子累积过头的信号需要把rho调大一点。5. 收敛性验证与lambda自动选择调参前先看这两条曲线最后一篇我讲一个很实用的技巧用残差曲线判断迭代是否已经跑到头再用 L 曲线自动选正则化参数。很多人在 ADMM 地里跑完后只看 PSNR 峰值但lambda本身是靠网格搜索调出来的离开测试环境后并不稳定。先把残差曲线画出来。r_prim是原始残差r_dual是对偶残差两个都下降到初始值的千分之一时就可以认为收敛% 记录并绘制收敛曲线 for k 1:maxIter x_old x; % ... ADMM核心迭代 ... r_prim(k) norm(x(:) - z(:), 2); r_dual(k) norm(rho * (z - z_old), 2); end % 同时观察两条曲线防止单条误判 semilogy(r_prim, LineWidth, 1.5); hold on; semilogy(r_dual, LineWidth, 1.5); legend(r\_prim (原始残差), r\_dual (对偶残差)); xlabel(迭代轮数); ylabel(残差);如果发现r_prim一直下降但r_dual持平说明rho太小乘子更新跟不上反之r_dual下降但r_prim停滞则说明rho太大约束被拉得过死。修正rho之后再来选lambda。lambda的自动选择用 L 曲线法把保真项残差和 TV 项幅度作为一个二轴曲线的横纵坐标用不同lambda跑完迭代后连成一条 L 形曲线最模糊的拐角就是平衡点% L曲线选lambda lambdas logspace(-4, -1, 15); for i 1:length(lambdas) x_r admm_tv_solve(A, y, lambdas(i), rho, maxIter); fidelity(i) norm(A * x_r - y, 2); tv_norm(i) sum(sqrt(diff(x_r,1,1).^2 diff(x_r,1,2).^2), all); end plot(log(fidelity), log(tv_norm), o-);拐角通常在lambda跨越两个数量级的位置取拐点附近的fidelity和tv_norm再去跑最终版本。这样选出来的lambda不依赖测试图像的 PSNR 先验直接应用到同噪声水平的临床数据上重建切片的边缘梯度会更接近手头 CT 数据本身的分布。本文还有配套的精品资源点击获取