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

密级： 公开 NANCHANG UNIVERSITY 学 士 学 位 论 文 THESIS OF BACHELOR（20092013年）题 目 基于Matlab实现的CT重建算法仿真比较研究 学 院： 信息工程 系 电子信息工程 专业班级： 学生姓名： 学号： 指导教师： 职称： 起讫日期： 2013.3.18 2013.6.7 南 昌 大 学学士学位论文原创性申明本人郑重申明：所呈交的论文是本人在导师的指导下独立进行研究所取得的研究成果。除了文中特别加以标注引用的内容外，本论文不包含任何其他个人或集体已经发表或撰写的成果。对本文的研究作出重要贡献的个人和集体，均已在文中以明确方式表明。本人完全意识到本申明的法律后果由本人承担。作者签名： 日期：学位论文版权使用授权书本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅。本人授权南昌大学可以将本论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。 保密，在 年解密后适用本授权书。本学位论文属于 不保密。（请在以上相应方框内打“”）作者签名： 日期：导师签名： 日期：摘要基于Matlab实现的CT重建算法仿真比较研究 专 业： 学 号： 学生姓名： 指导教师： 摘要计算机断层成像（即CT（Computerized Tomography）技术是一项近年来快速发展多学科交叉的先进技术，除了医学应用，还应用在工业无损检测、射电天文学、精密仪器反演等多个重要领域。CT图像重建是CT技术的核心，重建算法的优劣直接关系到对检测结果判断的准确性。医学CT普遍采用解析成像技术，这主要因为医用CT可以采集到完整的投影数据。然而在工业应用中，一般得不到完整的投影数据，因此要利用不完整的投影数据得到较清晰的图像，必须采用迭代重建法。针对以上情况，本课题作了以下两方面的分析比较。1.能获取完整投影数据时，用计算机模拟了解析法的图像重建。发现图像重建的速度快，且重建的图像与原图的视觉相似度接近1。2.不能获取完整投影数据时，用计算机模拟了迭代算法的图像重建。先用ART迭代算法引入松弛参数，分析了松弛参数重建图像的影响。接着用SIRT算法重建图像，这种算法所得的图像比ART算法较好。最后在ART算法基础上采用全差分对重建图像进行优化，图像质量有了明显提高。关键字：图像重建算法；解析重建；迭代重建AbstractBased on Matlab realize the CT reconstruction algorithm simulation comparative studyAbstractComputed tomography (CT (Computerized Tomography) technology is a rapid development in recent years, multidisciplinary advanced technology, in addition to medical applications, but also used in industrial nondestructive testing, radio astronomy, precision instruments and other important areas of the inversion. CT image reconstruction CT technology is the core of the reconstruction algorithm is directly related to the accuracy of the test results determine. Medical CT imaging commonly used analytical technique, which can be collected mainly due to a complete medical CT projection data. However, in industrial applications, the general lack of complete projection data, and therefore the use of incomplete projection data to be more clear image, an iterative reconstruction method must be used.For the above, this project made the following two aspects of the analysis and comparison.1. can get the full projection data, using a computer simulation of the analytical method of image reconstruction. Found fast image reconstruction and the reconstructed image with the original visual similarity is close to 1.2.can not get the full projection data, using a computer simulation of the iterative algorithm for image reconstruction. First with the introduction of ART relaxation parameter iterative algorithm analyzes the impact of relaxation parameters of the reconstructed image. Followed by SIRT algorithm for reconstruction of images, the resulting image of this algorithm is better than the ART algorithm. Finally ART algorithm based on a fully differential optimized for image reconstruction, image quality has improved significantly.Keywords:Image reconstruction algorithm；analytic reconstruction；iterative reconstructionI目录目录摘 要IAbstractII第一章 绪论11.1 引言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 Radon变换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 全差分迭代算法（TV）274.5 迭代重建算法的实现和对比分析27第五章 结论37参考文献38附录40谢辞60III第一章绪论第一章 绪论1.1 引言自从X射线发现后，医学上就开始用它来探测人体疾病。但是，由于人体内有些器官对X线的吸收差别极小，因此X射线对那些前后重叠的组织的病变就难以发现。于是，美国与英国的科学家开始了寻找一种新的东西来弥补用X线技术检查人体病变的不足。CT(电子计算机体层成像)是70年代初放射诊断的项重突破，CT不是X线摄影，而是用X线对人体扫描（图1.1所示），取得信息,经电子计算机处理而获得的重建图像它能使传统的X线检查难以显示的器官及其病变显示成像，而且图像逼真，解剖关系明确，从而扩大了人体的检查范围，大大提高了病变的早期检出率和诊断准确率。这种检查简便、安全、无痛苦、无创伤、无危险，它促进了医学影像诊断学的发展，CT的研制成功被誉为自伦琴发现X射线以后，放射诊断学上最重要的成就。发明者亨斯菲尔德和科马克共同获得了1979年的诺贝尔奖金。CT最初只用于头部检查，1974年又出现了全身的CT。在短短10余年间，CT已遍及全球，从第一代发展到第五代。我国各大城市医院所使用的CT属第三代。全身CT可以作头、胸、腹、骨盆的横断扫描，也可作甲状腺、脊柱、关节和软组织及五官等小部位的区域扫描。CT最适于查明占位性病变如肿瘤、囊肿、增大的淋巴结、血肿、脓肿和肉芽肿的大小、形态、目和侵犯范围，它可以决定某些器官癌肿的分期和是否能进行手术切除。在某些情况下，CT还能区别病变的病理特性如实性、囊性、血管性、炎性、钙性、脂肪等。CT成像技术的基本原理是根据人体中不同组织对X射线的衰减作用不同，用X线束对人体某部一定厚度的层面进行扫描，由探测器接收透过该层面的X线，转变为可见光后，由光电转换变为电信号，再经模拟/数字转换器（analog/digital converter）转为数字，输入计算机处理。图像形成的处理有如对选定层面分成若干个体积相同的长方体，称之为体素（voxel）。扫描所得信息经计算而获得每个体素的X线衰减系数或吸收系数，再排列成矩阵，即数字矩阵（digital matrix），数字矩阵可存贮于磁盘或光盘中。经数字/模拟转换器（digital/analog converter）把数字矩阵中的每个数字转为由黑到白不等灰度的小方块，即像素（pixel），并按矩阵排列，即构成CT图像。所以，CT图像是重建图像。每个体素的X线吸收系数可以通过不同的数学方法算出。如图1.1所示。 图1.1 计算机断层成像示意图1.2 CT的分类在投影重建的原理和方法的框架下，根据投影方式不同产生了不同种类的CT成像术。透射断层成像透射断层成像（Transmission Computed Tomography,TCT，简称CT）系统中，从发射源射出的X射线穿透物体到达接收器。射线在通过物体时被物体吸收一部分，余下部分被接收器接收。由于物体各部分对射线的吸收不同，所以接收器获得的射线强度实际上反映了物体各部分对射线的吸收情况。由于透射成像术中采用X射线，所以也叫X-CT。接收器接收到的模拟信号经模数转换器转换成数字信号后，把代表着不同角度下的投影数据送给计算机，它再运用复杂的数学方法重建出物体截面的图像。在医学上它主要针对骨质成像。XCT成像原理如图1.2所示：图1.2 X射线CT成像原理发射断层成像在发射断层成像（Emission Computed Tomography,ECT）系统中，发射源在物体内部，一般是将具有放射性的离子注入物体内部，从物体外检测放射出来的量。通过这种方法可以了解离子在物体内运动情况的分布，从而可以检测到与生理相关的状况。现在通常用的ECT主要有两种：正电子成像PET（Positron Emission Tomography）；单光子成像SPECT（Single Positron Emission Computed Tomography）。它们的区别是，采用单光子辐射器时，检测器所接收的光子计数便是投影测量数据；而采用正电子辐射器时，由于在正电子辐射点的毫米级范围内，一个正电子会俘获一个负电子，并产生沿相反方向运动的两个中子，因而测量数据是在一个很小时间间隔内由两个性质相反的检测器检测到的光子对的数目。ECT主要针对代谢过程的成像。磁共振成像磁共振成像（Magnetic Resonance Imaging,MRI）在早期称为核磁共振（Nuclear Magnetic Resonance,NMR）。它的发明者斯坦福大学的F.Bloch和哈佛大学的E.M.Purcell为此获得了1952年度的诺贝尔物理奖。它的工作原理简介如下：氢核以及其它具有奇数个质子和中子的核，包含具有一定磁动量和旋量的质子。如果把它们放在磁场中，它们会像陀螺在地球重力场中一样在磁场中进动。一般情况下质子在磁场中时随意排列着的，当一个适当强度和频率的共振场信号作用于物体时，质子吸收能量并转向与磁场相交的朝向。如果此时将共振场信号除去，质子吸收的能量释放并被接收器检测到。根据检测到的信号就可以确定质子的密度。通过控制所用的共振场信号和磁场强度，可每次检测到沿透过目标中一条线的信号。换句话说，检测到的信号是MRI信号沿直线的积分。所以检测目标的工作成了投影重建的问题。由这些信号可以得到极好的器官解剖图，在细节上远远超过由X射线CT产生的图像。由于它能够很好的保留丰富的软组织图像，所以在脑功能研究上，MRI胜过ECT。超声断层成像超声断层成像（Ultrasound CT,UCT）与上述的CT不同的地方是超声射线在人体组织中不一定按直线传播，尤其是在软组织和硬组织的分界面会产生折射效应。因此超声成像主要用于人体中没有硬组织（如骨骼）部分的检测，特别是妇女乳房癌的检测。反射断层成像反射断层成像（Reflection CT,RCT）也是利用投影重建的原理工作的。常见的一个例子是雷达系统，其中的雷达图是物体反射的回波所产生的。例如在前视雷达（FLR）中，雷达发射器从空中向地面发射无线电波；雷达接收器在特定角度所接收到的回波强度是地面发射量的一个扫描段的积分。通过发展硬件设备和改进软件算法，可以不断降低噪声和提高图像的清晰度。CT技术在临床的许多领域得到愈来愈广泛的应用。应该指出，上面介绍的各类CT的功能是相互补充而不是相互替代的。1.3 CT的组成结构以上几种成像技术尽管使用的物理系统各不相同，但所有系统所测量到的数据均具有物体某种感兴趣的物理特性在空间分布的积分形式。成像技术广义上讲大致可分为三个部分：信息的收集、重建以及处理和传递。数据采集，即获得携带有物体内部特征信息的数据。所采用的物理测量系统包括放射源、检测器及配套的机械、电路设备。为进行计算机处理，系统还需将测量得到的模拟量转化成数字量。图像重建，即利用各种反演算法从所测数据重建出物体内部某种物理量的分布图。这一步是成像技术的核心。图像后处理，包括图像去噪、锐化，特征提取以及图像的存储、压缩和传输等。1.4 论文的组织结构论文主要研究辐射成像技术中的透射断层成像理论。讲述了CT重建的基本原理第一章绪论部分，首先介绍了计算机断层成像术（CT）的发展历史，及其医学意义根据射线产生的物理方式对CT具体的分类。之后详细介绍了透射断层成像的研究现状，进而对于目前投影重建中存在的问题进行了阐述，从而引出了论文的工作目标和内容。第二章详细介绍了CT图像重建系统的基本原理，根据Beer定理对断层进行体素划分，通过联立方程组求解方程组的方式重建图像。但当断层被划分为更小更多的体素时，即使采用今天最先进的计算机技术解这一组方程组也不是件容易的事。所以这个方法在后来实际的图像重建中已不再使用，从而引出了现在CT图像重建中常用的方法解析法。第三章详细介绍了现在CT常用的图像重建算法解析法的算法原理，研究进展。先讲述了直接反投影重建算法，但由于直接反投影法存在图像伪迹，从而提出了滤波反投影重建算法。在运算的过程中根据卷积定理从而提出了卷积反投影重建算法，最后对上述三种解析算法重建的图像进行了对比分析。第四章迭代重建算法部分，上述解析算法是目前医用CT中最常用的重建算法。但它并不适用于所有投影重建场合。在实际应用中，有时无法测得大量的投影，有时投影不是均匀的分布在180度或360度范围内。例如，为了避免物体（如脏器）运动，或为了减小剂量，投影数据采集不足；又如某物体在一个方向上尺寸特别长，则在这一方向上的投影数据无法或较难测得。再如用成像法勘探地球资源时，用孔穴法采集的数据，数量不足又不均匀等上述滤波反投影法无能为力，于是提出了图像重建的迭代重建算法。本章介绍了迭代算法的成像模型及其成像原理。最后演示了一下迭代算法重建图像的效果。第五章总结部分，本章对论文中讲述的CT重建算法进行了认真的总结和比较。探讨了尚需进一步解决的问题及今后有待努力的方向。4第二章CT图像重建的基本原理第二章 CT图像重建的基本原理2.1 所需数据的扫描获取假设采用单射线源平移旋转方式来获取图像重建所需的数据，该方法需要一个射线发生器和一个射线探测器。用X射线发生器发出单色细束X射线，对成像对象进行投射，如图2.1所示。设X线源发出的X射线强度是，探测器接收到的X射线强度是一次扫射完成后将X射线与探测器平移一个小的步长，再次对成像对象进行投射，探测器又接收到一个新的透射强度，在该投射角度上反复小步长平移投射并记录检测的数据后，将X射线源与探测器绕被测对象旋转一个小的角度后在进行平行投射并平移，如此重复多次，直到旋转180度为止。取得这些不同投射位置的透射强度后，就可以用这些数据，采用一定的算法重建出X射线源环绕的断层的图像。图2.1 单射线源平移旋转扫描示意图2.2 图像重建算法CT与普通的X线摄影术之间有着非常重要的区别。在CT技术中，组织对X线的局部衰减特性被用于离散成像；而在通常的X线摄影术中，这种衰减信息则重叠在X线底片上。组织对X线的这种局部衰减特性是X线与被检测物体之间的若干相互作用过程的产物。例如光电吸收过程和康普顿散射过程。这些过程中的每一种都有它自己的发生几率。几率也是辐射能量的函数,因为从X线管产生的X线由全部能谱组成，具有“线衰减系数”的某种组织的衰减性质是一个复杂的函数，随着辐射情况的变化，它可以表现出不同的数值。对于能量为的单能射线，假设其投射的物体是均匀且为同一介质，被投射物体的线性衰减系数为，厚度为，则在入射线强度为，出射X射线强度为时有即（2.1）这就是著名的Lambert-Beer定律。实际被投射的物体是非均匀的，如果在X射线束通过的路径上，介质是非均匀的，则可以将X射线穿过的介质沿扫描路径人为划分为大小相等的个小方块，每个小块的厚度为，块内视为均匀介质。这样，每个小块有一个均匀的线性衰减系数，如图所示，设各个方块的线性衰减系数分别为。X射线通过每一个方块的衰减为，射线通过第二个方块的衰减为，通过第个方块的衰减为。将上述各式代入整理，有：图2.2 X射线通过非均匀介质时的体素划分（2.2）即（2.3）式中，是通过射线强度检测得到的检测值，是一个已知量，称为投影。这样，上式中，等号右边可视为一个常数，等待求解的是等式左边的每一个方块的衰减系数。由于未知量太多，这个方程没有唯一解。考虑到围绕受检体的一次扫描可以得到一个有关衰减系数的方程，如果将整个感兴趣断层都划分为大小相等的小方块，每个小方块由于一个固定的衰减系数，这样，上述单束射线环绕扫描得到的每个穿透强度值都可以形成一个关于衰减系数的方程。只要扫描线不重复且足够多，这些衰减系数的确定值肯定可以通过解线性方程组的方法解出唯一实数解。把每个方块的衰减系数值用灰度表示，就可以重建出以衰减系数为特征的断层图像。这种图像重建方法称为方程组法。值得说明的是，当上图中时，投影就是（2.4）式中，衰减系数是随路径连续变化的函数。上式说明，入射X射线强度与出射X射线强度之比经对数运算后，表示沿射线路径上衰减系数的线积分，而投影与射线穿越介质的路径长度成正比。然而，一般二维断层图像至少应划分出160160个像素，如要用方程组法重建这个断层图像，就要25600个独立的方程联立求解，这是一个相当费时的任务。尽管1967年世界首台CT试验机采用的就是这个图像重建方法，但当断层被划为更小更多的体素时，即使采用今天最先进的计算机技术，解一组这样的方程组也不是件容易的事。所以，这种方法在后来的实际应用中已不再使用，而是采用下一章将要介绍的反投影解析法。7第三章 CT图像解析重建算法第三章 CT图像解析重建算法通过上一章已经知道CT图像重建的基本原理，但由于使用最简单直接的方程组法不易求得断层的各个像素的衰减系数。这一章介绍现在CT最常用的图像重建算法反投影解析法。主要介绍了由投影重建图像的系统模型，引进Radon变换和傅立叶变换，它们对应于X 射线CT的二维重建问题；且给出了Radon变换在图像重建中的具体形式，即截面函数沿着直线的线积分等于其 Radon变换，而图像重建过程即是将截面函数沿许多不同角度下直线的线积分所产生的投影数据进行逆Radon变换从而得到截面函数。另外由逆Radon变换的推导给出了两种常用的图像解析重建法：滤波反投影法和卷积反投影法。3.1 反投影重建法（BP）反投影是一个求和的过程，它把所有角度上的数据都加在一起。它把投影域中的数据沿着原来的投影路径“涂抹”回去，但不改变数据的值。直接反投影重建是指，断层平面中某像素的线性衰减系数可由平面内所有经过该像素射线的投影和得出，而整幅图像重建可以看做为所有方向上扫描投影的累加。通俗地讲，所谓直接反投影，就是在对某体面一个方向的扫描完成后，以得到的投影为灰度值沿着扫描路径经过的体素回抹到体素对应的像素上。改变方向后的多次扫描形成多次回抹，同一像素上多次回抹的灰度累加即完成图像重建。图3.1 断层的像素和射线模型图3.2 反投影重建示例如图3.2所示，其中（a）为原图像像素，（b）为反投影重建后图像，（c）为（b）中像素值除以投影线数后得到的平均值。有上述例子可以看出原图中像素值为零的像素，经反投影重建后不再为零。即有伪迹。直接反投影图像产生星状伪迹的原因在于，该方法是把取自有限物体空间的投影均匀地回抹（反投影）到了射线所及的无限空间的各个像素上，包括原来像素值为0的点。图3.3 反投影重建的等效系统设为处于处的一个点源，即二维断面上的一个冲击函数，这时的输出就是系统的冲击响应。已知（3.1）可见，反投影重建过程后得到的图像不是原来的图像，系统存在星状伪迹。根据信号与系统理论，对于图3.3有：（3.2）式中，两个星号表示二维卷积。如果想要去除星状伪迹，可在系统的输出端再加上一个滤波器，设其时域函数为，传递函数为，为了使该加了滤波器的系统输出的图像等于系统输入的原图，也就是要求（3.4）即（3.5）对该式取二维傅立叶变换，得到：（3.6）即（3.7）这是一个二维滤波器。即使加上某种近似，这种滤波器也不容易实现。如果的值取无穷，则是一个不可能实现的理想滤波器。为此提出了滤波反投影图像重建法。3.2 滤波反投影图像重建法（FBP）为了消除直接反投影法产生的星状伪迹，提出了滤波反投影重建的方法，这种方法是在直接反投影重建方法的基础上增加一个滤波器。增加滤波器的思路有两种，如图3.4所示。图3.4 两种不同的去除星状伪迹反投影图像重建思路所谓滤波反投影重建法是把获得的投影函数先做卷积处理，即人为设计一种一维滤波器，用它对所得的投影函数进行改造，然后把这些改造过的投影函数进行反投影和累加等处理，就可以达到在一定程度上消除星状伪迹的目的。在介绍滤波反投影重建法之前先补充滤波反投影重建法要用到的数学知识和投影定理中心切片定理。3.3 Radon变换与傅立叶变换 为了有效和快速地对图像进行处理，常需要将原定义在图像空间的图像以某种形式转换到另外一些空间，并利用在这些空间的特有性质方便地进行一定的加工，最后再转换回图像空间以得到所需的效果。通常正变换：图像空间到其他空间反变换：其他空间到图像空间3.3.1 Radon变换Radon变换1917年由奥地利数学家Radon提出，是CT图像重建思想的数学表达。Radon从理论上证明了某物理量的二维分布函数可由该函数在其定义域内的所有方向上的线积分完全确定。Radon变换的意义在于，只要知道了一个未知二维分布函数的所有方向上的线积分，那么就能够求得该二维分布函数。所谓CT断层的图像重建，就是求取能够反应断层内部结构和组成的某物理量的二维分布。断层扫描获得的投影实际上可视为被测物体断层内部结构和组成的不同方向上的线积分，所以Radon变换的正变换和逆变换正好对应了CT技术的射线投影和图像二维分布函数的重建过程。对的Radon变换定义为沿由和定义的直线的线积分。若设表示物体的一个二维断层分布，通过Radon变换的到沿不同方向直线的线积分相当于得到物体不同方向的投影。图3.5 定义Radon变换的坐标系统上述线积分可写为：（3.8）如果借助Delta函数，上述线积分还可写为：（3.9）上式通常称为的Radon变换。由上述可知，关于某直线的Radon变换就是沿该直线上的一维投影。实际上当将也视作可变参数时，就是一个二元函数，此时Radon变换是空间域到变换域的映射。域中的每个点对应于空间域中的一条直线，这里实际上就是图像分析中直线的Hough变换。3.3.2 傅立叶变换的定义（3.10）其中（实自变量的复值函数）（3.11）称为的Fourier变换，记为。称为的Fourier逆变换，记为。3.4 中心切片定理中心切片定理的含义是平行投影的一维傅立叶变换等同于原始物体的二维傅立叶变换的一个切片。将已知投影数据通过一个简单的二维傅立叶反变换可以得到物体截面的一个评估。图3.6 中心切片定理示意图通过上图可知：某断层（或它对应的图像）在视角为时得到的平行投影（函数）的一维傅立叶变换，等于二维傅立叶变换过原点的一个垂直切片，且切片与轴相交成角。数学表达如下：首先定义代表截面的函数的二维傅里叶变换：（3.12）类似的定义某个角度下的一条投影的傅里叶变换：（3.13）举个最简单的例子：沿着=0 的直线方向，物体投影的傅里叶变换等于频域里面= 0的情形（3.13）因为相位因子不再依赖于，在此可以将积分分成两部分，（3.14）从平行投影的定义看，可以将上式括号中的部分看成是沿着常量的积分（3.15）将其代入(3.14)得到（3.16）等式右边表示的是投影的一维傅里叶变换；因此可以得到垂直方向投影和投影函数的二维变换的关系式：（3.17）这是傅里叶切片定理最简单的形式。很明显物体与坐标系统间的位置方向是独立的。例如，在图3.6中，如果,坐标系统为,坐标系统旋转，式(3.16)定义的投影的傅里叶变换等于物体沿着一旋转了角的直线的二维傅里叶变换。从而引出了成像系统中著名的中心（傅里叶）切片定理。即某图像在视角时的平行投影的傅里叶变换给出了图像二维傅里叶变换在与轴夹角为的一个切片，此切片通过原点。换句话说，的傅里叶变换在数值上对应于沿着图中所示直线BB的值。通过将坐标系统旋转得到坐标系统后，傅里叶切片定理可以得到更广泛的应用。在系统中沿着常量为的直线进行的投影可表示成：（3.18）其傅里叶变换（3.19）将以上两式合并，可以发现（3.20）借助旋转坐标定义可以得到关于坐标系统的变换形式 （3.21）等式右边表示一个空间域的傅里叶变换，则（3.22）这个等式是直接射线断层成像的基础，从而也证明了中心（傅里叶）切片定理。上面的结果显示：从取图像分别在与轴成角度的直线上的一系列投影，并对这些投影逐个进行傅里叶变换，可以得到图像的二维傅里叶变换在与轴成角的直线上的值。理论上，当投影线束无限多时，就能求得在频域上相应线束的无限多条直线的值，从而可以利用傅里叶逆变换得到原图。傅里叶逆变换可以采用直角坐标系公式也可以采用极坐标公式：（3.23）（3.24）其中（3.25） 但从多个一维傅立叶函数拟合二维傅立叶函数需要足够精确的傅立叶空间插值。由于理想的插值方法还在探索中，这阻碍了基于中心切片定理的二维傅立叶逆变换图像重建法的发展。3.5 平行投影重建算法3.5.1 算法原理傅里叶切片定理仅仅提供了一种断层成像的简单概念模型，实际应用中需要一种不同的算法。滤波反投影算法是目前广泛应用于所有直线透射断层成像的算法。此算法源于中心切片定理，不仅极其准确而且保证了快速执行的稳定性。它是通过转换极坐标中傅里叶逆变换和重新确定积分限来实现的。首先介绍平行束投影的反投影算法。先回顾一下傅里叶反变换的公式：（3.26）将频域内的直角坐标系换成极坐标系，有如下替换：（3.27）改变微分项（3.28）则傅里叶反变换的极坐标形式：（3.29）此积分可以分成两部分，从和,（3.30)利用特性，上式可以写成,（3.31）此处。利用(3.21)式用视角处的投影的傅里叶变换代替二维傅里叶变换，可以得到：（3.32）上式是滤波反投影重建的图像表达式。对上式进行不同的数学变形，将导致不同的物理解释，相应于不同的重建算法。由此得到的重建算法有滤波（卷积）反投影法(F（C）BP)和Radon反变换重建法。3.5.2 滤波（卷积）反投影有卷积定理：函数卷积的傅立叶变换是函数傅立叶变换的乘积，即：一个域中的卷积相当于另一个域中的乘积。例如：时域中的卷积等于频域中的乘积。设在某一旋转角时，采的投影，滤波函数为，则滤波后的投影为（3.33）由于数据在计算机实现时是离散的，故应进行离散卷积。反投影是滤波反投影算法的另一个核心，其实就是前面介绍的直接反投影算法，只不过是在反投影前将投影数据进行了卷积滤波预处理。由此可得极坐标下卷积反投影重建图像如下：（3.34）3.5.3 Radon反变换重建算法在CT成像原理中已经介绍过，对图像某一方向的投影表示沿射线路径方向上衰减系数的积分，而对图像进行某一视角上的Radon变换求得的的就是图像该方向上的投影，所以对投影进行Radon反变换就能实现图像重建。对滤波后的图像投影数据进行反Radon变换得到重建图像的方法就是Radon反变换重建法。3.6 解析重建算法的实现和对比分析针对本章介绍的几种重建算法，文本用matlab进行了模拟重建，下面采用的是调用系统函数phantom生成的标准Sheep Logon头模型。如图3.7。图3.7 Sheep Logon头模型用iradon函数用不同滤波器进行滤波反投影重建的仿真图像如图3.8所示： 图3.8 对标准头模型解析重建的图像对比下面是对几组图像处理中常用的图片直接读入再进行解析法模拟重建的效果图，重建过程中采用了各种不同的滤波器。图3.9 采用各种不同的滤波器进行滤波反投影重建效果对比在图像重建中，一般用视觉相似度SSIM对图像重建的效果进行评估。视觉相似度越接近1说明重建图像与原图越相进。采用不同的滤波器进行滤波反投影重建的图像与原始图像相比的视觉相似度及重建所用的时间如表3.1所示。表3.1采用不同滤波器滤波反投影重建图像与原图的视觉相似度对比滤波器名称SSIM时间（秒）直接反投影0.0328460.401583Ram-Lak0.999965 0.450677Shepp-Logan0.9999660.337064Hann0.9999680.349063Hamming0.9999680.286156Cosine0.9999670.293614通过图3.9各重建图像效果的对比和表3.1中视觉相似度的对比可知，直接反投影重建由于没有滤波，存在伪迹，图像不清晰，与原图的视觉相似度很小。采用了不同的滤波器，相对于没有滤波的情况，视觉相似度都很接近1，说明与原图比较接近，图像质量有了明显的提高。以上介绍的是平行束投影重建算法。CT扫描方式有平行束和扇束扫描两种扫描方式。扇束投影重建算法大致可分为两类：一类是重排算法，即把一个视图中采集到的扇形数据重新组合成平行的射线投影数据，再用上面所介绍的平行束算法重建。这种方法对扇束投影的视角与每一扇束投影各射线间的增角有一定的约束条件。另一种算法是扇束投影直接重建算法，这种算法不必先把数据重排，而是根据扇束投影数据自身的特点，发展出直接用扇形束投影数据重建图像的方法。扇束扫描自身的算法有两种，区别在于抽样的方式上，一种是等角抽样，即射线是按等角分布的，将检测器等角分布在以放射源为中心的圆弧或直线上；另一种是等间距抽样，检测器单元在直线上作等间距据分布。这里不作详细介绍。21第四章 迭代重建算法第四章 迭代重建算法4.1 迭代重建算法的思想原理迭代重建算法的概念与前面讲述的解析法最大的区别在于前者一开始就把连续的图像离散化。在解析算法中，假设图像是连续的，每个像素只是一个点，这些离散的点是以图像显示为目的的，这些点的选择可以是随意的，与图像重建无关。但在迭代算法中，每个像素是个小面积。这些小面积在计算当前图像的投影数据时要用到。像素模型对重建图像的质量好坏影响很大。图像重建的解析算法与图像重建的迭代算法之间的另一个区别是解析算法着力对一个积分方程求解，而迭代算法着力对一个线性方程组求解。迭代的主导思想是，假设断层截面由一个未知的矩阵组成，然后由测量投影数据建立一组未知向量的代数方程式，通过方程组求解未知图像向量。此方法虽然比变换重建法易于理解，但对于图像要求较高的应用来说，这一方法缺少精度和速度。然而在不可能获取大量投影数据或者投影不均匀分布于180或360之间情况下，迭代法是一个非常有效的方法。4.1.1 迭代重建的数学模型图4.1 断层截面与投影的几何布置如图4.1所示，图中将截面离散化，即在图像上叠加一个方格网。每一个小方格内像素值是一常量，令表示第个单元的常量值，N表示像素的总数。对于迭代法而言射线的定义稍有不同。具体地讲，此时穿过平面的射线具有一定的宽度，为了说明这一点在图中用阴影加重了第根射线，且每一根都是如此，绝大多数情况下射线宽度大约等于图像像素单元的宽度。此时线积分被称为射线和。如图4.1所示，投影数据也采用单索引表示，设是图中由第 i 根射线测量的射线和。则的关系可以表示如下，（4.1）式中为射线总数（或投影总数），而为加权因子，其代表第个像素对第根射线积分的贡献。其值等于第根射线穿过第个像素所占部分的面积，如图中所示。需要注意，大多数均为0，对于任意一条射线和，只有很少的像素单元有贡献。如果值很小的话，可以用传统的矩阵理论方法对式(4.1)进行求逆。实际中对于128 128 的图像阵，=16384，如果要重建256256的图像，值将更大。同时往往具有与同样量级的数。假设等于16384的情况，式(4.1)中的矩阵的大小等于16384 16384，而这将使任何直接求解矩阵的逆变得无法实现。当然，当测量数据混有噪声或时，即使值很小，直接矩阵求逆也是不允许的，这种情况可能不得不使用最小二乘法来求解。对于为很大值的情况，有一些很吸引人的迭代方法来求解式(4.1)。这些方法是基于“投影法”的，最早被Kaczmarz提出，此后由Tanabe进一步阐明。为了说明这些方法中的计算步骤，可以把(4.1)式写成如下形式，（4.2）令。因此，由表示的一幅图像可以看成是维空间中的一个点。在这一空间中，上面方程组中每一个方程式表示一个超平面。当这些方程组的唯一解存在的话，这些超平面的交点为一个点，这就是(4.2)式方程组的解。这一概念进一步用图4.2说明，为理解方便，这里仅考虑有两个像素的情况，它们满足下面的方程， （4.3）图4.2 Kaczmarz交替投影法求解线性方程组的示意图找寻图4.2中解的计算过程如下：先从一个初始图像开始，把这一初始向量投影到第一根射线，把投影点再投到第二根射线上，然后将投影点再投回到第一根射线上，依此类推。若唯一解存在，迭代终究会收敛到那一点。对于这一方法的计算机实现而言，首先在可行解中估计一个初始解。这一估计解是维空间的一个点，以表示。一般情况下可以设置。将这一初始点投影到式(4.2)中第一个等式所表示的超平面得到，图4.2所示为两个变量的情况。再将投影到式(4.3)中第二个等式所表示的超平面得到，依次下去。当投影到(4.2)式中的第个等式所表示的超平面时，得到的过程可以用下面的数学方式表示（4.4）式中，而与其自身的内积。 如前述，迭代重建的计算过程如下：给定一个初始解，将其投影到方程(4.2)表示的某个超平面上，得到另一个解，再将其投向下一个超平面上，依此循环，当所有方程所表示的超平面都被投影过后，得到。在下一次迭代中，作为初始解，再依次投影到(4.2)式中表示的各个超平面上，获得，如此重复。Tanabe证明，若方程(4.2)有唯一解，则有（4.5）当式(4.2)描绘的M个超平面彼此正交，且原方程有唯一解，经一次迭代循环便能收敛。若M个超平面彼此成很小的夹角，则解的收敛速度较慢。很明显，超平面间的夹角很大程度上影响了求解的收敛速度。由此可知，为了快速得到收敛解，有两种途径可以解决：采用Gram-Schmidt过程超将平面正交化，如Ramakrishnan等提出了一种两两正交的方案；仔细选择投影超平面的顺序，以增大两次投影超平面间的夹角。当MN，式(4.2)则不会收敛到唯一解，而在超平面交点附近产生振荡效应；当MN，式(4.2)式无唯一解，可能会有无限多个解。在(4.4)式基础上发展了一系列其他的迭代重建法，很多都是其的近似表达。为了在后面部分更好的说明，不妨将式(4.4)重新整理，（4.5）式中（4.6）式(4.5)说明将第次迭代解投影到第个超平上时，第个像素的灰度值（其当前值为）可以通过差值得到，（4.7）4.2 代数重建法ART上一节介绍了迭代算法的基本原理，我们知道迭代法的要旨是：从一个假设的初始图像出发，采用迭代的方法，将根据人为设定并经理论计算得到的投影值同实验测得的投影值比较，不断进行逼近，按照某种最优化准则寻找最优解。这种方法的优点是更容易把与空间几何相关的、与测量值条件相关的各种校正因子包括进去。例如容易进行空间分辩的不均匀性校正；散射衰减的校正；物体几何形状约束和平滑性约束等操作。现在讲述有上述原理提出的迭代算法。首先介绍ART方法。本文用三种算法实现了ART迭代重建。4.2.1 Kaczmarz 方法Kaczmarz算法是ART类算法中最经典的一种算法。这种方法是一种所谓的行的操作方法，因为每次迭代都要“扫过”矩阵A中所有的行。这种方法引入的松弛参数。它的实现公式为：（4.8）其中，松弛参数=1，为射线总条数。其基本步骤为：设定任一初始图像向量，一般情况下取；迭代步骤：式(4.8)，射线号循环选择。TanabeTanabe1971证明了，若测量方程(4.3)是相合的，上述迭代算法收敛到的一个解。否则，Eggermont及Censor引进了松弛系列，目的是有效的减小噪声。松弛系数通常是一个关于迭代次数的一个函数，随着的增大而减小，且满足。由此引出了下面两种ART方法Randomized Kaczmarz方法（随机Kaczmarz方法）和symmetric Kaczmarz方法（对称Kaczmarz方法）。通过改变从此系数和迭代次数来改善重建图像的质量。跟Kaczmarz相比，不同之处在于symmetric Kaczmarz方法的而Kaczmarz方法的只取而Randomized Kaczmarz方法（随机Kaczmarz方法）随机选择，但选中的概率正比于。4.