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]的重建过程。[!--empirenews.page--]
2.4 目标重建过程
医用CT的典型图像矩阵512×512。对于50cm重建FOV，每个像素尺寸大约1mm。根据Nyquist采样理论，这样的采样密度所支持的最高频率成分是5lp/cm(线对/厘米)。如果想检查具有更高空间辨率成分的解剖结构，图像像素间的采样距离必须减小。这可以通过增大重建图像矩阵尺寸或减小重建FOV来实现。增大的图像尺寸不仅影响重建速度(因为要被重建像素数与矩阵尺寸成平方关系)，而且增加存储量。
另一可选方案是减少重建FOV。因为大多数高分辨率应用只需要检查很小一个区域(如内听管或脊椎骨)，缩小的FOV不会产生限制。考虑到重建被定位到一个较小区域这个事实，该方法经常被称为“目标”或“缩放”重建。
目标重建过程类似于全FOV重建。一旦得到一个滤波投影(滤波过程和全FOV过程一致)，反投影将缩小的FOV映射到投影上去。例如，假设对一个 重建FOV进行像素驱动反投影，该FOV中心相对于系统旋转中心的坐标 。进一步假设图像矩阵尺寸为n×n，图像矩阵中心标记为 。一个位于 的图像像素可以映射到中心为旋转中心的原始坐标系统中的一个点，根据以下等式：
(2.20)
由于知道如何对位于(x.y)的点进行全FOV重建的反投影，滤波投影可以按照类似于全FOV重建的方式定位、插值，并加到重建图像中。
第三章 Matlab下模型的建立
3.1 Simulink简介
由于Matlab具有语法筒单、易学、好写以及有强大的运算及绘图能力和强大且多样化的各种工具箱可供使用的优点，我们决定在Matlab下面进行模型的建立，对比m文件和simulink的优缺点，我们采用较为直观的Simulink仿真形式进行仿真。
Simulink是一个用来对动态系统进行建模、仿真和分析的软件包，它支持连续、离散及两者混合的线性和非线性系统，也支持具有多种采样频率的系统。在Simulink环境中，利用鼠标就可以在模型窗口中直观地“画”出系统模型，然后直接进行仿真。它为用户提供了方框图进行建模的图形接口，采用这种结构画模型就像你用手和纸来画一样容易。它与传统的仿真软件包微分方程和差分方程建模相比，具有更直观、方便、灵活的优点。Simulink包含有SINKS(输入方式)、SOURCE(输入源)、LINEAR(线性环节)、NONLINEAR(非线性环节)、CONNECTIONS(连接与接口)和EXTRA(其他环节)子模型库，而且每个子模型库中包含有相应的功能模块。用户也可以定制和创建用户自己的模块。
用Simulink创建的模型可以具有递阶结构，因此用户可以采用从上到下或从下到上的结构创建模型。用户可以从最高级开始观看模型，然后用鼠标双击其中的子系统模块，来查看其下一级的内容，以此类推，从而可以看到整个模型的细节，帮助用户理解模型的结构和各模块之间的相互关系。在定义完一个模型后，用户可以通过Simulink的菜单或Matlab的命令窗口键入命令来对它进行仿真。菜单方式对于交互工作非常方便，而命令行方式对于运行一大类仿真非常有用。采用SCOPE模块和其他的画图模块，在仿真进行的同时，就可观看到仿真结果。除此之外，用户还可以在改变参数后来迅速观看系统中发生的变化情况。仿真的结果还可以存放到Matlab的工作空间里做事后处理。
模型分析工具包括线性化和平衡点分析工具、Matlab的许多工具及Matlab的应用工具箱。由于Matlab和Simulink的集成在一起的，因此用户可以在这两种环境下对自己的模型进行仿真、分析和修改。
3.2 模块划分
模型共分3部分，第一部分是counters，主要是确定图像的大小;第二部分是filtering，完成数据的滤波过程;最后是back projection对数据进行反投影运算。
滤波器模块采用斜滤波器的模型，如图3.2所示：
反投影算法的理论依据是“中心切片定理”，该定理指出：线性衰减系数函数在某一方向上的投影函数的一维傅立叶变换函数是此衰减系数函数的二维傅立叶变换函数或其在平面上沿同一方向且过原点的直线上的值，根据这个原理，我们很容易得出反投影算法的理论公式，并可以得出仿真模型，如图3.3所示：
3.3功能介绍与仿真设置
整个模型实现的就是FBP算法，其中的每一个模块采用的都是定点算法，对于运算过程中的数据精度需要提前设置。
定点模型数据按照功能的不同划分为七个部分，每部分的精度设置如下：
fpw_cos = 12 正弦与余弦函数的位宽;
fpw_prj = 18 投影数据的位宽;
fps_prj = fpw_prj – 1 – nextpow2(max(max(abs(R)))) 投影数据的比例;
fpw_fil = 16 滤波器的位宽;
fps_fil = fpw_fil – 1 – nextpow2(max(abs(h1)))) l滤波器的比例;
fpw_out = 18 输出数据的位宽;
fps_out = fps_prj –(fpw_prj – fpw_out) 输出数据的比例。
3.4仿真结果
仿真图像采用120像素点的S-L头像，如下图3.4所示：
因为采集的头像仅为120个像素点，并且投影间隔较短，所以重建效果并不是非常理想，但是当采集点数上升和加密投影间隔之后，重建效果就会趋于理想，此结果已经可以证明算法的正确性。
第四章 FPGA的实现
4.1Modelsim仿真
4.1.1 Modelsim简介
Modelsim仿真工具是Model公司开发的。它支持Verilog、VHDL以及他们的混合仿真，它可以将整个程序分步执行，使设计者直接看到他的程序下一步要执行的语句，而且在程序执行的任何步骤任何时刻都可以查看任意变量的当前值，可以在Dataflow窗口查看某一单元或模块的输入输出的连续变化等，比Quartus自带的仿真器功能强大的多，是目前业界最通用的仿真器之一。
4.1.2 仿真结果
通过Modelsim仿真工具，可以得到滤波反投影重建的仿真结果如图4.1所示：
4.2ISE仿真
4.2.1Xilinx ISE开发平台简介
Xilinx ISE是Xilinx公司的EDA软件开发系统，是一个集成化环境，主要由项目导航工具、设计输入工具、逻辑综合工具、设计实现工具、设计约束图形编辑接口等组成一个软件平台[13]。①项目导航工具是基本窗口界面，用来访问ISE软件系统的各种各种工具箱;②设计输入工具包括：电路逻辑输入工具—电路图编辑器、硬件描述语言输入工具—硬件描述语言编辑器(HDL Editor)、状态机编辑器、硬件描述语言测试生成器;③逻辑综合工具将硬件描述语言代码经过综合优化后输出EDIT格式电路逻辑连接(网表);④设计实现工具用于面向FPGA的设计实现中的布局布线，并且可以对网表反标注以便提供给仿真工具进行后仿真验证;⑤设计约束图像编辑接口包含图像化的设计约束编辑接口，实现控制逻辑块的位置约束和时间约束。ISE软件界面如图4.2所示：
4.2.2ISE的仿真
通过ISE的平台操作，可以得出以下Counters RTL级仿真和滤波器RTL级电路图：
4.3系统联调
系统联调时，首先在Matlab里面建立好Shepp-Logan的模型，并且通过雷登变换计算出头像的投影数据，通过USB口发送至开发板，数据在板上运行完毕后，再将数据发送回PC机，送到Matlab中进行图像的现实处理，即可完成整个过程，下图为最终所得到的重建图形，结果证实整个过程是可行的。
第五章 总结
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图像重建算法，同时要在算法本身上面实现理论上的突破，在浮点转定点的过程中，进一步合理的设置位宽，保证重建结果的精度。
结论语
本文针对CT图像重建方法进行了研究，提出了二维图像重建的滤波反投影FBP算法。该算法具有重建速度快、空间和密度分辨率高等优点，缺点是对投影数据的完备性要求较高，从数学上讲，只有获得被检试件所有的Radon变换数据(完全投影数据)后才能精确重建其切片图像。接着采用计算机仿真对算法进行了验证。仿真实验表明，当投影数据较完备时，滤波反投影算法能获得较好的重建质量。需要指出的是上述算法是基于平行束投影数据的，对于扇形束投影数据，可采用数据重排或直接扇形束重建方法。尽管滤波反投影算法与其他图像重建算法相比具有较高的运算速度，但由于CT重建的数据量庞大，使得图像重建的计算非常耗时，尤其是对于高分辨率图像重建。因此，提高该算法的重建速度仍是今后的研究方向。
下载文档
