DIY编程器网

标题: CT图像重建算法的FPGA实现 [打印本页]

作者: admin 时间: 2015-4-27 16:08
标题: 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射衰减值，再经模/数转换，使各体素不同的衰减值形成相应各像素的不同灰度，各像素所形成的矩阵图像即为该层面不同密度组织的灰度图像。

图1.1 医学螺旋CT机实物照片

图1.2 CT检查示意图

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

图1.3螺旋X-CT扫描机的系统结构示意图

     如图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.4 Radon变换原理

1.3.3 傅里叶切片定理

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

图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算法中，即使重建物体断面中一个小区域的图像，也得围绕整个断面采集投影数据。而局部重建算法，仅需围绕感兴趣区域及其邻域采集投影数据，即可重建感兴趣区域的图像。局部重建算法可减少数据采集时间和重建时间，降低人体(或生物体)的放射摄入量。

1.3.5国内外研究现状

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

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

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.3

函数（实线）和有限带宽斜变滤波器傅里叶变换（虚线）的比较

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

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

图2.4 平行投影重建过程的流程图

     2.4 目标重建过程
     医用CT的典型图像矩阵512×512。对于50cm重建FOV，每个像素尺寸大约1mm。根据Nyquist采样理论，这样的采样密度所支持的最高频率成分是5lp/cm(线对/厘米)。如果想检查具有更高空间辨率成分的解剖结构，图像像素间的采样距离必须减小。这可以通过增大重建图像矩阵尺寸或减小重建FOV来实现。增大的图像尺寸不仅影响重建速度(因为要被重建像素数与矩阵尺寸成平方关系)，而且增加存储量。

     另一可选方案是减少重建FOV。因为大多数高分辨率应用只需要检查很小一个区域(如内听管或脊椎骨)，缩小的FOV不会产生限制。考虑到重建被定位到一个较小区域这个事实，该方法经常被称为“目标”或“缩放”重建。

     目标重建过程类似于全FOV重建。一旦得到一个滤波投影(滤波过程和全FOV过程一致)，反投影将缩小的FOV映射到投影上去。例如，假设对一个

重建FOV进行像素驱动反投影，该FOV中心相对于系统旋转中心的坐标

。进一步假设图像矩阵尺寸为n×n，图像矩阵中心标记为

。一个位于的图像像素可以映射到中心为旋转中心的原始坐标系统中的一个点，根据以下等式：

(2.20)

（2.21）

     由于知道如何对位于(x.y)的点进行全FOV重建的反投影，滤波投影可以按照类似于全FOV重建的方式定位、插值，并加到重建图像中。
用Simulink创建的模型可以具有递阶结构，因此用户可以采用从上到下或从下到上的结构创建模型。用户可以从最高级开始观看模型，然后用鼠标双击其中的子系统模块，来查看其下一级的内容，以此类推，从而可以看到整个模型的细节，帮助用户理解模型的结构和各模块之间的相互关系。在定义完一个模型后，用户可以通过Simulink的菜单或Matlab的命令窗口键入命令来对它进行仿真。菜单方式对于交互工作非常方便，而命令行方式对于运行一大类仿真非常有用。采用SCOPE模块和其他的画图模块，在仿真进行的同时，就可观看到仿真结果。除此之外，用户还可以在改变参数后来迅速观看系统中发生的变化情况。仿真的结果还可以存放到Matlab的工作空间里做事后处理。

     模型分析工具包括线性化和平衡点分析工具、Matlab的许多工具及Matlab的应用工具箱。由于Matlab和Simulink的集成在一起的，因此用户可以在这两种环境下对自己的模型进行仿真、分析和修改。

         第四章 FPGA的实现
     4.1Modelsim仿真
     4.1.1 Modelsim简介
     Modelsim仿真工具是Model公司开发的。它支持Verilog、VHDL以及他们的混合仿真，它可以将整个程序分步执行，使设计者直接看到他的程序下一步要执行的语句，而且在程序执行的任何步骤任何时刻都可以查看任意变量的当前值，可以在Dataflow窗口查看某一单元或模块的输入输出的连续变化等，比Quartus自带的仿真器功能强大的多，是目前业界最通用的仿真器之一。

     4.1.2 仿真结果

通过Modelsim仿真工具，可以得到滤波反投影重建的仿真结果如图4.1所示：

图4.1 仿真结果

     4.2ISE仿真
     4.2.1Xilinx ISE开发平台简介
     Xilinx ISE是Xilinx公司的EDA软件开发系统，是一个集成化环境，主要由项目导航工具、设计输入工具、逻辑综合工具、设计实现工具、设计约束图形编辑接口等组成一个软件平台[13]。①项目导航工具是基本窗口界面，用来访问ISE软件系统的各种各种工具箱；②设计输入工具包括：电路逻辑输入工具—电路图编辑器、硬件描述语言输入工具—硬件描述语言编辑器（HDL Editor）、状态机编辑器、硬件描述语言测试生成器；③逻辑综合工具将硬件描述语言代码经过综合优化后输出EDIT格式电路逻辑连接（网表）；④设计实现工具用于面向FPGA的设计实现中的布局布线，并且可以对网表反标注以便提供给仿真工具进行后仿真验证；⑤设计约束图像编辑接口包含图像化的设计约束编辑接口，实现控制逻辑块的位置约束和时间约束。ISE软件界面如图4.2所示：

图4.2 ISE软件界面

4.2.2ISE的仿真
通过ISE的平台操作，可以得出以下Counters RTL级仿真和滤波器RTL级电路图：

图4.2 Counters RTL级仿真电路图

图4.3 滤波器RTL级电路仿真结果（部分图）

         第五章总结
     5.1 研究过程中所遇到的问题
        在整个项目过程中所遇到的问题都与滤波反投影FBP算法有关，一是，计算过程中的数据量较大；二是，对于要产生高分辨率的重建图像，那么就需要高像素的投影数据和较小的投影间隔，这些数据提升的同时会加大计算的数据量；还有就是高精度的重建工作是非常困难的，因为在浮点转定点运算的过程中，要考虑到量化的影响，所以我们要正确的选择位宽来实现定点运算，保证浮点转定点之后的误差最小。

     5.2 FBP算法所面临的挑战

FBP 算法的运算时间约99 %消耗在加权反投影阶段。因此,要提高FBP 重建的速度,须减少卷积反投影阶段的三角函数和浮点乘除的运算。国内外学者做过这方面的研究。如文献所提出的查表法,事先将固定的反投影的加权值和反投影的位置以参数表的形式储存起来,计算时再从表中查取。由于这两个值都与投影角成函数关系,因此,随着投影数量增大,参数表的规模将显著增大,既消耗大量内存,查表又需要用很多时间。

     5.3 课题的研究前景
     目前，市场上虽然已有商品化的CT三维重建软件，但是国内外各CT厂家的图像数据格式不同，并且格式不公开；另外文件格式经过加密处理，图像格式转换极为困难。西方主要发达国家的CT、MRI等医学成像设备和医学处理软件远远超过了我国。除了通用电气、西门子、东芝等医学成像设备制造商开发的配套软件产品外，比较成熟的医学图像处理软件还包括美国宾州大学开发的3D VIEWNIX系统、美国生物动力学研究中心开发的ANALYSE系统和德国汉堡大学开发的Voxel－Man系统等。这些系统各有长短，但大都在UNIX工作站上开发，其设备价格昂贵，不适合中国国情，并且其模型无法根据手术需要进行编辑，进而无法开展进一步手术计划和模拟等。目前，在我国各大医院以进口CT机型为主，而在软件的使用上，受到一定的限制。针对以上情况，从二维的算法研究入手，并在此基础上研制自己的三维重建系统，并为进一步的分析诊断、手术预测和模拟服务，具有一定的学术价值和实用意义。

     5.4 下一步计划
     目前绝大多数情况是通过观察人体的每一切片图像来进行诊断的。由于人体器官构造的复杂性和形态多样性以及病变位置与形态的不可预知性，没有相当的专业知识和实践经验，很难读懂这种二维切片图像，更难从这种二维切片构想出组织器官的立体型态和相互关系，所以仅仅从二维图像难以满足医疗诊断的要求。

     随着断层投影(CT)、核磁共振(MRI)、超声等医学成像技术的产生和发展，人们可以得到人体内部器官的二维数字断层图像序列或三维数据(称为医学体数据)。作为科学计算可视化的一个重要研究分支，医学体数据的三维重建是要在计算机上对这些离散数据进行拟合，将其转变成为具有直观三维效果的图像，利用人类视觉系统特性来展示人体器官的三维形态，从而提供若干用传统手段无法获得的解剖结构信息，并为进一步模拟操作提供可视交互手段。

     所以下一步计划实现三维空间的CT图像重建算法，同时要在算法本身上面实现理论上的突破，在浮点转定点的过程中，进一步合理的设置位宽，保证重建结果的精度。

欢迎光临 DIY编程器网 (http://diybcq.com/)