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

第一章 绪论

1.1 引言

计算机断层摄像技术CT(Computerized Tomography)是20世纪医学的重大成果之一，该成果将计算机应用于医学领域，使医学射线学发生了革命性的变化。自从CT问世以来，计算机科学与生物医学工程相结合，形成了计算机医学图像研究的新领域，并为生命科学的研究提供了新的方法，成为近年来世界科技界最活跃、最富有生机和成就的领域之一。

医学影像学将数字图像处理技术和计算机图形学技术广泛的应用于生物医学领域中，通过把人体的内部结构以图像或图形的方式显示出来，提高了医疗诊断的可靠性，使治疗能够准确和彻底。

1.2 医用CT的简介

CT是计算机X射线断层造影术(Computerized Tomography)的简写。CT的发明是20世纪后期最重大的科技成果之一，由Hounsfield于1969年设计成功，1972年公诸于世。

CT利用人体各种组织(包括正常和异常组织)对X射线的吸收不等这一特性，将人体某一选定层面分成许多立方体小块(这些立方体小块称为体素)，X射线通过人体测得每一体素的密度或灰度，即为CT图像上的基本单位，称为像素。它们排列成行列方阵，形成图像矩阵。当X射线球管从一方向发出X射束穿过选定层面时，沿该方向排列的各体素均在一定程度上吸收一部分X射线，使X射线衰减。当该X射线束穿透组织层面(包括许多体素)为对面探测器接收时X射线量已衰减很多，为该方向所有体素X射线衰减值的总和。然后X射线球管转动一定角度，再沿另一方向发出X射线束，则在其对面的探测器可测得沿第2次照射方向所有体素X射线衰减值的总和;以同样方法反复多次在不同方向对组织的选定层面进行X射线扫描，即可得到若干个X射线衰减值总和。在上述过程中，每扫描一次，即可得一方程。该方程中X射线衰减总量为已知值，而形成该总量的各体素X射线衰减值是未知值。经过若干次扫描，即可得一联立方程组，经过计算机运算可解出这一联立方程组，而求出每一体素的X射衰减值，再经模/数转换，使各体素不同的衰减值形成相应各像素的不同灰度，各像素所形成的矩阵图像即为该层面不同密度组织的灰度图像。

 

 

螺旋CT检查包括两方面的基本内容：一是X射管及探测器连续360°旋转;二是患者同时随检查床匀速推进，如图1.2所示。在扫描时间内，X射线焦点对病人作螺旋式运动，并同时收集这一范围的全部扫描数据，用线性内插法重建图像。

 

如图1.3所示，医用X-CT机的系统结构主要由六大部分组成，其各部分的作用如下：

(1)X射线源：产生用来检测被测物的X射线，X射线源包括X射线球管源(能量在450kV以下)和直线加速器(能量在2MeV以上)。射线源的能量，决定了穿透能力。

(2)探测系统：包括准直器、传感器、信号处理和信号传输等部分，是获取信号的关键部分，也是决定CT性能的关键部分之一。穿过被测物的X射线首先通过准直器准直并离散化，传感器先将射线转换成电信号，信号处理电路再将不规则的信号转换成标准的信号传输到计算机接口。

(3)计算机采集系统：主要由特殊的专用的多信道数据采集接口电路和计算机软硬件组成。完成数据采集、转换、校正、处理等。将采集的数据处理成标准的文件格式，供图像重建、处理使用。

(4)机械扫描系统：作为各部分的载体并提供CT扫描所需的多个自由度的高精度运动。

(5)自动控制系统：包括检测、驱动、控制器(计算机)，完成扫描运动控制、系统逻辑和程控、状态监控和安全保护，协调整机工作，并完成系统自检与数据诊断。

(6)图像处理系统：包括图像处理计算机硬件和软件，如用于图像重建与处理的高速计算机、大屏幕图像显示器、大容量数据存贮器、图像拷贝输出设备(打印机)、系统软件及专用软件。完成数据校正、图像重建、处理、分析、测量、图像输出、存贮、显示等。

我们所研究的CT图像重建部分处于图像处理系统中，是整个系统的瓶颈所在，也是决定系统整个过程所消耗时间的关键部分。

1.3 CT图像重建技术概述

1.3.1 CT图像重建的简介

我们试图重建的物体可被看作是某种函数的二维分布。对于CT，该函数代表物体线性衰减系数。关于断层重建问题的描述，我们可以假设采集了一组测量结果，每个测量结果代表沿着特定的射线路径，物体衰减系数的累加或线积分。这些测量结果是在不同角度和到旋转中心的不同距离上获取的。为避免数据采样的冗余，我们假设测量按以下次序进行。首先沿着彼此平行且等间距的路径进行一组测量。这些测量结果构成一次“观测”或一组“投影”。在略微改变的角度下重复同样的测量。持续该过程直到覆盖整个360°(理论上仅有180°平行投影是必要的)。在整个过程中，相邻两次观测之间的角度增量保持不变，并且被扫描物体在同一位置固定不动。CT重建的问题就是，我们如何基于这些测量结果来估计被扫描物体的衰减系数分布。

CT图像重建问题是一个有趣而复杂的课题。它的公示表达可以追溯到1917年，当时Radon(雷登)[2]首先找到了从函数线积分重建该函数的求解方法。随着20世纪70年代后期和80年代早期临床实用CT扫描机的发展，该领域的研究活动有了极大的发展。大量研究论文、会议论文汇编、书籍章节，甚至教科书都关注这个课题[3,4]。提出了许多技术，它们在计算复杂性、空间分辨率、时间分辨率、噪声、临床治疗方案、灵活性以及伪像各方面具有不同的折中平衡。

1.3.2 Radon(雷登)变换

CT的基本思想源于1917年奥地利数学家Radon提出的Radon变换。

Radon变换的内容可以表述为：若已知某函数，

 

如图1.4所示，其沿直线S的线积分为：

 

(1.1)

 

则

 

(1.2)

 

式(1.1)为Radon变换，实际上就是物体的投影，式(1.2)为Radon反变换，即根据投影数据 重建函数 。

 

 

1.3.3 傅里叶切片定理

傅里叶切片定理的含义是：平行投影的一维傅里叶变换等同于原始物体的二维傅里叶变换的一个切片。即是指出线性衰减系数函数f(x,y)在某一方向上的投影函数gθ(R)的一维傅立叶变换函数Gθ(ρ)是f(x,y)的二维傅立叶变换函数F(u,v)或F(ρ,θ)(极坐标形式)在(ρ,θ)平面上沿同一方向且过原点的直线上的值，如图1.5所示。

 

为此，我们在不同的角度下取得足够多的投影函数数据，并作它们的傅立叶变换，那么变换后的数据就将充满整个(u,v)平面。一旦频域函数F(u,v)或F(ρ,θ)的全部值得到后，将其作一次傅立叶反变换，就得到原始的衰减系数函数f(x,y),即

 

(1.3)

 

令u=ρcosθ，v=ρsinθ，则式(1)可进一步变形为

 

 

 

 

 

(1.4)

 

式中 ，表示对投影函数 的傅里叶变换函数进行滤波变换，其中 为滤波函数。

 

由傅立叶变换性质可知，频域中的滤波运算可等效地在空域中用卷积运算来完成，因此由(2)可得到

 

(1.5)

 

式中h(R)为滤波函数 的空域形式， ，因而这种方法也称为卷积反投影方法。

 

利用中心切片定理[5]及二维FFT反变换法重建图像，由于勿须反投影运算,因而速度快,但图像重建过程中,需要内插运算，因而重建图像精度相对较低。

首先求出各投影数据的一维傅里叶变换，在不同的投影角度下所到的一维变换函数可构成完整的二维傅里叶变换函数，将此二维函数作一次反傅里叶变换，就得到重建图像。为了在二维逆变换中采用快速傅里叶变换算法，通常在逆变换前要将极坐标转化为直角坐标的形式。

傅里叶变换法重建法的特点是变换速度快，但精度不如滤波反投影法。算法的关键是将弧形的的极坐标数据转换成直角坐标数据时，由于在边缘区高频数据减少，因而造成误差，但傅里叶变换重建法重建速度比滤波反投影可提高2-3倍一在弧形极坐标数据向直角坐标系转化时，最简单的是最邻近内插法，当然这种方法精度最低，双线性内插重建图像精度好于最邻近内插法，而且计算又不复杂。

解决的方法是扩大计算区域，通过外延数据附加上一些格外的点，即计算更多的像素点以减小边缘的误差。如重建图像为M×M，则可计算3M×3M区域内的FFT变换，当然这是以增加了计算量为代价的。傅里叶变换重建图像算法在内插网格点上进行一些适当的选择。如使径向点取在直角坐标网格的线上，这样只需一次内插，而重建图像精度有了较大的改进。

1.3.4 CT图像重建的几种算法

在实际重建当中所存在的问题是，虽然Radon给出了一个数学公式，但是我们需要一个有效的算法来解决它，图像重建的算法有很多，大致分为三类：精确算法、近似算法和迭代算法。近似算法中，以滤波反投影算法(Filter back projection，FBP)最具代表性，应用最为广泛。选代算法中，代数重建算法(Algebraic reconstruction technique，ART)是提出最早并最为人们熟悉的算法。迭代型算法(如代数重建算法等)具有许多优点，但由于计算量大、重建时间长.在很长一段时间内限制了其在医学和工业CT领域的应用。提高迭代型算法的计算速度一直是人们关注的问题。近年来人们提出了不少提高迭代算计算速度的方法，加上近年来计算机计算速度的迅速提高，迭代算法重新受到人们青睐。此外，由于应用的需要，局部重建算法(Local Reconstruction Algorithm， LocalRA)也在近十年中有了较大的发展。在传统全局CT算法中，即使重建物体断面中一个小区域的图像，也得围绕整个断面采集投影数据。而局部重建算法，仅需围绕感兴趣区域及其邻域采集投影数据，即可重建感兴趣区域的图像。局部重建算法可减少数据采集时间和重建时间，降低人体(或生物体)的放射摄入量。[!--empirenews.page--]

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]