医学成像 TP4：2D PET Plug-and-Play 重建

实验背景与目标

正电子发射断层扫描（PET）通过检测放射性示踪剂发出的射线来重建人体内部的代谢活动图像。在实际临床应用中，为了降低患者所受的辐射剂量，往往需要减少注射的放射性药物量或缩短扫描时间，这直接导致采集到的投影数据（sinogram）中光子计数减少，统计噪声显著增加。

本实验的核心目标是探索Plug-and-Play（PnP）方法在PET图像重建中的应用，即将传统优化算法中的近端算子（proximal operator）替换为预训练的去噪器（如深度神经网络），从而在保持数据保真度的同时引入更强的先验信息。

实验环境与数据

首先使用EM（Expectation-Maximization）算法对两组数据分别进行100次迭代重建，得到基准结果。

【图片：Phantom】

Phantom图像结构细节完整，这是重建算法需要恢复的目标。

【图片：EM low dose】

低剂量EM重建结果呈现出严重的噪声，边界模糊

【图片：EM normal dose】

正常剂量EM重建结果的噪声水平明显低于低剂量重建。灰白质对比度更好，结构边界更清晰，但仍存在一定程度的噪声。这一结果作为本实验的参考标准，目标是通过PnP方法使低剂量重建达到接近此质量水平。

Part I: ADMM重建

ADMM算法原理

PET图像重建可以表述为以下优化问题：

\hat{\theta} = \arg\min_{\theta} \Phi(\theta, y) + \lambda U(D\theta)

其中\Phi(\theta, y)是数据保真项，描述重建图像\theta与测量数据y之间的一致性

交替方向乘子法（ADMM）通过引入辅助变量将原问题分解为若干子问题交替求解。对于上述问题，ADMM迭代格式为：

\theta^{(k+1)} = \text{prox}_{\Phi/\rho}(z^{(k)} - u^{(k)}) \quad \text{(S1)}

z^{(k+1)} = \text{prox}_{\lambda U(D\cdot)/\rho}(\theta^{(k+1)} + u^{(k)}) \quad \text{(S2)}

u^{(k+1)} = u^{(k)} + \theta^{(k+1)} - z^{(k+1)} \quad \text{(S3)}

Question 1: ADMM实现

见.ipynb

Question 2: 参数影响分析

\lambda参数的影响

在固定\rho = 10^{-6}、迭代次数n_{it} = 50的条件下，测试了\lambda \in \{0.0001, 0.001, 0.01, 0.1\}四个值对重建质量的影响。

【图片：ADMM-TV, lambda=0.0001, rho=1e-06】

当\lambda = 0.0001时，正则化强度过弱，TV约束几乎不起作用。重建结果与原始低剂量EM重建效果一样，噪声过多

【图片：ADMM-TV, lambda=0.001, rho=1e-06】

当\lambda = 0.001时，TV正则化开始发挥作用，噪声被显著抑制。但正则化效果过于强，导致图像过于平滑，细节丢失

【图片：ADMM-TV, lambda=0.01, rho=1e-06】

【图片：ADMM-TV, lambda=0.1, rho=1e-06】

当\lambda = 0.1时，同上

当\lambda = 0.01 和 \lambda = 0.1 时，正则化过强，图像几乎完全平滑化，结构信息完全丢失

\rho参数的影响

后续为了选择最好的正则化参数，我们固定\lambda = 0.005、迭代次数n_{it} = 50的条件下，测试了\rho \in \{10^{-7}, 10^{-6}, 10^{-5}, 10^{-4}\}四个值。

【图片：ADMM-TV, lambda=0.005, rho=1e-07】

当\rho = 10^{-7}时，由于 \rho 过小，\lambda/\rho 过大，导致正则化强度过高，图像细节丢失

【图片：ADMM-TV, lambda=0.005, rho=1e-06】

同上

【图片：ADMM-TV, lambda=0.005, rho=1e-05】

当\rho = 10^{-5}时，重建质量变好，可以看到边界和结构信息，且噪声水平得到了抑制，但整体和理想情况差距仍然较大

【图片：ADMM-TV, lambda=0.005, rho=1e-04】

当\rho = 10^{-4}时，正则化效果过弱（\lambda/\rho = 50），图像有大量噪声。

因此，\rho和\lambda需要联合调节，以在去噪声(\lambda)和收敛稳定性(\rho)之间获得最好的平衡

Part II: 去噪器对比

Question 3: 后处理去噪器比较

本部分对低剂量EM重建结果应用多种后处理去噪方法，包括传统方法（TV近端算子、高斯滤波、BM3D）和深度学习方法（net0、net1、net2），比较它们的去噪效果、执行时间和模型复杂度。

TV近端算子去噪

【图片：TV denoised low dose】

TV去噪结果显示噪声得到一定程度的抑制，但仍有较多残留噪声

高斯滤波去噪

【图片：Gaussian denoised low dose】

高斯滤波结果显示噪声被有效平滑，图像整体更加连续。但高斯滤波是各向同性的，在抑制噪声的同时也模糊了边缘，导致边界的对比度降低。

BM3D去噪

【图片：BM3D denoised low dose】

BM3D通过在图像中搜索相似块并进行协同滤波，能够在保持边缘的同时有效抑制噪声。重建结果略优于TV和高斯滤波。

深度学习去噪器

【图片：Denoised with net1】

net1网络的去噪结果中噪声被有效抑制，结构细节保持良好，且图像没有被过度平滑，非常漂亮

【图片：Denoised with net2】

net2网络的去噪效果与net1类似，同样实现了良好的噪声抑制和细节保留。两个网络的输出在视觉上难以区分

Part III: PnP ADMM

Question 4: 神经网络替代TV

PnP（Plug-and-Play）方法的核心思想是将ADMM中的正则化近端算子替换为任意去噪器。在标准ADMM框架中，步骤(S2)计算TV正则化的近端算子，在PnP-ADMM中，这一步被替换为直接调用去噪网络：

z^{(k+1)} = D_\theta(\theta^{(k+1)} + u^{(k)})

其中D_\theta表示预训练的去噪神经网络。实验使用参数\rho = 8 \times 10^{-10}，迭代次数n_{it} = 200，分别测试net1和net2两个网络。

为监控算法收敛性，记录:

原始残差（primal residual）\|\theta^{(k)} - z^{(k)}\|
对偶残差（dual residual）\|z^{(k+1)} - z^{(k)}\|

当两个残差都趋向于零时，ADMM算法收敛。

net1的PnP ADMM结果

【图片：reconstruction with PnP ADMM - net1】

【图片：Ref（正常剂量EM重建）】

【图片：ordose（Phantom）】

net1的PnP ADMM重建结果质量良好，噪声被有效抑制，灰白质边界清晰，整体效果接近正常剂量重建和Phantom。图像没有明显的伪影或异常。

【图片：MSE - net1】

MSE曲线在前20次迭代内快速下降，随后保持稳定，收敛到一个稳定的解。

【图片：Log likelihood - net1】

对数似然值速上升，小幅波动后稳定。

【图片：Primal residual - net1】

【图片：Dual residual - net1】

原始残差和对偶残差都在持续下降，表面算法正在收敛

net2的PnP ADMM结果

【图片：reconstruction with PnP ADMM - net2】

net2的PnP ADMM重建结果很差，有明显的异常伪影，

【图片：MSE - net2】

MSE曲线表面算法不稳定，最后MSE急速上升，说明重建结果反而越来越偏离目标

【图片：Log likelihood - net2】

【图片：Primal residual - net2】

原始残差和对偶残差同样不稳定，🕛上升时而下降，说明算法无法收敛

【图片：Dual residual - net2】

收敛性分析

这一差异的根本原因在于神经网络是否满足ADMM收敛所需的条件。传统ADMM的收敛性依赖于近端算子的性质（如非扩张性），而任意的神经网络去噪器并不一定具备这些性质。net1恰好满足或接近满足收敛条件，而net2则严重违反这些条件，导致算法发散。后续问题将通过分析网络的Lipschitz常数来进一步解释这一现象。

Question 5: \rho敏感性分析

本部分固定使用net1，测试\rho \in \{8 \times 10^{-11}, 8 \times 10^{-10}, 8 \times 10^{-9}, 8 \times 10^{-8}\}四个值对PnP ADMM重建的影响，迭代次数n_{it} = 200。

\rho = 8 \times 10^{-11}

【图片：PnP ADMM net1, rho=8e-11】

【图片：MSE, rho=8e-11】

当\rho = 8 \times 10^{-11}时，重建图像中噪声过多，去噪效果不充分。MSE不断下降并稳定

\rho = 8 \times 10^{-10}

【图片：PnP ADMM net1, rho=8e-10】

【图片：MSE, rho=8e-10】

当\rho = 8 \times 10^{-10}时，重建质量最佳，图像噪声被有效抑制，结构细节保持良好。MSE快速下降并稳定

\rho = 8 \times 10^{-9}

【图片：PnP ADMM net1, rho=8e-09】

【图片：MSE, rho=8e-09】

\rho = 8 \times 10^{-8}

【图片：PnP ADMM net1, rho=8e-08】

【图片：MSE, rho=8e-08】

当\rho = 8 \times 10^{-9} 和\rho = 8 \times 10^{-8}时，图像过度平滑，MSE随后上升，即神经网络对图像的过度平滑反而使重建结果偏离目标。

Question 6: 谱范数关系

compute_norm_jac函数使用幂迭代法计算Jacobian矩阵J的谱范数。根据幂法原理，该方法适用于计算矩阵的主特征值（最大模特征值）。对于一般的非对称矩阵J，无法直接应用幂迭代；但矩阵J^T J是对称半正定的，满足幂法的适用条件。

幂迭代应用于J^T J时，迭代格式为：

\mathbf{v}^{(k+1)} = \frac{J^T J \mathbf{v}^{(k)}}{\|J^T J \mathbf{v}^{(k)}\|_2}

根据幂法收敛性理论，序列(\mathbf{v}^{(k)})^T (J^T J) \mathbf{v}^{(k)}收敛到J^T J的主特征值\|J^T J\|。由于J^T J的特征值恰为J的奇异值的平方，其最大特征值等于J的最大奇异值的平方，即：

\|J^T J\| = \sigma_{\max}^2(J) = \|J\|^2

因此，函数最后执行.sqrt()操作，得到\|J\| = \sqrt{\|J^T J\|}。

Question 7: 收敛时u的值

ADMM对偶变量的更新公式为：

u^{(k+1)} = u^{(k)} + \theta^{(k+1)} - z^{(k+1)}

当算法收敛时，原始残差\|\theta^{(k)} - z^{(k)}\| \to 0，即\theta^{(k)} \to z^{(k)}。此时更新量\theta^{(k+1)} - z^{(k+1)} \to 0，对偶变量u不再变化，稳定在某个固定值u^*。

从优化理论的角度，u是增广拉格朗日函数中对应等式约束\theta = z的拉格朗日乘子。收敛时u^*代表该约束在最优解处的对偶变量值。在PnP框架下，由于去噪器D不一定对应某个显式的正则化函数，u^*的具体物理意义可能不如传统ADMM那样明确，但其数学性质保持不变：当原始变量收敛时，对偶变量也收敛到一个稳定值。

Question 8: FNE条件检验

PnP ADMM的收敛性理论指出，当去噪器D满足FNE（Firmly Non-Expansive，坚定非扩张）条件时，算法收敛。FNE条件要求算子2D - I是1-Lipschitz的，即其Jacobian矩阵的谱范数不超过1：

\|J_{2D-I}\| \leq 1

使用compute_norm_jac函数在低剂量EM重建图像上计算net1和net2对应的2D - I的Jacobian谱范数，结果如下：

网络	\|J_{2D-I}\|
net1	0.998
net2	29.29

net1的谱范数约为0.998，小于1，满足FNE条件。这解释了为什么Q4中net1的PnP ADMM能够稳定收敛：MSE持续下降、残差趋于零、重建质量良好。

net2的谱范数约为29.29，远大于1，严重违反FNE条件，这直接导致了Q4中观察到的发散现象：迭代过程中误差不断累积放大，最终MSE上升、残差振荡、重建失败。

Question 9: 迭代过程中谱范数监控

本部分在PnP ADMM迭代过程中，每一步计算去噪器输入点处2D - I的Jacobian谱范数，监控其随迭代的变化。

【图片：Jacobian spectral norm along iterations】

图中蓝色曲线为net1，橙色曲线为net2，红色虚线为FNE阈值（\|J_{2D-I}\| = 1）。

net1的谱范数在迭代中一直低于FNE阈值，满足FNE条件

net2的谱范数从最初的23不断下降，直至收敛都无法超越FNE阈值，不满足收敛条件，这解释了Q4中观察到的发散行为

Question 10: FNE修正

当去噪器D对应的2D - I的Lipschitz常数为L > 1时，可以构造一个新的去噪器D'使其满足FNE条件：

D' = \left(1 - \frac{1}{2L}\right)I + \frac{1}{2L}(2D - I) = \left(1 - \frac{1}{2L}\right)I + \frac{1}{L}D - \frac{1}{2L}I = \left(1 - \frac{1}{L}\right)I + \frac{1}{L}D

对于net2，L \approx 29.3，因此修正后的去噪器为：

D'(x) = \left(1 - \frac{0.5}{29.3}\right)x + \frac{0.5}{29.3} \cdot \text{net2}(x) \approx 0.983x + 0.017 \cdot \text{net2}(x)

使用修正后的FNE去噪器运行PnP ADMM，参数\rho = 4 \times 10^{-8}，迭代次数n_{it} = 200。

【图片：reconstruction with PnP ADMM FNE】

【图片：MSE】

修正后的重建结果显示算法逐渐收敛，MSE不断下降，但图像中仍残留大量噪声，效果不如NET1（MSE值远高于net1）

【图片：Primal residual】

【图片：Dual residual】

原始残差，对偶残差逐渐下降，证明了修正后net2满足FNE条件，保证了PnP ADMM的收敛性。

Part IV: PnP Forward-Backward（可选）

Question A1: Backtracking线搜索

Forward-Backward算法中使用backtracking线搜索来自适应地选择步长\tau。该方法的优缺点如下：

优点

自适应步长选择可以让算法无需预先知道目标函数的Lipschitz常数，通过检验下降条件自动调整步长。如果步长过大，算法会减小步长，这保证了算法的收敛性，很适合像ET重建这类难以精确估计Lipschitz常数的问题

缺点

计算量大，且步长只能单调缩小（乘以因子\eta < 1），无法在后续迭代中恢复到较大步长，收敛速度慢

Question A2: L2-TV正则化FB算法

Forward-Backward算法用于求解以下优化问题：

\min_{\theta} \left\{ -L(\theta) + \lambda \|\nabla \theta\|_2^2 + \iota_{[0,+\infty)^N}(\theta) \right\}

FB迭代格式为：

\theta^{(k+1)} = \text{prox}_{\tau L + \iota_{[0,+\infty)^N}}\left(\theta^{(k)} - \tau \lambda \nabla r(\theta^{(k)})\right)

实验使用参数n_{it} = 50，\lambda = 10^{-9}。

【图片：reconstruction with FB - TV -l2 regularization】

FB算法的重建结果很好，噪声较小，没有过度平滑

Question A3: GSnet的PnP版本

GSnet是一种特殊设计的去噪器，其输出形式为：

\text{GSnet}(\theta) = N(\theta) + J_{N(\theta)}^\top (\theta - N(\theta))

经过一系列数学变换后，公式可以表达为：

\text{GSnet}(\theta) = \nabla \left( \frac{1}{2}\|\theta - N(\theta)\|_2^2 \right)

因此，将GSnet用于PnP Forward-Backward等价于使用正则化项r(\theta) = \frac{1}{2}\|\theta - N(\theta)\|_2^2。这个正则化项惩罚图像与其去噪版本之间的距离，鼓励重建结果接近去噪后的图像。

PnP FB的迭代格式变为：

u = \tau \lambda \cdot \text{GSnet}(\theta^{(k)}) + (1 - \tau \lambda) \theta^{(k)}

\theta^{(k+1)} = \text{prox}_{\tau L + \iota_{[0,+\infty)^N}}(u)

实验使用参数n_{it} = 30，\lambda = 10^{-8}。

【图片：reconstruction with PnP FB】

PnP FB的重建结果是最好的，非常漂亮，这说明通过深度学习所表达的隐式正则化比手动设计正则化参数更优

分享