CT重建 平行射线滤波反投影

计算机断层重建(CT)是一个比较热门的领域，这篇文章简单介绍了反投影方法的重建过程。

参考资料：冈萨雷斯，《数字图像处理》，电子工业出版社。

直接反投影

该方法是沿着射线来的方向把一维信号反投影回去，可以想象成把投影穿过图像区域反“涂抹”回去。注意到相隔180°的投影互为镜像，因此，为了产生重建所要求的所有投影，只需要按照圆周的一半角度增量来考虑即可。下图是不同投影角度数量下直接反投影法的结果

从图中可以看到，直接反投影法有明显的伪影，也被称之为“晕环”。伪影的亮度比物体第，但比背景高。越多的投影使得伪影亮度越低，但会带来模糊。在CT重建中，去除伪影是一个重要的问题。

投影与Radon变换

Radon变换描述了投影产生时的数学过程。考虑直线$y=ax+b$，将其用极坐标来表示可写为：
$$x\cos \theta +y\sin \theta =\rho$$
对于图像$f(x,y)$，沿该线积分为：
$$g(\rho ,\theta )={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)\delta (x\cos \theta +y\sin \theta -\rho )dxdy$$
其中$\delta$为示性函数，或者是和示性函数差不多的冲激函数。离散情况下可写为：
$$g(\rho ,\theta )=\sum _0^{M-1}\sum _0^{N-1}f(x,y)\delta (x\cos \theta +y\sin \theta -\rho )$$

进一步可以得到直接反投影的数学表示：
$$f(x,y)={\int }_0^{\pi }g(\rho ,\theta )d\theta$$

滤波反投影法(FBP)

傅里叶切片定理

关于$\rho$投影的一维傅里叶变换为
$$G(\omega ,\theta )={\int }_{-\infty }^{\infty }g(\rho ,\theta )e^{-j2\pi \omega \rho }d\rho$$
傅里叶变换可以看成从空间域到频率域上的转换，这里$\omega$为频率变量。将$g$带入上式中可得：
$$\begin{array}{rl}G(\omega ,\theta )& ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)\delta (x\cos \theta +y\sin \theta -\rho )e^{-j2\pi \omega \rho }dxdyd\rho \\ & ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)\left[{\int }_{-\infty }^{\infty }\delta (x\cos \theta +y\sin \theta -\rho )e^{-j2\pi \omega \rho }d\rho \right]dxdy\\ & ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)e^{-j2\pi \omega (x\cos \theta +y\sin \theta )}dxdy\\ & =F(\omega \cos \theta ,\omega \sin \theta )\end{array}$$
其中$F(u,v)$为$f(x,y)$的二维傅里叶变换，即：
$$F(u,v)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)e^{-j2\pi (ux+vy)}dxdy$$

这说明了，对投影$g$的一维傅里叶变换，等于对图像$f$二维傅里叶变换在某一条线下的切片。这就是傅里叶切片定理。

平行射线下的滤波反投影重建

傅里叶逆变换：
$$f(x,y)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }F(u,v)e^{j2\pi (ux+vy)}dudv$$
令$u=\omega \cos \theta ,v=\omega \sin \theta$，则有
$$f(x,y)={\int }_0^{2\pi }{\int }_{-\infty }^{\infty }F(\omega \cos \theta ,\omega \sin \theta )e^{j2\pi \omega (x\cos \theta +y\sin \theta )}\omega d\omega d\theta$$
应用傅里叶切片定理，以及$G(\omega ,\theta +\pi )=G(-\omega ,\theta )$可得
$$\begin{array}{rl}f(x,y)& ={\int }_0^{2\pi }{\int }_{-\infty }^{\infty }G(\omega ,\theta )e^{j2\pi \omega (x\cos \theta +y\sin \theta )}\omega d\omega d\theta \\ & ={\int }_0^{\pi }{\left[{\int }_{-\infty }^{\infty }|w|G(\omega ,\theta )e^{j2\pi \omega \rho }d\omega \right]}_{\rho =x\cos \theta +y\sin \theta }d\theta \end{array}$$

方括号里可以看成一维傅里叶逆变换，只是附加了$|\omega |$项，这一项可以看为一个斜坡滤波器。因为$|\omega |$在无穷处都到达$+\infty$,所以傅里叶反变换无定义。实践中常常通过加窗的方式，使它在定义的范围之外为0。

简单的方法是通过方波限制，但方波有不希望的振铃特性，这会为图像重建带来伪影。因此可以考虑平滑窗，比如汉明窗或者韩窗：
$$h(w)=\{\begin{array}{rlr}& c+(c-1)\cos \frac{2\pi \omega }{M-1}& 0⩽\omega ⩽(M-1)\\ & 0& 其他\end{array}$$

总结： 获得完整反投影图像是由如下步骤得到的：

• 计算每个投影的一维傅里叶变换。
• 用滤波函数$|\omega |$乘以每个傅里叶变换，即乘以一个合适的窗。
• 得到每个滤波后的变换的一维傅里叶反变换。
• 对所有的一维反变换积分（求和）。

卷积与傅里叶反变换

考虑傅里叶反变换可以写成卷积的形式，可以有
$$\begin{array}{rl}f(x,y)& ={\int }_0^{\pi }{\left[|\omega |G(\omega ,\theta )e^{j2\pi \omega \rho }d\omega \right]}_{\rho =x\cos \theta +y\sin \theta }d\theta \\ & ={\int }_0^{\pi }{\left[s(\rho )\star g(\rho ,\theta )\right]}_{\rho =x\cos \theta +y\sin \theta }d\theta \\ & ={\int }_0^{\pi }\left[{\int }_{-\infty }^{\infty }g(\rho ,\theta )s(x\cos \theta +y\sin \theta -\rho )d\rho \right]d\theta \end{array}$$
其中$\star$表示卷积，$s(\rho )$表示$|\omega |$的傅里叶逆变换。

注：

• 重建过程中不需要存储所有的反投影图像，单个求和运算仅被最后的反投影图像更新。
• 斜坡滤波器在频率域的直流项归零，故每幅反投影图像的均值将为零。这意味着每幅反投影图像都将有负像素和正像素。简单的解决办法是所有像素减去最小值再放缩。
  $$f^{\prime }=K\left[\frac{f-\min f}{\max (f-minf)}\right]$$