滤波反投影

滤波反投影法重建CT 图像实验指导书

一、 实验目的

1． 了解傅立叶变换法、直接反投影法重建CT 图像的原理；

2． 掌握滤波反投影法重建CT 图像的原理和基本方法。

二、 实验器材

装有MATLAB 程序的PC 机，滤波反投影法图像重建演示软件，投影数据。

三、 实验原理

CT 图像重建问题实际上就是如何从投影数据中解算出成像平面上各像素点的衰减系数。图像重建的算法有多种，如反投影法、傅立叶变换法、迭代法、滤波反投影法等。在介绍算法前，有必要先介绍从投影重建图像的重要依据，即中心切片定理。

1． 中心切片定理

密度函数(,)f x y 在某一方向上的投影函数()g R θ的一维傅立叶变换函数()g θρ是原密度函数(,)f x y 的二维傅立叶变换函数(,)F ρθ在(,)ρθ平面上沿同一方向且过原点的直线上的值。

图1 中心切片定理

2．傅立叶变换法

如果在不同角度下取得足够多的投影函数数据，并作傅立叶变换，根据中心切片定理，变换后的数据将充满整个(,)u v 平面。一旦频域函数(,)F u v 或(,)F ρβ的全部值都得到后，将其做傅立叶反变换，就能得到原始的密度函数(,)f x y ，即所要重建的图像。

上述图像重建算法称为傅立叶变换法，图2给出了傅立叶变换重建方法的流程图。图中指出，对于每次测得的投影数据先作一维傅立叶变换。根据中心切片定理，可将此变换结果看成二维频率域中同样角度下过原点的直线上的值。在不同投影角下所得的一维变换函数可在频域中构成完整的二维傅立叶变换函数，将此二维变换函数做一次逆变换，就得到了所要求的空间域中的密度函数。为了在二维逆变换中采用快速傅立叶变换算法，通常在逆变换前要将极坐标形式的频域函数变换成直角坐标形式的数据。

图2 傅立叶变换重建图像的过程

采用傅立叶变换法重建图像时，投影函数的一维傅立叶变换在频域中为极坐标形式，把极坐标形式的数据通过插补运算转换为直角坐标形式的数据时，计算工作量较大。此外，在极坐标形式的频域数据中，离原点较远的频率较高的部分数据比较稀疏，当这些位置上的数据转换到直角坐标下时，需经插补，这将引入一定程度的误差。即，在重建的图像中，高频分量可能有较大失真。

3．直接反投影法

直接反投影法把每次测得的投影数据“原路”反投影到投影线的各个像素上。即指定投影线上所有各点的值等于所测得的投影值。如图3所示的例子中，被探查物体只是在原点位置上的一个致密点。在一个角度上测量到其投影值时，就把这个值赋给投影线上的所有点。于是，从不同角度进行反投影后的重建图像是由以原点为中心的一系列辐射线。显然，原点位置上的分布密度最高；愈往四周，密度愈低。这当然可以说是粗略地把图像恢复出来了。但问题是除了密度最大的中心点，四周出现了逐渐变浅的云晕状阴影。

图3直接反投影重建法

研究发现，反投影后重建的密度函数(,)b f x y 与实际的密度函数(,)f x y 满足如下卷积关系 1(,)(,)(1)b f x y f x y r =** 上式说明了直接反投影法造成图像模糊的原因。一个以δ函数形式出现的密度函数用直接反投影法重建后成了一个带“长尾巴”的1/r 形式的图像，如图4所示。

图4直接反投影重建法造成的伪像

为了得到真实的重建密度函数，必须设法消除1/r 的影响。一种可能的方法是把式(1)转换到频域。因为时域中两个函数的卷积的傅立叶变换是这两个函数分别的傅立叶变换的乘积，因此式(1)相应的频域表达式为

1

(,)(,)(2)b F F ρβρβρ=

式中(,)b F ρβ和(,)F ρβ分别是(,)b f x y 与(,)f x y 的傅立叶变换，1/ρ为1/r 的傅立叶变换。从公式(2)可得

1(,)[(,)](3)b f x y F ρρβ-=

从式(3)可知，为了获得真实的密度函数(,)f x y ，可以先求出反投影函数(,)b f x y 的傅立叶变换(,)b F ρβ，在频域中对(,)b F ρβ加上权重ρ后求其逆傅立叶变换，就能得到所要的密度函数(,)f x y ，如图5所示。用这样的方法重建图像当然是可行的，但它还是没有避免计算二维傅立叶变换的问题。两次二维傅立叶变换所花费的时间还是相当可观的。

图5 直接反投影法的图像校正

4．滤波反投影法

直接反投影法重建的图像是模糊的。虽然这种模糊的图像可以经过修正后再现原始的密度函数，但修正的过程是很费时的。这就是先反投影，后修正重建方法存在的问题。滤波反投影法采用先修正、后反投影的做法，同样可得到原始的密度函数。其基本方法是：在某一投影角下取得了投影函数（一维函数）后，对此一维投影函数作滤波处理，得到一个经过修正的投影函数；然后再将此修正后的投影函数作反投影运算，得到所需的密度函数。这一方法中要解决的主要问题是如何修正投影函数才能使之在反投影后能重建原密度函数。图6给出了滤波反投影法重建图像的流程图。

图6 滤波反投影法

滤波反投影法重建图像有以下几个步骤：

（1）对某一角度下的投影函数作一维傅立叶变换；

（2）对（1）的变换结果乘上一维权重因子 ；

（3）对（2）的加权结果作一维逆傅立叶变换；

（4）用（3）中得出的修正过的投影函数做直接反投影；

（5）改变投影角度，重复（1）~（4）的过程，直到完成全部180度的反投影。

5．滤波函数

滤波函数的选取是滤波反投影法的关键问题。在实际系统中选择滤波函数时还要考虑许多其他因素，包括系统的带宽、信噪比与分辨率等。下面介绍两中不同类型的滤波函数。

（1）R-L 滤波函数

R-L 滤波函数的基本出发点是认为实际的二维图像函数总有一个频率上限，所以频域中的滤波函数ρ可表示为 0,()(4)0R L H ρρρρ-⎧≤=⎨⎩

其滤波函数如图7所示。

图7 R-L 滤波函数

连续的R-L 卷积函数为

22000()[2sin (2)sin ()](5)R L h R c R c R ρρρ-=-

离散化的R-L 卷积函数为

22221,

04()0,(6)1,R L n T h nT n n n T

π-⎧=⎪⎪=⎨⎪⎪-⎩为偶数

为奇数

图8（a ）（b ）分别为连续形式、离散形式的R-L 卷积函数。

图8 R-L 卷积函数，（a ）连续形式，（b ）离散形式