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

1.3.5国内外研究现状

同类课题所研究的技术基本上被国外所垄断，国内尚未有人提出，国内现在所使用的技术是利用PC机上软件来实现图像的重构，所需时间较长，如果用FPGA来实现的话，速度可以提高数十乃至上百倍。

1.4研究背景及意义

在当今社会大力发展医疗卫生条件的背景下，许多医院迫切需要先进的CT来为患者诊断病情，现在的CT技术被国外所垄断，设备也都在200万以上，只有极少数医院有能力配备，所以急需研发具有自主知识产权的产品，把价格控制在50万以内。CT的关键技术之一是快速断层图像重建技术，本课题的立足点就在于利用FPGA的高度并行性，实现CT断层图像重建算法，满足实际产品速度要求，为实现CT国产化准备，推动社会医疗卫生条件的发展。

第二章 滤波反投影算法

2.1 滤波反投影算法介绍

尽管傅里叶切片定理提供了断层成像重建的一个直接方案，在真正实现过程中，它提出了一些难题。首先，傅里叶空间中产生的采样模式不是笛卡儿坐标的。傅里叶切片定理说明一次投影的傅里叶变换是二维傅里叶空间中通过原点的一条直线。结果，不同投影采样落到极坐标栅格上。为了执行二维傅里叶变换，这些采样不得不被插值或重新栅格化到一个笛卡儿坐标中。二维频率域中的插值不像真实空间中的插值一样直接。在真实空间里，一个插值误差局限于像素所在的小区域。然而，对于频域插值，这个特性不再有效，因为二维傅里叶空间中每个采样表示某一个空间频域(在水平和垂直方向上)。于是，在傅里叶空间中一个单独采样点上产生的误差会影响整个图像(经过傅里叶反变换后)的外貌。为阐明傅里叶域插值的敏感性，进行下面的简单实验。扫描一个肩部模体，并在512×512矩阵中重建，矩阵用f(x,y)表示，其中x=0，1，…512，y= 0，1，…，512。下一步，执行图像的二维离散傅里叶变换，得到一个函数F(u,v)，其中u=0，1，…，511，v=0，1，…，511。注意F(u,v)是一个512×512复数矩阵。在该矩阵中，F(00)代表图像的直流成分。如果简单地进行函数F(u,v)的离散傅里叶反变换，将得到原始图像f(x,y)。注意函数F(u,v)是我们试图采用平行投影进行估计的量值(傅里叶切片定理)。

直接傅里叶域重建的另一缺点是进行目标重建的困难性。目标重建是在CT中常用的技术，用来检查物体中一个小区域的精密细节。如果能以某种方式把重建“聚焦”在感兴趣区，物体的细节就可以更好地显现。采用直接傅里叶重建方法，需要用大量的0填充F(u,v)，以进行必要的频率域插值。傅里叶反变换的大小和目标ROI的尺寸成反比。对非常小的ROI，矩阵尺寸庞大以至无法管理。尽管其他技术可以用来克服其中一些困难，这些技术的实现仍不直截了当。因此，必须研究傅里叶切片定理的替代实现方法。滤波反投影算法是目前得到广泛应用的基于变换法的图像重建算法，它具有重建速度快、空间和密度分辨率高等优点，缺点是对投影数据的完备性要求较高[7]，从数学上讲，只有获得被检试件所有的Radon变换数据(完全投影数据)后才能精确重建其切片图像。

2.2 滤波反投影算法公式的推导

我们从傅里叶变换和傅里叶反变换是共扼算子这一众所周知的事实开始。图像函数f(x,y)可以通过傅里叶反变换从它的傅里叶变换F(u,v)中恢复，

(2.1)

与推导傅里叶切片定理时进行的坐标变换类似，我们从笛卡儿直角坐标(u,v)转换到极坐标

。

坐标转换的目的是以更自然的数据采集形式表达数值F(u,v)。坐标转换如下：

，

， (2.2)

和

(2.3)

将等式(2.2)和(2.3)带入到(3.1)，得到

(2.4)

利用公式

中描述的傅里叶切片定理，我们用代替，建立如下关系：

(2.5)

.

对于平行采样几何束，在投影采样中存在一个微妙的对称性：

(2.6)

通过研究一组相差180°的平行束的采样几何，这个特性可以很容易理解。两组投影正好代表同一组射线路径。基于傅里叶变换的特性，对于相应的傅里叶变换对来说，存在一个简单关系：

(2.7)

将等式(2.7)代入等式(2.5)，我们得到下面等式：

(2.8)

通过在旋转坐标系(s,t)中表达上面的等式，并利用等式：

.

中指出的关系，我们得到下面等式：

. (2.9)

这里，是在角度投影的傅里叶变换。内部积分是数值的傅里叶反变换。在空间域，它代表一个经频域响应为的函数滤波后的投影。我们称之为“滤波投影”。

如果用标记等式(2.9)的内部积分所代表的角上的滤波投影：

. (2.10)

等式(2.8)可以下面形式重写

(2.11)

变量是从点(x,y)到一条通过坐标系原点，并与x轴成角的直线的距离。等式(2.11)说明，重建图像f(x,y)在位置(x,y)，是通过该点的所有滤波投影采样的累加。另外，我们还可以选择关注一个特定滤波投影采样，研究它对重建图像的贡献。因为代表与产生投影采样的射线路径重叠的一条直线，的强度沿着直线均匀地加到重建图像。结果，滤波投影采样的值沿着整个直线路径被“涂抹”或“叠加”。

我们还可以给出滤波反投影方法的一个直观解释。基于傅里叶切片定理，物体的二维傅里叶变换是通过将许多一维傅里叶变换拼起来得到的。理论上，如果假定一次投影的傅里叶变换形状像一个切成薄片的“派”，我们可以简单地把每个楔子插入适当位置，以得到物体的一个二维傅里叶变换。不幸的是，每个投影傅里叶变换形状类似在频率空间的一个长条。如果简单地计算所有投影傅里叶变换的和(假设在角度上等间隔)，中心区域被人为地增强，而外侧区域数值不足。为了用条形区域估计“派”状区域，我们可以给条形傅里叶变换乘以一个函数，该函数在靠近中心位置强度低，靠近边缘时强度高。例如，可以将投影的傅里叶变换与该频率处“派”状楔子的宽度相乘。如果假设N个投影在180°内均匀间隔，每个楔子的宽度在频率是。权重函数的最终作用是加权长条的累加与“派”状楔子的累加具有相同的“质量”。[8]