CT断层图像重建算法研究

1、CT断层图像重建算法研究专业：通信工程 姓名：刘明帅 指导教师：骆岩红摘 要CT技术是一种融合了射线光电子学、信息学、微电子学等学科的新兴技术，因为其先进的无损检测技术，所以被广泛地应用于医学、航天、生物等多个领域。随着科技的进步，图像重建技术开始应用于X射线中，这是数字图像处理的一个重大进步。如何能重建出高质量的图像，取决于所采用的重建算法。从图像重建的角度来看，主要分为解析法与迭代法。解析法是利用解析、变换重建公式来构建重建图像。它具有容易实现,速度较快,且能重建出高质量的图像的特点,但是对投影数据完备性要求高。迭代法是利用求解线性方程组来重建图像，它能够在投影数据信噪较低条件下，获得高质

2、量图像。本文将从原理、应用、与优缺点的角度来分析两种算法，重点对解析法中的滤波反投影算法从平行束与扇束投影方式进行研究，最后通过Visual C+与MATLAB软件相结合的方式对图像重建，并分析各参数对重建图像的影响。关键字：CT技术 图像 重建算法 滤波反投影算法AbstractCT technology is a emerging technology that blend of the Ray optoelectronics, microelectronics and informatics subject. Because of its advanced nondestructive

3、testing technology, it is widely used in medical, aerospace, biological and other fields. With the progress of science and technology, Image reconstruction technology is applied to the X ray, This is a major progress of digital image processing. How to rebuild the high quality images, depends on the

4、 reconstruction algorithm you adopt. From the perspective of image reconstruction, it mainly divided into the analytical method and iteration method. Analytical method use the analysis and transform formula to build image reconstruction.It has the characteristics of implementating easily and fa

5、st,and reconstructing out high quality images,but the demand of the projection data is high.Iterative method is used to solve the linear system of equations to reconstruction image, the projection data under the condition of low signal-to-noise, it can get high quality image.This article we will be

6、from the point of view of the principle ,application,and the advantages and disadvantages to analysis the two kinds of algorithms, focusing on studying the analytical method of filter back projection algorithm from the parallel beam and fan beam projection methods , finally, combining the softw

7、are of Visual c + + with MATLAB software to image reconstruction, and analyzes the influence of various parameters on the reconstruction imageKey words: CT technology image reconstruction algorithm Filtered Backprojection Algorithm目录第一章 绪论- 5 -1.1 CT技术与图像重建概述- 6 -1.2 CT和重建技术的发展及研究现状- 6 -1.3 研究的目的与意义

8、- 9 -第二章 CT成像原理和图像重建算法- 9 -2.1 CT成像原理与系统组成- 9 -2.2 CT成像系统扫描方式的发展- 10 -2.3 CT断层图像原理- 12 -2.4图像重建算法概述- 13 -解析类方法- 13 -传统迭代类方法- 14 -第三章 CT图像重建算法实现原理的研究- 14 -3.1图像重建系统中的数学概念及变换- 14 -3.1.1 投影与反投影- 14 -3.1.2 Radon变换及其反变换- 15 -傅里叶变换- 16 -中心切片定理- 17 -3.2解析类重建算法- 18 -直接傅里叶变换算法- 18 -反投影重建算法- 19 -3.3 迭代类重建算法-

9、21 -3.3.1 代数迭代重建算法- 22 -(1) ART算法- 22 -(2)同时代数重建算法- 22 -3.3.2 影响代数迭代重建算法的因素- 23 -3.3.3 ART重建算法与SART- 26 -3.3.4 统计迭代重建算法- 27 -（1）最小二乘图像重建算法- 27 -第四章 滤波反投影重建算法- 30 -4.1平行束滤波反投影重建算法- 31 -卷积反投影重建算法- 31 -滤波函数的选择- 32 -4.2 扇形束滤波反投影算法- 34 -等夹角扇形束滤波反投影算法- 35 -等间距扇形滤波反投影算法- 36 -4.3滤波反投影重建算法的软件实现- 37 -4.4滤波反投影

10、算法和ART算法的对比- 49 -第五章 总结与展望- 50 -5.1总结- 50 -5.2新兴的迭代算法- 51 -5.3 展望- 52 -参考文献- 53 -第一章 绪论1.1 CT技术与图像重建概述所谓的CT技术就是X射线计算机断层成像，它是一种新兴技术，发展于20世纪80年代，它将射线光电子学、信息科学、微电子学、计算机科学等学科结合在一起。我们知道，当X射线照射不同物体时，每个物体对这种射线的吸收和透射率不同，而重建正是利用这种原理，射线照射后，利用探测器进行接收，这样我们就可以根据衰减数得到其分布图像，这就是CT成像技术的基础。CT技术在对物体进行检测时，不用破坏物体内部结构,正是

11、由于这种无损检测技术，所以被广泛应用于医学、生物、航天、航空等多个领域。图像重建是由物体的截断面向该平面做投影，根据投影所得的函数来重建截断面的过程。随着时间的推移，人类在科学技术上也有了重大进步，尤其是计算机技术的高速发展，也推动了图像重建的发展，在医学领域应用最为显著，它大大的丰富了对于人体内脏器官检查的手段，为更加准确地诊断疾病提供了强有力的依据。根据原始数据获取方法和重建原理不同，可分为透射断层重建成像、发射断层重建成像、反射断层重建成像1。1.2 CT和重建技术的发展及研究现状早在1895年，伦琴发现了X射线，这就是CT技术发展的早期萌芽。自此之后，人们也意识到CT技术在成像上有很好

12、地发展前景。果如其然，不就之后人们就看到了诸如旋转阳极X射线管、影像增强管等东西。然而，图像重建技术应用于X放射线医学却是一个重大突破，对身体内脏疾病诊断有着非凡的意义。1917年，奥地利数学家J.Radon发表了图像投影理论，可是由于当时科学技术受限，所以没能将其实现。到了二十世纪六十年代，由于计算机技术的迅猛发展，图像重建技术被不少学者用来创造性地探索研究，在这其间，英国物理学家Cormack提出了一种方案使得X射线投影在医学领域中的应用变成可能。1967年，在英国EMI实验中心从事计算机和重建技术研究的Hounsfield，他发现X射线从不同方向透过物体后，对这些衰减量做测试就能得到物体

13、内部结构，通过不懈的研究，终于在1971年，他研制出了世界上第一台临床使用的计算机断层扫描装置，即CT机。1972年，这台CT成功为一名妇女诊断。颅脑CT的安装与使用，标志CT时代的到来。图1-1和1-2分别表示CT扫描成像设备与临床图像。图1-1 CT扫描成像设备 图l-2 CT临床图像20世纪80年代，可以不间断采集投影数据的螺旋式CT的诞生与使用，能获得高质量的三维图像。这期间，我国分别于1983年和1985年研制出了第一台颅脑扫描装置和第二代CT设备。1989年，我国又成功研制开发出了第一台全身常规CT-2000。1992年，我国生产出了第一台螺旋扫描设备。随着人类的研究探索，多层螺旋

14、在单层的基础上发展起来，它具有覆盖面积更广，扫描速度更快，得到更高质量三维图像，辐射低等特点。然而，被业界誉为CT技术的一重大突破的双能CT于2005年11月诞生了，它能以较低的剂量重建出具有很高时间分辨率的图像。1.3 研究的目的与意义图像重建是利用物体截断面的投影来重建它本身的图像，它不需要将物体解剖，而正是由于这种先进的无损检测技术，所以被广泛地应用于检测和观察中。现实生活中，有些物体在构建其本身图像时不能破坏它的物理性，因此图像重建尤其在医学领域中，它作为一种先进的检测技术，能更准确地检测人体内部的器官。图像重建其实可以看做是一类特殊的图像复原技术。现实中，有多种重建算法，不同的重建算

15、法得出不同质量的重建图像。目前，滤波反投影算法在CT断层图像重建中应用最为广泛，因为它具有很高的精度，能够快速实现，且它的基本算法容易在软件和硬件上实现。第二章 CT成像原理和图像重建算法2.1 CT成像原理与系统组成 一般来说，CT成像过程大致可归纳为:用仪器扫描需要检测的物体,由探测器接收透过物体的X射线,经过模数转换后成为数字信号,即原始数据。然后经过射束硬化、解除零点漂移等预处理,以便获得更精确投影数据；再经过图像重建与处理后,可得到断层图像输出显示或存储至设备。简单地说，CT成像系统组成如下：(1)扫描设备，包括滑环和(平板)检测器，X射线管，模数信号转换系统等2；(2)原始数据处理

16、与重建设备，即电子计算机系统；(3)图像的存储设备以及用来显示的设备。 CT成像过程以及其系统组成如图2-1所示。2-1 CT系统组成与成像过程2.2 CT成像系统扫描方式的发展自1971年9月世界上第一台CT装置诞生以后，CT装置很快推广到其他领域，由于各领域的不同需要，系统的扫描方式也发生很大变化。目前为止，通用的扫描方式有四种：第一代CT扫描方式最简单，就是平行束扫描。它由单一的射线管和在测试区相对应的探测器。扫描中，射线源和探测器同步旋转和平移运动时，探测器记录X射线衰减情况，在不同角度下，会产生不同投影数据，但该光片不能得出被检测物体的衰减分布情况，需要分度旋转和重复平移。如图2-2

17、所示。 2-2 第一代CT系统扫描方式 2-3 第二代CT系统扫描方式第二代与第一代CT扫描方式不同的是平移的次数减少了。该系统包括一个产生窄角扇束射线源和小型的探测器阵列，由于扇束的扇角比较小，因此需要将射线源和探测器旋转前进行一定笔束平移。如图2-3所示。第三代CT系统使用扇束扫描方式，它是使用产生大角度的X射线源，且使用探测器数量较多的阵列。这样，射线源与探测器组成的扇形可以覆盖整个被检测物体，不需平移，只需旋转就可以。如图2-4所示。第四代CT系统同样也是扇形束扫描方式，探测器阵列使用固定环形探测器阵列，如图2-5所示。它仅仅是射线源进行旋转。被检测物体放在环形探测器的中心，扇形束扇角

18、的大小可以根据检测区域来适当调整，射线源每次进行一定分度的旋转，最后获得能够重建被检测物体内部图像的投影数据。 2-4 第三代CT系统扫描方式 2-5 第四代CT系统扫描方式2.3 CT断层图像原理根据投影图像的重建原理以及方法，可以分为不同的类型。我们主要研究透射成像原理，即断层成像技术。断层摄影成像的方式通常为发射接收，即用平行的X光线从不同的方向对物体进行照射，然后记录每个方向的透射场，如图2-6所示。 2-6断层图像获取示意图 2-7图像在角下投影示意图其中，X射线源产生平行的X射线对物体照射，设入射光强度为 ，接收器阵列得到的透射光度为。依次类推，X源和X射线接收器沿中心转一个角度到

19、，这样就可以一组投影数据为，其中的范围是到，每隔一定间隔给予改变，这样我们可以得到投影向量。在数学上，图像的投影可以描述如下。如图2-7所示，图像函数为，穿过该线的一条线称为射线。如果将函数在某条射线上的积分值化成一个集合，这个集合所表示的含义就是投影值。从原点向射线作一条垂线，此垂线作为新坐标的一个轴，并构成一个新坐标系，这样可以得出坐标系仅是坐标系旋转角的结果，二者存在下列变换关系3：(2-1)因而射线积分表达为：=(2-2) 断层成像技术就是从不同射线角度，不同检测器的位置的许多投影值重建原始图像的过程，反应了处的密度。2.4图像重建算法概述图像重建在CT技术中发挥着重要作用。本质上说，

20、它是按照采集后的数据，利用电子计算机来求解图像矩阵中的像素，然后重新构建图像的过程。从CT图像重建的角度看,主要分为两大类方法，即解析类方法与迭代类方法。2.4.1解析类方法解析类重建方法是利用解析、变换重建公式来构建重建图像。目前，常用的解析类重建方法主要有:(1)滤波反投影算法(FBP)；(2)直接傅立叶变换算法(FBP)；(3)卷积反投影算法( CBP)。解析类算法因为其容易实现,速度较快,且能重建出高质量的图像的特点,因此CT系统被较为广泛地应用起来，尤其是在医用CT系统中。然而滤波反投影算法应用最为广泛。但是，其最大的不足之处在于对投影数据完备性的要求偏高,如果能够在投影数据输入给解

21、析法之前，把不利于投影数据精确性的因素给与纠正，这样就可以得到满意的重建图像。 2.4.2传统迭代类方法迭代方法于1970年第一次被引入图像重建的领域,它是利用数学级数迭代的原理来完成图像重建的技术。迭代重建算法分别为代数迭代算法与统计迭代算法。代数迭代算法可分为:ART和SIRT型。统计迭代重建算法(SIR)是对投影数据进行准确性分析,通过精细的统计估计,逐步提高重建图像的质量。迭代图像重建技术的最大优势,在于它能够在投影数据信噪比很低的情况下，甚至不完全的条件下，依然能够获得较高质量的重建图像。然而迭代算法是逼近原始图像，因此在实现时，运行时间相对较长且数据存储量大。第三章 CT图像重建算

22、法实现原理的研究3.1图像重建系统中的数学概念及变换 投影与反投影我们知道，射线穿过物体之后，由于物体的吸收或散射作用，致使我们在检测时会发现射线的强度会发生衰减。通常我们用衰减系数表示衰减程度。设一物质是分均匀的，一个面上衰减系数为，射线穿过该物质后，入射强度由变为,射线在面内的路径如下图3-1。 3-1 射线的衰减由Beer定律得出关系式： (3-1)在CT系统中，我们可以把看作是射线在空气中扫描时探测器测得的数据，是射线通过物体后衰减后数据，由于它们都是测量值，可以得： (3-2) 这样我们就把上面的积分集合称之为投影数据。将其推广得出，如果射线从不同方向照射时，它所对应的路径上投影数据

23、都能得到，这样构成一个投影数据集合。图像重建就是利用投影数据集合来计算的过程。 反投影的定义取决于投影的定义。但它不是投影运算的逆运算。数学语言来说，反投影算子不是投影算子的逆算子。 Radon变换及其反变换1917年，Radon变换的提出使CT的逆问题得到了解决，它也是CT图像重建的理论基础。如图3-2所示。 3-2 Radon变换示意图设直线L的方程为： (3-3)其中为直线到源点的距离，为轴正方向与直线垂线的角度，则与平面上直线L互相对应。记为要重建的图像函数。有变换式如下，(3-4) 称该式为的Radon变换。其实这相当于一个算子，记为，它将空间中函数与空间中函数联系在一起。因为是沿的

24、积分，可以认为空间上的一个点相应于空间上的一条直线L，用向量记，或者 (3-5) 经过变换，上式变为 (3-6) 这样，相应的重建过程就是Radon反变换。3.1.3傅里叶变换一元连续函数傅里叶变换式定义为：(3-7) 其中，。若给定的为连续函数，则反函数也存在，用傅里叶逆变换得到其原函数，即：(3-8) 上俩式构成傅里叶变换对。上述式子推广到二维中，设二元连续函数为，其傅里叶变换及其反变换定义式下： (3-9) (3-10) 对于傅里叶变换，有一个重要定理，如下： (3-11)以上式子就是我们常用的卷积定理。3.1.4中心切片定理 中心切片定理是断层成像的理论基础。二维图像的中心切片定理指出

25、：二维函数图像在视角为时的投影的一维傅里叶变换,的二维傅里叶变换与探测器平行的方向，而且过原点的的一个切片。是切片与轴所成的夹角。原理如下图3-3所示。3-3 二维中心切片原理图 若探测器绕物体旋转至少180度，物体的二维傅里叶变换所对应的探测器方向的中心片段就可以将整个傅里叶空间进行覆盖4。3.2解析类重建算法解析类算法因其利用变换、重建公式，速度较快，实现容易，因此被广泛应用。直接傅里叶变换算法直接傅里叶变换算法是利用二维傅里叶变换来连接极变量函数与，它是直接利用投影定理。如下图3-4。3-4 投影参数说明图将极坐标形式Radon变换改写成直角坐标系下的形式，即：(3-12)其中表示平行束

26、投影。进一步转化有：(3-13)使用该方法进行图像重建的过程是，先采集投影，再傅里叶变换，将投影放在傅里叶空间上，处理完所有投影后，在进行二维傅里叶逆变换，最后得到最终结果。但该算法一般不用于重建，仅作为理论研究，因为它的逆变换要得到星状分布的变换插值为均匀网格形式，这是很困难的，同时它的计算量大，不能即刻成像，且重建后图像质量比较差。反投影重建算法反投影算法是最基本、最简单的算法，它的基本思想是断层内该点的密度值就是经过该平面上的一点的射线投影和。我们将“取投影”“反投影重建”“重建后图像”看作一个系统，得出系统原理图如下图3-5所示。 取投影反投影重建原理重建后图像 3-5 反投影系统原理

27、图设一位于坐标原点上的点源为断层图像上的唯一像点。当系统扫描方式为平移旋转时，即射线先平移，再旋转一定角度，知道累计转角为为止。其中，为原坐标系与旋转坐标系的夹角，这样，投影位置可由来确定。如下图3-6所示。 3-6 平移/旋转扫描方式所用系统为固定坐标，为旋转坐标系，为极坐标设为离散取值，当时，相应投影为：(3-14)若，相应投影为：(3-15)根据反投影公式的定义，点的图像在坐标系表示为： (3-16)式中，是投影数。上式反映的物理意义是：把经过某点所有的投影值做平均之后就是该点的密度值。如果有限区间将射线扩增至无限条，这样就得到连续投影。这样我们可以得到投影表达式为：(3-17) 式中，

28、积分区间为，因为忽略射线硬化条件下，它与内投影值等效。在输入图像为点源条件下，我们得(3-18)由式看出，反投影重建算法的系统扩展函数不是函数，因此系统不完美，图像密度为零的点，重建后不一定为零，可能致图像失真5。对于最为常用的滤波反投影算法，我将着重在下一章介绍。3.3 迭代类重建算法由于计算机运算能力越来越强，使得迭代算法在应用中也越来越多地得到重视。现实中，迭代法重建图像效果好，空间分辨率高，尤其是当采集的数据不完全时，这一点体现更为突出。然而因为它是渐渐靠近原始图像的，所以它在运算时间上会比较长，同样带来的是数据存储量大。迭代法不是找到一个准确的解析式，这是与解析法最大的不同，一般来说

29、，迭代算法的步骤是：先给断层图像赋一个初始估计值，根据此值算出理论投影值，将理论投影值与实际的相比较，按照一定的原则对原始图像进行修正之后与理论值比较，在修正，如此循环，直到达到满意的效果。其实概括地讲，就是假设，比较，修正。目前，迭代图像重建算法有两类，即代数迭代重建算法与统计迭代重建算法。代数迭代重建算法是以代数方程理论为基础，主要有一般ART算法和同时代数重建算法。而统计迭代重建算法是基于各种统计准则，主要有最小均方误差、最大似然估计等重建算法。3.3.1 代数迭代重建算法(1) ART算法ART算法有时也称为Kaczmarz算法，它的主要思路是让当前所有估算的图像在每一次更新中满足一个

30、方程，在迭代修正过程中，每次只考虑一个投影单元的投影值。原理如下图3-7所示。比较并计算修正量理论投影值实际投影值断层图像估计值迭代循环环环环修改计算3-7 迭代算法计算过程示意图最常用的ART算法是基于交替投影法进行迭代修正的，它的图像更新式为：(3-19)为松弛参数。 (2)同时代数重建算法由于ART算法对图像值进行修正时只依赖一条投影带上数据，因此人们又提出了同时代数重建算法（SART算法），它是在校正像素单元的图像值之前，计算出像素单元上所有投影估计值与实际值的差别，并求出来，再利用平均值对图像进行修正。SART算法的迭代公式如下：(3-20)这样，通过各条投影带上的平均值，可以减小误

31、差，避免对重建结果带来过大影响，同时它又抑制图像重建过程中的噪声。3.3.2 影响代数迭代重建算法的因素ART重建算法之所以不能广泛应用，是因为它是一个反复迭代与修正的过程，计算量大。通过长期探索与研究，得到了影响ART算法的很多因素。如下：(1)基函数的选择；(2)松弛参数的选择；(3)迭代次数的最优设计；(4)数据访问方式。基函数的选择直接影响重建算法的计算量。如果基函数选择过于复杂，重建的图像更加逼近真实的，但是计算量太庞大，如果选择很简单，虽然计算量小，但重建图像质量不高。 ART重建算法就是解线性方程并伴随松弛参数的迭代过程，这样选择合理的松弛参数会很大程度地影响运算速度。现如今，松

32、弛参数的选择的基本标准是。只有在这个区间迭代重建算法才能按照标准并保证它的收敛性6。虽然ART迭代重建是不断重建的过程，但是它的重建质量不与迭代次数成正比，而是在达到某个迭代次数后，质量就开始下降。因此选择最优的迭代次数是必要的，不仅节省时间又有高质量的重建图像。下面我们用一个实例来说明这种情况。我们选用一个圆形工件，该圆形工件内部含有碳、气泡、铅等杂质，工件的一些杂质参数7如下表3.1所示 。 名称半径(cm)圆心 坐标(x,y)吸收系数(1/cm) 铁 5.0(0.0, 0.0) 0.435 碳 2.0 (0.0, 2.5) 0.115 空气 0.5 (0.0, -1.0) 0.000 铅

33、 1.0 (0.0, -3.0) 0.800表3.1 工件的杂质参数表根据参数，我们到模型的二维灰度分布图及y轴上一维分布图如下图3.8所示。3-8 模型的二维灰度分布图及y轴上一维分布图选择松弛因子为0.1,抽样间距为0.1cm，重建的像素为200200的重建图像，得到二维灰度分布与y轴上一维分布图如下图3-9至3-13所示。 3-8 迭代次数为1次 3-9 迭代次数为3次 3-10 迭代次数为5次 3-11 迭代次数为10次 3-12 迭代次数为15次 通过上述图片，我们可以看出第1到第3次迭代图像的质量有明显的改善，第5到第10次的迭代图像质量缓慢变好，从第10次到第15次迭代图像质量明

34、显下降。因此，并不是一味的增加迭代次数就可以改善重建图像质量。 ART重建算法与SARTART算法每次迭代时只有一条投影值，如果在这条射线在投影时引入了噪声，这样对重建图像也会带来误差。而SART重建算法的每个像素校正值是所有射线的误差累计和，这样它能有效地抑制测量数据中的噪声。SART算法的迭代速度比ART算法快很多，大约差一个数量级。ART在选择松弛因子方面范围较大，在(0,2)之间，而SART算法在(0,0.001)之间，选择性较窄。尽管SART算法抗噪性较好，但是由于在选择松弛因子方面比较小，同样使得重建图像很模糊，如下图3-13和3-14所示，在高斯随机噪声为情况下，迭代次数为5次的

35、条件下，ART与SART算法的重建图像。 3-13 松弛因子为0.1迭代次数为5次的ART重建图像 3-14 松弛因子为0.0001迭代次数为5次的SIRT重建图像从上图看出，虽然SART算法重建图像比较光滑，但是图像比较模糊。3.3.4 统计迭代重建算法(1)OS-EM算法在OS-EM（分成有序的子集来求期望值的极大值）算法中，数据被分成不同的子集，按照给定的次序走访不同的子集，每访问一个子集图，像就更新一次8。子集的分割有很大的自由度，可以根据需要进行调整，不必按照投影角度来划分。一般来说，若使用N个子集，收敛速度大约提升N倍。但是根据人们的经验，通常让子集分割达到收敛速度提升10倍左右，

36、在这样条件下噪声破坏不明显。第四章 滤波反投影重建算法滤波反投影算法（FBP）是CT图像重建中的一种经典算法，因为它的计算效率高，需求资源少，所以被广泛地应用在各个领域上，特别是医学领域。反投影重建算法是直接对投影进行重建并得出图像，这种方法重建得出的图像质量不高，容易失真。与反投影重建算法不同的是，滤波反投影算法是先对投影数据进行滤波，再对数据进行反投影重建，这样得出的图像更加清晰准确。根据成像系统采集数据方式来分，又可分为平行束滤波反投影重建算法和扇形束滤波反投影算法9。FBP重建算法有两种计算方法，一种与中心切片定理有密切关系，叫做卷积反投影重建算法。另一种是Radon反变换，我们常用第

37、一种，下面主要对第一种做详细介绍。4.1平行束滤波反投影重建算法4.1.1卷积反投影重建算法卷积反投影算法是以中心切片定理为基础产生的。由中心切片定理知，重建图像的二维傅里叶变换可由在不同视角下的投影的，FBP重建算法坐标示意图如下图4-1所示。 4-1 FBP算法坐标系需要重建图像为：(4-1)如图所示，可以得到(4-2)通过上式变换可得到， (4-3)式中，,。事实上，图像重建的过程可以看成是由一系列坐标变换得到的。由上列推倒知，FBP算法是在一定视角下投影，然后进行滤波投影，再做反投影，把这些反投影值累加就可以得到重建图像。滤波函数的选择滤波函数的选择在FBP算法中扮演一个重要的角色。理

38、想的滤波函数可以使重建图像更加准确清晰。上面式子中的滤波函数是理想化的，不是平方可积的，因此滤波函数无法实现。合理地选择滤波函数是必要的，常用的主要有R-L和S-L两种滤波函数。(1) R-L滤波函数系统函数表示为：(4-4)其中，代表采样间隔，为最高不失真频率。系统冲激函数表为 (4-5)我们在对图像滤波时，常常用离散化的函数，这样将带入上述式子，这样离散化的冲激函数表示为：(4-6)函数的连续空域特性如下图4-2所示。 4-2 冲激函数空域特性通过上式，我们可以看到它形式简单，实用性强，用它来重建的图像轮廓上清楚。但是也有缺点，有吉布斯效应，表现为明显的振荡响应。此外，重建的图像容易受到噪

39、声的影响。(2) S-L滤波函数 此函数缓解了R-L滤波函数的不足，选择了不同的窗函数,所对应的滤波函数也不相同。该滤波函数所对应的系统函数为：(4-7)这里的以及的含义与R-L滤波器中的一样。它的冲激响应函数为：(4-8)同上述R-L滤波函数一样，可得到冲激函数离散化形式为： (4-9) 函数的空域特性如下图4-3所示。 4-3 函数空域特性利用S-L滤波函数来重建图像的优点是，振荡响应小，而且在投影数据有噪声的条件下，它的重建质量也较R-L好，但是在低频段上， R-L滤波函数的重建质量相对较高10。4.2 扇形束滤波反投影算法由于扇形束CT系统探测器阵列方式有圆弧形和直线型，这样扇形射线就

40、分为等角型和等间距型。相应的重建算法就有等夹角型和等间距型扇形束滤波反投影算法。如下图所示4-4和4-5所示。 4-4 等夹角型FBP算法 4-5 等间距型FBP算法下面将对两种方法分别讨论。等夹角扇形束滤波反投影算法等夹角扇形束的参数关系图如下图4-6所示。 4-6 等夹角扇形束参数关系图 其中，表示为扇形区域内部物体，为射线源，弧线为探测器所在位置， 为射线源到中心的距离，为中心射线，扇形的位置由和夹角决定，任一条射线由决定，由于坐标系固定，因此射线位置就取决于。设投影函数为，通过投影函数重建图像。扇形束的FBP重建算法可由平行束的FBP重建算法推导出来。首先，将公式变换得到：(4-10)

41、其中，,为物体断层内任一点的极坐标。通过上式可以看出，扇形束的FBP重建算法是对平行束FBP重建算法的加权与修正。等间距扇形滤波反投影算法 等间距扇形束参数关系图如图4-7所示。 4-7 等间距扇形束参数关系图同样的，是扇形区域内部的物体，线段为探测器所在位置，为任意射线，为射线源与中心距离，由于坐标系固定，这样由确定。同样是虚拟探测器。设投影函数为，现利用关系求重建图像。同等夹角扇束FBP重建表达式推导过程一样，由于平行束的数据与扇形束的数据不同，只需将一些变量修改，经过变换，最后得到的等间距扇束FBP重建算法的表达式为：(4-11)其中，是等效投影，为断层内任一点极坐标。4.3滤波反投影重

42、建算法的软件实现这里，我们选择MATLAB与Visual C+相结合的方式对图像进行重建。我们知道，Visual C+软件更注重类的创建，而且通过在每个类中编写相应类的程序能更加深入地了解图像重建的过程，然而对于MATLAB软件来说，因为它本身具有相应的函数，因此利用它能简洁地分析参数对重建图像的影响。首先我们选用经典的shepp-logan头部模型来研究，该模型是由11个不同密度的椭圆叠加而成的，这种模型被用来进行图像的仿真与评估算法的性能。在Visual C+软件中，根据面向对象的思想，在photostar平台中创建CCTEmulate类来实现仿真。实验中，读入图像为灰度图像即可，在CCT

43、Emulate类中设置模型函数的菜单消息，要确定像素在哪个椭圆中。Shepp-logan函数的思路是：检测图像的每个像素是否依次在某一个椭圆中，若在椭圆旋转角度外，置为背景色，若在某个椭圆内，置为其灰度值。这些操作都在归一化坐标系中完成。在相应类中，编写程序如下。BOOL CCTEmulate:SheepLogan()ASSERT(m_pimageproject != NULL); int nWidth = pImageObject->GetWidth(); int nHeight = pImageObject->GetHeight();int nNumBits = pImageO

44、bject->GetNumBits();int nWidthBytes; char *pBuffer = (char *) pImageObject->GetDIBPointer( &nWidthBytes );if( pBuffer = NULL ) return(FALSE);BITMAPFILEHEADER *pBFH; BITMAPINFOHEADER *pBIH; RGBQUAD *pRGBPalette; unsigned char *pBits; int nNumColors = pImageObject->GetNumColors();pBFH = (

45、BITMAPFILEHEADER *) pBuffer;pBIH = (BITMAPINFOHEADER *) &pBuffersizeof(BITMAPFILEHEADER); pRGBPalette = (RGBQUAD *) &pBuffersizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER);pBits = (unsigned char *) &pBuffersizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER)+nNumColors*sizeof(RGBQUAD);for

46、 (int i=0;i<nWidth;i+) for (int j=0;j<nHeight;j+) int k=0; switch(k) case 0: /椭圆aif (EllipseCaculate(0,0,0.92,0.69,0,1,i,j)>1) pBits(nWidth-1-i)*nWidthBytes +j = 0*255; break;else pBits(nWidth-1-i)*nWidthBytes +j = 1*255; case 1:/椭圆bif (EllipseCaculate(0,-0.0184,0.874,0.6624,0,1,i,j)>1)

47、break;else pBits(nWidth-1-i)*nWidthBytes +j = 0.5*255;case 2: /椭圆c if (EllipseCaculate(0.22,0,0.31,0.11,0.3090,0.9511,i,j)<=1) pBits(nWidth-1-i)*nWidthBytes +j = int(0.3*255); break; 以下椭圆本分省略。 在产生投影数据时，要确定合理的投影数，即在归一化坐标系中选择理想的角度进行投影，这里我们将图像的指定直线像素置为255，每10度画线一次，根据要求编写程序。int nProjectionHalfLength

48、= (nWidth>nHeight)? nWidth:nHeight;nProjectionHalfLength = (int)(nProjectionHalfLength+1)/2.0); int nProjectionLength = 2*nProjectionHalfLength -1; for (int i = 0;i < 180; i+) /进行投影，每一度一次 float r1 = i/180.0*3.1415; /把角度转换为弧度 float r2 = (i+180)/180.0*3.1415;for (int k = 0; k <nProjectionLeng

49、th; k+)/法线长度 float L = float(k - nProjectionHalfLength + 1); /原图像的列,从0开始计算 int n = 0; while (k<nProjectionHalfLength-1 && n<nWidth) /r2 float x = n-X; float y = (L - x*cos(r2)/sin(r2); int m = int(Y - y); n+; if (m<0 | m>=nHeight) continue; pProjDai*nProjectionLength + k = pProjD

50、ai*nProjectionLength + k + pBits(nHeight-1-m)*nWidthBytes + n; n = 0; /对每一条直线x轴都要扫描一次 while (k>=nProjectionHalfLength-1 && n<nWidth) /r1 float x = n-X; float y = (L - x*cos(r1)/sin(r1); int m = int(Y - y);开始计算 n+; if (m<0 | m>=nHeight) continue; pProjDai*nProjectionLength + k = p

51、ProjDai*nProjectionLength + k + pBits(nHeight-1-m)*nWidthBytes + n; 选择滤波反投影重建图像时，需要解决滤波器的采样频率，从空间坐标系与归一化坐标系的关系看，滤波器采样频率在图像采样频率的2倍以上即可。我们接下来进行构造滤波器，编写相应的程序。根据以上基本思想，我们编写程序图，进入界面如下图4-8所示。4-8 编写程序运行图运行之后如下图4-9所示，4-9生成TEST文本图选择仿真按钮，生成shepp-logan模型图片如下4-10所示， 4-10shepp-logan模型图 4-11CT扫描过程图对shepp-logan模型图

52、片进行重建仿真，得到过程图片如上图4-11所示。重建仿真之后，选择滤波反投影成像之后，我们可以看到有两种滤波函数，一种是R-L滤波器，另一种是S-L滤波器。如下图4-12所示。4-12两种滤波函数界面选择R-L滤波器进行仿真，得到图像如下图4-13所示。 4-13 R-L滤波重建图像 4-14 S-L滤波重建图像选择S-L滤波器进行仿真，得到图像如上图4-14所示。 通过滤波反投影重建图像之后，我们可以看到，图像能大致反应原来图像的特性，但是分辨率比较低，原始图像中的下三个椭圆几乎无法形成，图像的上下边缘有大面积的白色伪迹。相比于R-L滤波器，S-L滤波器重建效果要好些，白色杂点的状况相对较少

53、11。其实我们可以用MATLAB仿真软件进行，因为MATLAB里面含有自带的函数，只需调整相应参数即可。下面用MATLAB来实现图像重建，并分析投影数对重建图像的影响。用MATLAB中图像处理工具箱中的phantom生成shepp-logan模型。程序：P=phantom(256); Imshow(P)用MATLAB中的radon函数获得shepp-logan模型的投影数据。程序：thetal1=0:10:170;R1,xp=radon(P,thetal); thetal2=0:5:175;R2,xp=radon(P,theta2); thetal3=0:2:178;R3,xp=radon(P,theta3);显示18、36、90个角度投影数据程序如下：figure,imagesc(theta1,xp,R1);xlabel('the