CT断层图像重建算法研究.pdf

《CT断层图像重建算法研究.pdf》由会员分享，可在线阅读，更多相关《CT断层图像重建算法研究.pdf（38页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

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

2、对投影数据完备性要求高。迭代法是利用求解线性方程组来重建图像，它能够在投影数据信噪较低条件下，获得高质量图像。本文将从原理、应用、与优缺点的角度来分析两种算法，重点对解析法中的滤波反投影算法从平行束与扇束投影方式进行研究，最后通过 Visual C+与MATLAB 软件相结合的方式对图像重建，并分析各参数对重建图像的影响。关键字：CT 技术 图像 重建算法 滤波反投影算法 Abstract CT technology is a emerging technology that blend of the Ray optoelectronics,microelectronics and infor

3、matics subject.Because of its advanced nondestructive 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

4、rebuild the high quality images,depends on the 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 has the characteristics of implemen

5、tating easily and fast,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 article we w

6、ill be 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 software

7、 of Visual c+with MATLAB software to image reconstruction,and analyzes the influence of various parameters on the reconstruction image Key words:CT technology image reconstruction algorithm Filtered Backprojection Algorithm 目录 第一章 绪论 CT 技术与图像重建概述 所谓的 CT 技术就是 X 射线计算机断层成像，它是一种新兴技术，发展于 20 世纪 80 年代，它将射线

8、光电子学、信息科学、微电子学、计算机科学等学科结合在一起。我们知道，当 X 射线照射不同物体时，每个物体对这种射线的吸收和透射率不同，而重建正是利用这种原理，射线照射后，利用探测器进行接收，这样我们就可以根据衰减数得到其分布图像，这就是 CT 成像技术的基础。CT 技术在对物体进行检测时，不用破坏物体内部结构,正是由于这种无损检测技术，所以被广泛应用于医学、生物、航天、航空等多个领域。图像重建是由物体的截断面向该平面做投影，根据投影所得的函数来重建截断面的过程。随着时间的推移，人类在科学技术上也有了重大进步，尤其是计算机技术的高速发展，也推动了图像重建的发展，在医学领域应用最为显着，它大大的丰

9、富了对于人体内脏器官检查的手段，为更加准确地诊断疾病提供了强有力的依据。根据原始数据获取方法和重建原理不同，可分为透射断层重建成像、发射断层重建成像、反射断层重建成像1。CT 和重建技术的发展及研究现状 早在 1895 年，伦琴发现了 X 射线，这就是 CT 技术发展的早期萌芽。自此之后，人们也意识到 CT 技术在成像上有很好地发展前景。果如其然，不就之后人们就看到了诸如旋转阳极 X 射线管、影像增强管等东西。然而，图像重建技术应用于 X 放射线医学却是一个重大突破，对身体内脏疾病诊断有着非凡的意义。1917 年，奥地利数学家发表了图像投影理论，可是由于当时科学技术受限，所以没能将其实现。到了

10、二十世纪六十年代，由于计算机技术的迅猛发展，图像重建技术被不少学者用来创造性地探索研究，在这其间，英国物理学家 Cormack 提出了一种方案使得 X 射线投影在医学领域中的应用变成可能。1967 年，在英国 EMI 实验中心从事计算机和重建技术研究的Hounsfield，他发现 X 射线从不同方向透过物体后，对这些衰减量做测试就能得到物体内部结构，通过不懈的研究，终于在 1971 年，他研制出了世界上第一台临床使用的计算机断层扫描装置，即 CT机。1972 年，这台 CT 成功为一名妇女诊断。颅脑 CT 的安装与使用，标志 CT 时代的到来。图 1-1 和 1-2 分别表示 CT 扫描成像设

11、备与临床图像。图 1-1 CT 扫描成像设备 图 l-2 CT 临床图像 20 世纪 80 年代，可以不间断采集投影数据的螺旋式 CT 的诞生与使用，能获得高质量的三维图像。这期间，我国分别于 1983 年和1985 年研制出了第一台颅脑扫描装置和第二代 CT 设备。1989 年，我国又成功研制开发出了第一台全身常规 CT-2000。1992 年，我国生产出了第一台螺旋扫描设备。随着人类的研究探索，多层螺旋在单层的基础上发展起来，它具有覆盖面积更广，扫描速度更快，得到更高质量三维图像，辐射低等特点。然而，被业界誉为 CT 技术的一重大突破的双能 CT 于 2005 年 11 月诞生了，它能以较

12、低的剂量重建出具有很高时间分辨率的图像。研究的目的与意义 图像重建是利用物体截断面的投影来重建它本身的图像，它不需要将物体解剖，而正是由于这种先进的无损检测技术，所以被广泛地应用于检测和观察中。现实生活中，有些物体在构建其本身图像时不能破坏它的物理性，因此图像重建尤其在医学领域中，它作为一种先进的检测技术，能更准确地检测人体内部的器官。图像重建其实可以看做是一类特殊的图像复原技术。现实中，有多种重建算法，不同的重建算法得出不同质量的重建图像。目前，滤波反投影算法在 CT 断层图像重建中应用最为广泛，因为它具有很高的精度，能够快速实现，且它的基本算法容易在软件和硬件上实现。第二章 CT 成像原理

13、和图像重建算法 CT 成像原理与系统组成 一般来说，CT 成像过程大致可归纳为:用仪器扫描需要检测的物体,由探测器接收透过物体的 X 射线,经过模数转换后成为数字信号,即原始数据。然后经过射束硬化、解除零点漂移等预处理,以便 获得更精确投影数据；再经过图像重建与处理后,可得到断层图像输出显示或存储至设备。简单地说，CT 成像系统组成如下：(1)扫描设备，包括滑环和(平板)检测器，X 射线管，模数信号转换系统等2；(2)原始数据处理与重建设备，即电子计算机系统；(3)图像的存储设备以及用来显示的设备。CT 成像过程以及其系统组成如图 2-1 所示。2-1 CT 系统组成与成像过程 CT 成像系统

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

15、系统包括一个产生窄角扇束射线源和小型的探测器阵列，由于扇束的 扇角比较小，因此需要将射线源和探测器旋转前进行一定笔束平移。如图 2-3 所示。第三代 CT 系统使用扇束扫描方式，它是使用产生大角度的 X 射线源，且使用探测器数量较多的阵列。这样，射线源与探测器组成的扇形可以覆盖整个被检测物体，不需平移，只需旋转就可以。如图 2-4 所示。第四代 CT 系统同样也是扇形束扫描方式，探测器阵列使用固定环形探测器阵列，如图 2-5 所示。它仅仅是射线源进行旋转。被检测物体放在环形探测器的中心，扇形束扇角的大小可以根据检测区域来适当调整，射线源每次进行一定分度的旋转，最后获得能够重建被检测物体内部图像

16、的投影数据。2-4 第三代 CT 系统扫描方式 2-5 第四代 CT 系统扫描方式 CT 断层图像原理 根据投影图像的重建原理以及方法，可以分为不同的类型。我们主要研究透射成像原理，即断层成像技术。断层摄影成像的方式通常为发射接收，即用平行的 X 光线从不同的方向对物体进行照射，然后记录每个方向的透射场，如图 2-6所示。2-6 断层图像获取示意图 2-7 图像在角下投影示意图 其中，X 射线源产生平行的 X 射线对物体照射，设入射光强度为 0I，接收器阵列得到的透射光度为I1,t。依次类推，X 源和 X射线接收器沿中心转一个角度到2，这样就可以一组投影数据为 2,tI，其中的范围是0到180

17、，每隔一定间隔给予改变，这样我们可以得到投影向量。在数学上，图像的投影可以描述如下。如图 2-7 所示，图像函数为yxf,，穿过该线的一条线称为射线。如果将函数在某条射线上的积分值化成一个集合，这个集合所表示的含义就是投影值。从原点向射线作一条垂线，此垂线作为新坐标的一个轴t，并构成一个新坐标系st,，这样可以得出st,坐标系仅是yx,坐标系旋转角的结果，二者存在下列变换关系3：cossinsincosstyx(2-1)因而射线积分表达为：P1t=AByxf,dxdy=smsmstf,1ds(2-2)断层成像技术就是从不同射线角度，不同检测器的位置的许多投影值Pt重建原始图像yxf,的过程，y

18、xf,反应了yx,处的密度。图像重建算法概述 图像重建在 CT 技术中发挥着重要作用。本质上说，它是按照采集后的数据，利用电子计算机来求解图像矩阵中的像素，然后重新构建图像的过程。从 CT 图像重建的角度看,主要分为两大类方法，即解析类方法与迭代类方法。2.4.1 解析类方法 解析类重建方法是利用解析、变换重建公式来构建重建图像。目前，常用的解析类重建方法主要有:(1)滤波反投影算法(FBP)；(2)直接傅立叶变换算法(FBP)；(3)卷积反投影算法(CBP)。解析类算法因为其容易实现,速度较快,且能重建出高质量的图像的特点,因此CT 系统被较为广泛地应用起来，尤其是在医用 CT 系统中。然而

19、滤波反投影算法应用最为广泛。但是，其最大的不足之处在于对投影数据完备性的要求偏高,如果能够在投影数据输入给解析法之前，把不利于投影数据精确性的因素给与纠正，这样就可以得到满意的重建图像。2.4.2 传统迭代类方法 迭代方法于 1970 年第一次被引入图像重建的领域,它是利用数学级数迭代的原理来完成图像重建的技术。迭代重建算法分别为代数迭代算法与统计迭代算法。代数迭代算法可分为:ART 和 SIRT型。统计迭代重建算法(SIR)是对投影数据进行准确性分析,通过精细的统计估计,逐步提高重建图像的质量。迭代图像重建技术的最大优势,在于它能够在投影数据信噪比很低的情况下，甚至不完全的条件下，依然能够获

20、得较高质量的重建图像。然而迭代算法是逼近原始图像，因此在实现时，运行时间相对较长且数据存储量大。第三章 CT 图像重建算法实现原理的研究 图像重建系统中的数学概念及变换 3.1.1 投影与反投影 我们知道，射线穿过物体之后，由于物体的吸收或散射作用，致使我们在检测时会发现射线的强度会发生衰减。通常我们用衰减系数表示衰减程度。设一物质是分均匀的，一个面上衰减系数为yx,，射线穿过该物质后，入射强度由0I变为I,射线在面内的路径如下图 3-1。3-1 射线的衰减 由 Beer 定律得出关系式：LdlyxeII,0 (3-1)在 CT 系统中，我们可以把0I看作是射线在空气中扫描时探测器测得的数据，

21、I是射线通过物体后衰减后数据，由于它们都是测量值，可以得：IIdlyxlpL0ln,(3-2)这样我们就把上面的积分集合称之为投影数据。将其推广得出，如果射线从不同方向照射时，它所对应的路径上投影数据lp都能得到，这样构成一个投影数据集合。图像重建就是利用投影数据lp集合来计算yx,的过程。反投影的定义取决于投影的定义。但它不是投影运算的逆运算。数学语言来说，反投影算子不是投影算子的逆算子。3.1.2 Radon 变换及其反变换 1917 年，Radon 变换的提出使 CT 的逆问题得到了解决，它也是CT 图像重建的理论基础。如图 3-2 所示。3-2 Radon 变换示意图 设直线 L 的方

22、程为：sincosyxp (3-3)其中p为直线到源点的距离，为x轴正方向与直线垂线的角度，则,p与oxy平面上直线 L 互相对应。记yxf,为要重建的图像函数。有变换式如下，dlyxfpLf,(3-4)称该式为yxf,的 Radon 变换。其实这相当于一个算子，记为R，它将yx,空间中函数与,p空间中函数联系在一起。因为,p是f沿L的积分，可以认为,p空间上的一个点相应于yx,空间上的一条直线 L，用向量记yxX,，或者 cossinsincostpytpx (3-5)经过变换，上式变为 dttptpfpfcossin,sincos,(3-6)这样，相应的重建过程就是 Radon 反变换。3

23、.1.3 傅里叶变换 一元连续函数傅里叶变换式定义为：dxexfFxi2(3-7)其中，1i。若给定的 F为连续函数，则反函数也存在，用傅里叶逆变换得到其原函数 xf，即：dueFFFxfxi211(3-8)上俩式构成傅里叶变换对。上述式子推广到二维中，设二元连续函数为yxf,，其傅里叶变换及其反变换定义式下：dxdyeyxfFyxi21221,(3-9)21221211221,ddeFFFyxfyxi(3-10)对于傅里叶变换，有一个重要定理，如下：vuHvuFyxhyxfF,(3-11)以上式子就是我们常用的卷积定理。3.1.4 中心切片定理 中心切片定理是断层成像的理论基础。二维图像的中

24、心切片定理指出：二维函数图像yxf,在视角为时的投影 rxP的一维傅里叶变换,yxf,的二维傅里叶变换,21 FF与探测器平行的方向，而且过原点的的一个切片。是切片与1轴所成的夹角。原理如下图 3-3 所示。3-3 二维中心切片原理图 若探测器绕物体旋转至少 180 度，物体的二维傅里叶变换21,F所对应的探测器方向的中心片段就可以将整个傅里叶空间进行覆盖4。解析类重建算法 解析类算法因其利用变换、重建公式，速度较快，实现容易，因此被广泛应用。3.2.1 直接傅里叶变换算法 直接傅里叶变换算法是利用二维傅里叶变换来连接极变量函数,rf与,F，它是直接利用投影定理。如下图 3-4。3-4 投影参

25、数说明图 将极坐标形式 Radon 变换改写成直角坐标系下的形式，即：dxdytxyyxftPP sincos,(3-12)其中 PP表示平行束投影。进一步转化有：yxyxjyxddefFyxFyxfyxFF 22212,(3-13)使用该方法进行图像重建的过程是，先采集投影，再傅里叶变换，将投影放在傅里叶空间上，处理完所有投影后，在进行二维傅里叶逆变换，最后得到最终结果。但该算法一般不用于重建，仅作为理论研究，因为它的逆变换要得到星状分布的变换插值为均匀网格形式，这是很困难的，同时它的计算量大，不能即刻成像，且重建后图像质量比较差。3.2.2 反投影重建算法 反投影算法是最基本、最简单的算法

26、，它的基本思想是断层内该点的密度值就是经过该平面上的一点的射线投影和。我们将“取投影”“反投影重建”“重建后图像”看作一个系统，得出系统原理图如下图 3-5 所示。3-5 反投影系统原理图 取投影 反投影重建 原理 重建后图像 设一位于坐标原点上的点源为断层图像上的唯一像点。当系统扫描方式为平移旋转时，即射线先平移，再旋转一定角度，知道累计转角为180为止。其中，为原坐标系与旋转坐标系的夹角，这样，投影位置可由,rx来确定。如下图 3-6 所示。3-6 平移/旋转扫描方式所用系统 yx,为固定坐标，rryx,为旋转坐标系，,r为极坐标 设为离散取值，当1时，相应投影为：11cos|,11rxd

27、yyxydyxfxpxprrrrrrrrrr(3-14)若n，相应投影为：nrpncos(3-15)根据反投影公式的定义，点rryx,的图像在坐标系表示为：Niirp1cos11(3-16)式中，N，N是投影数。上式反映的物理意义是：把经过某点所有的投影值做平均之后就是该点的密度值。如果有限区间将射线扩增至无限条，这样就得到连续投影。这样我们可以得到投影表达式为：drrf0cos1,(3-17)式中，积分区间为,0，因为忽略射线硬化条件下，它与2,0内投影值等效。在输入图像为点源条件下，我们得 22111yxr(3-18)由式看出，反投影重建算法的系统扩展函数不是函数，因此系统不完美，图像密度

28、为零的点，重建后不一定为零，可能致图像失真5。对于最为常用的滤波反投影算法，我将着重在下一章介绍。迭代类重建算法 由于计算机运算能力越来越强，使得迭代算法在应用中也越来越多地得到重视。现实中，迭代法重建图像效果好，空间分辨率高，尤其是当采集的数据不完全时，这一点体现更为突出。然而因为它是渐渐靠近原始图像的，所以它在运算时间上会比较长，同样带来的是数据存储量大。迭代法不是找到一个准确的解析式，这是与解析法最大的不同，一般来说，迭代算法的步骤是：先给断层图像赋一个初始估计值，根据此值算出理论投影值，将理论投影值与实际的相比较，按照一定的原则对原始图像进行修正之后与理论值比较，在修正，如此循环，直到

29、达到满意的效果。其实概括地讲，就是假设，比较，修正。目前，迭代图像重建算法有两类，即代数迭代重建算法与统计迭代重建算法。代数迭代重建算法是以代数方程理论为基础，主要有一般 ART 算法和同时代数重建算法。而统计迭代重建算法是基于各种统计准则，主要有最小均方误差、最大似然估计等重建算法。3.3.1 代数迭代重建算法(1)ART 算法 ART 算法有时也称为 Kaczmarz 算法，它的主要思路是让当前所有估算的图像在每一次更新中满足一个方程，在迭代修正过程中，每次只考虑一个投影单元的投影值。原理如下图 3-7 所示。3-7 迭代算法计算过程示意图 最常用的 ART 算法是基于交替投影法进行迭代修

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

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

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

33、上一维分布图 选择松弛因子为,抽样间距为 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 次迭代图像质量明显下降。因此，并不是一味的增加迭代次数就可以改善重建图像质量。3.3.3 ART 重建算法与 SART ART 算法每

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

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

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

37、束滤波反投影重建算法 4.1.1 卷积反投影重建算法 卷积反投影算法是以中心切片定理为基础产生的。由中心切片定理知，重建图像的二维傅里叶变换可由yxf,在不同视角下的投影的,|,121PxpFFFr固定，FBP 重建算法坐标示意图如下图 4-1 所示。4-1 FBP 算法坐标系 需要重建图像为：21221211221,ddeFFFyxfrfyxi(4-1)和如图所示，cos21，sin22，cosrxr，可以得到 cos221ryx(4-2)通过上式变换可得到，2112,FFyxfrf drg0,cos(4-3)式中，11 Fxhr,11PFxpr,cosrrxpxhrg。事实上，图像重建的过

38、程可以看成是由一系列坐标变换得到的。由上列推倒知，FBP 算法是在一定视角下投影，然后进行滤波投影，再做反投影，把这些反投影值累加就可以得到重建图像。4.1.2 滤波函数的选择 滤波函数的选择在 FBP 算法中扮演一个重要的角色。理想的滤波函数可以使重建图像更加准确清晰。上面式子中的滤波函数是理想化的，不是平方可积的，因此滤波函数无法实现。合理地选择滤波函数是必要的，常用的主要有 R-L 和 S-L 两种滤波函数。(1)R-L 滤波函数 系统函数表示为：BrectWHLR2(4-4)其中 其他,021,12dBBrect，d代表采样间隔，dB21为最高不失真频率。系统冲激函数表为 BxcBBx

39、cBrr222sin2sin2(4-5)我们在对图像滤波时，常常用离散化的函数，这样将ndxr带入上述式子，这样离散化的冲激函数表示为：奇数，偶数ndnnndndhLR22221,00,41(4-6)函数的连续空域特性如下图 4-2 所示。4-2 冲激函数空域特性 通过上式，我们可以看到它形式简单，实用性强，用它来重建的图像轮廓上清楚。但是也有缺点，有吉布斯效应，表现为明显的振荡响应。此外，重建的图像容易受到噪声的影响。(2)S-L 滤波函数 此函数缓解了 R-L 滤波函数的不足，选择了不同的窗函数 W,所对应的滤波函数也不相同。该滤波函数所对应的系统函数为：BrectBBBrectBcHLS

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

41、射线就分为等角型和等间距型。相应的重建算法就有等夹角型和等间距型扇形束滤波反投影算法。如下图所示 4-4 和 4-5 所示。4-4 等夹角型 FBP 算法 4-5 等间距型 FBP 算法 下面将对两种方法分别讨论。4.2.1 等夹角扇形束滤波反投影算法 等夹角扇形束的参数关系图如下图 4-6 所示。4-6 等夹角扇形束参数关系图 其中，,rf表示为扇形区域内部物体，0S为射线源，弧线21DD为探测器所在位置，D为射线源到中心的距离，00DS为中心射线，扇形的位置由00DS和夹角决定，任一条射线ES0由决定，由于坐标系固定，因此射线位置就取决于,。设投影函数为,p，通过投影函数重建图像,rf。扇

42、形束的 FBP 重建算法可由平行束的 FBP 重建算法推导出来。首先，将公式变换得到：dhDpLrf020222sin2cos,1,(4-10)其中，sin222DrrDL，Lrcosarcsin0,r为物体断层内任一点的极坐标。通过上式可以看出，扇形束的 FBP 重建算法是对平行束 FBP 重建算法的加权与修正。4.2.2 等间距扇形滤波反投影算法 等间距扇形束参数关系图如图 4-7 所示。4-7 等间距扇形束参数关系图 同样的，,rf是扇形区域内部的物体，线段21DD为探测器所在位置，BS0为任意射线，D为射线源与中心距离，由于坐标系固定，这样BS0由s,确定。21DD同样是虚拟探测器。设

43、投影函数为,sp，现利用关系求重建图像,rf。同等夹角扇束 FBP 重建表达式推导过程一样，由于平行束的数据与扇形束的数据不同，只需将一些变量修改，经过变换，最后得到的等间距扇束 FBP 重建算法的表达式为：dsgspUrfe202,1,(4-11)其中DrDUsin，22,sDDspspe是等效投影，shsg21，,r为断层内任一点极坐标。滤波反投影重建算法的软件实现 这里，我们选择 MATLAB 与 Visual C+相结合的方式对图像进行重建。我们知道，Visual C+软件更注重类的创建，而且通过在每个类中编写相应类的程序能更加深入地了解图像重建的过程，然而对于 MATLAB 软件来说

44、，因为它本身具有相应的函数，因此利用它能简洁地分析参数对重建图像的影响。首先我们选用经典的 shepp-logan 头部模型来研究，该模型是由 11 个不同密度的椭圆叠加而成的，这种模型被用来进行图像的仿真与评估算法的性能。在 Visual C+软件中，根据面向对象的思想，在 photostar 平台中创建 CCTEmulate 类来实现仿真。实验中，读入图像为灰度图像即可，在 CCTEmulate 类中设置模型函数的菜单消息，要确定像素在哪个椭圆中。Shepp-logan 函数的思路是：检测图像的每个像素是否依次在某一个椭圆中，若在椭圆旋转角度外，置为背景色，若在某个椭圆内，置为其灰度值。这

45、些操作都在归一化坐标系中完成。在相应类中，编写程序如下。BOOL CCTEmulate:SheepLogan()ASSERT(m_pimageproject!=NULL);int nWidth=pImageObject-GetWidth();int nHeight=pImageObject-GetHeight();int nNumBits=pImageObject-GetNumBits();int nWidthBytes;char*pBuffer=(char*)pImageObject-GetDIBPointer(&nWidthBytes);if(pBuffer=NULL)return(FAL

46、SE);BITMAPFILEHEADER*pBFH;BITMAPINFOHEADER*pBIH;RGBQUAD*pRGBPalette;unsigned char*pBits;int nNumColors=pImageObject-GetNumColors();pBFH=(BITMAPFILEHEADER*)pBuffer;pBIH=(BITMAPINFOHEADER*)&pBuffersizeof(BITMAPFILEHEADER);pRGBPalette=(RGBQUAD*)&pBuffersizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER);

47、pBits=(unsigned char*)&pBuffersizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER)+nNumColors*sizeof(RGBQUAD);for(int i=0;inWidth;i+)for(int j=0;j1)pBits(nWidth-1-i)*nWidthBytes+j=0*255;break;else pBits(nWidth-1-i)*nWidthBytes+j=1*255;case 1:/椭圆b if(EllipseCaculate(0,0,1,i,j)1)break;else pBits(nWidth-1

48、-i)*nWidthBytes+j=*255;case 2:/椭圆c if(EllipseCaculate,0,i,j)nHeight)nWidth:nHeight;nProjectionHalfLength =(int)(nProjectionHalfLength+1)/;int nProjectionLength=2*nProjectionHalfLength-1;for(int i=0;i 180;i+)/进行投影，每一度一次 float r1=i/*;/把角度转换为弧度 float r2=(i+180)/*;for(int k=0;k nProjectionLength;k+)/法线长

49、度 float L=float(k-nProjectionHalfLength+1);/原图像的列,从0开始计算 int n=0;while(knProjectionHalfLength-1&nnWidth)/r2 float x=n-X;float y=(L-x*cos(r2)/sin(r2);int m=int(Y-y);n+;if(m=nHeight)continue;pProjDai*nProjectionLength+k=pProjDai*nProjectionLength+k+pBits(nHeight-1-m)*nWidthBytes+n;n=0;/对每一条直线x轴都要扫描一次

50、while(k=nProjectionHalfLength-1&nnWidth)/r1 float x=n-X;float y=(L-x*cos(r1)/sin(r1);int m=int(Y-y);开始计算 n+;if(m=nHeight)continue;pProjDai*nProjectionLength+k=pProjDai*nProjectionLength+k+pBits(nHeight-1-m)*nWidthBytes+n;选择滤波反投影重建图像时，需要解决滤波器的采样频率，从空间坐标系与归一化坐标系的关系看，滤波器采样频率在图像采样频率的 2 倍以上即可。我们接下来进行构造滤波