CT重建看起来是个硬件问题实际上是个数学问题。做医学图像处理或者工业无损检测的朋友应该都有体会投影数据转一圈采回来算法才是决定图像质量的关键。传统滤波反投影FBP在数据完整时又快又稳但一旦投影角度减少、剂量降低星状伪影和噪声放大能把图像变得没法看。这时候就要转到优化视角——把“重建”写成“求解一个带正则项的优化问题”用约束把先验信息塞进去。这篇文章就是围绕这套思路展开的实战笔记主角是原始对偶算法同时配套对比LASSO、TV正则化、ADMM、FISTA这几种常见解法所有代码均已用Matlab实现并整理成可直接运行的源码包。适合刚接触压缩感知重建、想搞懂算法原理与代码对应关系的同学也适合已经有了FBP基础、想进阶到迭代重建的工程师。我最早接触这个方向时最困惑的不是公式推导而是不知道这些优化算法到底谁该用在什么位置。后来把CTPD压缩感知CT原始对偶重建完整跑通又对比了ADMM和FISTA在一组真实投影数据上的表现才慢慢理清楚。这篇就把整个过程写出来从数学模型、算法迭代机制、Matlab代码骨架到参数调节和踩坑经验尽量让读者照着走一遍就能上手。1. 为什么CT图像重建要写成优化问题1.1 从“反投影”到“解方程”把CT重建想象成一个解方程的过程是最直观的理解方式。离散化之后成像过程可以写成y A x e其中 x 是待重建的图像向量长度是 N²假设图像是 N×N 像素y 是探测器采集到的投影数据长度是“探测器单元数×投影角度数”A 是系统矩阵每一行描述一条射线穿过图像时的加权积分关系e 是测量噪声。在Matlab里A 和 A 的转置操作可以用 radon 和 iradon 模拟但要注意 iradon 默认带滤波后面我会专门说这个坑。如果显式地构造 A那么 A 是一个巨型稀疏矩阵尺寸大约是 “投影数据量×像素数” 的量级。2048×2048 图像配上几百个投影角度时A 的行数轻松上百万直接存满矩阵几乎不可能所以工程上必须用函数句柄或者稀疏矩阵来操作。FBP的直接思想是用滤波反投影来近似 A 的转置然后叠加各个角度。它本质上是把反投影当成重建结果这种方式在投影数据完备时没问题但角度不足时A 是列数远大于行数的欠定矩阵方程有无穷多个解反投影给出的只是其中一个带有大量伪影的解。迭代重建的思路就是在这些解里挑一个符合“图像先验”的。1.2 欠定问题与稀疏先验CT重建最常见的场景是稀疏角度或者低剂量扫描。前者是故意减少投影帧数以降低扫描时间或辐射剂量后者是每个角度下剂量变低导致噪声变大。两种场景下采集到的信息都不足以唯一确定像素值这就必须引入正则化项。压缩感知理论给了我们一个强有力的工具如果信号在某个变换域里是稀疏的那么从远少于奈奎斯特采样数的测量中可以精确重建该信号。CT图像的梯度域天然稀疏——图像内部大范围区域灰度变化平缓只有器官边界、骨骼边缘这些地方梯度值较大。这就直接指向了全变分Total Variation, TV正则化。TV正则化的数学形式分两种各向异性TV(x) Σ |∇x_i|₁ Σ (|D_x x| |D_y x|)各向同性TV(x) Σ sqrt((D_x x)² (D_y x)²)各向同性的几何意义更符合“边缘保持”它在梯度幅值上做 L1 惩罚会同时惩罚水平和垂直方向的梯度重建出的边缘更平滑、没有方向性偏差。各向异性实现起来简单一点但对角方向的边缘会有些失真。实际工程里我默认用各向同性后面代码也按这个写。1.3 从L2到L1再到TVLASSO问题可以写成min_x 0.5 ||A x - y||₂² λ ||x||₁这里的 L1 范数促使人 x 本身稀疏。如果图像本身不稀疏就要先做变换。医学图像在小波域、曲波域里确实更稀疏但变换会引入额外的计算开销而且在TV已经能给出不错结果时单纯L1变换域并不总是最优。L2正则Tikhonov做约束会出现过度平滑因为 L2 对梯度接近于零的小值不敏感不会强制稀疏性图像边缘被均化。L1 反而保留了尖锐边缘。但从实际效果看直接对像素值做 L1 并不适合大多数CT图像——CT图像的像素值本身并不稀疏稀疏的是梯度。这也是为什么 TV 正则化在CT重建里地位这么高它把“梯度稀疏”这个物理性很强的先验直接编进了目标函数。所以最终我们面对的问题形态通常是min_x 0.5 ||A x - y||₂² λ TV(x)接下来问题就是这个非光滑、不可导的目标函数怎么高效求解FISTA、ADMM、原始对偶算法给出了三条不同的路径。2. 四类核心算法的迭代机制拆解2.1 FISTA加速近端梯度的正确打开方式FISTAFast Iterative Shrinkage-Thresholding Algorithm是 ISTA 加上 Nesterov 动量加速的版本。ISTA 的思路很朴素因为 TV 项不可导但它有近端算子proximal operator所以把目标函数拆成两部分——光滑的数据保真项 f(x) 0.5||Ax-y||²和非光滑的正则项 g(x) λ TV(x)。每一步迭代就是x_{k1} prox_{tg}(x_k - t ∇f(x_k))其中 t 是步长∇f(x_k) A^T(Ax_k - y)。prox 运算可以理解成“在正则项约束下做一次近端投影”。对于 L1 正则prox 就是软阈值函数对于 TVprox 是一个去噪子问题没有直接的闭式解需要内层迭代。FISTA 的改进是在每一步加入外推z_k x_k ((k-1)/(k2)) (x_k - x_{k-1}) x_{k1} prox_{tg}(z_k - t ∇f(z_k))这个外推系数是 Nesterov 动量的一种形式它让收敛率从 ISTA 的 O(1/k) 提升到 O(1/k²)。直观感受是同样迭代100步FISTA 的重建质量明显好于 ISTA尤其在迭代初期差距特别大。FISTA 的关键参数是步长 t必须满足 t ≤ 1/L其中 L 是 ∇f 的 Lipschitz 常数也就是 A^TA 的最大特征值。L 估计小了算法会震荡甚至发散估计大了收敛变慢。工程上用幂迭代法对 A^TA 估计 L20次迭代就能得到足够准确的值。2.2 ADMM把大问题拆成三个小问题ADMMAlternating Direction Method of Multipliers的出发点完全不同。它引入辅助变量把原来纠缠在一起的数据项和正则项分离开。以 TV 为例写成min_x 0.5 ||A x - y||₂² λ ||z||_TV约束 x z然后构造增广拉格朗日函数L_ρ(x, z, u) 0.5||Ax-y||² λ TV(z) (ρ/2)||x - z u||²迭代就变成了三个子问题x 更新求解一个 L2 正则最小二乘问题对图像重建来说这步是计算大头的部分z 更新对 x u 做一次 TV 去噪等价于 prox_{λ/ρ} TV 运算u 更新对偶变量累加残差 u u (x - z)x 更新这步需要解 (A^TA ρI)x A^Ty ρ(z - u)。当 A 是大规模矩阵时不能直接求逆要借助共轭梯度CG迭代求解。如果 A 是显式稀疏矩阵也可以预分解。ADMM 的迷人之处在于三个子问题各自独立、交替求解天然适合分块并行。这也是为什么“并行 ADMM”会成为热词——大数据规模下图像可以分块每一块的子问题同时解再同步对偶变量。工程实现时要注意 ρ 的选取ρ 太小对偶收敛慢ρ 太大x 和 z 的耦合太强迭代容易出现锯齿状震荡。实际中 ρ 取 0.5~2 之间比较稳后面细说。2.3 原始对偶算法从鞍点视角看重建原始对偶算法以 Chambolle-Pock 为代表是处理 TV 这类非光滑正则最优雅的方案之一。核心思路是把原问题改写成鞍点问题min_x max_y ⟨Ax, y⟩ g(x) - f*(y)对 CT 重建来说f(v) 0.5||v - y||²其 Fenchel 共轭为 f*(y) 0.5||y||² ⟨y, y_obs⟩。Chambolle-Pock 的迭代格式是y_{k1} prox_{σf*}(y_k σ A x̄_k)x_{k1} prox_{τg}(x_k - τ A^T y_{k1})x̄_{k1} x_{k1} θ(x_{k1} - x_k)这里的 x̄ 是外推项θ 一般取1。σ 和 τ 是步长只要满足 στ L² 1L 是 ||A||₂算法就能保证收敛。数据项的对偶 prox 有一个闭式解prox_{σf*}(z) (z - σ y_obs) / (1 σ)计算极其简单。正则项的 prox 同样是 TV 去噪子问题需要调用内层算法。因为整个算法不需要像 FISTA 那样显式计算 A^TA 的范数只需求 ||A||₂且外推公式天然内建所以对 CT 重建这类问题实现起来非常顺滑。这也是标题里 CTPD 这套组合能成立的原因之一原始对偶框架对 TV 这类非光滑正则特别友好。2.4 选型逻辑先别急着写代码四类算法里FISTA 实现简单、参数最少但对 Lipschitz 常数敏感ADMM 拆解清晰、方便并行但需要调 ρ原始对偶框架适合 TV 类非光滑正则不需要显式 Lipschitz 常数但需要内层去噪子过程配合普通 ISTA 基本用于教学验证工程上已经被 FISTA 取代。我的建议是如果你只是想快速验证一个重建思路用 FISTA如果未来要扩展到大规模并行用 ADMM如果确定用 TV 正则且想要稳定收敛直接上原始对偶。下面一节给出三套代码骨架可对照复制。3. Matlab实现要点与核心代码骨架3.1 系统矩阵A的两个实现陷阱Matlab 里做 CT 模拟最顺手的自然是 radon/iradon。但这里有两个大坑第一iradon 默认是带斜坡滤波Ram-Lak的那是 FBP 的转置近似不是严格的 A^T。所有迭代重建算法里A^T 必须是 A 的严格伴随算子否则梯度方向错误迭代一定发散。正确的做法是ATfun (y) iradon(y, angles, linear, none)把滤波去掉只保留反投影。第二不要轻易显式构造 A 大矩阵。我见过不少同学在 256×256 图像上尝试用 sparse 存 A结果内存直接爆掉。我的做法是把投影函数封装成两个函数句柄N 128; % 图像尺寸 angles 0:3:177; % 60个投影角度 Afun (x) radon(reshape(x, N, N), angles); ATfun (y) reshape(iradon(y, angles, linear, none), [], 1);这样 A 的“乘法”和“转置乘法”都变成了函数调用内存开销被压到最低。代价是无法使用预分解等矩阵分解技巧但配合 CG 或原始对偶迭代是足够的。3.2 FISTA-TV与ADMM-TV代码骨架FISTA-TV 的主循环function x fista_tv(y, Afun, ATfun, N, lambda, maxIter) x zeros(N*N, 1); z x; t 1; % 步长由Lipschitz常数决定 for k 1:maxIter grad ATfun(Afun(z) - y); x_new z - t * grad; x_new prox_tv(reshape(x_new, N, N), t*lambda, 20); x_new x_new(:); alpha (k - 1) / (k 2); z x_new alpha * (x_new - x); x x_new; end end其中 prox_tv 是 TV 去噪子问题可以用 Chambolle 投影算法迭代也可以直接调用现成实现。内层迭代次数 20 次左右足够太少会导致外部计算不稳定太多则拖慢整体速度。ADMM-TV 的主循环function x admm_tv(y, Afun, ATfun, N, lambda, rho, maxIter) x zeros(N*N, 1); z x; u x; for k 1:maxIter rhs ATfun(y) rho * (z - u); x pcg((v) ATfun(Afun(v)) rho*v, rhs, 1e-6, 100); z prox_tv(reshape(x u, N, N), lambda/rho, 20); z z(:); u u x - z; end end这里的 x 更新用 pcg 解正规方程是个稳定选择。如果 A 是显式稀疏矩阵可以提前分解 A^TA ρI。3.3 原始对偶CT重建的代码骨架原始对偶算法的代码骨架更简洁function x primal_dual_ct(y, Afun, ATfun, N, lambda, sigma, tau, maxIter) x zeros(N*N, 1); xbar x; p zeros(size(y)); theta 1; for k 1:maxIter p (p sigma * (Afun(xbar) - y)) / (1 sigma); x_new x - tau * ATfun(p); x_new prox_tv(reshape(x_new, N, N), tau*lambda, 20); x_new x_new(:); xbar x_new theta * (x_new - x); x x_new; end end关键点是步长选择σ 和 τ 要满足 στ ||A||₂² 1。通常在估计 ||A||₂ 时可以做个近似比如 σ τ 0.1 在多数模拟CT问题上都能稳定收敛。如果想更精确可以做几次幂迭代求 ||A||₂。注意 p 的对偶更新里直接包含了 y 的信息这就是前面说的f*的prox闭式解。3.4 迭代停止与结果评估迭代重建最怕跑着跑着看不到收敛信号。我一般用两个判据一个是相邻两次迭代的图像变化量 norm(x_new - x)/norm(x) 小于某个阈值比如1e-5另一个是数据残差 norm(Afun(x) - y) 的下降曲线不再明显变化。两者结合能有效判断是否已经收敛。重建质量的定量评估也不能只看视觉效果。常用的指标是 MAE平均绝对误差、PSNR峰值信噪比和 SSIM结构相似性。在模拟数据里因为有无噪声的 Ground TruthPSNR 和 SSIM 能相对客观地反映算法好坏。实际使用中我建议对同一组数据跑完 FISTA、ADMM、原始对偶三套然后保存每步的残差曲线和最终重建图做成一张对比图放在结果报告里说服力最强。4. 实测对比收敛速度、重建质量与参数敏感性4.1 同一组投影数据上的表现我在 Shepp-Logan 幻影上做过一组实测图像尺寸 128×128投影角度 60 个角度间隔 3 度加入了1%的高斯噪声。用 FISTA、ADMM、原始对偶三套算法各跑100次迭代结果记录如下算法100次迭代耗时PSNR(dB)SSIM残差下降趋势FISTA-TV约2.1s28.30.86前40次显著后程趋缓ADMM-TV约3.8s29.10.88前20次很快之后平稳原始对偶约3.2s29.50.89全程平稳无震荡从表里可以看到FISTA 虽然单次迭代快但收敛到接近同样质量需要更多轮次ADMM 前期下降猛烈但最后的精细修正略慢原始对偶在整体质量和稳定性上最均衡。这个结果不意外——原始对偶算法对 TV 项的耦合方式最自然不会因为外推过度而产生伪影振荡。4.2 参数敏感性分析lambda 是 TV 权重的核心它对重建结果的影响最大。lambda 太小正则化不足重建图像噪声和伪影明显lambda 太大图像被过度平滑细小结构丢失。根据实测128×128 图像、数据归一化到 [0,1] 的情况下lambda 取 0.01~0.05 是安全区间。最好先跑一个快速的小迭代数实验观察残差和图像平均梯度随 lambda 的变化曲线在拐点附近取值。ADMM 的 rho 对收敛速度影响大但对最终结果影响相对小。rho 调太小0.01以下时对偶变量 u 累积缓慢收敛特别慢rho 调太大10以上时图像容易出现“棋盘格”式振荡。我的经验是 rho1 起步效果不满意再按 10 倍步进调整。原始对偶的 σ、τ 是成对出现的。当 τ 固定时σ 过大直接导致震荡发散。我习惯固定 τ0.1σ 从 0.1 开始下调观察前10次迭代的残差是否单调下降。也有一个更聪明的做法用 SART 或者 CFBP 的范数估计算法估算 ||A||₂然后 σ τ 1/||A||₂这时收敛最快且不会发散。4.3 什么场景该用哪个算法总结几条我在实际项目里总结的选型经验投影角度少稀疏角度但图像本身结构规则比如工业零件检测首选原始对偶TV因为 TV 能有效抑制有限角度伪影。投影数据噪声大比如低剂量CTFISTA 配合小 lambda 更稳健因为它对噪声的响应最“温和”。图像尺寸大、需要实时性ADMM 配合 GPU 并行化有天然优势因为三个子问题可以分别并行。只需要快速实验、不想调太多参数FISTA 的默认参数通常就能跑出一个可用的结果。5. 踩坑记录这些坑比算法本身更值得注意5.1 iradon默认滤波导致A^T错误的经典灾难这是我见过最多人踩的坑没有之一。很多人刚上手迭代重建直接用 iradon(y, angles) 作为 ATfun然后就发现迭代发散或者图像边缘出现奇怪的振铃。原因之前说过默认的斜坡滤波是 FBP 的组成部分不是严格的 A^T。A^T 应该只做反投影不滤波。改成 none 之后算法立刻恢复正常。这个坑我建议所有读者遇到第一个奇怪现象时先检查这行代码。5.2 Lipschitz常数估计不准的连锁反应FISTA 对 Lipschitz常数 L 很敏感。L 估计偏大会导致步长过小收敛速度被严重拖慢表面上看像是“算法跑得慢”实际上是步长被误杀L 估计偏小则直接造成震荡。幂迭代法估计 L 的代码很简单但注意要作用在 A^TA 上而不是 A 上。如果原始问题有条件数很大的 A幂迭代的收敛速度也偏慢可能需要 50 次以上迭代。此时我干脆换用原始对偶算法省去这个烦恼。5.3 lambda、rho、sigma的调参顺序调参顺序非常重要乱调会让人完全失去方向。我的顺序是先固定 rho1sigmatau0.1只调 lambda找到能接受的图像质量。固定 lambda再调 rho 或 sigma/tau优化收敛速度。最后微调迭代次数和内层 TV 去噪迭代次数。先调正则强度再调算法参数。反过来调很容易陷入“质量不好→疯狂调参数→参数互相干扰→质量更差”的怪圈。5.4 初始化与停机准则的隐性影响初始化用零向量最保险。有人觉得用 FBP 结果初始化能让算法更快收敛实际试下来反而容易出问题FBP 的伪影分布和 TV 正则的梯度结构耦合会让算法陷入一个“伪影看起来像边缘”的局部结构导致最终结果比零初始化更差。停机准则方面只看残差容易误判因为残差下降到一个平台后图像可能还在缓慢变化。要同时监视图像变化量当两者都进入平台时再停。对于在线重建场景还建议加一个最大迭代次数的硬限制防止单帧耗时过长。6. 最后再分享一个小技巧我一般会在实际项目里对同一组数据同时跑 FISTA、ADMM、原始对偶三套算法但只让其中一套作为线上默认其余作为对照。这么做的好处是换数据、换扫描参数时可以通过对比快速发现是不是算法本身出了问题还是数据、参数的问题。有一次我在一批工业扫描数据上发现 ADMM 结果带条纹对照原始对偶结果正常排查了半天发现是 rho 忘记随数据噪声水平调整了。如果没有对照组这类问题会浪费大量时间。如果是从零开始的项目建议先把原始对偶算法跑通因为它对参数最宽容、对 TV 正则的兼容性最好。跑通之后再扩展 FISTA 或 ADMM这时候你对问题本身已经有手感后面的路会顺很多。