CT图像重建算法的FPGA实现 (二)

2.3 计算机实现的理论研究

在程序中，滤波反投影算法的步骤为：

投影数据采集

对投影数据做FFT变换

滤波

反投影数据

逆FFT变换

等式(2.8)不能以它现有形式直接实现，只要考虑公式(2.10)的解释，就很容易理解这一点。基于傅里叶变换的特性，我们知道在傅里叶域中两个函数相乘等价于两个相应空间域函数的卷积。 在空间域中的对应函数是被测平行投影 。对应滤波函数 的空间领域(或冲激响应) ，就是该函数的傅里叶反变换，

(2.12)

并不存在。必须研究一个代替方法。一个这样的方法是把有限带宽函数引入公式中。例如在上式中设置t=0，让我们考虑 的值。 代表在曲线 下面的面积。当 。因此，等式(2.8)不能以它现有形式实现。必须研究一个代替方法。一个这样的方法是把有限带宽函数引入公式中。

假设投影的傅里叶变换是有限带宽的。换句话说，在频率间隔 以外能量为0.在这个假设下，等式(2.10)可以按下面形式表示：

(2.13)

等式(2.13)指出，要计算滤波的投影 ，只需要进行投影 的傅里叶变换以得到 ，在 范围内乘以 ，并进行傅里叶反变换。不幸的是，有两个因素使这个看似简单的问题变得复杂：被截断的滤波核的离散化以及环状卷积的性质。要彻底理解滤波核问题，让我们首先在空间域中推导理想滤波核。为保证无混叠采样，投影带宽T必须满足Nyquist(奈奎斯特)采样准则：

(2.14)

其中 是投影采样间隔(单位为 )。在该条件下，初始的斜变函数 实际上是与窗函数 相乘：

(2.15)

其中

滤波函数 在图2.1中描绘。现在，滤波器冲激响应可以描述如下

. (2.16)

注意由于 的 的一个实偶函数，相应的冲激响应 也是t的一个实偶函数。

图2.1 有限带宽斜变滤波器的频率表示

注意，投影以间隔 采样。根据卷积理论，等式(2.9)可以写为

(2.17)

其中 是满足条件

的 值。这里，我们利用被扫描物体具有有限空间紧支集这一事实。在滤波投影的离散实现时，我们只对在 整数倍处的滤波数值感兴趣。把 代入等式(2.16)中，得到

(2.18)

滤波函数的冲激响应在图2.2中画出。在该图中，我们设 。如果用

表示在角度 下投影的离散采样，等式(2.10)中描述的滤波投影可以表达为一个空间域卷积：

(2.19)

图2.2 斜边滤波器的冲激响应

这里，我们利用了每个投影在空间上具有有限紧支集的事实。即在下标范围以外， 为0.这意味着，要确定 ，我们只需利用在范围 内的 。

尽管等式(2.19)的离散卷积实现可以直接得到被滤波的投影，当N很大时，往往在频率域中执行运算效率更高[使用快速傅里叶变换(FFT)运算]。对于目前一台典型的CT扫描机，一次单独投影的采样数N接近1000.因此，我们希望得到 序列的频率域形式。在有限范围内 的离散傅里叶变换 与等式(2.15)描述的 不同，如图2.3所示。二者之间主要差别是直流成分。尽管差相当小，它对重建图像CT数准确度的影响不能忽略。

现在我们考虑循环卷积[9]的问题。等式(2.10)中描述的原始滤波运算需要一个非周期性的卷积。当这个运算在频率域中执行时，只能是周期性或循环卷积。如果直接实线前面所述的运算序列，可能产生干涉伪像。这就是所谓的卷绕(warp-around)效应，或者周期间干涉。为了避免伪像，需要在傅里叶变换和滤波运算之前为每个投影填补足够数量的零。零的最少数目必须等于初始投影采样数减1(即N-1)。

图2.3中所示斜变滤波器

的特性表明，相对于低频成分， 更突出强调高频成分。事实上，斜变滤波器表现有些好像微分运算符。因此，可以把滤波运算想象为一个反卷积过程，去掉了反投影产生的模糊。

函数(实线)和有限带宽斜变滤波器傅里叶变换(虚线)的比较在等式(2.15)中，我们采用了一个简单的矩形窗函数来限制滤波核。可以另外修改窗函数，以改变滤波器的频率响应。实际应用中，窗函数经常被作为一个工具，用来修改重建图像的噪声特性，以实现空间分辨率和图像噪声之间的折中。

在许多用于数值计算和图像的高级语言软件系统中，如Matlab(The MathWorks,Natick,MA)或IDL(Research Systems,Inc,Boulder,CO),矢量或矩阵可以直接表示成变量。还可以针对矢量定义不同运算符。在这样的环境中，平行反投影的实现变得相当容易。对于每个被测投影(在数据预处理或预调理后)，投影被填补足够多的0以避免“周期间”干扰。对补零后的投影进行傅里叶变换，并且被变换的投影乘以一个滤波函数[10]。然后，对结果进行傅里叶反变换，得到一个被滤波的投影。该投影被反投影(通过“像素驱动”或“射线驱动”)到图像矩阵。为了提高空间分辨率，滤波投影经常在反投影过程之前进行预插值。在投影数据集合中对每次投影重复整个过程。图2.4显示一个流程图，描述了对于平行束投影[11]的重建过程。