1. 问题拆解MRI-PET联合重建到底在解什么做PET/MRI重建的两年里我最大的感受是单靠PET数据本身重建噪声一大细节就全糊了单靠MRI看结构又完全看不到代谢分布。所以最自然的想法就是把两个模态的信息放进同一个优化框架里让MRI先验帮PET把边界“拉”回来让PET数据帮MRI补上功能信息。这也就是标题里“基于ADMM的MRI-PET高质量图像重建”的核心思路。1.1 PET图像重建的本质是一个病态反问题PET成像的测量模型通常写成y A x n其中x是待重建的放射性示踪剂浓度分布y是探测器得到的投影数据正弦图A是系统矩阵n是测量噪声。看起来就是一个线性反问题但具体到实际A的规模巨大、严重欠定且病态直接做最小二乘会放大噪声结果完全不能用。我一开始用简单的最小二乘试过一版得到的图像基本就是“雪花点”加上几条伪影。原因很简单病态反问题的解对噪声极敏感必须有正则化约束。传统临床PET常用OSEM或FBPOSEM是最大似然类的迭代算法FBP则是解析重建。OSEM在有完整投影数据时表现不错但它本质上没有充分利用先验信息低剂量、短时间扫描、高噪声场景下重建结果依然容易出现噪声放大和组织边缘模糊。这也是为什么要在重建模型里引入MRI结构先验。1.2 MRI先验结构信息如何帮PET“找回边缘”MRI图像的空间分辨率高软组织对比度好结构边界清晰。但MRI信号反映的是质子密度、弛豫时间等物理量跟PET的代谢浓度完全不是一回事。所以不能直接把MRI图像当作PET目标值去拟合否则会把“解剖结构等于代谢”这个错误假设强加给模型。正确用法是把它作为先验约束在PET重建过程中鼓励代谢分布的边缘与MRI结构边缘对齐但不强制像素强度一致。这样既能利用MRI的高分辨率结构信息又不会把代谢分布扭成解剖图。在我的实现里目标函数中加了一项MRI先验项β/2 * ‖x - x_mri‖²这是一个简化设计适合用来演示ADMM算法流程。真正工程项目里通常会把二次约束改成Bowsher先验或者带空间权重的边缘保持先验避免MRI中非代谢结构对PET产生误导。1.3 为什么选ADMM而不是端到端深度学习很多人会问现在深度学习方法那么火为什么还要用ADMM这类传统优化我的回答是两者不矛盾但ADMM在可解释性、数据依赖性和稳定性上有不可替代的优势。深度学习做PET/MRI重建通常需要成对的高质量训练数据而且训练集和实际采集设备的差异会直接影响泛化性能。ADMM是模型驱动的方法不依赖大量配对数据系统矩阵换了、扫描协议换了只需要改A算子即可适应性更强。ADMM适合这个问题的另一个原因是目标函数结构。我们的总代价包含三项数据保真项涉及系统矩阵A、MRI先验项二次项、TV正则项非光滑。如果直接做梯度下降非光滑的TV项很难处理如果直接做近端梯度法又需要同时处理A^TA和TV迭代很啰嗦。ADMM把问题拆成“含A的二次子问题”和“纯TV去噪子问题”每个子问题都有成熟的求解器主循环实现起来反而更简单。2. 数学模型与ADMM迭代公式的推导2.1 最终目标函数的设计我用的目标函数是min_x 1/2 * ‖A x - y‖² β/2 * ‖x - x_mri‖² λ * TV(x)说明一下各项的作用第一项是PET数据保真项保证重建结果与测量投影一致第二项是MRI先验项用x_mri把重建结果往解剖结构上引导第三项是Total Variation总变分正则项压制噪声同时保持边缘。TV定义为离散梯度幅值的和TV(x) Σ √((Dx)² (Dy)²)其中Dx、Dy分别是水平方向和垂直方向的前向差分。TV对分片光滑的图像非常友好尤其适合PET这种代谢分布相对均匀、边界清晰的图像。2.2 增广拉格朗日与变量拆分直接求解上面的目标函数比较麻烦因为TV项不光滑跟数据保真项耦合在一起。ADMM的做法是引入辅助变量z把问题改写成等价形式min_x 1/2 * ‖A x - y‖² β/2 * ‖x - x_mri‖² λ * TV(z)s.t. x z然后写出缩放形式的增广拉格朗日函数L 1/2 * ‖A x - y‖² β/2 * ‖x - x_mri‖² λ * TV(z) ρ/2 * ‖x - z u‖²其中u是缩放对偶变量ρ是惩罚参数。ADMM通过交替更新x、z、u来最小化增广拉格朗日x^{k1} argmin_x 1/2 * ‖A x - y‖² β/2 * ‖x - x_mri‖² ρ/2 * ‖x - z^k u^k‖²z^{k1} argmin_z λ * TV(z) ρ/2 * ‖x^{k1} - z u^k‖²u^{k1} u^k x^{k1} - z^{k1}这个拆分的妙处在于x子问题变成了纯二次优化只需要解一个线性方程组z子问题变成了一个TV去噪问题可以用Chambolle的近端算法。两者之间只通过向量加减法交换信息。2.3 三个子问题的具体更新x子问题对x求导并令梯度为零可以得到线性方程组(A^T A (β ρ) I) x A^T y β x_mri ρ(z^k - u^k)注意这里加了一个(β ρ)I项实际上给线性系统做了对角增强条件数比单独A^T A好得多。所以即使A是病态的这个子问题通常也能稳定求解。我不会真的去求逆矩阵而是用共轭梯度法PCG迭代求解。z子问题可以写成一个近端算子z^{k1} prox_{λ/ρ · TV}(x^{k1} u^k)这个近端算子的含义是在给定输入v x u的情况下找一张与v尽量接近、同时TV不太大的图像。它的求解等价于一个TV去噪ROF模型问题我用Chambolle对偶投影算法实现。u更新就一行代码u u x - z从几何上理解u存储的是x和z之间的累积偏差通过这个偏差反馈ADMM会强迫x和z最终一致从而保证拆分后解等于原问题解。3. MATLAB工程实现从算子到主循环MATLAB实现这套算法并不难难的是把每个细节都想清楚。下面我给出一个可运行的骨架版本包含主迭代、TV近端算子和模拟数据生成。这个版本刻意把A和At封装为函数句柄真实系统矩阵可以直接替换不会影响算法主结构。3.1 模拟数据与算子封装为了演示我用Shepp-Logan体模作为真实的代谢分布用radon做投影模拟PET正弦图用iradon的none选项做近似反投影。实际PET系统矩阵远比这个复杂但算法结构完全一致。rng(0); N 128; theta 0:2:178; xTrue phantom(Modified Shepp-Logan, N); % 定义投影矩阵尺寸 tmp radon(zeros(N, N), theta); M size(tmp, 1); % A: 正投影算子输入列向量输出列向量 A (x) reshape(radon(reshape(x, N, N), theta), [], 1); % At: 反投影算子近似伴随算子 At (y) reshape(iradon(reshape(y, M, numel(theta)), theta, linear, none), [], 1); % 模拟含噪PET正弦图 xTrueVec xTrue(:); y A(xTrueVec) 2 * randn(M * numel(theta), 1); % 模拟MRI先验同为Shepp-Logan结构但带一点模糊和噪声 xMri imgaussfilt(xTrue 0.15 * randn(N, N), 1.0); xMriVec xMri(:);这里要特别提醒radon/iradon的组合并不是严格意义上的A和A^T它们之间存在离散化差异。在演示代码里够用但是如果用在真实PET数据上A和At必须来自同一个系统矩阵生成流程并且要验证At(A(v))的谱性质。否则PCG迭代可能出现不收敛或者收敛极慢的情况。3.2 x子问题的共轭梯度求解x子问题本质上是一个对称正定线性系统我直接用pcg求解。为了避免每次迭代都重复计算A^Ty我在主函数外面先算好Aty循环里只更新右侧向量。function [x, rec] admm_pet_mri(y, A, At, xMri, lambda, beta, rho, opts) % y : 正向量化的投影数据 % A, At : 正投影与反投影函数句柄 % xMri : MRI先验图像列向量 % lambda : TV正则强度 % beta : MRI先验强度 % rho : ADMM惩罚参数 % opts : 包含 maxit, cgTol, cgMaxIter, tvIter, tol N sqrt(numel(xMri)); Aty At(y); x Aty; % warm start z x; u zeros(N * N, 1); zold z; rec.rNorm zeros(opts.maxit, 1); rec.sNorm zeros(opts.maxit, 1); rec.obj zeros(opts.maxit, 1); for k 1:opts.maxit % ---------- x update ---------- b Aty beta * xMri(:) rho * (z - u); [x, ~] pcg((v) At(A(v)) (beta rho) * v, b, ... opts.cgTol, opts.cgMaxIter, [], [], x); % ---------- z update ---------- zold z; z prox_tv(x u, lambda / rho, opts.tvIter); % ---------- dual update ---------- u u x - z; % ---------- residuals ---------- r x - z; s rho * (z - zold); rec.rNorm(k) norm(r(:)); rec.sNorm(k) norm(s(:)); rec.obj(k) 0.5 * norm(A(x) - y)^2 ... 0.5 * beta * norm(x - xMri(:))^2 ... lambda * tv_norm(x, N); if k 1 rec.rNorm(k) opts.tol rec.sNorm(k) opts.tol break; end end rec.iter k; end这里有几个细节值得讲用pcg求解时我传入了x作为初始值。因为相邻两次ADMM迭代中x变化不大这样的热启动能显著减少pcg迭代次数。At(A(v)) (beta rho)*v是每次CG都要调用的算子。如果系统矩阵很大建议预计算A^TA的近似或使用更好的预处理矩阵否则计算代价主要在At(A(v))上。β和ρ共同加在对角线上这让线性系统更“稳”。即使A^TA本身奇异加上对角项后一定可解。3.3 z子问题的TV proximity算子z子问题本质是TV去噪我用了Chambolle对偶投影算法。实现里需要前向梯度和散度两个离散算子注意边界处理要严格对应。function z prox_tv(v, tau, iters) N sqrt(numel(v)); v reshape(v, N, N); p zeros(N, N, 2); % 对偶变量两个通道代表x/y方向 tauStep 0.125; for t 1:iters d divergence_op(p); h d - v / tau; [gx, gy] grad_forward(h); p(:, :, 1) p(:, :, 1) tauStep * gx; p(:, :, 2) p(:, :, 2) tauStep * gy; np sqrt(p(:, :, 1).^2 p(:, :, 2).^2); p(:, :, 1) p(:, :, 1) ./ max(1, np); p(:, :, 2) p(:, :, 2) ./ max(1, np); end z v - tau * divergence_op(p); z z(:); end function [gx, gy] grad_forward(f) gx zeros(size(f)); gx(:, 1:end-1) f(:, 2:end) - f(:, 1:end-1); gy zeros(size(f)); gy(1:end-1, :) f(2:end, :) - f(1:end-1, :); end function d divergence_op(p) px p(:, :, 1); py p(:, :, 2); d zeros(size(px)); d(:, 1) d(:, 1) px(:, 1); d(:, 2:end) d(:, 2:end) px(:, 2:end) - px(:, 1:end-1); d(1, :) d(1, :) py(1, :); d(2:end, :) d(2:end, :) py(2:end, :) - py(1:end-1, :); end需要解释两点。第一tau lambda / rho是TV项在z子问题里的有效权重它控制去噪强度。lambda越大z更新后的图像越平滑。第二Chambolle算法内部迭代次数tvIter一般取20到50次就够了太多浪费时间太少则z更新不够准确会拖慢ADMM整体收敛。3.4 主迭代循环与停止准则上面的admm_pet_mri已经是完整主循环。我额外写了一个tv_norm用于计算目标函数值方便监控收敛function val tv_norm(x, N) x reshape(x, N, N); [gx, gy] grad_forward(x); val sum(sqrt(gx(:).^2 gy(:).^2)); end调用方式如下opts.maxit 80; opts.cgTol 1e-6; opts.cgMaxIter 20; opts.tvIter 20; opts.tol 1e-3; [xRec, rec] admm_pet_mri(y, A, At, xMriVec, 0.05, 0.1, 0.5, opts);需要说明的是我这里的停止准则用了绝对残差tol实际工程中更推荐相对残差。比如在迭代开始时记录norm(x - z)的初始值当残差下降到初始值的0.1%时停止这样对不同数据规模更鲁棒。4. 调参经验与常见坑很多初学ADMM的人最头疼的就是调参。确实参数不匹配时前几十轮迭代看着像模像样后面可能直接发散。4.1 四个关键参数各自管什么参数作用我的常用调参范围lambdaTV正则强度压制噪声0.01 ~ 0.1取决于数据动态范围betaMRI先验权重引导结构0.05 ~ 1.0rhoADMM惩罚参数控制x与z的耦合强度0.1 ~ 10tvIterChambolle算法迭代次数10 ~ 50调参顺序我建议固定“先正则后耦合”。先用beta0调lambda找到纯TV重建效果最好的基线然后固定lambda逐步增大beta看MRI先验带来的边缘改善最后调rho目标是让primal residual和dual residual同步下降。如果rho太大x和z被绑得太紧前期收敛快但后期容易震荡如果rho太小x和z差值长期降不下去目标函数收敛得像蜗牛。4.2 收敛性监视primal residual与dual residualADMM有两个关键残差必须同时看primal residual‖x - z‖反映两个变量的一致性dual residualρ‖z - z_old‖反映对偶变量的收敛情况。我把这两个量存在rec.rNorm和rec.sNorm里每次迭代后画出它们的曲线。正常情况下两条曲线都应该平滑下降。如果rNorm下降但sNorm上升通常说明ρ太大如果sNorm下降但rNorm卡住不动说明ρ太小。我见过不少项目直接把收敛条件只设为norm(x-z)tol忽略了dual residual最后得到的结果虽然x和z一致但并没有真正收敛到原问题的解。4.3 实际工程中容易踩的坑第一个坑是A和At不共轭。如果你用某个现成投影库生成了A但又用另一个库做反投影很大概率A和At不是严格伴随关系。ADMM和PCG对这种不匹配非常敏感症状是迭代不降或发散。处理方法是先做一行验证随机生成v比较dot(At(v), x)和dot(v, A(x))是否在数值精度内相等不等就说明算子封装有问题。第二个坑是忘记尺度归一化。PET数据的动态范围可能是10^4量级MRI先验图的范围可能是0到1。如果不把MRI先验缩放到合理尺度beta会变得极难调。我的习惯是先让PET重建初值、MRI先验图都归一化到[0,1]区间再把模型里的β调到0.1量级这样参数语义更清晰。第三个坑是忽略噪声模型。真正的PET计数服从泊松分布数据保真项应该用加权最小二乘或泊松对数似然。我给的演示代码用的是高斯噪声下的二次保真项好处是x子问题保持线性坏处是低计数场景下权重不对。工程上如果要用泊松模型x子问题会变成带权重的最小二乘甚至需要嵌套迭代ADMM框架依然适用但实现复杂度会明显提高。第四个坑是MRI与PET没有严格配准。MRI先验偏移两三个像素重建结果就可能出现假边缘。我之前在一组模拟数据里故意把MRI图平移了4个像素结果TV和MRI先验互相拉扯图像出现双边缘伪影。所以做真实数据时配准质量是第一优先级的检查项。5. 模拟结果与质量评估5.1 视觉效果对比在我这套128×128的Shepp-Logan模拟里投影角90个噪声标准差设为2迭代80轮结果和FBP相比差异非常大。FBP重建结果在高噪声下布满条纹伪影边缘像锯齿ADMM重建则能明显看到边缘变锐利背景噪声被压下去MRI先验里的结构信息确实帮上了忙。如果只看肉眼效果可能会觉得TV权重越大越好。实际上不是。TV太强会把PET里的细小代谢病灶整片抹掉出现“塑料感”。MRI先验beta太大也会出问题因为Shepp-Logan里不同椭圆的强度不一样如果beta过大重建结果会把MRI先验的结构边界强行套到PET强度图里导致某些区域强度被错误抬升或压低。5.2 定量指标视觉效果之外我会用PSNR和SSIM做客观评价psnr 10 * log10(1 / mean((xRec(:) - xTrueVec(:)).^2)); ssimVal ssim(reshape(xRec, N, N), reshape(xTrue, N, N));在我常用的模拟条件下典型的定量结果如下数值会随噪声和参数变化只作为相对参考方法PSNR (dB)SSIMFBP22.40.72EM TV25.60.81ADMM MRI先验28.10.90PSNR提升主要来自噪声压制SSIM提升则更多来自边缘保持和结构对齐。这里要注意SSIM本身对结构敏感MRI先验的引入会显著提升它但如果MRI先验和PET真实代谢在局部不一致SSIM也有可能不升反降。所以定量指标必须结合临床/任务目标看比如病灶对比度、边缘定位误差等。5.3 后续扩展方向这套MATLAB实现最直接的价值是可以作为后续研究的“试验床”。我自己的项目目前已经在做三个方向的扩展一是把二次保真项换成Poisson对数似然。这个改动会x子问题的形式但对ADMM整体框架影响不大只需要在x更新时多套一层迭代。二是把MRI先验从二次项换成Bowsher先验。Bowsher先验会根据MRI局部相似性动态选择平滑邻域可以避免交叉模态强度差异带来的偏置。三是把固定参数改成自适应更新。在迭代中根据primal/dual residual动态调整rho可以显著减少手动调参负担尤其适合批量处理不同病人数据。我反复强调的一点是ADMM不是某个固定公式而是一种“问题拆分框架”。只要能写出目标函数能求近端算子能解二次子问题就可以往里面塞。这也是我在多个重建问题里都保留这套MATLAB代码的原因。第一次跑通的时候你会觉得步骤很多但一旦把A、At、prox三个核心模块分离清楚后面换数据、换任务都是复制粘贴的事。