滤波反投影

从傅里叶切片定理到 FBP

反投影、傅里叶切片定理、Ramp 滤波和离散 FBP。

1 简单反投影

反投影算子把投影值沿原射线路径均匀“涂回”图像：

\[ \mathcal{B}p(x,y) =\int_0^\pi p_\theta(x\cos\theta+y\sin\theta) \,\mathrm{d}\theta. \tag{1}\]

简单反投影能够形成物体轮廓，但结果会出现明显的 \(1/r\) 型模糊。原因是低频成分在极坐标频域中被过度累积。

2 一维投影傅里叶变换

定义投影关于 \(s\) 的一维傅里叶变换

\[ P_\theta(\omega) =\int_{-\infty}^{\infty} p_\theta(s)e^{-i2\pi\omega s}\,\mathrm{d}s. \tag{2}\]

二维物体的傅里叶变换为

\[ F(u,v) =\int_{\mathbb{R}^2} f(x,y)e^{-i2\pi(ux+vy)}\,\mathrm{d}x\,\mathrm{d}y. \tag{3}\]

定理 1 (傅里叶切片定理) 投影 \(p_\theta(s)\) 的一维傅里叶变换，等于二维物体频谱 \(F(u,v)\) 沿角度 \(\theta\)、通过原点的一条径向切片：

\[ P_\theta(\omega) =F(\omega\cos\theta,\omega\sin\theta). \tag{4}\]

3 逆变换与 Ramp 滤波

将二维逆傅里叶积分转换到极坐标频域，面积元满足

\[ \mathrm{d}u\,\mathrm{d}v =|\omega|\,\mathrm{d}\omega\,\mathrm{d}\theta. \tag{5}\]

其中雅可比因子 \(|\omega|\) 正是 Ramp 滤波器。定义滤波投影

\[ q_\theta(s) =\mathcal{F}^{-1} \left\{ |\omega|P_\theta(\omega) \right\}(s), \tag{6}\]

则二维物体可写为

\[ f(x,y) =\int_0^\pi q_\theta(x\cos\theta+y\sin\theta) \,\mathrm{d}\theta. \tag{7}\]

式 7 就是平行束滤波反投影公式。

4 窗函数滤波器

理想 Ramp 滤波器增强高频，也会放大噪声。实践中常使用

\[ H(\omega)=|\omega|W(\omega), \tag{8}\]

其中 \(W(\omega)\) 可以是 Shepp–Logan、Cosine、Hamming 或 Hann 窗。窗口越强，高频噪声越低，但空间分辨率也会下降。

5 离散算法

对每个角度 \(\theta_k\)：

1. 对离散投影 \(p_k[m]\) 做一维 FFT；
2. 乘以离散 Ramp 或窗化 Ramp；
3. 做逆 FFT 得到 \(q_k[m]\)；
4. 对每个重构像素计算 \(s=x\cos\theta_k+y\sin\theta_k\)；
5. 在探测器坐标上插值 \(q_k(s)\) 并对所有角度累加。

离散近似通常写为

\[ f[i,j] \approx \Delta\theta \sum_{k=0}^{N_\theta-1} q_k\!\left(x_i\cos\theta_k+y_j\sin\theta_k\right). \tag{9}\]

下一节通过可运行的 Python 数值实验 比较反投影与 FBP。