实验目的
本实验旨在探索PET图像重建中不同经典先验的效果。通过在模拟的2D数据和真实的3D数据上进行实验，深入理解最大似然期望最大化（MLEM）算法的特性，以及最大后验概率（MAP）重建中不同势函数和解剖学先验对图像质量的影响。实验的核心是观察信号恢复（通过对比度恢复系数CRC衡量）与噪声之间的权衡关系。
体模设置
体模为一个半径100 mm的圆柱形背景，背景活度值设为1.0。在背景中均匀分布着8个半径为10 mm的对比度圆盘，它们的相对活度值分别为4.0、3.0、2.5、2.0、1.5、1.25、0.5和0.0。其中前6个圆盘为正对比度（活度高于背景），最后两个为负对比度（活度低于背景，其中0.0表示冷区域）。
Q1：有噪声数据的事件数选择
实验设置
选择一个合适的真符合事件数 N，使得重建图像具有适当的噪声水平——既能清晰展示体模结构，又有足够明显的噪声以便后续观察正则化的效果。重建采用MLEM算法，设置4次迭代和14个子集。实验测试了6个不同的事件数：10^4、10^5、10^6、10^7、10^8、10^9。
重建结果
每组实验生成4幅图像：原始体模（Phantom）、单个切片的重建结果（Slice 0）、45个切片的均值图像（Mean across 45 slices）、以及45个切片的标准差图像（Std-dev across 45 slices）。标准差图像反映了图像的噪声水平。
【这里有一张图片：10^4 事件数的重建结果】
【这里有一张图片：10^5 事件数的重建结果】
【这里有一张图片：10^6 事件数的重建结果】
【这里有一张图片：10^7 事件数的重建结果】
【这里有一张图片：10^8 事件数的重建结果】
【这里有一张图片：10^9 事件数的重建结果】
结果分析
当事件数为 10^4 时，图像几乎完全被噪声淹没，无法分辨任何对比度圆盘，噪声极其严重。
当事件数增加到 10^5 时，图像质量变好了一点，但图像整体仍被噪声主导
随着事件数不断提升，图像质量不断提升，所有对比度圆盘都能被清晰分辨，在事件数为 10^8 时，图像已经非常清晰，噪声很低，事件数为 10^9 时，图像质量接近理想状态
参数选择
综合以上分析，本实验选择 N = 10^7作为后续实验的事件数，这个事件数在计算效率和图像质量之间取得了良好的平衡，适合进行多组对比实验。。因为 N = 10^6无法看到上面两个点，而 N = 10^9 的噪声太小，无法提现后续正则化效果
Q2：MLEM迭代次数的定性分析
实验设置
本问题旨在定性评估MLEM算法中迭代次数对重建图像质量的影响。实验参数设置为24次迭代、7个子集，优化算法为MLEM。每次迭代后保存重建结果，生成包含原始体模、单切片重建、均值图像和标准差图像的对比图。
在实验过程中，最初采用 10^9 事件数进行重建，但发现该事件数下图像噪声极低，标准差图像的数值范围仅为0至0.04左右，迭代过程中噪声变化不明显，难以观察到MLEM算法在噪声与信号之间的权衡特性。因此，将事件数调整为 10^7，此时图像具有适中的噪声水平，能够更清晰地展示迭代次数对图像质量的影响。
重建结果
【这里有一张图片：10^9 事件数下1-24次迭代的重建图像对比】
【这里有一张图片：10^7 事件数下1-24次迭代的重建图像对比】
结果分析
从 10^9 事件数的重建结果来看，由于初始噪声水平极低，即使在第1次迭代时图像就已经相当清晰，所有对比度圆盘都能被辨认。随着迭代次数增加，图像变化非常微小，标准差图像始终保持在很低的水平，这使得观察MLEM算法的收敛特性变得困难。
相比之下，10^7 事件数的重建结果展示了更为明显的迭代演化过程。在第1次迭代时，图像整体较为模糊，对比度圆盘的边界不够清晰，均值图像显示出较低的对比度恢复程度。随着迭代次数的增加，可以观察到以下变化趋势：
高对比度圆盘（活度值4.0、3.0、2.5）最先变得清晰，大约在第5-8次迭代后其亮度值接近真实值。中等对比度圆盘（活度值2.0、1.5、1.25）的恢复速度稍慢，需要更多迭代才能达到稳定。负对比度圆盘（活度值0.5、0.0）的恢复最为缓慢，即使在24次迭代后，冷区域的对比度仍未完全恢复到理论值。
在噪声演化方面，噪声随迭代次数单调递增。第1次迭代时标准差范围较小，图像相对平滑；到第24次迭代时，标准差范围显著增大，图像中出现明显的噪声纹理。
Q3：MLEM迭代次数的定量分析
实验设置
本问题通过对比度恢复系数（CRC）和噪声指标对MLEM算法进行定量分析。CRC定义为测量对比度与真实对比度的比值，其中对比度为（圆盘活度 - 背景活度）/ 背景活度。CRC值为1表示完全恢复，小于1表示对比度低估，大于1表示过估计。噪声指标采用背景中心区域标准差图像的均值与均值图像的比值（即变异系数）。
CRC与噪声曲线
【这里有一张图片：10^9 事件数下的CRC vs Noise曲线】
【这里有一张图片：10^7 事件数下的CRC vs Noise曲线】
结果分析
关于收敛速度的差异，高对比度圆盘的CRC收敛最快。以活度值4.0的圆盘为例，其CRC在噪声约为0.4时就已接近1.0，而此时低对比度圆盘的CRC仍在0.7-0.8左右。活度值为0.5和0.0的负对比度圆盘收敛最慢，最终未能达到完全恢复。
Q4：MAP先验的定性分析
实验设置
本问题旨在定性评估MAP重建中不同先验对图像质量的影响。实验采用BSREM（Block Sequential Regularized Expectation Maximization）优化算法，设置24次迭代和7个子集，事件数为 10^7。正则化强度参数 \beta 取值为 [0, 0.125, 0.5, 2.0, 9.0]，其中 \beta = 0 等价于无正则化的MLEM算法，\beta = 9.0 选取的依据是使该条件下的噪声水平接近MLEM第1次迭代的噪声水平。
实验测试了四种先验配置：
- Quadratic（二次势函数），对所有像素差异施加相同的二次惩罚
- Huber势函数，对小差异采用二次惩罚而对大差异采用线性惩罚，具有边缘保持特性
- Nuyts相对差分势函数，惩罚强度与像素值成比例
- Quadratic-Bowsher解剖先验，利用MRI解剖图像信息，仅在相似组织之间进行平滑
Quadratic先验结果
【这里有一张图片：BSREM-quadratic不同beta值的重建结果】
Huber先验结果
【这里有一张图片：BSREM-huber不同beta值的重建结果】
Nuyts先验结果
【这里有一张图片：BSREM-nuyts不同beta值的重建结果】
Bowsher解剖先验结果
【这里有一张图片：BSREM-quadratic-bowsher不同beta值的重建结果】
结果分析
从四组重建图像中可以观察到 \beta 值对图像质量的系统性影响。当 \beta = 0 时，所有四种配置的结果完全相同，等价于标准MLEM重建，图像中存在明显的噪声纹理
随着 \beta 值增大，图像逐渐变得平滑。在 \beta = 0.125 时，噪声略有降低但变化不明显，\beta = 9.0 时，图像达到最强的平滑效果，噪声大幅降低，但同时可以观察到对比度有所损失，尤其是低对比度圆盘变得更加难以分辨
比较四种先验的差异:
Quadratic先验产生全局均匀的平滑效果，在大 \beta 值时会导致边缘模糊。
Huber先验在保持类似降噪效果的同时，对比度圆盘的边缘相对更加锐利，这是因为Huber势函数对大梯度（即边缘）的惩罚从二次变为线性，减少了对边缘的过度平滑。
Nuyts相对差分先验的特点是平滑强度与局部活度成比例，因此低活度区域（如背景）的平滑效果更强，而高活度区域保持相对较好的对比度。
Bowsher解剖先验中的平滑仅发生在MRI图像中相似的组织之间，对比度圆盘的边缘保持最好，即使在 \beta = 9.0 时，圆盘边界仍然清晰，在噪声被抑制的同时，没有被过度正则化（平滑）
Q5：MAP先验的定量分析
实验设置
本问题通过CRC与噪声曲线对MLEM和四种MAP先验进行定量比较。与Q3相同，CRC衡量对比度恢复程度，噪声指标为背景中心区域的变异系数。每条曲线上的点对应不同的 \beta 值（从 \beta = 0 到 \beta = 9.0），\beta 越大对应的点噪声越低、位于曲线左侧。
各先验的CRC与噪声曲线
【这里有一张图片：BSREM-quadratic的CRC vs Noise曲线】
【这里有一张图片：BSREM-huber的CRC vs Noise曲线】
【这里有一张图片：BSREM-nuyts的CRC vs Noise曲线】
【这里有一张图片：BSREM-quadratic-bowsher的CRC vs Noise曲线】
各方法综合对比
【这里有一张图片：所有方法在8个对比度圆盘上的CRC vs Noise曲线汇总对比图】
结果分析
在相同噪声水平下，Bowsher解剖先验（紫色曲线）的CRC普遍最高，尤其在低噪声区域优势明显。这是因为Bowsher先验利用了MRI提供的解剖边界信息，能够在边界处避免过度平滑，从而在降噪的同时更好地保持对比度。三种非解剖先验（quadratic、huber、nuyts）的表现较为接近，曲线基本重叠，其中Huber和Nuyts先验在某些对比度水平下略优于Quadratic先验，但差异不大。
但该四种正则化方法都比标准的最大似然估计MLEM方法结果更好，代表着正则化的必要性
Q6：MLEM算法实现
算法原理
MLEM（Maximum Likelihood Expectation Maximization）算法是PET图像重建中最基础的迭代算法。其目标是在泊松噪声模型下最大化测量数据的似然函数。MLEM的迭代更新公式为：
其中 \theta 表示待重建的图像，m 表示测量数据（正弦图），b 表示加性校正项（包括散射和随机符合），A 表示正向投影算子（从图像空间到正弦图空间），A^T 表示反向投影算子（从正弦图空间到图像空间），A^T 1 表示灵敏度图像（对全1正弦图进行反向投影的结果）。
正弦图可视化
【这里有一张图片：事件正弦图】
重建结果
【这里有一张图片：30次迭代后的脑部PET重建图像】
结果分析
由于采用了无子集的标准MLEM算法，图像中可见明显的噪声纹理，这与Q2和Q3的模拟实验结果一致：MLEM算法在提高对比度的同时会累积噪声。若要获得更平滑的图像，可以减少迭代次数（牺牲对比度），或者采用MAP方法引入正则化约束。
Q7：真实脑部数据重建对比
实验设置
本问题将MLEM和MAP+Bowsher先验应用于真实的脑部FDG-PET数据，对比两种方法在实际临床数据上的表现。
MLEM重建参数设置为6次迭代、28个子集，投影算子采用Joseph方法，并使用高斯点扩散函数（PSF）进行建模校正。
MAP重建采用OSL（One Step Late）算法，这是一种适用于list-mode数据的MAP优化方法。正则化项采用Markov随机场（MRF）先验，邻域定义为6 mm半径的球形区域，并结合Bowsher解剖加权。解剖图像为同机采集的T1加权MRI（StudyID1622_0_0_MRI_3DT1.hdr）。正则化强度 \beta = 0.0007，迭代次数为8次，子集数为28。
MLEM重建结果
【这里有一张图片：MLEM重建的脑部PET图像（轴向、冠状、矢状三视图）】
MAP+Bowsher重建结果
【这里有一张图片：OSL+Bowsher重建的脑部PET图像（轴向、冠状、矢状三视图），以及对应的参数配置截图】
结果分析
MAP+Bowsher重建结果在保持相似对比度的同时，噪声得到明显抑制，其结果比MLEM的结果更清晰，颗粒状噪声消失，边界判断更加准确，整体看起来既保持了边界锐度，又整体平滑，而且消除了噪音的影响