基于Matlab实现的CT重建算法仿真比较研究毕业设计论文

1、密级：公开南昌大孽NANCHANGUNIVERSITY学士学位论文THESISOFBACHELOR题目基于Matlab实现的CT重建算法仿真比较研究学院：信息工程系电子信息工程专业班级：南昌大学学士学位论文原创性申明本人郑重中明：所呈交的论文是本人在导师的指导下独立进行研究所取得的研究成果。除了文中特别加以标注引用的内容外，本论文不包含任何其他个人或集体已经发表或撰写的成果。对本文的研究作出重要贡献的个人和集体，均已在文中以明确方式表明。本人完全意识到本申明的法律后果由本人承担。作者签名：日期：学位论文版权使用授权书本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保留并向国家有

2、关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅。本人授权南昌大学可以将本论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。保密口，在年解密后适用本授权书。本学位论文属于不保密口。（请在以上相应方框内打）作者签名：日期：日期:导师签名:摘要基于Matlab实现的CT重建算法仿真比较研究摘要计算机断层成像(即CT(ComputerizedTomography)技术是一项近年来快速发展多学科交叉的先进技术，除了医学应用，还应用在工业无损检测、射电天文学、精密仪器反演等多个重要领域。CT图像重建是CT技术的核心，重建算法的优劣直接关系到对检

3、测结果判断的准确性。医学CT普遍采用解析成像技术，这主要因为医用CT可以采集到完整的投影数据。然而在工业应用中，一般得不到完整的投影数据，因此要利用不完整的投影数据得到较清晰的图像，必须采用迭代重建法。针对以上情况，本课题作了以下两方面的分析比较。1,能获取完整投影数据时，用计算机模拟了解析法的图像重建。发现图像重建的速度快，且重建的图像与原图的视觉相似度接近1。2,不能获取完整投影数据时，用计算机模拟了迭代算法的图像重建。先用ART迭代算法引入松弛参数，分析了松弛参数重建图像的影响。接着用SIRT算法重建图像，这种算法所得的图像比ART法较好。最后在ARTS法基础上采用全差分对重建图像进行优

4、化，图像质量有了明显提高。关键字：图像重建算法；解析重建；迭代重建AbstractBasedonMatlabrealizetheCTreconstructionalgorithmsimulationcomparativestudyAbstractComputedtomography(CT(ComputerizedTomography)technologyisarapiddevelopmentinrecentyears,multidisciplinaryadvancedtechnology,inadditiontomedicalapplications,butalsousedinindustri

5、alnondestructivetesting,radioastronomy,precisioninstrumentsandotherimportantareasoftheinversion.CTimagereconstructionCTtechnologyisthecoreofthereconstructionalgorithmisdirectlyrelatedtotheaccuracyofthetestresultsdetermine.MedicalCTimagingcommonlyusedanalyticaltechnique,whichcanbecollectedmainlydueto

6、acompletemedicalCTprojectiondata.However,inindustrialapplications,thegenerallackofcompleteprojectiondata,andthereforetheuseofincompleteprojectiondatatobemoreclearimage,aniterativereconstructionmethodmustbeused.Fortheabove,thisprojectmadethefollowingtwoaspectsoftheanalysisandcomparison.1 .cangetthefu

7、llprojectiondata,usingacomputersimulationoftheanalyticalmethodofimagereconstruction.Foundfastimagereconstructionandthereconstructedimagewiththeoriginalvisualsimilarityiscloseto1.2 .cannotgetthefullprojectiondata,usingacomputersimulationoftheiterativealgorithmforimagereconstruction.Firstwiththeintrod

8、uctionofARTrelaxationparameteriterativealgorithmanalyzestheimpactofrelaxationparametersofthereconstructedimage.FollowedbySIRTalgorithmforreconstructionofimages,theresultingimageofthisalgorithmisbetterthantheARTalgorithm.FinallyARTalgorithmbasedonafullydifferentialoptimizedforimagereconstruction,imag

9、equalityhasimprovedsignificantly.Keywords:Imagereconstructionalgorithm；analyticreconstruction；iterativereconstructionII目录目录摘要IAbstractII第一章绪论11.1 弓I言11.2 CT的分类21.3 CT的组成结构41.4 论文的组织结构4第二章CT图像重建的基本原理62.1 所需数据的扫描获取62.2 图像重建算法6第三章CT图像解析重建算法93.1 反投影重建法（BP）93.2 滤波反投影图像重建法（FBP）113.3 Radon变换与傅立叶变换123.3.1 R

10、adon变换123.3.2 傅立叶变换的定义133.4 中心切片定理143.5 平行投影重建算法173.5.1 算法原理173.5.2 滤波（卷积）反投影183.5.3 Radon反变换重建算法183.6 解析重建算法的实现和对比分析19第四章迭代重建算法224.1 迭代重建算法的思想原理224.1.1 迭代重建的数学模型224.2 代数重建法（ART）254.3 联立迭代重建法（SIRT）264.4 全差分迭彳t算法（TV）274.5 迭代重建算法的实现和对比分析27第五章结论37参考文献38III目录附录41谢辞63IV第一章绪论第一章绪论1.1引言自从X射线发现后，医学上就开始用它来探测

11、人体疾病。但是，由于人体内有些器官对X线的吸收差别极小，因此X射线对那些前后重叠的组织的病变就难以发现。于是，美国与英国的科学家开始了寻找一种新的东西来弥补用X线技术检查人体病变的不足。CT（电子计算机体层成像）是70年代初放射诊断的一项重大突破，CT不是X线摄影，而是用X线对人体扫描（图1.1所示），取得信息，经电子计算机处理而获得的重建图像。它能使传统的X线检查难以显示的器官及其病变显示成像，而且图像逼真，解剖关系明确，从而扩大了人体的检查范围，大大提高了病变的早期检出率和诊断准确率。这种检查简便、安全、无痛苦、无创伤、无危险，它促进了医学影像诊断学的发展，CT的研制成功被誉为自伦琴发现X

12、射线以后，放射诊断学上最重要的成就。发明者亨斯菲尔德和科马克共同获得了1979年的诺贝尔奖金。CT最初只用于头部检查，1974年又出现了全身的C工在短短10余年间，CT已遍及全球，从第一代发展到第五代。我国各大城市医院所使用的C嘉属第三代5。全身CT可以作头、胸、腹、骨盆的横断扫描，也可作甲状腺、脊柱、关节和软组织及五官等小部位的区域扫描。CT最适于查明占位性病变如月中瘤、囊月中、增大的淋巴结、血月中、脓月中和肉芽月中的大小、形态、数目和侵犯范围，它可以决定某些器官癌月中的分期和是否能进行手术切除。在某些情况下，CT还能区别病变的病理特性如实性、囊性、血管性、炎性、钙性、脂肪等。CT成像技术的

13、基本原理是根据人体中不同组织对X射线的衰减作用不同，用X线束对人体某部一定厚度的层面进行扫描，由探测器接收透过该层面的X线,转变为可见光后，由光电转换变为电信号，再经模拟/数字转换器（analog/digitalconverter）转为数字，输入计算机处理。图像形成的处理有如对选定层面分成若干个体积相同的长方体，称之为体素（voxel）1。扫描所得信息经计算而获得每个体素的X线衰减系数或吸收系数，再排列成矩阵，即数字矩阵（digitalmatrix）,数字矩阵可存贮于磁盘或光盘中。经数字/模拟转换器（digital/analogconverter）把数字矩阵中的每个数字转为由黑到白不等灰度的小

14、方块，即像素（pixel）,并按矩阵排列，即构成CT图像。所以，CT图像是重建图像。每个体素的X线吸收系数可以通过不同的数学方法算出。第一章绪论如图1.1所示图1.1计算机断层成像示意图1.2CT的分类在投影重建的原理和方法的框架下，根据投影方式不同产生了不同种类的CT成像术。透射断层成像透射断层成像(TransmissionComputedTomography,TCT,简称CT)系统中，从发射源射出的X射线穿透物体到达接收器。射线在通过物体时被物体吸收一部分，余下部分被接收器接收。由于物体各部分对射线的吸收不同，所以接收器获得的射线强度实际上反映了物体各部分对射线的吸收情况。由于透射成像术中

15、采用X射线，所以也叫X-CT。接收器接收到的模拟信号经模数转换器转换成数字信号后，把代表着不同角度下的投影数据送给计算机，它再运用复杂的数学方法重建出物体截面的图像4。在医学上它主要针对骨质成像。XCT成像原理如图1.2所示：图1.2X射线CT成像原理第一章绪论发射断层成像在发射断层成像(EmissionComputedTomography,ECT系统中，发射源在物体内部，一般是将具有放射性的离子注入物体内部，从物体外检测放射出来的量。通过这种方法可以了解离子在物体内运动情况的分布，从而可以检测到与生理相关的状况。现在通常用的ECT主要有两种：正电子成像PET(PositronEmission

16、Tomography)；单光子成像SPECT(SinglePositronEmissionComputedTomography)。它们的区别是，采用单光子辐射器时，检测器所接收的光子计数便是投影测量数据；而采用正电子辐射器时，由于在正电子辐射点的毫米级范围内，一个正电子会俘获一个负电子，并产生沿相反方向运动的两个中子，因而测量数据是在一个很小时间问隔内由两个性质相反的检测器检测到的光子对的数目。ECT主要针对代谢过程的成像。磁共振成像磁共振成像(MagneticResonanceImaging,MRI)在早期称为核磁共振(NuclearMagneticResonance,NMR)。它的发明者斯

17、坦福大学的F.Bloch和哈佛大学的E.M.Purcell为此获得了1952年度的诺贝尔物理奖。它的工作原理简介如下：氢核以及其它具有奇数个质子和中子的核，包含具有一定磁动量和旋量的质子。如果把它们放在磁场中，它们会像陀螺在地球重力场中一样在磁场中进动。一般情况下质子在磁场中时随意排列着的，当一个适当强度和频率的共振场信号作用于物体时，质子吸收能量并转向与磁场相交的朝向。如果此时将共振场信号除去，质子吸收的能量释放并被接收器检测到。根据检测到的信号就可以确定质子的密度。通过控制所用的共振场信号和磁场强度，可每次检测到沿透过目标中一条线的信号。换句话说，检测到的信号是MRI信号沿直线的积分。所以

18、检测目标的工作成了投影重建的问题。由这些信号可以得到极好的器官解剖图，在细节上远远超过由X射线CT产生的图像。由于它能够很好的保留丰富的软组织图像，所以在脑功能研究上，MRI胜过ECT超声断层成像超声断层成像(UltrasoundCT,UCT与上述的CT不同的地方是超声射线在人体组织中不一定按直线传播，尤其是在软组织和硬组织的分界面会产生折射效应。因此超声成像主要用于人体中没有硬组织(如骨骼)部分的检测，特别是妇女乳房癌的检测。反射断层成像反射断层成像(ReflectionCT,RCT)也是利用投影重建的原理工作的。常见的一个例子是雷达系统，其中的雷达图是物体反射的回波所产生的。例如在前第一章

19、绪论视雷达（FLR中，雷达发射器从空中向地面发射无线电波；雷达接收器在特定角度所接收到的回波强度是地面发射量的一个扫描段的积分。通过发展硬件设备和改进软件算法，可以不断降低噪声和提高图像的清晰度。CT技术在临床的许多领域得到愈来愈广泛的应用。应该指出，上面介绍的各类CT的功能是相互补充而不是相互替代的601.3 CT的组成结构以上几种成像技术尽管使用的物理系统各不相同，但所有系统所测量到的数据均具有物体某种感兴趣的物理特性在空间分布的积分形式。成像技术广义上讲大致可分为三个部分：信息的收集、重建以及处理和传递。数据采集，即获得携带有物体内部特征信息的数据。所采用的物理测量系统包括放射源、检测器

20、及配套的机械、电路设备。为进行计算机处理，系统还需将测量得到的模拟量转化成数字量。图像重建，即利用各种反演算法从所测数据重建出物体内部某种物理量的分布图。这一步是成像技术的核心。图像后处理，包括图像去噪、锐化，特征提取以及图像的存储、压缩和传输等。1.4 论文的组织结构论文主要研究辐射成像技术中的透射断层成像理论。讲述了CT重建的基本原理第一章绪论部分，首先介绍了计算机断层成像术（CT）的发展历史，及其医学意义根据射线产生的物理方式对CT具体的分类。之后详细介绍了透射断层成像的研究现状，进而对于目前投影重建中存在的问题进行了阐述，从而引出了论文的工作目标和内容。第二章详细介绍了CT图像重建系统

21、的基本原理，根据Beer定理对断层进行体素划分，通过联立方程组求解方程组的方式重建图像。但当断层被划分为更小更多的体素时，即使采用今天最先进的计算机技术解这一组方程组也不是件容易的事。所以这个方法在后来实际的图像重建中已不再使用，从而引出了现在CT图像重建中常用的方法解析法。第三章详细介绍了现在CT常用的图像重建算法解析法的算法原理，研究进第一章绪论展。先讲述了直接反投影重建算法，但由于直接反投影法存在图像伪迹，从而提出了滤波反投影重建算法。在运算的过程中根据卷积定理从而提出了卷积反投影重建算法，最后对上述三种解析算法重建的图像进行了对比分析。第四章迭代重建算法部分，上述解析算法是目前医用CT

22、中最常用的重建算法。但它并不适用于所有投影重建场合。在实际应用中，有时无法测得大量的投影，有时投影不是均匀的分布在180度或360度范围内。例如，为了避免物体（如脏器）运动，或为了减小剂量，投影数据采集不足；又如某物体在一个方向上尺寸特别长，则在这一方向上的投影数据无法或较难测得。再如用成像法勘探地球资源时，用孔穴法采集的数据，数量不足又不均匀等上述滤波反投影法无能为力，于是提出了图像重建的迭代重建算法。本章介绍了迭代算法的成像模型及其成像原理。最后演示了一下迭代算法重建图像的效果。第五章总结部分，本章对论文中讲述的CT重建算法进行了认真的总结和比较。探讨了尚需进一步解决的问题及今后有待努力的

23、方向。第二章CT图像重建的基本原理第二章CT图像重建的基本原理2.1 所需数据的扫描获取假设采用单射线源平移一旋转方式来获取图像重建所需的数据，该方法需要一个射线发生器和一个射线探测器。用X射线发生器发出单色细束X射线，对成像对象进行投射，如图2.1所示。设X线源发出的X射线强度是I。，探测器接收到的X射线强度是I。一次扫射完成后将X射线与探测器平移一个小的步长，再次对成像对象进行投射，探测器又接收到一个新的透射强度，在该投射角度上反复小步长平移投射并记录检测的数据后，将X射线源与探测器绕被测对象旋转一个小的角度后在进行平行投射并平移，如此重复多次，直到旋转18。度为止。取得这些不同投射位置的

24、透射强度后，就可以用这些数据，采用一定的算法重建出X射线源环绕的断层的图像。1M痔图2.1单射线源平移旋转扫描示意图2.2 图像重建算法CT与普通的X线摄影术之间有着非常重要的区别。在CT技术中，组织对X线的局部衰减特性被用于离散成像；而在通常的X线摄影术中，这种衰减信息则重叠在X线底片上。组织对X线的这种局部衰减特性是X线与被检测物体之间的若干相互作用过程的产物。例如光电吸收过程和康普顿散射过程。这些过程中的每一种都有它自己的发生几率。几率也是辐射能量的函数，因为从X线管产生的X线由全部能谱组成，具有“线衰减系数u”的某种组织的衰减性质是一个复杂第二章CT图像重建的基本原理的函数，随着辐射情

25、况的变化，它可以表现出不同的数值。对于能量为E的单能射线，假设其投射的物体是均匀且为同一介质，被投射物体的线性衰减系数为u,厚度为x,则在入射线强度为I。，出射X射线强度为I时有II0eux即uxln(I0/I)(2.1)这就是著名的Lambert-Beer定律。实际被投射的物体是非均匀的，如果在X射线束通过的路径l上，介质是非均匀的，则可以将X射线穿过的介质沿扫描路径l人为划分为大小相等的n个小方块，每个小块的厚度为d,块内视为均匀介质。这样，每个小块有一个均匀的线性衰减系数，如图所示，设各个方块的线性衰减系数分别为Ui,U2,U3,.,Un。X射线通过每一个方块的衰减为IiIeu1d，射线

26、通过第二个方块的衰减为I2Iieu2d，通过第n个方块的衰减为IInIneund。将上述各式代入整理，有：图2.2X射线通过非均匀介质时的体素划分(2.2)uiu2u3unp/d(2.3)(uiu2u3.un)dIIoe式中，pln(J)是通过射线强度检测得到的检测值，是一个已知量，称为I投影。这样，上式中，等号右边可视为一个常数，等待求解的是等式左边的每一个方块的衰减系数。由于未知量太多，这个方程没有唯一解。考虑到围绕受检体的一次扫描可以得到一个有关衰减系数的方程，如果将整个感兴趣断层都划分为大小相等的小方块，每个小方块由于一个固定的衰减系数，这样，上述单束射线环绕扫描得到的每个穿透强度值都

27、可以形成一个关于衰减系数的方程。只要扫描线不重复且足够多，这些衰减系数的确定值肯定可以通过解线性方程组的方法解第二章CT图像重建的基本原理出唯一实数解。把每个方块的衰减系数值用灰度表示,就可以重建出以衰减系数为特征的断层图像。这种图像重建方法称为方程组法。值得说明的是，当上图中d0时，投影p就是(2.4)pln(，)u(l)dl式中，衰减系数u(l)是随路径l连续变化的函数。上式说明，入射X射线强度与出射X射线强度之比经对数运算后，表示沿射线路径上衰减系数的线积分,而投影p与射线穿越介质的路径长度成正比。然而，一般二维断层图像至少应划分出160X160个像素，如要用方程组法重建这个断层图像，就

28、要25600个独立的方程联立求解，这是一个相当费时的任务。尽管1967年世界首台CT试验机采用的就是这个图像重建方法，但当断层被划为更小更多的体素时，即使采用今天最先进的计算机技术，解一组这样的方程组也不是件容易的事。所以，这种方法在后来的实际应用中已不再使用，而是采用下一章将要介绍的反投影解析法。第三章CT图像解析重建算法第三章CT图像解析重建算法通过上一章已经知道CT图像重建的基本原理，但由于使用最简单直接的方程组法不易求得断层的各个像素的衰减系数。这一章介绍现在CT最常用的图像重建算法一一反投影解析法。主要介绍了由投影重建图像的系统模型，引进Radon变换和傅立叶变换，它们对应于X射线C

29、T的二维重建问题；且给出了Radon变换在图像重建中的具体形式，即截面函数沿着直线的线积分等于其Radon变换，而图像重建过程即是将截面函数沿许多不同角度下直线的线积分所产生的投影数据进行逆Radon变换从而得到截面函数。另外由逆Radon变换的推导给出了两种常用的图像解析重建法：滤波反投影法和卷积反投影法。3.1 反投影重建法（B?反投影是一个求和的过程，它把所有角度上的数据都加在一起。它把投影域中的数据沿着原来的投影路径“涂抹”回去，但不改变数据的值网。直接反投影重建是指，断层平面中某像素的线性衰减系数可由平面内所有经过该像素射线的投影和得出，而整幅图像重建可以看做为所有方向上扫描投影的累

30、加。通俗地讲，所谓直接反投影，就是在对某体面一个方向的扫描完成后，以得到的投影为灰度值沿着扫描路径经过的体素回抹到体素对应的像素上。改变方向后的多次扫描形成多次回抹，同一像素上多次回抹的灰度累加即完成图像重建。图3.1断层的像素和射线模型第三章CT图像解析重建算法图3.2反投影重建示例如图3.2所示，其中(a)为原图像像素，(b)为反投影重建后图像，(c)为(b)中像素值除以投影线数后得到的平均值。有上述例子可以看出原图中像素值为零的像素，经反投影重建后不再为零。即有伪迹。直接反投影图像产生星状伪迹的原因在于，该方法是把取自有限物体空间的投影均匀地回抹(反投影)到了射线所及的无限空间的各个像素

31、上，包括原来像素值为0的点。口(工扫描后取投崽反投意图像重建图3.3反投影重建的等效系统设u(x,y)为处于x0,y0处的一个点源，即二维断面(x,y)上的一个冲击函数(x,y),这时的输出f(x,y)就是系统的冲击响应h(x,y)。1已知h(x,y)&飞(3.1)可见，反投影重建过程后得到的图像不是原来的图像(x,y)，系统存在星状伪迹。根据信号与系统理论，对于图3.3有：f(x,y)u(x,y)*h(x,y)(3.2)式中，两个星号表示二维卷积410第三章CT图像解析重建算法如果想要去除星状伪迹，可在系统的输出端再加上一个滤波器，设其时域函数为q(x,y)，传递函数为Q(1，2),为了使该

32、加了滤波器的系统输出的图像等于系统输入的原图u(x,y),也就是要求q(x,y)*h(x,y)(x,y)(3.4)即1q(x,y)*；22(x，y)(3.5).xy对该式取二维傅立叶变换，得到：1,Q(1,2):221(3.6)xy即q(1.2)T77(3.7)这是一个二维滤波器。即使加上某种近似，这种滤波器也不容易实现。如果1和2的值取无穷，则是一个不可能实现的理想滤波器。为此提出了滤波反投影图像重建法。3.2 滤波反投影图像重建法(FBP为了消除直接反投影法产生的星状伪迹，提出了滤波反投影重建的方法，这种方法是在直接反投影重建方法的基础上增加一个滤波器。增加滤波器的思路有两种，如图3.4所

33、示。第三章CT图像解析重建算法物体新射t原始图像先用护影:后：维涯液输出图像（同质图像）划体断层原姆图像）输出图像同悚图像）1皿先对投影函数一摊滤波再反投影帝建图3.4两种不同的去除星状伪迹反投影图像重建思路所谓滤波反投影重建法是把获得的投影函数先做卷积处理，即人为设计一种一维滤波器，用它对所得的投影函数进行改造，然后把这些改造过的投影函数进行反投影和累加等处理，就可以达到在一定程度上消除星状伪迹的目的。在介绍滤波反投影重建法之前先补充滤波反投影重建法要用到的数学知识和投影定理一一中心切片定理。3.3 Radon变换与傅立叶变换为了有效和快速地对图像进行处理，常需要将原定义在图像空间的图像以某

34、种形式转换到另外一些空间，并利用在这些空间的特有性质方便地进行一定的加工，最后再转换回图像空间以得到所需的效果。通常正变换：图像空间到其他空间反变换：其他空间到图像空间3.3.1 Radon变换Radon变换1917年由奥地利数学家Radon提出，是CT图像重建思想的数学表达。Radon从理论上证明了某物理量的二维分布函数可由该函数在其定义域内的所有方向上的线积分完全确定。Radon变换的意义在于，只要知道了一个未知二维分布函数的所有方向上的线积分，那么就能够求得该二维分布函数。所谓CT断层的图像重建，就是求取能够反应断层内部结构和组成的某物理量的二维分布。断层扫描获得的投影实际上可视为被测物

35、体断层内部结构和组成的不同方向上的线积分，所以Radon变换的正变换和逆变换正好对应了CT技术的射线投影和图像二维分布函数的重建过程。12第三章CT图像解析重建算法对f(x,y)的Radon变换Rf(p,)定义为沿由p和定义的直线l的线积分。若设f(x,y)表示物体的一个二维断层分布，通过Radon变换的到f沿不同方向直线的线积分相当于得到物体不同方向的投影。图3.5定义Radon变换的坐标系统上述线积分可写为:Rf(p,)f(x,y)dl(3.8)如果借助Delta函数，上述线积分还可写为：Rf(p,)f(x,y)(pxcosysin)dxdy(3.9)上式通常称为f(x,y)的Radon变

36、换。由上述可知，f(x,y)关于某直线的Radon变换就是f(x,y)沿该直线上的一维投影。实际上当将也视作可变参数时，Rf(p,)就是一个二元函数，此时Radon变换是空间域(x,y)到变换域(p,)的映射。域(p,)中的每个点对应于空间域(x,y)中的一条直线，这里(p,)实际上就是图像分析中直线的Hough变换10。3.3.2傅立叶变换的定义一1一：(3.10)f2-f()eideid其中13第三章CT图像解析重建算法(3.11)F()f(t)eitdt(实自变量的复值函数)称为f的Fourier变换，记为Ff(t)0；F()eitd称为F()的Fourier逆变换，记为F1f()。3.

37、4中心切片定理中心切片定理的含义是平行投影的一维傅立叶变换等同于原始物体的二维傅立叶变换的一个切片。将已知投影数据通过一个简单的二维傅立叶反变换可以得到物体截面的一个评估。通过上图可知：某断层(或它对应的图像)f(x,y)在视角为时得到的平行投影(函数)的一维傅立叶变换，等于f(x,y)二维傅立叶变换F(1，2)过原点的一个垂直切片，且切片与轴1相交成角。数学表达如下：首先定义代表截面的函数的二维傅里叶变换：(3.12)F(u,v)f(x,y)ej2(uxvy)dxdy类似的定义某个角度下的一条投影P(t)的傅里叶变换:第三章CT图像解析重建算法S()P(t)ej2tdt(3.13)举个最简单

38、的例子：沿着=0的直线方向，物体投影的傅里叶变换等于频域里面=0的情形F(u,0)f(x,y)ej2uxdxdy(3.(13)因为相位因子不再依赖于y,在此可以将积分分成两部分，F(u,0)f(x,y)dyej2uxdx(3.(14)从平行投影的定义看，可以将上式括号中的部分看成是沿着常量x的积分P0(x)f(x,y)dy(3.(15) 将其代入(3.14)得到F(u,0)P0(x)ej2uxdx(3.(16)等式右边表示的是投影P0的一维傅里叶变换；因此可以得到垂直方向投影和投影函数的二维变换的关系式：F(uQ)S0(u)(3.17)这是傅里叶切片定理最简单的形式。很明显物体与坐标系统间的位

39、置方向是独立的。例如，在图3.6中，如果(t,s),坐标系统为(x,y),坐标系统旋转，式(3.16)定义的投影的傅里叶变换等于物体沿着一旋转了角的直线的二维傅里叶变换。从而引出了成像系统中著名的中心(傅里叶)切片定理20即某图像f(x,y)在视角时的平行投影的傅里叶变换给出了图像二维傅里叶变换F(u,v)在与u轴夹角为的一个切片，此切片通过原点。换句话说，的P(t)傅里叶变换在数值上对应于F(u,v)沿着图中所示直线BB的值。通过将(x,y)坐标系统旋转得到(t,s)坐标系统后，傅里叶切片定理可以得到更广泛的应用。在(t,s)系统中沿着常量为t的直线进行的投影可表示成：P(t)f(t,s)d

40、s(3.18)第三章CT图像解析重建算法其傅里叶变换S()P(t)ej2tdt(3.19)将以上两式合并，可以发现S()f(t,s)dsej2tdt(3.20)借助旋转坐标定义可以得到关于(x,y)坐标系统的变换形式S()f(x,y)ej2(xcos()ysin(dxdy(3.21)等式右边表示一个空间域的傅里叶变换，ucos(),vsin()则S()F(,)F(cos(),sin()(3.22)这个等式是直接射线断层成像的基础，从而也证明了中心(傅里叶)切片定理。上面的结果显示：从取图像f(x,y)分别在与x轴成1,2,.,k角度的直线上的一系列投影，并对这些投影逐个进行傅里叶变换，可以得到

41、图像f(x,y)的二维傅里叶变换F(u,v)在与u轴成1,2,.,k角的直线上的值。理论上，当投影线束无限多时，就能求得F(u,v)在频域上相应线束的无限多条直线的值，从而可以利用傅里叶逆变换得到原图f(x,y)。傅里叶逆变换可以采用直角坐标系公式也可以采用极坐标公式：f(x,y)F(u,v)ej2(uxvy)dudv(3.23)f(x,y)F(,)ej2I|dd0(3.24)其中(3.25)txcosysin但从多个一维傅立叶函数拟合二维傅立叶函数需要足够精确的傅立叶空间插值。由于理想的插值方法还在探索中，这阻碍了基于中心切片定理的二维傅立叶逆变换图像重建法的发展。第三章CT图像解析重建算法

42、3.5 平行投影重建算法3.5.1 算法原理傅里叶切片定理仅仅提供了一种断层成像的简单概念模型，实际应用中需要一种不同的算法。滤波反投影算法是目前广泛应用于所有直线透射断层成像的算法。此算法源于中心切片定理，不仅极其准确而且保证了快速执行的稳定性。它是通过转换极坐标中傅里叶逆变换和重新确定积分限来实现的。首先介绍平行束投影的反投影算法。先回顾一下傅里叶反变换的公式：f(x,y)j2(uxvy)F(u,v)edxdy(3.26)将频域内的直角坐标系(u,v)换成极坐标系(，)，有如下替换:ucos(),vsin()(3.27)(3.28)(3.29)此积分可以分成两部分,从0180和180360

43、改变微分项dudvdd则傅里叶反变换的极坐标形式:2j2(xcosysin)f(x,y)F(,)edd002f(x,y)F(,)ej2(xc0s00ysin)ddF()e，2(xcos()ysin()dd0(3.30)利用特性F(,)F(17第三章CT图像解析重建算法f(x,y)F(,)ej2I|dd00(3.31)此处txcosysin。利用(3.21)式用视角处的投影的傅里叶变换S()代替二维傅里叶变换F(,),可以得到:f(x,y)S()ej2=|dd00(3.32)上式是滤波反投影重建的图像表达式。对上式进行不同的数学变形，将导致不同的物理解释，相应于不同的重建算法。由此得到的重建算法

44、有滤波(卷积)反投影法(F(C)BP)和Radon反变换重建法。3.5.2滤波(卷积)反投影有卷积定理：函数卷积的傅立叶变换是函数傅立叶变换的乘积，即：一个域中的卷积相当于另一个域中的乘积。例如：时域中的卷积等于频域中的乘积。设在某一旋转角时，采的投影p(t,),滤波函数为h(t),则滤波后的投影为(3.33)由于数据在计算机实现时是离散的，故应进行离散卷积。反投影是滤波反投影算法的另一个核心，其实就是前面介绍的直接反投影算法，只不过是在反投影前将投影数据进行了卷积滤波预处理卷积反投影重建图像如下：由此可得极坐标下f(r,)g(t,)|trcos()d0(3.34)3.5.3Radon反变换重

45、建算法在CT成像原理中已经介绍过，对图像某一方向的投影表示沿射线路径方向上衰减系数的积分，而对图像进行某一视角上的Radon变换求得的的就是图像该第三章CT图像解析重建算法方向上的投影，所以对投影进行Radon反变换就能实现图像重建。对滤波后的图像投影数据进行反Radon变换得到重建图像的方法就是Radon反变换重建法。3.6 解析重建算法的实现和对比分析针对本章介绍的几种重建算法，文本用matlab进行了模拟重建，下面采用的是调用系统函数phantom生成的标准SheepLogon头模型。如图3.7。K-edtInsertToolsDesktopMidcwH守p口二iljI片晨，。4X聂|口

46、门dJL.b,2；J怛拉回博图3.7SheepLogon头模型用iradon函数用不同滤波器进行滤波反投影重建的仿真图像如图3.8所示:图3.8对标准头模型解析重建的图像对比第三章CT图像解析重建算法下面是对几组图像处理中常用的图片直接读入再进行解析法模拟重建的效果图，重建过程中采用了各种不同的滤波器。间厚始图像。直接反投影重建依JR白讣匕归漉渡反投身重建图像(d)Shepp-ogan滤波反投劈重建国候原始图像Hann滤波反投量重建(f)H白mming漉港反投寄重建图像小/向他滤波反投登重建图像图3.9采用各种不同的滤波器进行滤波反投影重建效果对比20第三章CT图像解析重建算法在图像重建中，一

47、般用视觉相似度SSIM8对图像重建的效果进行评估。视觉相似度越接近1说明重建图像与原图越相进。采用不同的滤波器进行滤波反投影重建的图像与原始图像相比的视觉相似度及重建所用的时间如表3.1所示。表3.1采用不同滤波器滤波反投影重建图像与原图的视觉相似度对比滤波器名称SSIM时间（秒）直接反投影0.0328460.401583Ram-Lak0.9999650.450677Shepp-Logan0.9999660.337064Hann0.9999680.349063Hamming0.9999680.286156Cosine0.9999670.293614通过图3.9各重建图像效果的对比和表3.1中视

48、觉相似度的对比可知，直接反投影重建由于没有滤波，存在伪迹，图像不清晰，与原图的视觉相似度很小。采用了不同的滤波器，相对于没有滤波的情况，视觉相似度都很接近1,说明与原图比较接近，图像质量有了明显的提高。以上介绍的是平行束投影重建算法。CT扫描方式有平行束和扇束扫描两种扫描方式。扇束投影重建算法大致可分为两类：一类是重排算法，即把一个视图中采集到的扇形数据重新组合成平行的射线投影数据，再用上面所介绍的平行束算法重建。这种方法对扇束投影的视角与每一扇束投影各射线间的增角有一定的约束条件。另一种算法是扇束投影直接重建算法，这种算法不必先把数据重排，而是根据扇束投影数据自身的特点，发展出直接用扇形束投

49、影数据重建图像的方法。扇束扫描自身的算法有两种，区别在于抽样的方式上，一种是等角抽样，即射线是按等角分布的，将检测器等角分布在以放射源为中心的圆弧或直线上；另一种是等间距抽样，检测器单元在直线上作等间距据分布。这里不作详细介绍。21第四章迭代重建算法第四章迭代重建算法4.1 迭代重建算法的思想原理迭代重建算法的概念与前面讲述的解析法最大的区别在于前者一开始就把连续的图像离散化。在解析算法中，假设图像是连续的，每个像素只是一个点，这些离散的点是以图像显示为目的的，这些点的选择可以是随意的，与图像重建无关。但在迭代算法中，每个像素是个小面积。这些小面积在计算当前图像的投影数据时要用到。像素模型对重

50、建图像的质量好坏影响很大。图像重建的解析算法与图像重建的迭代算法之间的另一个区别是解析算法着力对一个积分方程求解，而迭代算法着力对一个线性方程组求解。迭代的主导思想是，假设断层截面由一个未知的矩阵组成，然后由测量投影数据建立一组未知向量的代数方程式，通过方程组求解未知图像向量。此方法虽然比变换重建法易于理解，但对于图像要求较高的应用来说，这一方法缺少精度和速度。然而在不可能获取大量投影数据或者投影不均匀分布于180或360之间情况下，迭代法是一个非常有效的方法。4.1.1 迭代重建的数学模型!就S海为凭空Sssslr万方今隹r一n黛缪滕0IIRZ勿钦经图4.1断层截面与投影的几何布置如图4.1

51、所示，图中将截面离散化，即在图像f(x,y)上叠加一个方格网。每一个小方格内像素值f(x,y)是一常量，令Xj表示第j个单元的常量值，N表示像素的总数。对于迭代法而言射线的定义稍有不同。具体地讲，此时穿过(x,y)平面的射线具有一定的宽度，为了说明这一点在图中用阴影加重了第i根射线，22第四章迭代重建算法且每一根都是如此，绝大多数情况下射线宽度大约等于图像像素单元的宽度。此时线积分被称为射线和。如图4.1所示，投影数据也采用单索引表示，设y是图中由第i根射线测量的射线和。则Xj和y的关系可以表示如下，NajXjVii1,2,.,M(4.1)ji式中M为射线总数(或投影总数)，而a。为加权因子，

52、其代表第j个像素对第i根射线积分的贡献。其值等于第i根射线穿过第j个像素所占部分的面积，如图中所示。需要注意，大多数a。均为0,对于任意一条射线和，只有很少的像素单元有贡献。如果M和N值很小的话，可以用传统的矩阵理论方法对式(4.1)进行求逆。实际中对于128X128的图像阵，N=16384,如果要重建256X256的图像，N值将更大。同时M往往具有与N同样量级的数。假设M和N等于16384的情况，式(4.1)中的矩阵aj的大小等于16384X16384,而这将使任何直接求解矩阵的逆变得无法实现。当然,当测量数据混有噪声或MN时，即使N值很小，直接矩阵求逆也是不允许的，这种情况可能不得不使用最

53、小二乘法来求解。对于M和N为很大值的情况，有一些很吸引人的迭代方法来求解式(4.1)。这些方法是基于“投影法”的，最早被Kaczmarz提出14,此后由Tanabe15进一步阐明。为了说明这些方法中的计算步骤，可以把(4.1)式写成如下形式，(a11x1a12X2.a1NXNY1a21X1a22X2.a2NXNy2iaM1X1aM2X2.aMNXNyM,、(4.(2)令X(X1,X2,.,Xn)T。因此，由X表示的一幅图像可以看成是N维空间中的一个点。在这一空间中，上面方程组中每一个方程式表示一个超平面。当这些方程组的唯一解存在的话，这些超平面的交点为一个点，这就是(4.2)式方程组的解。这一

54、概念进一步用图4.2说明，为理解方便，这里仅考虑有两个像素xX2的情况，它们满足下面的方程，23第四章迭代重建算法七为02X2V1(4.(3)a21x1a2&y2图4.2Kaczmarz交替投影法求解线性方程组的示意图找寻图4.2中解的计算过程如下：先从一个初始图像X0开始，把这一初始向量投影到第一根射线，把投影点再投到第二根射线上，然后将投影点再投回到第一根射线上，依此类推。若唯一解存在，迭代终究会收敛到那一点。对于这一方法的计算机实现而言，首先在可行解中估计一个初始解。这一估计解是N维空间的一个点，以X0(x0，x；,.,xN)表示。一般情况下可以设置X0(0,0,.0)。将这一初始点投影

55、到式(4.2)中第一个等式所表示的超平面得到Xii1AXv(4.4)X1,图4.2所示为两个变量的情况。再将X1投影到式(4.3)中第二个等式所表示的超平面得到X2,依次下去。当Xk1投影到(4.2)式中的第i个等式所表示的超平面时，得到Xk的过程可以用下面的数学方式表示式中A(七82，.为)，而aA为A与其自身的内积。如前述，迭代重建的计算过程如下：给定一个初始解，将其投影到方程(4.2)表示的某个超平面上，得到另一个解，再将其投向下一个超平面上，依此循环，当所有方程所表示的超平面都被投影过后，得到XM。在下一次迭代中，XM作为初始解，再依次投影到(4.2)式中表示的各个超平面上，获得X2M,如此重复。Tanabe15证