基于OS-EM算法的图像迭代重建系统

资 源 简 介

本项目实现了一种高性能的OS-EM图像重建算法。该算法作为极大似然估计方法的重要演进，主要用于解决医学影像（如PET、SPECT和CT）及工业CT中的逆问题求解。算法的核心思路是通过将所有的投影切片数据划分为多个有序子集，并在每一个子步中利用单一子集的投影误差来更新当前估算的图像像素值。这种方法在数学上近似于全数据的梯度下降，但在实际操作中极大地缩短了达到收敛标准所需的时间成本。该项目详细实现了图像的正向投影匹配（Forward Projection）、投影数据的比值归一化、反向投影校正（Back Pro

详 情 说 明

基于有序子集极大似然期望最大化 (OS-EM) 的图像重建系统

项目介绍

本项目实现了一种高性能的医学与工业图像迭代重建算法——OS-EM（Ordered Subsets Expectation Maximization）。该算法是对经典极大似然期望最大化（MLEM）算法的改进，通过将投影数据划分为多个子集并进行增量式更新，在保持图像重建质量的同时，显著提升了收敛速度。本项目模拟了从原始断层扫描（如CT、PET）到最终图像恢复的全过程，旨在解决低信噪比环境下的图像逆问题求解。

功能特性

1. 加速收敛：通过子集化策略（Ordered Subsets），在单次全量迭代中多次更新图像，比传统MLEM快数倍。
2. 强大的抗噪性：在正向投影中引入高斯噪声，验证算法在低信噪比环境下的鲁棒性。
3. 伪影抑制：相比传统的滤波反投影（FBP）算法，能有效抑制条纹伪影并保留更多边缘细节。
4. 实时监控：系统自动计算并绘制迭代过程中的均方根误差（RMSE）与峰值信噪比（PSNR）变化曲线。
5. 灵活的几何配置：支持自定义图像分辨率、视角数量、子集数量和噪声水平。
6. 对比评估：内置标准FBP算法作为基准，并提供局部区域放大对比功能。

系统要求

1. 运行环境：MATLAB R2016b 或更高版本。
2. 必备工具箱：Image Processing Toolbox（用于执行Radon和Inverse Radon变换）。
3. 硬件建议：4GB以上内存，支持图形显示。

核心功能实现逻辑

项目逻辑严格遵循医学物理重建的标准流程，主要分为五个核心模块：

1. 参数初始化与子集划分

系统设定默认分辨率为256像素，总计180个投影角度，并划分为10个有序子集。为了保证子集之间具有较好的正交性和采样代表性，算法采用等间隔采样分配策略，即每个子集间隔固定的步长选取角度数据。

1. 物理过程模拟（正向投影）

使用Shepp-Logan头模型作为重建目标。通过Radon变换模拟X射线穿透物体的物理过程，生成原始正弦图。系统向投影数据中注入最大值5%的高斯随机噪声，模拟真实的探测器物理噪声，并强制执行非负性约束。

1. 规范化因子预计算

这是OS-EM实现的关键步骤。算法对全1的投影数据进行反向投影，为每个子集计算独立的规范化因子（Normalization Factors）。这一步骤补偿了反向投影过程中由于几何权重不平衡带来的偏差，确保护度更新的平稳性。

1. 核心迭代引擎

重建引擎利用双层嵌套循环实现：外层控制总迭代次数，内层遍历各个有序子集。

• 正向投影匹配：将当前估算的图像投影至当前子集的对应角度。
• 误差比值计算：计算测量投影值与估算投影值的比值。
• 反向投影校正：将投影域的误差比值映射回图像空间。
• 图像更新：基于乘法更新准则，利用校正因子修改当前像素值，随后进行非负性约束处理。

1. 性能度量与可视化

每次迭代后，系统会实时对比重建图像与原始模型，计算RMSE和PSNR。最后生成的交互式图表包括原始模型、含噪正弦图、FBP重建基准、OS-EM最终结果以及性能指标演化曲线。

关键函数与技术分析

• radon / iradon：作为仿真的核心，分别代表了物理检测的正向过程和反向投影的数学基础。在OS-EM中，iradon通过禁用内部滤波器（'none'参数）来实现散点式的反向更新。
• 非负性约束：根据物理意义，物质密度不能为负。代码在投影数据生成和图像更新步骤中均实施了非负滤波，有效防止了数值发散。
• 正交子集策略：代码通过 indices = s:subsets_count:angles_count 的逻辑分配角度，确保了每一个子集都能覆盖尽可能广的视角范围，这是OS-EM收敛稳定性的数学保障。
• 细节放大对比：通过对中心区域的局部切片对比，直观展示了迭代算法在去除噪声干扰和伪影方面的优势。

使用方法

1. 在MATLAB中打开主程序。
2. 直接运行主函数。
3. 观察控制台输出的迭代信息，包括当前步骤、RMSE和PSNR。
4. 迭代完成后，会自动弹出两个分析窗口：一个展示全局重建效果，另一个展示局部细节对比。
5. 如需测试不同工况，可修改程序顶部的 image_size, subsets_count 或 noise_level 等参数。