毕业论文-基于matlab实现的ct重建算法仿真比较研究

密级公开NANCHANGUNIVERSITY学士学位论文THESISOFBACHELOR（20092013年）题目基于MATLAB实现的CT重建算法仿真比较研究学院信息工程系电子信息工程专业班级学生姓名学号指导教师职称起讫日期2013318201367基于MATLAB实现的CT重建算法仿真比较研究摘要计算机断层成像（即CT（COMPUTERIZEDTOMOGRAPHY）技术是一项近年来快速发展多学科交叉的先进技术，除了医学应用，还应用在工业无损检测、射电天文学、精密仪器反演等多个重要领域。CT图像重建是CT技术的核心，重建算法的优劣直接关系到对检测结果判断的准确性。医学CT普遍采用解析成像技术，这主要因为医用CT可以采集到完整的投影数据。然而在工业应用中，一般得不到完整的投影数据，因此要利用不完整的投影数据得到较清晰的图像，必须采用迭代重建法。针对以上情况，本课题作了以下两方面的分析比较。1能获取完整投影数据时，用计算机模拟了解析法的图像重建。发现图像重建的速度快，且重建的图像与原图的视觉相似度接近1。2不能获取完整投影数据时，用计算机模拟了迭代算法的图像重建。先用ART迭代算法引入松弛参数，分析了松弛参数重建图像的影响。接着用SIRT算法重建图像，这种算法所得的图像比ART算法较好。最后在ART算法基础上采用全差分对重建图像进行优化，图像质量有了明显提高。关键字图像重建算法；解析重建；迭代重建BASEDONMATLABREALIZETHECTRECONSTRUCTIONALGORITHMSIMULATIONCOMPARATIVESTUDYABSTRACTCOMPUTEDTOMOGRAPHYCTCOMPUTERIZEDTOMOGRAPHYTECHNOLOGYISARAPIDDEVELOPMENTINRECENTYEARS,MULTIDISCIPLINARYADVANCEDTECHNOLOGY,INADDITIONTOMEDICALAPPLICATIONS,BUTALSOUSEDININDUSTRIALNONDESTRUCTIVETESTING,RADIOASTRONOMY,PRECISIONINSTRUMENTSANDOTHERIMPORTANTAREASOFTHEINVERSIONCTIMAGERECONSTRUCTIONCTTECHNOLOGYISTHECOREOFTHERECONSTRUCTIONALGORITHMISDIRECTLYRELATEDTOTHEACCURACYOFTHETESTRESULTSDETERMINEMEDICALCTIMAGINGCOMMONLYUSEDANALYTICALTECHNIQUE,WHICHCANBECOLLECTEDMAINLYDUETOACOMPLETEMEDICALCTPROJECTIONDATAHOWEVER,ININDUSTRIALAPPLICATIONS,THEGENERALLACKOFCOMPLETEPROJECTIONDATA,ANDTHEREFORETHEUSEOFINCOMPLETEPROJECTIONDATATOBEMORECLEARIMAGE,ANITERATIVERECONSTRUCTIONMETHODMUSTBEUSEDFORTHEABOVE,THISPROJECTMADETHEFOLLOWINGTWOASPECTSOFTHEANALYSISANDCOMPARISON1CANGETTHEFULLPROJECTIONDATA,USINGACOMPUTERSIMULATIONOFTHEANALYTICALMETHODOFIMAGERECONSTRUCTIONFOUNDFASTIMAGERECONSTRUCTIONANDTHERECONSTRUCTEDIMAGEWITHTHEORIGINALVISUALSIMILARITYISCLOSETO12CANNOTGETTHEFULLPROJECTIONDATA,USINGACOMPUTERSIMULATIONOFTHEITERATIVEALGORITHMFORIMAGERECONSTRUCTIONFIRSTWITHTHEINTRODUCTIONOFARTRELAXATIONPARAMETERITERATIVEALGORITHMANALYZESTHEIMPACTOFRELAXATIONPARAMETERSOFTHERECONSTRUCTEDIMAGEFOLLOWEDBYSIRTALGORITHMFORRECONSTRUCTIONOFIMAGES,THERESULTINGIMAGEOFTHISALGORITHMISBETTERTHANTHEARTALGORITHMFINALLYARTALGORITHMBASEDONAFULLYDIFFERENTIALOPTIMIZEDFORIMAGERECONSTRUCTION,IMAGEQUALITYHASIMPROVEDSIGNIFICANTLYKEYWORDSIMAGERECONSTRUCTIONALGORITHM；ANALYTICRECONSTRUCTION；ITERATIVERECONSTRUCTION目录摘要IABSTRACTII第一章绪论111引言112CT的分类213CT的组成结构414论文的组织结构4第二章CT图像重建的基本原理621所需数据的扫描获取622图像重建算法6第三章CT图像解析重建算法931反投影重建法（BP）932滤波反投影图像重建法（FBP）1133RADON变换与傅立叶变换12331RADON变换12332傅立叶变换的定义1334中心切片定理1435平行投影重建算法17351算法原理17352滤波（卷积）反投影18353RADON反变换重建算法1836解析重建算法的实现和对比分析19第四章迭代重建算法2241迭代重建算法的思想原理22411迭代重建的数学模型2242代数重建法（ART）2543联立迭代重建法（SIRT）2644全差分迭代算法（TV）2745迭代重建算法的实现和对比分析27第五章结论37参考文献38附录40谢辞60第一章绪论11引言自从X射线发现后，医学上就开始用它来探测人体疾病。但是，由于人体内有些器官对X线的吸收差别极小，因此X射线对那些前后重叠的组织的病变就难以发现。于是，美国与英国的科学家开始了寻找一种新的东西来弥补用X线技术检查人体病变的不足。CT电子计算机体层成像是70年代初放射诊断的项重突破，CT不是X线摄影，而是用X线对人体扫描（图11所示），取得信息,经电子计算机处理而获得的重建图像它能使传统的X线检查难以显示的器官及其病变显示成像，而且图像逼真，解剖关系明确，从而扩大了人体的检查范围，大大提高了病变的早期检出率和诊断准确率。这种检查简便、安全、无痛苦、无创伤、无危险，它促进了医学影像诊断学的发展，CT的研制成功被誉为自伦琴发现X射线以后，放射诊断学上最重要的成就。发明者亨斯菲尔德和科马克共同获得了1979年的诺贝尔奖金。CT最初只用于头部检查，1974年又出现了全身的CT。在短短10余年间，CT已遍及全球，从第一代发展到第五代。我国各大城市医院所使用的CT属第三代。全身CT可以作头、胸、腹、骨盆的横断扫描，5也可作甲状腺、脊柱、关节和软组织及五官等小部位的区域扫描。CT最适于查明占位性病变如肿瘤、囊肿、增大的淋巴结、血肿、脓肿和肉芽肿的大小、形态、目和侵犯范围，它可以决定某些器官癌肿的分期和是否能进行手术切除。在某些情况下，CT还能区别病变的病理特性如实性、囊性、血管性、炎性、钙性、脂肪等。CT成像技术的基本原理是根据人体中不同组织对X射线的衰减作用不同，用X线束对人体某部一定厚度的层面进行扫描，由探测器接收透过该层面的X线，转变为可见光后，由光电转换变为电信号，再经模拟/数字转换器（ANALOG/DIGITALCONVERTER）转为数字，输入计算机处理。图像形成的处理有如对选定层面分成若干个体积相同的长方体，称之为体素（VOXEL）。扫描所得信1息经计算而获得每个体素的X线衰减系数或吸收系数，再排列成矩阵，即数字矩阵（DIGITALMATRIX），数字矩阵可存贮于磁盘或光盘中。经数字/模拟转换器（DIGITAL/ANALOGCONVERTER）把数字矩阵中的每个数字转为由黑到白不等灰度的小方块，即像素（PIXEL），并按矩阵排列，即构成CT图像。所以，CT图像是重建图像。每个体素的X线吸收系数可以通过不同的数学方法算出。如图11所示。图11计算机断层成像示意图12CT的分类在投影重建的原理和方法的框架下，根据投影方式不同产生了不同种类的CT成像术。透射断层成像透射断层成像（TRANSMISSIONCOMPUTEDTOMOGRAPHY,TCT，简称CT）系统中，从发射源射出的X射线穿透物体到达接收器。射线在通过物体时被物体吸收一部分，余下部分被接收器接收。由于物体各部分对射线的吸收不同，所以接收器获得的射线强度实际上反映了物体各部分对射线的吸收情况。由于透射成像术中采用X射线，所以也叫XCT。接收器接收到的模拟信号经模数转换器转换成数字信号后，把代表着不同角度下的投影数据送给计算机，它再运用复杂的数学方法重建出物体截面的图像。在医学上它主要针对骨质成像。4XCT成像原理如图12所示图12X射线CT成像原理发射断层成像在发射断层成像（EMISSIONCOMPUTEDTOMOGRAPHY,ECT）系统中，发射源在物体内部，一般是将具有放射性的离子注入物体内部，从物体外检测放射出来的量。通过这种方法可以了解离子在物体内运动情况的分布，从而可以检测到与生理相关的状况。现在通常用的ECT主要有两种正电子成像PET（POSITRONEMISSIONTOMOGRAPHY）；单光子成像SPECT（SINGLEPOSITRONEMISSIONCOMPUTEDTOMOGRAPHY）。它们的区别是，采用单光子辐射器时，检测器所接收的光子计数便是投影测量数据；而采用正电子辐射器时，由于在正电子辐射点的毫米级范围内，一个正电子会俘获一个负电子，并产生沿相反方向运动的两个中子，因而测量数据是在一个很小时间间隔内由两个性质相反的检测器检测到的光子对的数目。ECT主要针对代谢过程的成像。磁共振成像磁共振成像（MAGNETICRESONANCEIMAGING,MRI）在早期称为核磁共振（NUCLEARMAGNETICRESONANCE,NMR）。它的发明者斯坦福大学的FBLOCH和哈佛大学的EMPURCELL为此获得了1952年度的诺贝尔物理奖。它的工作原理简介如下氢核以及其它具有奇数个质子和中子的核，包含具有一定磁动量和旋量的质子。如果把它们放在磁场中，它们会像陀螺在地球重力场中一样在磁场中进动。一般情况下质子在磁场中时随意排列着的，当一个适当强度和频率的共振场信号作用于物体时，质子吸收能量并转向与磁场相交的朝向。如果此时将共振场信号除去，质子吸收的能量释放并被接收器检测到。根据检测到的信号就可以确定质子的密度。通过控制所用的共振场信号和磁场强度，可每次检测到沿透过目标中一条线的信号。换句话说，检测到的信号是MRI信号沿直线的积分。所以检测目标的工作成了投影重建的问题。由这些信号可以得到极好的器官解剖图，在细节上远远超过由X射线CT产生的图像。由于它能够很好的保留丰富的软组织图像，所以在脑功能研究上，MRI胜过ECT。超声断层成像超声断层成像（ULTRASOUNDCT,UCT）与上述的CT不同的地方是超声射线在人体组织中不一定按直线传播，尤其是在软组织和硬组织的分界面会产生折射效应。因此超声成像主要用于人体中没有硬组织（如骨骼）部分的检测，特别是妇女乳房癌的检测。反射断层成像反射断层成像（REFLECTIONCT,RCT）也是利用投影重建的原理工作的。常见的一个例子是雷达系统，其中的雷达图是物体反射的回波所产生的。例如在前视雷达（FLR）中，雷达发射器从空中向地面发射无线电波；雷达接收器在特定角度所接收到的回波强度是地面发射量的一个扫描段的积分。通过发展硬件设备和改进软件算法，可以不断降低噪声和提高图像的清晰度。CT技术在临床的许多领域得到愈来愈广泛的应用。应该指出，上面介绍的各类CT的功能是相互补充而不是相互替代的。613CT的组成结构以上几种成像技术尽管使用的物理系统各不相同，但所有系统所测量到的数据均具有物体某种感兴趣的物理特性在空间分布的积分形式。成像技术广义上讲大致可分为三个部分信息的收集、重建以及处理和传递。数据采集，即获得携带有物体内部特征信息的数据。所采用的物理测量系统包括放射源、检测器及配套的机械、电路设备。为进行计算机处理，系统还需将测量得到的模拟量转化成数字量。图像重建，即利用各种反演算法从所测数据重建出物体内部某种物理量的分布图。这一步是成像技术的核心。图像后处理，包括图像去噪、锐化，特征提取以及图像的存储、压缩和传输等。14论文的组织结构论文主要研究辐射成像技术中的透射断层成像理论。讲述了CT重建的基本原理第一章绪论部分，首先介绍了计算机断层成像术（CT）的发展历史，及其医学意义根据射线产生的物理方式对CT具体的分类。之后详细介绍了透射断层成像的研究现状，进而对于目前投影重建中存在的问题进行了阐述，从而引出了论文的工作目标和内容。第二章详细介绍了CT图像重建系统的基本原理，根据BEER定理对断层进行体素划分，通过联立方程组求解方程组的方式重建图像。但当断层被划分为更小更多的体素时，即使采用今天最先进的计算机技术解这一组方程组也不是件容易的事。所以这个方法在后来实际的图像重建中已不再使用，从而引出了现在CT图像重建中常用的方法解析法。第三章详细介绍了现在CT常用的图像重建算法解析法的算法原理，研究进展。先讲述了直接反投影重建算法，但由于直接反投影法存在图像伪迹，从而提出了滤波反投影重建算法。在运算的过程中根据卷积定理从而提出了卷积反投影重建算法，最后对上述三种解析算法重建的图像进行了对比分析。第四章迭代重建算法部分，上述解析算法是目前医用CT中最常用的重建算法。但它并不适用于所有投影重建场合。在实际应用中，有时无法测得大量的投影，有时投影不是均匀的分布在180度或360度范围内。例如，为了避免物体（如脏器）运动，或为了减小剂量，投影数据采集不足；又如某物体在一个方向上尺寸特别长，则在这一方向上的投影数据无法或较难测得。再如用成像法勘探地球资源时，用孔穴法采集的数据，数量不足又不均匀等上述滤波反投影法无能为力，于是提出了图像重建的迭代重建算法。本章介绍了迭代算法的成像模型及其成像原理。最后演示了一下迭代算法重建图像的效果。第五章总结部分，本章对论文中讲述的CT重建算法进行了认真的总结和比较。探讨了尚需进一步解决的问题及今后有待努力的方向。第二章CT图像重建的基本原理21所需数据的扫描获取假设采用单射线源平移旋转方式来获取图像重建所需的数据，该方法需要一个射线发生器和一个射线探测器。用X射线发生器发出单色细束X射线，对成像对象进行投射，如图21所示。设X线源发出的X射线强度是，探测0I器接收到的X射线强度是一次扫射完成后将X射线与探测器平移一个小的步0I长，再次对成像对象进行投射，探测器又接收到一个新的透射强度，在该投射角度上反复小步长平移投射并记录检测的数据后，将X射线源与探测器绕被测对象旋转一个小的角度后在进行平行投射并平移，如此重复多次，直到旋转180度为止。取得这些不同投射位置的透射强度后，就可以用这些数据，采用一定的算法重建出X射线源环绕的断层的图像。图21单射线源平移旋转扫描示意图22图像重建算法CT与普通的X线摄影术之间有着非常重要的区别。在CT技术中，组织对X线的局部衰减特性被用于离散成像；而在通常的X线摄影术中，这种衰减信息则重叠在X线底片上。组织对X线的这种局部衰减特性是X线与被检测物体之间的若干相互作用过程的产物。例如光电吸收过程和康普顿散射过程。这些过程中的每一种都有它自己的发生几率。几率也是辐射能量的函数,因为从X线管产生的X线由全部能谱组成，具有“线衰减系数”的某种组织的衰减性U质是一个复杂的函数，随着辐射情况的变化，它可以表现出不同的数值。对于能量为的单能射线，假设其投射的物体是均匀且为同一介质，被投E射物体的线性衰减系数为，厚度为，则在入射线强度为，出射X射线强UX0I度为时有I即（21）UXEI0/LN0I这就是著名的LAMBERTBEER定律。实际被投射的物体是非均匀的，如果在X射线束通过的路径上，介质是L非均匀的，则可以将X射线穿过的介质沿扫描路径人为划分为大小相等的个LN小方块，每个小块的厚度为，块内视为均匀介质。这样，每个小块有一个均D匀的线性衰减系数，如图所示，设各个方块的线性衰减系数分别为。X射线通过每一个方块的衰减为，射线通过第二个方NU,321DUEI101块的衰减为，通过第个方块的衰减为。将上述各式代DUEI212NN入整理，有图22X射线通过非均匀介质时的体素划分（22）DUUNEI0321即（23）PN/321式中，是通过射线强度检测得到的检测值，是一个已知量，称为LN0IP投影。这样，上式中，等号右边可视为一个常数，等待求解的是等式左边的每一个方块的衰减系数。由于未知量太多，这个方程没有唯一解。考虑到围绕受检体的一次扫描可以得到一个有关衰减系数的方程，如果将整个感兴趣断层都划分为大小相等的小方块，每个小方块由于一个固定的衰减系数，这样，上述单束射线环绕扫描得到的每个穿透强度值都可以形成一个关于衰减系数的方程。只要扫描线不重复且足够多，这些衰减系数的确定值肯定可以通过解线性方程组的方法解出唯一实数解。把每个方块的衰减系数值用灰度表示，就可以重建出以衰减系数为特征的断层图像。这种图像重建方法称为方程组法。值得说明的是，当上图中时，投影就是0DP（24）DLUILN0式中，衰减系数是随路径连续变化的函数。上式说明，入射X射线强LUL度与出射X射线强度之比经对数运算后，表示沿射线路径上衰减系数的线积分，而投影与射线穿越介质的路径长度成正比。P然而，一般二维断层图像至少应划分出160160个像素，如要用方程组法重建这个断层图像，就要25600个独立的方程联立求解，这是一个相当费时的任务。尽管1967年世界首台CT试验机采用的就是这个图像重建方法，但当断层被划为更小更多的体素时，即使采用今天最先进的计算机技术，解一组这样的方程组也不是件容易的事。所以，这种方法在后来的实际应用中已不再使用，而是采用下一章将要介绍的反投影解析法。第三章CT图像解析重建算法通过上一章已经知道CT图像重建的基本原理，但由于使用最简单直接的方程组法不易求得断层的各个像素的衰减系数。这一章介绍现在CT最常用的图像重建算法反投影解析法。主要介绍了由投影重建图像的系统模型，引进RADON变换和傅立叶变换，它们对应于X射线CT的二维重建问题；且给出了RADON变换在图像重建中的具体形式，即截面函数沿着直线的线积分等于其RADON变换，而图像重建过程即是将截面函数沿许多不同角度下直线的线积分所产生的投影数据进行逆RADON变换从而得到截面函数。另外由逆RADON变换的推导给出了两种常用的图像解析重建法滤波反投影法和卷积反投影法。31反投影重建法（BP）反投影是一个求和的过程，它把所有角度上的数据都加在一起。它把投影域中的数据沿着原来的投影路径“涂抹”回去，但不改变数据的值。3直接反投影重建是指，断层平面中某像素的线性衰减系数可由平面内所有经过该像素射线的投影和得出，而整幅图像重建可以看做为所有方向上扫描投影的累加。通俗地讲，所谓直接反投影，就是在对某体面一个方向的扫描完成后，以得到的投影为灰度值沿着扫描路径经过的体素回抹到体素对应的像素上。改变方向后的多次扫描形成多次回抹，同一像素上多次回抹的灰度累加即完成图像重建。图31断层的像素和射线模型图32反投影重建示例如图32所示，其中（A）为原图像像素，（B）为反投影重建后图像，（C）为（B）中像素值除以投影线数后得到的平均值。有上述例子可以看出原图中像素值为零的像素，经反投影重建后不再为零。即有伪迹。直接反投影图像产生星状伪迹的原因在于，该方法是把取自有限物体空间的投影均匀地回抹（反投影）到了射线所及的无限空间的各个像素上，包括原来像素值为0的点。图33反投影重建的等效系统设为处于处的一个点源，即二维断面上的一个冲击,YXU0,YX,YX函数，这时的输出就是系统的冲击响应。FH已知（31）21,YXYH可见，反投影重建过程后得到的图像不是原来的图像，系统存在星状伪,迹。根据信号与系统理论，对于图33有（32）,YXHUYXF式中，两个星号表示二维卷积。4如果想要去除星状伪迹，可在系统的输出端再加上一个滤波器，设其时域函数为，传递函数为，为了使该加了滤波器的系统输出的图像,YXQ,21Q等于系统输入的原图，也就是要求,YXU（34）,YXHQ即（35）,1,2YXXYQ对该式取二维傅立叶变换，得到（36）1,221YXQ即（37）2121,Q这是一个二维滤波器。即使加上某种近似，这种滤波器也不容易实现。如果的值取无穷，则是一个不可能实现的理想滤波器。为此提出了滤波反21和投影图像重建法。32滤波反投影图像重建法（FBP）为了消除直接反投影法产生的星状伪迹，提出了滤波反投影重建的方法，这种方法是在直接反投影重建方法的基础上增加一个滤波器。增加滤波器的思路有两种，如图34所示。图34两种不同的去除星状伪迹反投影图像重建思路所谓滤波反投影重建法是把获得的投影函数先做卷积处理，即人为设计一种一维滤波器，用它对所得的投影函数进行改造，然后把这些改造过的投影函数进行反投影和累加等处理，就可以达到在一定程度上消除星状伪迹的目的。在介绍滤波反投影重建法之前先补充滤波反投影重建法要用到的数学知识和投影定理中心切片定理。33RADON变换与傅立叶变换为了有效和快速地对图像进行处理，常需要将原定义在图像空间的图像以某种形式转换到另外一些空间，并利用在这些空间的特有性质方便地进行一定的加工，最后再转换回图像空间以得到所需的效果。通常7正变换图像空间到其他空间反变换其他空间到图像空间331RADON变换RADON变换1917年由奥地利数学家RADON提出，是CT图像重建思想的数学表达。RADON从理论上证明了某物理量的二维分布函数可由该函数在其定义域内的所有方向上的线积分完全确定。RADON变换的意义在于，只要知道了一个未知二维分布函数的所有方向上的线积分，那么就能够求得该二维分布函数。所谓CT断层的图像重建，就是求取能够反应断层内部结构和组成的某物理量的二维分布。断层扫描获得的投影实际上可视为被测物体断层内部结构和组成的不同方向上的线积分，所以RADON变换的正变换和逆变换正好对应了CT技术的射线投影和图像二维分布函数的重建过程。对的RADON变换定义为沿由和定义的直线的线积分。,YXF,PRFPL若设表示物体的一个二维断层分布，通过RADON变换的到沿不同方向,F直线的线积分相当于得到物体不同方向的投影。图35定义RADON变换的坐标系统上述线积分可写为DLYXFPRF,（38）如果借助DELTA函数，上述线积分还可写为DXYXPYFPRFSINCO,（39）上式通常称为的RADON变换。,YXF由上述可知，关于某直线的RADON变换就是沿该直线上的一,YXF维投影。实际上当将也视作可变参数时，就是一个二元函数，此时RADON,PRF变换是空间域到变换域的映射。域中的每个点对应于空间域,YX,P中的一条直线，这里实际上就是图像分析中直线的HOUGH变换,YX。10332傅立叶变换的定义DEFTFII21（310）其中DTEFFI（实自变量的复值函数）（311）称为TF的FOURIER变换，记为。TFDEFTI21称为的FOURIER逆变换，记为。F1F34中心切片定理中心切片定理的含义是平行投影的一维傅立叶变换等同于原始物体的二维傅立叶变换的一个切片。将已知投影数据通过一个简单的二维傅立叶反变换可以得到物体截面的一个评估。图36中心切片定理示意图通过上图可知某断层（或它对应的图像）在视角为时得到的平,YXF行投影（函数）的一维傅立叶变换，等于二维傅立叶变换过原,F,21F点的一个垂直切片，且切片与轴相交成角。19数学表达如下首先定义代表截面的函数的二维傅里叶变换DXYEYXFVUFVUJ2,（312）类似的定义某个角度下的一条投影的傅里叶变换TPDTETPSJ2（313）举个最简单的例子沿着0的直线方向，物体投影的傅里叶变换等于频域里面0的情形DXYEFUFUJ2,0,（313）因为相位因子不再依赖于，在此可以将积分分成两部分，YDXEYXFUUJ2,0,（314）从平行投影的定义看，可以将上式括号中的部分看成是沿着常量的积分DYXFP,0（315）将其代入314得到DXEUFUJ20,（316）等式右边表示的是投影的一维傅里叶变换；因此可以得到垂直方向投影和0P投影函数的二维变换的关系式0,USF（317）这是傅里叶切片定理最简单的形式。很明显物体与坐标系统间的位置方向是独立的。例如，在图36中，如果,坐标系统为,坐标系统旋转，STYX式316定义的投影的傅里叶变换等于物体沿着一旋转了角的直线的二维傅里叶变换。从而引出了成像系统中著名的中心（傅里叶）切片定理。2即某图像在视角时的平行投影的傅里叶变换给出了图像二维傅里,YXF叶变换在与轴夹角为的一个切片，此切片通过原点。换句话说，的,VUF傅里叶变换在数值上对应于沿着图中所示直线BB的值。TP,VUF通过将坐标系统旋转得到坐标系统后，傅里叶切片定理可以得到,YXST更广泛的应用。在系统中沿着常量为的直线进行的投影可表示成,STDSTFTP,（318）其傅里叶变换DTETPSJ2（319）将以上两式合并，可以发现DTESTFSJ2,（320）借助旋转坐标定义可以得到关于坐标系统的变换形式,YXDXYEFSYJSINCO2,（321）等式右边表示一个空间域的傅里叶变换，则SI,VUINCOS,，FS（322）这个等式是直接射线断层成像的基础，从而也证明了中心（傅里叶）切片定理。上面的结果显示从取图像分别在与轴成角度的直线上,YXFXK,21的一系列投影，并对这些投影逐个进行傅里叶变换，可以得到图像的二,YXF维傅里叶变换在与轴成角的直线上的值。理论上，当投影线,VUFK,21束无限多时，就能求得在频域上相应线束的无限多条直线的值，从而可,以利用傅里叶逆变换得到原图。傅里叶逆变换可以采用直角坐标系公式,YXF也可以采用极坐标公式DUVEVUFYXFYXJ2,（323）FTJ|,20（324）其中SINCOSYXT（325）但从多个一维傅立叶函数拟合二维傅立叶函数需要足够精确的傅立叶空间插值。由于理想的插值方法还在探索中，这阻碍了基于中心切片定理的二维傅立叶逆变换图像重建法的发展。35平行投影重建算法351算法原理傅里叶切片定理仅仅提供了一种断层成像的简单概念模型，实际应用中需要一种不同的算法。滤波反投影算法是目前广泛应用于所有直线透射断层成像的算法。此算法源于中心切片定理，不仅极其准确而且保证了快速执行的稳定性。它是通过转换极坐标中傅里叶逆变换和重新确定积分限来实现的。首先介绍平行束投影的反投影算法。先回顾一下傅里叶反变换的公式DXYEVUFYXFVUJ2,（326）将频域内的直角坐标系换成极坐标系，有如下替换,VU,SIN,COSV（327）改变微分项DUV（328）则傅里叶反变换的极坐标形式DEFYXFYXJSINCO220,（329）此积分可以分成两部分，从18和360,DEFDEFYXFYXJYXJSINCOS220SINCO20,）（）（330利用特性，上式可以写成,FDEFYXFTJ|,20（331）此处。利用321式用视角处的投影的傅里叶变换代SINCOXTS替二维傅里叶变换，可以得到,FDESYXFTJ|,20（332）上式是滤波反投影重建的图像表达式。对上式进行不同的数学变形，将导致不同的物理解释，相应于不同的重建算法。由此得到的重建算法有滤波（卷积）反投影法F（C）BP和RADON反变换重建法。352滤波（卷积）反投影有卷积定理函数卷积的傅立叶变换是函数傅立叶变换的乘积，即一个域中的卷积相当于另一个域中的乘积。例如时域中的卷积等于频域中的乘积。设在某一旋转角时，采的投影，滤波函数为，则滤波后的投影,TPTH为（333）,TTTG由于数据在计算机实现时是离散的，故应进行离散卷积。反投影是滤波反投影算法的另一个核心，其实就是前面介绍的直接反投影算法，只不过是在反投影前将投影数据进行了卷积滤波预处理。由此可得极坐标下卷积反投影重建图像如下（334）0COS|,DTGRFRT353RADON反变换重建算法在CT成像原理中已经介绍过，对图像某一方向的投影表示沿射线路径方向上衰减系数的积分，而对图像进行某一视角上的RADON变换求得的的就是图像该方向上的投影，所以对投影进行RADON反变换就能实现图像重建。对滤波后的图像投影数据进行反RADON变换得到重建图像的方法就是RADON反变换重建法。36解析重建算法的实现和对比分析针对本章介绍的几种重建算法，文本用MATLAB进行了模拟重建，下面采用的是调用系统函数PHANTOM生成的标准SHEEPLOGON头模型。如图37。图37SHEEPLOGON头模型用IRADON函数用不同滤波器进行滤波反投影重建的仿真图像如图38所示图38对标准头模型解析重建的图像对比下面是对几组图像处理中常用的图片直接读入再进行解析法模拟重建的效果图，重建过程中采用了各种不同的滤波器。图39采用各种不同的滤波器进行滤波反投影重建效果对比在图像重建中，一般用视觉相似度SSIM对图像重建的效果进行评估。视8觉相似度越接近1说明重建图像与原图越相进。采用不同的滤波器进行滤波反投影重建的图像与原始图像相比的视觉相似度及重建所用的时间如表31所示。表31采用不同滤波器滤波反投影重建图像与原图的视觉相似度对比滤波器名称SSIM时间（秒）直接反投影00328460401583RAMLAK09999650450677SHEPPLOGAN09999660337064HANN09999680349063HAMMING09999680286156COSINE09999670293614通过图39各重建图像效果的对比和表31中视觉相似度的对比可知，直接反投影重建由于没有滤波，存在伪迹，图像不清晰，与原图的视觉相似度很小。采用了不同的滤波器，相对于没有滤波的情况，视觉相似度都很接近1，说明与原图比较接近，图像质量有了明显的提高。以上介绍的是平行束投影重建算法。CT扫描方式有平行束和扇束扫描两种扫描方式。扇束投影重建算法大致可分为两类一类是重排算法，即把一个视图中采集到的扇形数据重新组合成平行的射线投影数据，再用上面所介绍的平行束算法重建。这种方法对扇束投影的视角与每一扇束投影各射线间的增角有一定的约束条件。另一种算法是扇束投影直接重建算法，这种算法不必先把数据重排，而是根据扇束投影数据自身的特点，发展出直接用扇形束投影数据重建图像的方法。扇束扫描自身的算法有两种，区别在于抽样的方式上，一种是等角抽样，即射线是按等角分布的，将检测器等角分布在以放射源为中心的圆弧或直线上；另一种是等间距抽样，检测器单元在直线上作等间距据分布。这里不作详细介绍。第四章迭代重建算法41迭代重建算法的思想原理迭代重建算法的概念与前面讲述的解析法最大的区别在于前者一开始就把连续的图像离散化。在解析算法中，假设图像是连续的，每个像素只是一个点，这些离散的点是以图像显示为目的的，这些点的选择可以是随意的，与图像重建无关。但在迭代算法中，每个像素是个小面积。这些小面积在计算当前图像的投影数据时要用到。像素模型对重建图像的质量好坏影响很大。图像重建的解析算法与图像重建的迭代算法之间的另一个区别是解析算法着力对一个积分方程求解，而迭代算法着力对一个线性方程组求解。3迭代的主导思想是，假设断层截面由一个未知的矩阵组成，然后由测量投影数据建立一组未知向量的代数方程式，通过方程组求解未知图像向量。此方法虽然比变换重建法易于理解，但对于图像要求较高的应用来说，这一方法缺少精度和速度。然而在不可能获取大量投影数据或者投影不均匀分布于180或360之间情况下，迭代法是一个非常有效的方法。411迭代重建的数学模型图41断层截面与投影的几何布置如图41所示，图中将截面离散化，即在图像上叠加一个方格网。,YXF每一个小方格内像素值是一常量，令表示第个单元的常量值，N表,YXFJJ示像素的总数。对于迭代法而言射线的定义稍有不同。具体地讲，此时穿过平面的射线具有一定的宽度，为了说明这一点在图中用阴影加重了第,YX根射线，且每一根都是如此，绝大多数情况下射线宽度大约等于图像像素单元I的宽度。此时线积分被称为射线和。如图41所示，投影数据也采用单索引表示，设是图中由第I根射线测IY量的射线和。则的关系可以表示如下，IJYX和（41）NJIJIYXA1MI,21式中为射线总数（或投影总数），而为加权因子，其代表第个像素MIJAJ对第根射线积分的贡献。其值等于第根射线穿过第个像素所占部分的面积，IJ如图中所示。需要注意，大多数均为0，对于任意一条射线和，只有很少的IJ像素单元有贡献。如果值很小的话，可以用传统的矩阵理论方法对式N和41进行求逆。实际中对于128128的图像阵，16384，如果要重建N256256的图像，值将更大。同时往往具有与同样量级的数。假设M等于16384的情况，NM和式41中的矩阵的大小等于1638416384，而这将使任何直接求解IJA矩阵的逆变得无法实现。当然，当测量数据混有噪声或时，即使值很N小，直接矩阵求逆也是不允许的，这种情况可能不得不使用最小二乘法来求解。对于为很大值的情况，有一些很吸引人的迭代方法来求解式41。NM和这些方法是基于“投影法”的，最早被KACZMARZ提出，此后由TANABE进1415一步阐明。为了说明这些方法中的计算步骤，可以把41式写成如下形式，（42）MNMMNYXAXAXAYXAXAX212221211令。因此，由表示的一幅图像可以看成是维空间中TNX,21XN的一个点。在这一空间中，上面方程组中每一个方程式表示一个超平面。当这些方程组的唯一解存在的话，这些超平面的交点为一个点，这就是42式方程组的解。这一概念进一步用图42说明，为理解方便，这里仅考虑有两个像素的情况，它们满足下面的方程，21,X（43）112122AXY图42KACZMARZ交替投影法求解线性方程组的示意图找寻图42中解的计算过程如下先从一个初始图像开始，把这一初始0X向量投影到第一根射线，把投影点再投到第二根射线上，然后将投影点再投回到第一根射线上，依此类推。若唯一解存在，迭代终究会收敛到那一点。对于这一方法的计算机实现而言，首先在可行解中估计一个初始解。这一估计解是维空间的一个点，以表示。一般情况下可以设置N,0201NXX。将这一初始点投影到式42中第一个等式所表示的超平面得到0,0X，图42所示为两个变量的情况。再将投影到式43中第二个等式所表11示的超平面得到，依次下去。当投影到42式中的第个等式所表示2XKI的超平面时，得到的过程可以用下面的数学方式表示K（44）IIIIIAYX11式中，而与其自身的内积。,21INIIAAIIA为如前述，迭代重建的计算过程如下给定一个初始解，将其投影到方程42表示的某个超平面上，得到另一个解，再将其投向下一个超平面上，依此循环，当所有方程所表示的超平面都被投影过后，得到。在下一次迭代中，MX作为初始解，再依次投影到42式中表示的各个超平面上，获得，如MXM2此重复。TANABE证明，若方程42有唯一解，则有15（45）XKMKLIM当式42描绘的M个超平面彼此正交，且原方程有唯一解，经一次迭代循环便能收敛。若M个超平面彼此成很小的夹角，则解的收敛速度较慢。很明显，超平面间的夹角很大程度上影响了求解的收敛速度。由此可知，为了快速得到收敛解，有两种途径可以解决采用GRAMSCHMIDT过程超将平面正交化，如RAMAKRISHNAN等提出了一种两两正交的方案；仔细选择投影超平面的顺1序，以增大两次投影超平面间的夹角。123当MN，式42则不会收敛到唯一解，而在超平面交点附近产生振荡效应；当MN，式42式无唯一解，可能会有无限多个解。在44式基础上发展了一系列其他的迭代重建法，很多都是其的近似表达。为了在后面部分更好的说明，不妨将式44重新整理，（45）IKNKIIJIJAZYX12式中（46）NKIIIXAXAZ1式45说明将第次迭代解投影到第个超平上时，第个像素的灰度值1IIJ（其当前值为）可以通过差值得到，JXJ（47）IJXIKNKIIIJIJAZY1242代数重建法ART上一节介绍了迭代算法的基本原理，我们知道迭代法的要旨是从一个假设的初始图像出发，采用迭代的方法，将根据人为设定并经理论计算得到的投影值同实验测得的投影值比较，不断进行逼近，按照某种最优化准则寻找最优解。这种方法的优点是更容易把与空间几何相关的、与测量值条件相关的各种校正因子包括进去。例如容易进行空间分辩的不均匀性校正；散射衰减的校正；物体几何形状约束和平滑性约束等操作。现在讲述有上述原理提出的迭代算法。首先介绍ART方法。本文用三种算法实现了ART迭代重建。421KACZMARZ方法KACZMARZ算法是ART类算法中最经典的一种算法。这种方法是一种16所谓的行的操作方法，因为每次迭代都要“扫过”矩阵A中所有的行。这种方法引入的松弛参数。K它的实现公式为（48）,2,1|,121,0,MKIIKIIXIKMAXB其中，松弛参数1，为射线总条数。为迭代次数KI其基本步骤为设定任一初始图像向量，一般情况下取；NRX0TX0,0迭代步骤式48，射线号循环选择。1ODMITANABETANABE1971证明了，若测量方程43是相合的，上述迭代算法收敛到的一个解。否则，EGGERMONT及CENSOR引进了松弛系列，目BAX1718K的是有效的减小噪声。松弛系数通常是一个关于迭代次数的一个函数，随着I的增大而减小，且满足。由此引出了下面两种ART方法RANDOMIZEDI1IKACZMARZ方法（随机KACZMARZ方法）和SYMMETRICKACZMARZ方法（对称KACZMARZ方法）。通过改变从此系数和迭代次数来改善重建图像的质量。跟KACZMARZ相比，不同之处在于SYMMETRICKACZMARZ方法的而KACZMARZ方法的只取而RANDOMIZED2,31,21MII,21MKACZMARZ方法（随机KACZMARZ方法）随机选择，但选中的概率正比于。2|IA43联立迭代重建法SIRT对于44式方程组，还可以采用数值计算中的JACOBI迭代法来求解，与19ART迭代法不同，JACOBI迭代法采用的是并行迭代，即仅当所有投影数据都计算完以后，才对图像值进行更新。具体来说，采用47式计算，但的值并IJXJ不立即改变，在经历完所有的等式计算后，才修正的值，由此可以看出这种JX修正是每个像素对应于所有迭代变化的平均值。这种迭代重建算法被称为同步迭代重建法SIRT。具体步骤设定任一初始图像向量，一般情况下取；NRX0TX0,0迭代步骤411IMIIIIAY11此算法以减慢收敛速度换取比ART较好的图像质量。本论题模拟图像重建时应用的类似的方法有联合代数重建法SART（SIMULTANEOUSALGEBRAICRECONSTRUCTIONTECHNIQUE）、平均分量20231法（COMPONENTAVERAGINGMETHOD）、CIMMINO投影法（CIMMINOSPROJECTION25METHOD）、对角松弛征正交投影法（DIAGONALLYRELAXEDORTHOGONAL26PROJECTIONS）和经典LANDWEBER法（THECLASSICALLANDWEBERMETHOD）4。可以参考相关文献，在此不一一赘述。298744全差分迭代算法（TV）在上面ART的算法基础上，再加一些约束条件形成的另一种算法TV（TOTALVARIATION算法。待重建的图像表示为一个含有维的向量IMAGENF它的运算公式如下（412）TSTSTSTTSTTTVTSFFFF,21,2,1,|其中表示二维图像宽和高；21,WH和。IMAGENW45迭代重建算法的实现和对比分析针对42和43介绍的几种方法，本文进行了计算机模拟重建，采用的是自定义的SHEEPLOGAN头模型，大小为128128，如图43。图43SHEEPLOGAN头模型利用旋转角度为05180，平行射线束为75的平行扫描投影数据，进行KACZMARZ算法、SYMMETRICKACZMARZ算法、RANDOMIZEDKACZMARZ算法的重建图像如下。图44ART三种算法分别迭代10次的重建效果对于迭代法，迭代次数，扫描的角度和每个角度下射线的条数都会影响图像重建的效果。下面将从重建图像的视觉效果上用控制变量法分析比较迭代次数，扫描的角度以及每个角度下射线的条数对重建图像的具体影响。图45和图46是其他条件一样但每个角度的射线条数不同的情况下采用ART迭代算法重建的图像。图46和图47是比较扫描的角度对迭代重建图像的影响。图47和图48是比较分析迭代次数对迭代重建效果的影响。图45LENA图像ART迭代重建图像图46CORTEX图像ART迭代重建图像图47CENTAUR2图像ART迭代重建图像图48MRI图像ART迭代重建图像利用旋转角度为05180，平行射线束为75的平行扫描投影数据，进行SIRT中的联合代数重建法SART（SIMULTANEOUSALGEBRAICRECONSTRUCTIONTECHNIQUE）、平均分量法（COMPONENTAVERAGINGMETHOD）、2023125CIMMINO投影法（CIMMINOSPROJECTIONMETHOD）、对角松弛征正交投影法26（DIAGONALLYRELAXEDORTHOGONALPROJECTIONS）和经典LANDWEBER法（THE4CLASSICALLANDWEBERMETHOD）的重建图像如图49所示。2987图49SIRT各种算法的图像对比图410BARB图像SIRT迭代重建图411MRI图像SIRT迭代重建图412CAMERAMAN图像SIRT迭代重建图413CORTEX图像SIRT迭代重建如图41