滤波反投影（Filtered Back-Projection, FBP）算法是CT图像重建中广泛应用的一种方法，其核心思想是基于傅里叶切片定理，通过将投影数据在频域进行滤波处理，再进行反投影操作，从而得到高质量的重建图像。相比于直接反投影，FBP通过引入滤波器对投影数据进行预处理，可以显著减少图像模糊和噪声干扰，提升重建图像的质量[^2]。 下面是一个使用Python实现滤波反投影算法的示例。该实现基于NumPy和SciPy库，适用于模拟的Radon变换数据。 ### 代码实现 ```python import numpy as np import matplotlib.pyplot as plt from scipy.fftpack import fft, ifft, fftshift def radon_transform(image, theta=None, num_detectors=100): if theta is None: theta = np.linspace(0, np.pi, num=180, endpoint=False) projections = np.zeros((len(theta), num_detectors)) for i, angle in enumerate(theta): rotated = ndimage.rotate(image, np.degrees(angle), reshape=False) projections[i, :] = np.sum(rotated, axis=0) return projections, theta def filter_projection(projections, filter_type='ram-lak'): n_angles, n_detectors = projections.shape freqs = fftshift(np.fft.fftfreq(n_detectors)) if filter_type == 'ram-lak': filt = np.abs(freqs) elif filter_type == 'shepp-logan': filt = np.abs(freqs) filt[1:] *= np.sin(np.pi * freqs[1:]) / (np.pi * freqs[1:]) elif filter_type == 'cosine': filt = np.abs(freqs) * np.cos(np.pi * freqs / 2) filt[np.abs(freqs) > 0.5] = 0 elif filter_type == 'hamming': filt = np.abs(freqs) * (0.54 + 0.46 * np.cos(np.pi * freqs)) elif filter_type == 'hann': filt = np.abs(freqs) * (1 + np.cos(np.pi * freqs)) / 2 else: raise ValueError(f"Unknown filter type: {filter_type}") filtered_projections = np.zeros_like(projections) for i in range(n_angles): proj_fft = fftshift(fft(projections[i, :])) filtered_fft = proj_fft * filt filtered_projections[i, :] = np.real(ifft(ifftshift(filtered_fft))) return filtered_projections def back_project(filtered_projections, theta, image_size=100): reconstruction = np.zeros((image_size, image_size)) for i, angle in enumerate(theta): back_proj = np.zeros((image_size, image_size)) for x in range(image_size): for y in range(image_size): s = int((x - image_size // 2) * np.cos(angle) + (y - image_size // 2) * np.sin(angle)) s += image_size // 2 if 0 <= s < image_size: back_proj[x, y] += filtered_projections[i, s] reconstruction += back_proj return reconstruction # 模拟原始图像 image_size = 100 phantom = np.zeros((image_size, image_size)) phantom[30:70, 30:70] = 1 # 简单的方形物体 # 生成投影数据 theta = np.linspace(0, np.pi, num=180, endpoint=False) projections, _ = radon_transform(phantom, theta=theta) # 滤波处理 filtered_projections = filter_projection(projections, filter_type='ram-lak') # 反投影重建 reconstructed_image = back_project(filtered_projections, theta, image_size=image_size) # 显示结果 plt.figure(figsize=(12, 6)) plt.subplot(1, 2, 1) plt.imshow(phantom, cmap='gray') plt.title('Original Phantom') plt.subplot(1, 2, 2) plt.imshow(reconstructed_image, cmap='gray') plt.title('Reconstructed Image') plt.show() ``` ### 代码说明 1. **Radon变换**：`radon_transform` 函数用于生成模拟的Radon变换数据，即投影数据。它通过旋转图像并沿探测器方向求和来模拟实际的CT扫描过程。 2. **滤波器选择**：`filter_projection` 函数实现了多种常见的滤波器，如Ram-Lak、Shepp-Logan、Cosine、Hamming和Hann。不同的滤波器对重建图像的质量有不同的影响，用户可以根据需求选择合适的滤波器。 3. **反投影**：`back_project` 函数将滤波后的投影数据进行反投影操作，最终得到重建图像。 4. **图像显示**：最后使用 `matplotlib` 显示原始图像和重建图像，便于对比分析。 ### 注意事项 - 该代码使用了 `scipy.ndimage.rotate` 来模拟图像的旋转，实际应用中可能需要更高效的实现。 - 滤波器的选择会影响重建图像的质量，建议根据具体应用场景选择合适的滤波器。 - 为了简化代码，假设图像大小为 100x100，用户可以根据需要调整参数。 ### 总结 滤波反投影算法在CT图像重建中具有重要的地位，其结合了快速的傅里叶变换和反投影操作，能够在保证图像质量的同时实现高效的重建。上述代码提供了一个基础的实现框架，用户可以根据具体需求进一步优化和扩展。