CT图像重建之自编写滤波反投影逆Radon变换算法

资 源 简 介

本项目主要用于计算机断层扫描（CT）图像重建领域，旨在不依赖MATLAB图像处理工具箱内置iradon函数的情况下，依据数学原理底层实现逆Radon变换。核心功能围绕滤波反投影（Filtered Back Projection, FBP）算法展开，具体流程包括：首先，对输入的弦图（Sinogram）或投影数据进行一维快速傅里叶变换（FFT），将其从空域转换至频域；其次，在频域内设计并应用各种加权滤波器（如经典的Ram-Lak滤波器、Shepp-Logan滤波器或Cosine滤波器），以补偿直接反投影导致的低频放大和图像模糊现象，恢复图像的高频细节；再次，通过傅里叶逆变换（IFFT）将滤波后的数据还原至时域/空域；最后，编写反投影循环逻辑，利用坐标旋转矩阵和双线性插值等数值计算方法，将各角度的投影值精确映射回二维图像矩阵的对应像素点上，累加形成最终的重建图像。该项目不仅能输出重建后的图像，通常还包含与标准iradon函数结果的误差分析对比，以验证算法的准确性。

详 情 说 明

自编写滤波反投影算法实现的逆Radon变换项目

项目简介

本项目主要用于计算机断层扫描（CT）图像重建领域的教学与研究。项目的核心目标是在不依赖MATLAB图像处理工具箱内置 iradon 函数的前提下，依据数学原理底层实现逆Radon变换。通过编写滤波反投影（Filtered Back Projection, FBP）算法，项目成功模拟了CT成像的重建过程。

该实现涵盖了从正弦图（Sinogram）数据的频域滤波到空域反投影的全过程，并提供了多种滤波器选择以应对不同的图像重建需求。此外，项目还包含完整的算法性能评估模块，通过量化指标和可视化图表验证重建效果。

主要功能特性

• 底层算法实现：完全独立编写 FBP 算法逻辑，包括 FFT 滤波、滤波器频域设计及反投影插值，加深对 CT 重建原理的理解。
• 多滤波器支持：内置三种经典的频域滤波器，用户可对比不同滤波器对图像锐度与噪声的平衡效果：

* Ram-Lak (Ramp)：基础斜坡滤波器，保留高频信息。 * Shepp-Logan：结合 Sinc 窗函数，平滑高频噪声。 * Cosine：结合余弦窗函数，进一步抑制高频干扰。 * 频域转换：利用快速傅里叶变换（FFT）将投影数据转换至频域。 * 滤波处理：在频域内将投影数据与设计好的滤波器频响函数相乘。

• 数据仿真与处理：

* 支持生成 Shepp-Logan 头部模型作为原始图像（若环境不支持则自动降级为圆形几何体）。 * 模拟 Radon 变换过程生成正弦图（投影数据）。 * 为了提高 FFT 效率及防止混叠，算法自动对投影数据进行 2 的幂次零填充。

• 反投影重建：

* 基于坐标旋转矩阵计算像素点在不同角度下的投影位置。 * 利用线性插值（Linear Interpolation）获取精确的投影值。 * 实现图像矩阵的累加与归一化处理。

• 评估与可视化：

* 计算均方误差（MSE）和峰值信噪比（PSNR）作为客观评价指标。 * 生成包含原始图像、正弦图、三种滤波器重建结果及中心行剖面线对比的综合图表。

系统要求

• MATLAB (推荐 R2016b 及以上版本)
• 无需 Image Processing Toolbox 的 iradon 函数支持，但依赖基础数学库及 radon 函数用于生成仿真数据。

使用方法

直接运行主脚本 main 即可启动程序。程序将自动执行以下流程：

1. 初始化参数并生成仿真图像数据。
2. 执行 Radon 变换生成投影数据。
3. 依次调用自定义 FBP 算法，应用 Ram-Lak、Shepp-Logan 和 Cosine 滤波器进行重建。
4. 在控制台输出各算法的 MSE、PSNR 及耗时信息。
5. 弹出图形窗口展示详细的对比结果。

算法实现细节

本项目采取模块化设计，核心逻辑封装于自定义的 FBP 函数中，具体实现步骤如下：

1. 频域滤波器设计

算法首先在频域 [-1, 1] 范围内构建基础的斜坡（Ramp）响应 |w|，随后根据选择的滤波器类型应用不同的窗函数：

• Ram-Lak: 直接使用 |w|，无额外窗口。
• Shepp-Logan: 将 |w| 乘以 sinc 函数 sin(pi*f/2)/(pi*f/2)。
• Cosine: 将 |w| 乘以余弦函数 cos(pi*f/2)。
• 利用 ifftshift 调整滤波器响应顺序，确保零频分量位于数组首位，以匹配 FFT 的输出格式。

2. 投影数据滤波

• 对输入的正弦图数据在每一列（即每个投影角度）方向上进行 FFT 变换。此处采用了零填充策略（Padding），将长度扩展至最优的 2 的幂次，以消除时域混叠并提高计算速度。
• 将频域数据与预先设计好的滤波器进行点乘。
• 执行 IFFT（逆快速傅里叶变换）并取实部，随后截取有效数据长度，完成空域下的滤波操作。

3. 反投影（Back Projection）

这是重建过程中计算量最大的部分，逻辑如下：

• 坐标系建立：定义图像中心的笛卡尔坐标系，将图像像素索引映射为物理坐标。
• 旋转变换：对于每一个投影角度 $theta$，利用旋转公式 $t = x costheta + y sintheta$ 计算图像上每个像素点对应在探测器阵列上的位置 $t$。
• 数值插值：由于计算出的 $t$ 通常不是整数，算法使用 interp1 函数进行线性插值，从滤波后的投影数据中获取精确的密度值。
• 累加与归一化：将所有角度的投影值累加到图像矩阵中，并乘以系数 $frac{pi}{2 times text{角度总数}}$ 进行数值归一化，同时将负值截断为 0 以符合 CT 值的物理意义。

4. 误差分析

系统通过 calculate_metrics 模块计算重建图像与原始图像的差异：

• 预处理：将两幅图像通过 mat2gray 归一化至 [0, 1] 区间。
• MSE：计算差值图像的均方值。
• PSNR：基于 MSE 计算峰值信噪比，用于量化图像质量。

结果分析

程序运行结束后，会显示一个包含六个子图的窗口：

1. 原始图像：作为以及基准的 Shepp-Logan 模型。
2. 正弦图：显示 Radon 变换后的投影数据（横轴为角度，纵轴为传感器位置）。
3. Ram-Lak 重建：图像最清晰，保留了最多的高频细节，但同时也引入了最多的高频噪声，容易出现吉布斯现象。
4. Shepp-Logan 重建：通过 Sinc 窗平滑了高频分量，图像比 Ram-Lak 稍柔和，信噪比通常更高。
5. Cosine 重建：高频抑制最强，图像最为平滑，适合噪声较大的场景。
6. 剖面线对比：取图像中心一行的像素值绘制曲线，直观对比原始数据与 Ram-Lak 重建数据的吻合程度。