第三节 滤波──逆投影法

  滤波──逆投影法是当前用得较多的一种图像重建方法,在当代 射线CT系统中几乎都用这种方法构成系统。它的特点是精度高,能快速实现。滤波──逆投影法又叫卷积──逆投影法。这是因为频域上的滤波相当于空间域上的卷积运算。本节首先讨论平行投影方式的滤波──逆投影法的基本原理及其重建滤波器的设计方法,然后讨论扇束投影的滤波──逆投影法。

  一、基本原理

  在式(10.2.17)中,当用FFT计算投影数据的付立叶变换 [见(10.2.13)式]时,投影数据 总被有限截断。当 的取样间隔为 时,在频率 的变化范围将是 。于是式(10.2.17)又可近似写为

   (10.3.1)

  若记

   (10.3.2)

  因为 ,上式又可写成

   (10.3.3)

  那么由式(10.2.18)和(10.2.13)得

  (10.3.4)

  而重建图像 (10.3.5)

  所以,只要使投影数据 先和脉冲响应为式(10.3.2)表示的滤波器进行卷积,然后再由式(10.3.5)对不同旋转角 求和便可重建图像。事实上,式(10.3.2)表示的恰好是频率响应为 的滤波器,通常称为重建滤波器,(即 滤波器)。

  二、重建滤波器之设计

  从上述可知滤波──逆投影法重建图像的关键是如何设计重建滤波器。由式(10.3.2)定义的重建滤波器的频率响应为:

  

图 10-3-1

  重建滤波器是一个非因果(物理不可实现)的理想滤波器,其二维空间特性和 与R的关系如图10-3-1所示。对应的重建滤波器的脉冲响应由式(10.3.2)可求得

   (10.3.6)

  在离散情况,该脉冲响应将被有限截断。设在 的变化范围内,对其取 点,取样间隔为 ,即

   ,其中

  那么根据式(10.3.6)可求得

   (10.3.7)

  在 之间的滤波系数可用如下线性插值求得

  

  式中

  由式(10.3.7)、(10.3.8)给出的加权函数称为雷英欣俊──拉克西米纳内亚南(Ramachandran-Lakshminacrayanan)滤波器,简称R-L滤波器。其图形如图10-3-2实线所示。

  如用级坐标 以及用离散值代替连续积分时,式(10.3.4)可写为

   (10.3.9)

  式中 。再用极坐标和直角坐标的关系,求出对应于 值,从而求得 或用离散值表示的 ;最后根据式(10.2.2)求出重建之图像。

   (10.3.10)

  式(10.3.9)和(10.3.10)是一组便于用数字计算机和各种快速算法实现的公式。

图 10-3-2 其它重建滤波器

  类似于第六章中阐述的道理,R-L重建滤波器频率特性在截止频率 处突然截断,将空间域引起脉冲响应有较大的振荡,从而引起较大的噪声。希坡──洛干(shepp-Logan)对式(10.3.5)和(10.3.6)表示的R-L滤波器进行了修正,修正后的滤波器称S-L滤波器,它的频率响应函数为

   (10.3.11)

  其脉冲响应可由 取付立叶反变换求得,离散形式为

   (10.3.12)

  S-L滤波器的图形如图10-3-2虚线所示,显示它改善了脉冲响应特性。

  滤波──逆投影法进行图像重建的关键是设计一个具有良好响应特性的重建滤波器,除上述R-L,S-L滤波器外,还可构成其它滤波器。

  三、用于讨论的是由平行投影(如图10-1-2所示)数据的滤波──逆投影法算法

  这种算法所对应的装置是采用单束扫描方式,由一个探测器跟随射线源作同步旋转某一角度间隔,再重新收集新的投影数据,如图10-3-3所示。


图 10-3-3 单束扫描方式

  用单束扫描方式收集全部投影数据一般需要几分钟时间。现在CT采用的是更快的扇束扫描方式,如图10-3-4所示。其扇形张角为30度左右,探测器约250到350个。旋转扫描一周,收集全部投影数据仅需2.5秒左右,大大地缩短了工作时间。


图 10-3-4 扇束扫描方式

  扇束扫描一般有两种类型,主要取决于探测器安放的方式。图10-3-5表示了这两类投影的区别。图10-3-5(a)表示等角间隔射线组成扇束,其探测器安置在以扇束顶点为中心的圆弧上,探测器以等弧间隔排列,因而扇束是由一组等角间隔的射线组成。图10-3-5(b)表示等间距直线排列探测器,因而扇束是由一组角间隔不等的射线组成。这两种算法在原理上略有差别,主要表现在取得投影的方法上,其重建算法上也略有不同。本节主要介绍等角间隔射线组成的扇束投影的滤波──逆投影算法。

图 10-3-5 扇束投影的两种类型

图 10-3-6 由等角间隔扇束投影推导重建算法的各种算法

  设 表示某一扇束投影,如图10-3-6(a)所示,该扇束是由一组等角间隔的射线组成。这里 代表扇束的中心射线 轴的夹角, 代表扇束内任一射线与中心射线的夹角。对于任一射线 点的投影值若用前面讨论过的平行射线积分取得的话,那么该积分值应是平行投影 在图(a)中所标的 位置上的值。我们的目的在于求得 之间的关系,从而借助于平行投影的重建算法以导出扇束投影的算法。

  参看图10-3-6(a),根据平行投影 的求法,可知 的垂直于 的,且 长度就等于 。故有

   (10.3.13)

  其中 为射线源 到坐标原点 的距离。根据式(10.2.19)和(10.3.3)可知平行投影 的重建图像为(注意式(10.3.3中的参数 即为 ):

   (10.3.14)

  式中 是平行投影中 所能取的最大值,即当 时, 。式(10.3.14)只要求 在180度范围内的投影值,如果要求使 在360度范围内的投影值,其式(10.3.14)可改写成

   (10.3.15)

  利用极坐标关系[见图10-3-6(b)] ,式(10.3.15)可表示为

   (10.3.16)

  由式(10.3.13)的关系,可改写为

   (10.3.17)

  上式中 的积分从 ,覆盖了整个360度范围。考虑到以 为变量的函数都是周期为 的函数,所以积分限可以改成以0到 。同时,由图10-3-6(a)可见,最外侧的射线 对应的 值等于 ,故 的积分限可改为 。另外,式中项 就是平行投影 中沿 的射线积分值,它与扇束投影中沿 的射线积分值是完全一样的。因此可得到

   (10.3.18)

  式中 为扇束投影。为了得到计算机容易实现的重建算式,我们将函数 的变量改写成

   (10.3.19)

  令 为从射线源到点 的距离(见图10-3-6(b)),显然, 的函数。再令 是通过点 的射线与中心射线的夹角,很容易证明下式

   (10.3.20 a)

   (10.3.20 b)

  因而由 三个参数可完全确定

   (10.3.21 a)

   (10.3.21 b)

  利用式(10.3.20)可将式 的变量表示成

   (10.3.22)

  将式(10.3.22)代入(10.3.18)可得

   (10.3.23)

  由式(10.3.2)表示的重建滤波器为

   (10.3.24)

  因此

   (10.3.25)

  现用下面变换

   (10.3.26)

  就可将(10.3.25)式改写成

   (10.3.27)

  于是式(10.3.23)可表示成

   (10.3.28)

  其中 (10.3.29)

  计算机实现(10.3.28)式时,可看成是一种加权滤波──逆投影算法。为了解释此点,可将式(10.3.28)写成

   (10.3.30)

  其中

   (10.3.31)

   (10.3.32)

  所以,该算法要分三个步骤来实现:

  1.假设每一投影 以等角间隔 取样得到的数据为 ,其中 取整数值。 时对应于通过投影中心的射线, 是投影所取的角度/然后按式(10.3.32)求得相应的 值为

   (10.3.33)

  2.计算每一修正投影 的卷积得到相应的滤波投影

   (10.3.34)

  其中 是(10.3.29)式的取样值。

   (10.3.35)

  若把(10.3.7)式 的值代入(10.3.35)式,可得重建滤波器的离散冲激响应

   (10.3.36)

  这里应注意的是:当采用FFT在频域上实现(10.3.34)卷积运算时,需将运算的序列加以扩展,以避免产生重叠误差。显然为改善图像质量,实际使用时需以某种平滑滤波器进行滤波。此时

   (10.3.37)

  其中 是平滑滤波器的冲激响应。在频域上实现平滑滤波时,可与汉明窗函数等进行简单的相乘运算。

  3.对每个扇束滤波投影进行加权逆投影运算。此时,它与平行投影中的情况有很大的不同,其差别可用图10-3-7加以说明。在扇束的情况下,逆投影是沿着扇束中射线进行的,这可从(10.3.30)式结构看出:

图 10-3-7 扇束投影的滤波逆投影

  

  式中 是通过点 和射线与扇束中心射线的夹角。为了求得 对图(10.3.7 b)上点 的贡献值,先要找出通过该点的射线 所对应的 角。如果此 值与已知 的各 值并不重合,就需要用插值的方法求出 。最后将 值除以 ,即为 对此点的贡献值。