2. 平行束与扇束重建算法实现

滤波反投影（Filtered Backprojection, FBP）是CT图像重建的经典算法，也是理解现代迭代与深度学习重建方法的基础。本节从平行束FBP的数学推导出发，逐步过渡到扇束几何的重排与直接重建，并给出可运行的代码实现与体模验证方法。

一、平行束FBP的推导逻辑

平行束投影的物理模型是：X射线沿一组平行路径穿过物体，探测器记录每条路径上的线积分。设物体衰减系数为 f(x,y)，投影角为 θ，探测器位置为 t，则投影数据为：

p(θ, t) = ∫∫ f(x,y) · δ(x·cosθ + y·sinθ − t) dx dy

FBP的核心思想是：先对每个角度的投影做一维滤波，再将滤波后的投影沿原路径反投影回图像空间。其推导依赖傅里叶切片定理——投影的一维傅里叶变换等于物体二维傅里叶变换在对应角度上的切片。由此可得重建公式：

f(x,y) = ∫₀^π [ p(θ, t) * h(t) ] |_{t = x·cosθ + y·sinθ} dθ

其中 h(t) 是滤波核，* 表示一维卷积。反投影本身不是逆变换，必须配合滤波才能恢复高频细节；若直接反投影，图像会被 1/r 模糊核严重平滑。

二、Ram-Lak与Shepp-Logan滤波器

理想滤波核在频域为 |ω|，称为斜坡滤波器。但理想斜坡在奈奎斯特频率处突然截断，会在空域产生振荡，导致重建图像出现振铃伪影。工程上常用两种加窗修正：

• Ram-Lak滤波器：直接使用 |ω| 并在截止频率处硬截断，空间分辨率最高，但对噪声敏感，适合高剂量、低噪声的投影数据。
• Shepp-Logan滤波器：在斜坡上乘以 sinc 窗，即 |ω| · sin(πω/ω_c)/(πω/ω_c)，高频被平滑，振铃明显减弱，噪声抑制更好，但空间分辨率略降。

实际实现时，滤波通常在频域完成：对投影做FFT，乘以滤波器的离散频域响应，再IFFT回到空域。离散化时需注意补零以避免循环卷积混叠，并保证滤波器长度与投影长度匹配。

三、扇束重排与等距探测器重建

实际CT系统多采用扇束几何：X射线源发出锥形/扇形束，探测器阵列接收。扇束数据不能直接套用平行束FBP，需先处理几何差异。常见路线有两种：

重排法：将扇束投影数据插值重排为等效平行束数据，再调用平行束FBP。对于等距探测器，扇束中探测器位置 u 与平行束坐标 t 的关系为 t = D·sin(γ)，其中 γ = arctan(u/D)，D 为源到旋转中心距离。重排需要二维插值，会引入插值误差和分辨率损失，但实现简单、复用性好。

直接扇束FBP：不重排，直接在扇束几何下推导滤波反投影公式。其核心区别在于：滤波前需对投影乘以 D/√(D²+u²) 的余弦权重，反投影时距离权重为 1/L²，其中 L 是源到重建点的距离。直接法避免了重排插值，理论上分辨率更高，但实现更复杂。

四、代码实现要点

以下给出平行束FBP的Python核心实现，使用Shepp-Logan滤波器，并说明关键参数：

import numpy as np

def shepp_logan_filter(n, d=1.0):
    """生成Shepp-Logan滤波器频域响应"""
    omega = np.fft.fftfreq(n, d)
    filt = np.abs(omega)
    # sinc窗
    sinc = np.sinc(omega / (2 * omega.max()))
    return filt * sinc

def fbp_parallel(sinogram, theta, n):
    """sinogram: (angles, detectors)"""
    n_angles, n_det = sinogram.shape
    # 补零至2N
    pad = np.zeros((n_angles, 2 * n_det))
    pad[:, :n_det] = sinogram
    # 频域滤波
    filt = shepp_logan_filter(2 * n_det)
    filtered = np.real(np.fft.ifft(np.fft.fft(pad, axis=1) * filt, axis=1))
    filtered = filtered[:, :n_det]
    # 反投影
    recon = np.zeros((n, n))
    for i, th in enumerate(theta):
        for x in range(n):
            for y in range(n):
                t = (x - n/2) * np.cos(th) + (y - n/2) * np.sin(th) + n_det/2
                if 0 <= t < n_det - 1:
                    t0 = int(t)
                    w = t - t0
                    recon[y, x] += (1 - w) * filtered[i, t0] + w * filtered[i, t0+1]
    return recon * np.pi / n_angles

实现时需注意：投影角度通常覆盖 [0, π)，反投影时角度步长 Δθ 要参与归一化；探测器坐标原点应与旋转中心对齐；插值建议使用线性或三次插值，避免最近邻带来的块状伪影。

五、体模验证与结果判读

验证FBP实现最直接的方式是使用Shepp-Logan体模。该体模由多个椭圆组成，包含不同对比度和尺寸的结构，能同时检验空间分辨率、对比度恢复和伪影水平。验证流程如下：

1. 生成 256×256 的Shepp-Logan体模图像，像素值范围 [0,1]。
2. 用Radon变换生成 180 个角度、367 个探测器的平行束投影。
3. 分别用Ram-Lak和Shepp-Logan滤波器重建，比较图像质量。
4. 计算重建图像与原始体模的RMSE和SSIM，定量评价。

典型观察结果：Ram-Lak重建的边缘更锐利，但背景区域噪声和振铃更明显；Shepp-Logan重建的振铃被抑制，但小椭圆结构的边界略模糊。若重建图像出现明显条纹或中心亮斑，通常说明滤波器归一化错误或反投影角度权重未正确累加。若图像整体模糊，应检查投影角度是否覆盖完整 π 范围，以及探测器间距与重建像素间距是否匹配。