基于不同窗函数的滤波反投影(FBP)图像重建系统
项目介绍
本项目是一个基于MATLAB开发的计算机断层扫描(CT)图像重建仿真系统。其核心目标是实现并对比不同频域窗函数在滤波反投影(Filtered Back Projection, FBP)算法中的应用效果。通过模拟从二维物体生成投影数据(正弦图),再利用频域滤波和空间域反投影技术还原图像,本系统直观展示了医学影像重建的基本物理过程和数学逻辑。
功能特性
- 多种滤波器选择:系统集成了五种主流的滤波器窗函数,包括Ram-Lak(斜坡滤波器)、Shepp-Logan、Cosine(余弦)、Hamming(汉明)以及Hann(汉宁)窗。
- 全流程仿真:涵盖了从原始Shepp-Logan头模型生成、Radon变换正弦图获取、频域滤波处理到空间域反投影重建的完整闭环。
- 自定义算法实现:核心反投影逻辑未直接调用MATLAB内置的高级重建函数,而是通过手动实现线性插值和坐标映射,展示了算法底层原理。
- 数理质量评估:系统自动计算重建图像与原始图像之间的均方误差(MSE),并利用柱状图进行定量对比。
- 可视化展示:通过多子图布局,同步对比原始图、正弦图、五种不同滤波结果以及精度评价指标。
系统要求
- MATLAB R2016a 或更高版本
- Image Processing Toolbox(用于调用radon和phantom函数)
---
核心实现逻辑方案
1. 实验控制中心
系统首先定义图像尺寸为256x256像素,并设置投影角度为0至179度的全覆盖扫描。通过生成标准的Shepp-Logan模型作为目标物体,利用Radon变换模拟探测器接收到的投影信号(正弦图)。在得到重建结果后,系统执行归一化处理(0-1范围),以确保不同窗函数下的图像在同一亮度标准下进行MSE误差评估。
2. 频域滤波器的构建
在自定义重建函数中,系统首先将一维投影信号通过N轴点FFT(补零至2的幂次)转换至频域。
- 基础滤波器:以Ram-Lak(斜坡)滤波器为核心,其响应函数为频率的绝对值 |w|,用于补偿反投影带来的低频叠加效应。
- 窗函数修正:为了抑制吉布斯振铃效应和噪声,系统在斜坡滤波器的基础上乘以下列窗函数:
*
Shepp-Logan:利用sinc函数对高频进行平滑。
*
Cosine/Hamming/Hann:通过不同的余弦加权函数降低高频增益,平衡图像锐度与噪声。
3. 空间域反投影过程
这是本系统的核心底层逻辑,其步骤如下:
- 坐标预处理:利用meshgrid创建居中的二维像素网格,计算每个像素点在旋转坐标系下的位置。
- 坐标映射:基于中心切片定理,将像素点的(X, Y)坐标投影到对应旋转角度的探测器轴线(t轴)上。
- 手动线性插值:
* 系统通过计算投影位置的底数(floor)获取相邻索引。
* 手动计算插值权重(距离比例)。
* 对超出物理探测器范围的点设置合理掩码(valid_mask)以防止数组越界。
- 贡献累加:遍历所有扫描角度,将滤波后的贡献值回填至重建矩阵。
- 常数校正:最后根据采样角度总数,乘以 $pi / (2 times text{角度数})$ 进行离散化积分补偿。
---
关键算法分析
滤波器对图像质量的影响
- Ram-Lak:重建图像边缘最锐利,但对噪声极其敏感,容易出现明显的条纹伪影和吉布斯现象。
- Shepp-Logan:在保持分辨率的同时,通过sinc函数的低通特性微弱地抑制了高频噪声。
- Hamming/Hann/Cosine:这三类窗函数显著降低了图像的边缘锐度,但有效地清除了高频噪声,产生的图像视觉上更加平滑。
手动实现反投影的意义
代码中舍弃了效率更高的内置插值函数,改用手动计算
idx_low、
idx_high 和
dist,并直接操作内存数组。这种方法虽然增加了代码量,但清晰地揭示了CT成像中“像素回填”的数学本质,即如何将一维的滤波信号按照几何比例重新分配给二维空间中的每一个像素点。
使用方法
- 启动MATLAB,将当前工作目录切换至本项目文件夹。
- 直接在命令行窗口运行主程序函数。
- 系统将自动弹出仿真结果窗口,并同步在命令行输出每种窗函数对应的精确MSE数值。