工业CT窄角扇束的迭代法图像重建算法研究

《工业CT窄角扇束的迭代法图像重建算法研究》由会员上传分享，免费在线阅读，更多相关内容在学术论文-天天文库。

中文摘要摘要本文对工业CT窄角扇束扫描方式下的迭代图像重建算法进行了研究。为了提高重建图像质量，解决普通迭代法重建图像所产生的“盐和胡椒”现象。论文提出新的代数迭代重建算法(ARTti)，并且比较了新的重建算法和原有的迭代重建算法的原理和性能，重点研究了在不完全数据重建时，ARTII型算法的性能。迭代法是图像重建算法中的一大重要类别，尤其是在处理投影数据不完全及噪声干扰较为严重的情况下具有不可代替的优势。但是普通迭代法(ARTI型)在处理密度变化较大的工件的投影数据时，重建图像基础出现“盐和胡椒”现象，而且为了得到一个较高质量的图像，其迭代次数过多，导致重建效率降低。为了解决这一问题，论文首先分析了迭代法的过程，给出了其几何和物理性质。在此基础上，论文提出了新的迭代算法(ARTII型)。该算法改进已有的迭代过程，加入了新的迭代参数：收敛因子，通过改变迭代进程和调整收敛因子来提高重建图像的质量。此外论文还比较了该法和同时迭代法的性能。作者在微机上成功的实现了从数据采集仿真到图像重建的全过程，并且模拟实现了各种噪声源。给出了各种应用情况下的重建结果，并对结果进行了对比分析和评估，证明了舢盯II型算法在投影数据不完全时和噪声污染较大时可以有效的提高重建图像的质量。关键字：图像重建，投影，迭代，收敛因子 英文摘要AbstractInthispaper,theauthordiscussesindetailanefficientiterativeimagereconstructionalgorithminindustrialcomputerizedtomography(ICT)withthenarrowfan-beamprojection．Inordertoimprovethequalityofreconstructionimageandsolvetheproblemthat“saltandpepper'’phenomenaOCCurSinthereconstructionimage,thepaperproposesanovelAlgorithmReconstructionTechnique(ARTII)，comparethetheoryandcapabilityofbetweenthem．ThepaperfocusesontheARTII’Scapabilitywhenitcomestoreconstructionwi山incompleteandnoisyprojectiondata．Iterativereconstructionmethodisanimportantkindofreconstructionmethods，whichhassuperioradvantageswhenreconstructingimageswithincompleteandnoisyprojectiondata．However,inprocessingtheprojectiondataofinspectedworkpieceswithlargerangeofdensity,therewillbe‘‘saltandpepper'’inthereconstructionimageofARTI．Otherwise，inordertogetthehighqualityreconstructionimagesARTIneedadditionaliterativetimes，whichleadsintodecre嬲i119ofiterativeefficiency．Forthesakeofsolvingthisproblem，firstthepaperanalyzestheprocessofiterativemethodandpresentsthegeometryandphysicspropertyofit．Consideringthis，theauthorproposeanoveliterativealgorithm(ARTII)．Thenewmethodamelioratestheprocessofolditerativemethodandaddanewiterativeparameter：convergencefactor,whichcarlimprovethequalityofreconstructionimagebychangingtheiterativeprocessandadjustingconvergencefactor．Theauthorcompletestheentireprocedurefromdataacquisitiontoimagereconstructiononcompeer．TheefficiencyofARTIIisevaluatedincontrasttothatofitstraditionalserialcounterpart．ExperimentresultsdearlydemonstratethatARTIIcanefficientlyimprovequalityofreconstructedimageswhenprocessingtheincompleteprojectiondataornoisyprojectiondata．KeyWords：projectionimage，reconstruction,itemtivealgorithm，Convergencefactor,ART,SIRTⅡ 1引言1引言1．1工业CT技术发展状况工业cT，却工业计算机断层成像(IndustrialComputerizedTomography)，简称ICT，是计算机技术与放射探测技术相结合而形成的一种崭新的数值成像技术。工业CT在无损检测(NOT)与无损评价(NDE)领域得到了广泛的应用，它在无损的状态下获得被检断面的二维灰度图像，然后以图像的形式清晰、准确、直观的展现被检物体的内部特征、装配情况、材料密度和缺陷的位置、大小以及性质。工业CT在本质上是一种射线检测技术，但与常规的射线检测相比，有其独特的优点“1。工业CT可以给出被测工件的断层扫描图像，从图像上可以直观地看到目标细节的空间位置、形状、大小，感兴趣的目标不会受到周围其他细节的阻挡，图像容易识别和理解：它没有一般透视照相法普遍存在的影像重叠和模糊．空f司分辨率和密度分辨率比常规的射线检测技术高两个数量级。工业CT技术被国际无损检测界公认为最佳的无损检测手段。20世纪的70年代由于计算机技术和放射学的密切接台，产生了图像技术的新一代，其结果人们称之为是医学上的一次革命，这就是计算机断层成像技术“1(或称断层扫描)。CT的闷世存放射学界引起了爆炸性的轰动，被认为是继伦琴发现x射线后，工程界对放射学诊断的叉一划时代的贡献。1979年的诺贝尔生理和医学奖破例地授予两位没有专门医学经历的科学家：Oodfrey．N．Housfiold与AllanMacleodCoormack。从此，放射诊断学进入了CT时代。CT装置很快推广，并在很短期间内获得极大发展。1971年世界上第一台CT系统诞生，1973年美国MayoClinic和麻省总医院即相继安继了颅脑CT。1974年美国制成第一代全身CT。此后20多年来tCT装置很快得到推广，并获得了极大的发展。从第一台CT机问世至今，CT设备己从第1代发展至第5代，其各项性能和速度都有了很大提高。扫描时间从4分钟到5分钟缩短到1秒钟甚至更短，最新式的CT机扫描速度可达每秒24层(横断面图像)，可以跟上血液在器官和组织中的流动，因而可以对心脏作动态检查。后来，莱德利设计成功全身CT装置，进一步扩大了CT的检查范围，取得了更大的效益。除并射线CT外，其他型号的CT也相继问世，如单光子发射CT、核磁共振CT等均已付诸临床应用。超声CT、微波CT的研究也取得了极大的进展。毫无疑问，CT已成为影像诊断学领域中不可缺少的检查手段。计算机断层成像术作为一种高性能的无创伤诊断技术，显然已在医学成像领域确定了其不可动摇的地位。在医学应用中取得成功后，CT技术很快被应用到了工业检测上．发展成为工 重庆大学硕士学位论文业计算机断层扫描成像技术(ICT)。进入八十年代以来，发达的工业国家已经把盖射线或y射线的ICT用于航天、航空、军事、冶金、机械、石油、电力、地质等产业部门的无损检测及无损评价。从使用的射线源来分，工业CT主要分为x射线工业CT和T射线工业CT，Y射线源体积小，性能稳定而且穿透能力强：但J射线源的能量范围宽，检测效率高。工业CT技术发展非常迅速，1986年美国科学与公司推出由15Mev直线加速器提供的高能Z射线，可检测直径2ra的断面，密度分辨率达O．1％，最小可测裂纹宽为0．08ram。1992年美国研究与应用公司推出由60Mev直线加速器提供的高能射线，可检测直径为3．9m的固体火箭发动机。我国的CT技术起步较晚，于八十年代束才开始研究，到九十年代己逐步把ICT用于无损检测。不仅如此，CT作为一种技术，既有坚实的数学理论为依托，又有现代微电子与计算技术相支撑，必然在其他领域也会得到广泛应用。事实上，CT在工业生产上，在地球物理研究上，甚至在农业、林业和环境保护方面都已取得了令人瞩目的成果并展示了美好的前景，如用于反应堆组件的无损评估，火箭发动机、导弹等部件及钢板焊缝的无损检测，以及水泥制品的质量检查等，CT扫描仪均能明察秋毫；而常规的无损探伤技术对此则往往不能胜任。1988年9月，清华大学通过了“Y射线工业CT实验系统研究”的博士论文(系统采用C060，NaI单探头，平行束扫描，旋转极坐标重建法)。1991年10月，在熏庆大学通过了Y射线四探头实验系统鉴定。重庆大学ICT研究中心，是全国最早从事ICT研究的科研单位之一，在全体科研人员的共同努力下，作为火炬计划预备项目的我国第一台实用型Y射线工业CT机：XN-1300于1993年5月诞生。重庆大学ICT研究中心经多年研究，现已成功开发出一系列的y射线第一、第三代ICT系统，其性能获得业界的好评。2001年通过国家验收的CD-．400CG的主要技术指标达到了90年代未国际水平，标志着我国工业CT技术已步入先进国家的行列。为了更好进行工件检测(尤其是大直径、高密度工件的检测)以满足我国建设的需要，重庆大学ICT研究中心决定探索新的检测方法，研制第二代工业CT机。本论文课题就是在这样的情形下提出的，对原有CT机图像重建方法进行研究找出适合二代CT图像重建的方法，为重庆大学ICT研究中心的下一步研究作了铺垫和基础。1．2课题意义和国内外研究现状1．2．1课题的实际意义工业CT在无损检测与无损评价方面具有不可代替的重要性，而检测和评价的结果取决于工业CT重建图像(ReconstructionImage)的质量。应用工业CT设备最2 终获得被检测工件的断层图像，“图像重建技术”是ICT的关键技术之一。按照一定的算法将检测所得的投影数据进行数学运算处理，从而获得其断层图像的过程称为图像重建。工业CT所重建出来的图像质量主要是由运动控制、数据采集过程和重建算法等因素决定的。按照数据采集方式划分，工业CT可以分为五代：第一代CT：单源、单探测器的平行束扫描；第二代CT：单源、多探测器的窄角扇束扫描；第三代CT：单源、多探测器的广角扇柬扫描；第四代CT：单源、全周探测器的大扇束全包容扫描：第五代CT：多源、多探测器、扇形束扫描，无相对运动。目前应用最广泛的是广角扇束扫描方式，即三代扫描方式。但随着现代工业的发展需要，尤其是对于检测大直径、高密度工件检测的要求，传统的三代检测方式已不能满足发展的需求，故而需要开发单源、多探测器窄角扇束扫描方式的ICT系统，其原因在于：①传统的三代扫描的ICT机，多采用的是天然射线源，其能量较低无法穿透大型工件，而经过直线加速器产生的高能x射线，射线的出射扇束角较小(一般小于150)射线无法包容整个被检测的断面。②窄角扇束扫描方式是使用单放射源小角度多探测器，它介于第一代与第三代扫描方式之间，扫描效率明显高于第一代扫描方式。窄角扇束扫描方式下的图像重建算法是本课题研究的重点。现有文献中的图像重建算法主要是针对平行束和广角扇束扫描方式的，不能直接应用于窄角扇束扫描方式下所采集的投影数据。在重建算法的选择上，迄今为止，目前实际使用的CT中绝大多数采用的都是卷积反投影(ConventionBackVrojeetion)算法，尤其是在医用CT中。该方法图像质量较高，同时主要由加法和乘法运算组成，计算开销不大，也易于用硬件设计出专用的图像处理系统。但卷积反投影算法的一大弱点在于，当射线投影数I<(1／3)J(，为像素数目)、或物体密度变化剧烈(高对比度)、或具有强噪声场合，重建图像质量会急剧下降133．当数据完全时，代数重建方法在图像质量至少与卷积反投影算法相当，但当数据不完全或噪声污染较为严重的时候，代数重建算法具有不可代替的优势。对于大型工业CT机，大多采用高速赢线加速器来产生z射线源，这种源较天然的x或Y射线源而言，稳定性较差，而且是经高压产生射线，电磁干扰更大，由此带来的噪声污染更加严重。为了能够重建出质量较好的重建图像，我们拟采用代数熏建方法实现图像重建。但是由于重建速度较慢严重限制了它的实际应用。在此我们提出新的迭代算法来加快迭代过程，提高重建速度。此外随着计算技术的发展，计算机速度的提高，代数重建算法会得到更广泛的应用。因此，结合窄角扇束扫描方式实现代数重建算法重建图像，无疑具有重要的学术和实际意义。l_2-2国内外研究现状 重庆大学硕士学位论文迭代重建的概念首先由Gordon、Bender、HermanHl等人于1970年引入图像重建领域。当时遭到了不少学者的指责。1971年Herman等人曾同Cmwther、klug等在JournalofTheoreticalBiology上展开了一场辩论。Cmwther等认为迭代重建的理论根据不足。嗣后，Herman等证明迭代重建法是以“估计理论”作为其坚实的数学基础的”1。S．L．Wood(1978)29．用随机滤波型论解释了迭代重建法，并将迭代重建法与变换法作了统一的处理。经过多年的理论研究与实践，证明迭代重建法在某些方面有着变换法所不及的许多优点。滤波反投影重建算法无能为力的场合，均可借迭代重建法解决(或利用迭代法结合义其他算法解决)。迭代法的最大缺点是速度较慢，不能用硬件实现。但Nassi等将迭代法与卷积反投影算法结合应用解决心脏x—cT要求快速采集数与快速重建的矛盾。从现有文献资料来看，提高迭代法的迭代效率一直是CT图像重建研究的热点，Herman和Lorraine．B．Meyer‘61提出了分块迭代(Block．Iteractive)修正，H．Malcon．HudsontT]提出了子集排序(OrderingSubset)算法，其核心思想在于利用数据分组排列提高迭代速度，但仅适用于广角扇柬扫描方式，而且在提高收敛效率，减少迭代次数方面效果不很明显，无法满足课题要求；DidierLaRard和GuyMazaret引提出了心动阵列(SystolicArray)上利用并行结构实现图像的实时重建，并取得了较好的效果，但它需要专用的并行处理机，无法满足课题的应用要求：CaivirkA．Johnson”1提出了修改迭代算法的计算结构使之适用于并行处理(Parallel-Processing)，以达到缩短重建时间的目的，并收到了很好的效果，这是近年来图像重建方法的研究热点。1．3研究内容为了配合ICT中心的产品开发，本课题旨在研究在窄角扇束扫描方式下如何提高重建图像的质量，并且兼顾重建时间。工业CT和医学CT有很多相同之处，但由于各自的检测对象的不同，又有其各自的特点。医学CT的检测对象是活的生物体，这就要求用于检测的源剂量要小而且检测时间要尽可能的短，此外为了实时检测生物体变化的情况(如心脏的跳动)这就要求重建速度一定要快，这也是为什么医学CT中普遍采用卷积法重建图像的一个原因；工业CT检测的对象都是非生命体，为了达到提高检测分辨率的弱的，人们一般采用大剂量的辐射源，并且适当的加长检测扫描时间，大部分情况下不要求实时检测，这就是为什么可以在工业CT中采用迭代法进行图像重建的原因。迭代重建法又称直接重建法，按照修正顺序包含的是整个模型(矩阵)、一条射线还是一个单点，迭代重建法又可分为：同时重显技术(ILST)、代数重建法(ART)及同时迭代重建法(SIRT)。迭代法按照修正方法性质的不同又可分为“加”型和“乘”型算法。Housfield的第一台CT实际上用的是ART算法。ART类的迭代算法在重4 建效率上优势较为明显，但存在迭代不收敛的问题，而且在数据受到污染较为严重的时候，抑制噪声能力较弱，重建图像质量下降较快；SIRT类算法总是收敛的，在数据噪声特别大的时候或数据丢失较为严重之时，重建质量明显好于ART类算法，但其迭代次数较ART算法多，放重建速度较慢。本课题针对二代扫描所使用源的情况(高能x射线源)和实际检测要求，同时使用这两类重建算法进行模拟实验，找出最适合的图像重建方法。鉴于实验条件所限，不能获得真实物体的扫描数据。因此我们将软件仿其方式获得采集数据，人为加入噪声污染，比较这两种方法各自的优缺点及适用范围。 2图像重建理论分析2．1CT过程我们以描述)卜tT(窄角扇束扫描方式)为例简要说明一下CT的原理。HT机是由以下几个主要部分组成(见图2．1所示)：1。放射源2．检测器3．计算机4．显示器以及相关的辐射防护设备、运动控制设备和监视设备组成。图2．1j。cT机原理图Fig2．1theoryfigofX-Ctmachine从图2．1可见，工射线源与检测器同步地平行移动，被测工件作旋转运动。由Ⅳ射线源发出的均匀的x射线束穿过工件后，由于工件不同组织对Ⅳ射线的吸收率不同，因而检测器在接收到的不同强度的Ⅳ射线数据中，包含了反映工件组织的信息。将从检测器得到信号，通过模数转换器作模拟信号和数字信号的转换，再交由计算机作图像重建，最后以不同的灰阶形式在监视器上显示，或以数字形式存人计算机硬盘。这就是人们所希望的，不经过解剖或其他破坏性手段而得到的剖视图。计算机在)(-℃T机中的作用，是计算从检测器测量到的一系列数据，确定在每个射柬通路上的、不同点的工件组织的特性值，用这些数据重显工件组织特性的图像。 重庆大学硕士学位论文2．2图像重建的定义2．2．1图像的定义与数字化简单地说，一副图像包含了对所表示的客观存在的物体的一种描述信息。从物理上讲，图像是物质或能量的实际分布：从数学上讲，图像是由连续函数和离散函数所组成的集合。本文在讨论图像重建问题的时候，我们通常假设图像具有下面的性质：①图像区域是一个正方形，其中心位置与坐标系原点相重合；②图像由二维函数讹∥表示，其值在图像区域外为零，其中@力为空间点的坐标，任意点向力处函数幅值永，y)正比于图像在该点的光强度(灰度级)。③函数图像佩纠具有非负有界性，即满足0三讹叫乱。通常，在不引起混淆的情况下，将②中的函数称为“图像”。如果图像的区域不同，即使是相同的函数也可以得到不同的图像。在直角坐标系中，图像，是变量X，Y的函数，则，在点@力的值就称为图像在点@力的密度。～幅数值图像就是图像f酝，)0在二维空问坐标上和亮度上都已离散化了的图像，可以把一幅数字图像考虑为一个矩阵，其行和列标出了图像各个点在二维空间的位置，而矩阵元素的值，标出这些点的灰度等级，这样的数字阵列的元素叫做图像元素或象素【101。一般把图像区域分成nXn个相等的正方形，每一个小正方形就是一个象素。在CT中，图像区域裁是重建区域，被澳4物体断层重建图像包括在重建区域内，并且在点(x，y)上的图像密度值对应于点(x，y)在某种固定能量下的相对线性衰减系数值。2．2．2图像重建问题由投影进行图像重建的问题，过去30年中反复地出现在很多科学、医学和技术领域中。应用的范围之广令人吃惊。一方面，用电子显微镜得来的数据重建出噬菌体的分子结构：另一方面，利用送往地球大气层外的火箭所采集的数据重建出超新残星的x射线结构。这些应用看来似乎是各不相同的，但他们有着共同的数学基础和计算基础H”。在所有这些应用中，世界上影响最大的应用也许是在诊断医学方面：计算机断层扫描(CD的出现，已经引起了放射医学上的革命。穿过人体横截面沿着许多直线由测量到的x射线衰减的数据重建出人体横截面上的图像来。根据投影数据重建图像，是一个产生具有某种物理性质的二维分布图像过程，它由已知位置的穿过物体的有限条射线的线积分确定【l21。图像重建就是指根据物体横剖面的投影数据，经过计算机处理后得到该剖面的图像，它是一种由投影数据得到图像的处理技术。 2图像重建理论分析2．3图像重建原理2．3．1CT的物理基础当射线通过物体时，由于吸收和散射，射线将产生衰减，射线被衰减的程度取决于物体的密度、原子的组成以及射线的能量谱。射线通常为工y，0射线，医用CT采用的是低能J射线，工业CT对低能X射线、中能Y射线及高能工射线均有采用。在真空中，按照探测器方向离开辐射源的射线的光子将全部到达该探测器。当一种物体放在辐射源和探测器之间时，在探测器方向上离开辐射源的某些光子将会从射线束中消失，即被吸收和散射。射线通过物体，在任何一点上的衰减量取决于该点上的物质密度和射线束的能量分布。物体的密度分布往往是非均匀的，密度大的物质对射线的衰减大，反之则小。当仅用能量谱的一个能量级，即单色时，每一点将有一个唯一确定的衰减系数。确定断层衰减系数的分布将是断层图像重建的最终目的。图2．2射线在物体内部的衰减Fig2．2theradial’sattenuationintheworkpieee如图2．2，设射线通过物体内部的路径长度为L，将分为Ⅳ个微段：△11，△12⋯。△lIl；设各微段上各点的衰减系数相同，且相应的衰减系数为ul，u2⋯。!-tn，则沿射线在L路径上的总衰减值M为：M=△，l×■l+△忍XA2+．，．+△胁×12nⅣ=∑AliX,uii=1⋯2⋯．．n(2．1)1=1当n一一时，9 重庆大学硕士学位论文怍1im∑△ff×且f～一∞lll。I∥(，)甜(2．2)：其中：^和～沿L的总的衰减函数对单色射线，射线按比尔定律衰减，即：No=MxP“(2．3)其中：Ⅳ，——为射线的入射强度Ⅳ0一为射线穿透物体后的强度将式(2．2)式代入比尔定律有：No=NJxcxp[一f∥(f)讲】(2．4)L对(2．4)式两边取对数，整理后得：J∥(Oel=In(Nl／No)￡显然对确定的射线而言，P=In口n／No)为一常数，故可令：l∥(，)脚2，(2．5)：由此得到的P值，就称为投影值，反映了物体被检测面内射线的衰减特性。对式(2．5)式进行离散处理即为式(2．1)，表明了沿射线的物体衰减系数的累加和即为投影值。因此图像重建的问题就可以归结为：从投影数据P来计算函数值的反演问题。2．3．2图像重建的数学原理CT成像的目的是根据被检物体断面的二维分布。由于物质的原子序数和密度的差异，因而不同物质对z射线的吸收系数p不同，所以可以由“值的微小变化来区别不同的物质结构，并通过计算机图像处理获得高密度分辨率的图像。为了描述两个变量的函数，我们常用直角坐标(笛卡尔坐标)。因此，我们使用p(x，y)表示点(X’y)上的相对线性衰减系数，其中(X，y)被认为是直角坐标系统的，参看图2．3。可是，在更为数学化的工作中可以得出结论，使用极坐标(r，纠更加方便。极坐标(r，函)由如下公式和直角坐标扛，．y)相联系，即：r=√x2+y2工=rCOS庐10 2图像重建理论分析毋=培。1@／x)Y=rsin≯(2．6)＼。E＼⋯。㈣tia。⋯．E＼<：入E＼’．E图2．3CT中的射线积分Fig2．3radialintegralinCT我们使用这种两个极变量的函数来表示函数，它的值Ar,西)代表极坐标f}；纠的几何点上某个物理参数的值(如相对线性衰减系数值)。在CT图像重建问题中，图像区域即重建区域。在前面我们己将一幅图像定义为两个变量的函数，这个函数在图像区域的外面取值为零，这个区域是一个正方形(比如说大小为2EX2E)，它的中心在坐标系的原点．在下文中，我们用f表示两个极变量r和曲的函数，它被用来定义被重建的图像，我们知道：如果1rcos妒陋√2点或1rsiia≯险√2E／阮圳=O(2．7)现在，我们定义两个极变量的因数的雷当变换．首先，我们引进一个符号约定。雷当变换是一个算子，当它作用在一个函数上的时候，就产生另外一个函数。我们用大写的R表示Radon变化算子，如果f是一个函数，它的雷当变换的函数用可表示，J2，在它的域内点(f，口)上的值用Rf(1,即来表示。对于实数对(‘口)，f的雷当变换定义如下：当l；e：O时，[RUN，口)=广，(√瓦≯，O+tan一，(z／1))dz(2．8) 重庆大学硕士学位论文图2．4极坐标示意图Fi92，4ooIarcoordinatessketchm《p姗孙脯眭l=r赫ros(O-。(2．9)当l=0时则有m舶[彤】(o，口)2Ls(一z，目一万／2)出+J：f(z，e+u／2)dz(2·10)在公式(2．8)中，z=O对应于由原点出发并垂直于L的射线和L的交点。同时我们可以观察到：[Rf]q，臼)=[8，1】(一，目+石)=【Rf]q，口+2万)(2．11)又由式(2．7)可得．当m≥,／2E时，有：[Rf]q，口)=0(2．12)由以上方程可以看出，函数彤完全由它的各个点(‘0)(．√2E<，<√2E和D≤口O，l勺<了的平均值为工，并使下面的最大妒∞(工)，一．芘(x)=一∑(x，／Jx)tn(x／Jx)(5．18)』ll式中的平均值工可根据测得的投影数据进行估计，因为某一射线经过的像素是可知的，它的射线和是测得的，从而这一射线上的平均像素值是可知的。考虑覆盖图像区域的一组射线和就可求得足够数目的射线上的平均像素值，据此即可算出工。最大熵准则特别适用于投影数据不完全场合。参考一般的估计问题，我们可将从投影重建图像问题归结为：根据测量矢量只利用式(5．10)估计图像矢量X。故图像矢量的估计问题，可按下面三步解决：(1)选取基本图像；(2)选取最优准则；(3)求解相应方程得出合适的最佳图像。5．1．2ARTI型算法如前所述，我们应根据式(5．10)结合最优准则来求图像矢量x。如果不用最优准则，那么，或者得不到解，或者所得的解不是最优。ART是迭代重建法的一种形式[251，它的特点是：先假设一初始图像x∞’，然后根据x∞’，求一次近似图像xn’，再根据x“’求二次近似图像z幢’，如此继续，直至满足预定条件而后止。在根据X咄’求J‘k．‘’需加～校正值A工咄’。△x咄’只考虑一条射线(例如“号射线)的射线投影的影响，所修正的像素值也限于矗号射线经过的那些像素。下一次则考虑下一条射线即“，号射线。总之，每次校正依次考虑一条射线的射线和并校正该射线经过的像素。它的解是：工(“，任意：工(I)∥)+锭铲‰||‰112“‰木《=0：(5．19)其他； 5ART与SIRT算法‘=k(modI)+1=[k—Int@)I+1】(5．20)』故而ARTI型算法实现的一般步骤如下：(1)首先设定被重建图像的初始估计值x(”，在没有任何先验知识的情况下，一般加形ART算法把图像的所有像素值置为T；乘形ART算法把图像的所有像素值置为1。(2)使用校正公式校正图像各个像素的灰度值。(3)重复步骤(2)，依次使用每个投影去校正图像，一直到该式中的1个方程结束为止。至此，迭代运算的一次循环就算完成了。(4)检验图像的质量，如果图像中各个像素值的校正量可以忽略不记或是达到了预定的标准，整个迭代过程结束。否则返回步骤(2)，进入下一次的迭代循环。在检验图像质量确定是否再次迭代时，使用的就是我们前面叙述的最优准则，目前较多使用的是最小二乘方准则。如果使用最大熵准则就构成了我们前面所述的最大熵重建法(ME)。5．1．4ARTII型算法文献中有多种形式的ART，对应于一定的最优准则“1。但也有某些形式的ART并不确知它对应于哪个最优准则，然而行之有效，我们下面的ART就属于这一类n酏，我们把这类型算法称为ARTII型算法。它将求解方程式(5．10)转化为求解不等式：触酃(5．21)如何求解这一不等式方程?为此我们采用求解线性不等式方程组的松弛法。将式(5．21)展开有：‘1x≤P．i=1，2,3”：，。(5．22)它的解是：工(”，任意；工(★+1)=X(々)∥"×％笋‰‰T∥≤P腩；(5．23)其他；‘=k(m。∽+1=旷州亨)，+1](5．24)上面各式中七代表迭代次数，Int代表取整数，∥‘’称为松弛参数。可以证明，当0<岛s∥‘’兰s2<2(5．25)时，则用这个松弛法解不等式所生成的序列x(”，x‘”，⋯将收敛于R中的一个矢量37 重庆大学硕士学位论文“1(假定R非空)。现在我们来看解不等式松弛法求解的几何意义：H，=体∽x*P。}(5．26)即研是满足式(5．23)q=第i个方程的那些矢量的集合。显然日是R的子集。每个％代表一个超平面(Hyperplane)。若x的维数卢3，超平面为三维空间中的一个平面。如果工的维数J=2，超平面为二维空间(即平面)中的直线。将式(5．23)等号右边的第二项(校正项)改写为：等喾k：J|‰||2“如二尘∑』L“0(5．27)—々为‰的单位向量。故式(5．24)既描述了校正值的大小，又指出了它的方向。Ilo¨具体说来，P4一《z‘‘’p*，则校正值指向r,k的反方向，企图缩短x‘o’至量研的距离；否则，不会改变我们的估计值。校正移动量的大小依赖于∥’的大小，并具有以下几何性质：如果搿”2，工《“’将在衄的另一边，而且x‘“’比z‘‘’更远离岛。上述校正过程是一根射线一根射线依次进行的，故ART算法又称为逐线校正(Ray—by—RayCorrection)算法。文献【5】指出利用实际的z射线投影，根据式(5．23)进行重建图像实际上满足最小范数准则。APT算法根据迭代过程中校正方法的不同又分为加型ART算法和乘形ART算法。加型ART算法解的形式为：工(‘+1)=工(’’+△x(‘’f5．28)我们前面所述的ART算法就是加型ART算法。这里我们给出乘形ART算法，其解的形式为H1：x”’初始矢量，‘，维，常取全1矢量；、工卜斛、”工尹州’2，⋯／。歹§’2∞L／。其中 5ART与SIRT算法‘=k(modI)+1；五(‘)—松弛系数D<￡≤五(‘’≤7：工，>0一图像矢量x中的一个元素；r／k,)．一定义见式(3．2)，取值1或0。其他各量定义同前。可以证明：在一定条件下，式(5．19)所得解的序列工(”，工‘”，⋯收敛于：P=Rx的最大熵解。至此，我们可以将又ARTII型算法重建图像的实现概括为：(1)首先设定被重建图像的初始估计值石(0)，在没有任何先验知识的情况下，一般加形ART算法把图像的所有像素值置为0；乘形ART算法把图像的所有像素值置为1，选择迭代收敛系数∥)．(2)使用式(5．21)校正图像各个像素的灰度值。(3)重复步骤(2)，依次使用每个投影去校正图像，一直到该式中的1个方程结束为止。至此，迭代运算的一次循环就算完成了。(4)判断是否达到预定迭代次数，若是则迭代过程结束。否则重新计算i。返回步骤(2)，进入下一次的迭代循环。5．2用同时迭代重建法(SIRT)实现图像重建SIRT算法(SimultaneousIterafiveReconstructionTechniques)是与ART算法并行的另一种迭代算法。SIRT的提出托71旨在使重建图像对测量误差不敏感。众所周知，ART算法只用到—条射线的射线和(或射线投影)。如果这一射线投影包含误差，则所得的解跟着也引入误差。要减小误差，办法之一是使P=Rx+P的解满足最小二乘方准则，即：氟c功=喜[p，一芸吩三，]2=|lp—R三旷。；．：。，f-1＼』-‘／L)·juJ=(P-Rx)1(p一尺曲我们来看一下式(5·28)的物理意义。根据反投影定理可知R7P为p的反投影，p为测量矢量．因为jh表示(伪)射线的和，故R7b表示(伪)和的射线反投影。如要求解满足最小二乘方法则，那么，这两种情况下的反投影重建应该相等。或者说：估计图像所得射线和的反投影值等于实测投影的反投影值。这就是式(5．22)的物理涵义。 重庆大学硕士学位论文若月为非奇异阵，则其最小二乘方解为工=(R7五)。R7P(5．31)月7p代表反投影运算，如把(∽7尺)’1看作一二维滤波器，则式(5．21)就是反投影加二维滤波。式(5．29)又可化成x=R1(胜7-1P(5．32)式中，(尼R’)一1为一维p滤波器[271，施于p。这时，式(5．32)表示将一维投影滤波后进反投影，即滤波反投影算法。但一般不用直接解法，而用迭代法求解[281fx(o)：R7P{工(。+1’=x耻’+∥”(RrP—Rr。Rx耻’)(5．33)I=工仕’+∥”尺r(P一冗x仕’)式(5．33)的意义是：取测量值P的反投影作为初如图像，自属合理。在求(k+1)次估计z‘“1’时，利用k次估计工啪加上校正图像。校正图像正比于第k次估计的误差矢量的反投影R7(p一触忙’)，。因而每个像素的校正值实际就是通过该像素的所有射线和的误差值之累加，而不是只与一条射线有关。这是与ART算法的根本区别，也是SIRT能有效地抑制测量数据中噪声的根本原因。由于每—个像素的校正值是所合过该像素射线的共同贡献，一些随机误差被平均掉了。为便于迭代计算，我们将式(5．33)改写为：f石(o)=R7P似+1)纠～∥杰o，一‘r川‘(5‘34)Lf—I式(5．34)揭示的物理意义更加明显：只一r：x‘姊表示i号射线的投影值误差，n提供了i号射线经过的像素的信息。q≠o表示，号像素在f号射线的路径上。所以o；一r：x‘‘’)‘表示i号射线的投影误差对所经过的像素的反投影。例如对．，号像素的反投影为(A—r：x‘”)珞，勺为矢量n中的一个元素。于是，∑(Pl-rTx(k))‘f-I表示所有经过_，号像素的诸射线的投影误差对该像素值的校正。因此SIRT的校正过程被称为逐点校正(PointbyPointCorrection)。SIRT算法根据迭代过程中校正方法的不同也分为加型SIRT算法和乘形SIRT算法。以上所述的就是加型SIRT算法，乘形SIRT算法在实际中很少使用．这里我们将不作介绍了。由以上分析结合ART算法的实现过程，我们可以SIRT的实现步骤和ART的 5ART与SIRT算法实现步骤大致相同，只是迭代初始值和迭代中校正方法不同，故在实现中我们完成可以参照ART的实现步骤，修改迭代初始值和校正方法，实现SIRT算法重建图像。5．3技巧的使用在实践中已经发现，图像重建迭代算法的效率常常能够提高，方法是在各步迭代之间对图像向量作某种处理，这些处理在文献中被称为技巧。这些技巧的使用有时被用来加速收敛于满足规定的最优化准则的图像向量，有时它们可以使得移向一个与最优化规定函数不同的图像向量，但是这个图像向量仍然是被重建图像的一个较好的近似值t29}。在文献中较多被使用的技巧“’19’221有：有选择性的平滑、约束和归一化。这些技巧都是建立在重建图像的基图像是一Xn的数字化图像。仔细观察式(5．21)、式(5．27)和式(5．31)，我们可以发现影响迭代的过程主要有基图的选取、松弛参数、各投影数据之间的关系。而实际中大部分的技巧的提出都是从这三方面出发的。下面我们就这三个方面分别给予讨论。(1)基图的选取：初始值工”’的选择对于迭代过程的结果有很大的影响，特别是由于时间和费用的约束，迭代步数可能很受限制。在实践中，由于检测对象的多样性，找寻反映实际情况的多元随机变量的均值是很困难的。因而我们常常用其他方法(诸如卷积法或反投影法)的输出结果作为ART方法的初值x。’，甚至在实践中更是经常用均匀灰度的图像作为初值，可以将所估计的平均密度(甚至是空密度)用于每个象索。此外重建图像的灰度值总在0～255范围内以及线性衰减系数(在任何能量下)总是非负的，我们可以将这两个信息加入到迭代的每一步中去，以加快收敛速度。(2)松弛参数：在ART方法中，我们可利用的基本工具是松弛参数。在前面我们已经提到，如果松弛参数的值选得较低，就会发现使用ART型算法能获得好的重建结果，甚至是实验得到的数据。低松弛参数，能减少方程中不精确性的影响，防止“盐和胡椒”(saltandpepper)的产生，这些东西在使用高松弛参数的ART型重建中是经常看到的。但在某些情况下，有限制地使用高松弛参数是适当的。当求解一个不等式系时，如果仅是一个不等式更动一下，松弛参数取2，那末，迭代过程能够大为缩短，使得在与这个不等式相联系的超平面中，在迭代之前，造成了估计值的镜面反射。有选择性的使用反射，在某些情况下，能够确保有限收敛性。(3)各投影数据之间的关系：在系统中，不等式(或等式)次序的重排对于算41 重庆大学硕士学位论文法的实际性能也可能有很大的影响。利用标准几何结构的数据采集，我们已经发现，依次先处理同一方位图上所有射线的方程，然后处理与前一光源位置夹角为大角度(例如600)方位图上的射线方程，这样比较优越。此外，根据迭代过程本身有作者提出了一些技巧。根据以上分析以上影响迭代过程的主要因素，要加快迭代过程可以从两方面着手进行：一是利用数据约束条件，根据窄角扇束采集数据的特征，可以使用正交优化迭代Bo’或是Schmidt正交化[311方法来加快迭代过程；二是利用数据约束条件，从迭代过程本身着手，使用并行迭代算法”1和分块迭代算法以加快收敛过程。5．3．1正交优化的层次迭代算法我们在分析ART过程的几何意义时就已指出，ART迭代过程实际上是将一估计矢量X，在各个超平面上不停的作投影直至其足够接近实际的图像矢量(见式(5．24))。由此我们可知如果各个超平面之间如果两两正交，则只需一次迭代便可完成迭代过程。由此正交优化的层次迭代算法的一般过程就是：(1)将图像的估计矢量J作为初值，在任意超平丽∥上进行投影；(2)选择一个与∥线性相关度最小的剖面∥进行再投影：(3)在余下的M-2个超平面中再选取一个与^∥线性相关度最小的剖面∥进行再投影，直至使用完全部M个超平面；(4)检验图像是否达到要求，满足则停止迭代；否则，将迭代得到的向量z作为下一次迭代的初值重新进行迭代。＼fXlklf＼黾冯nJ』l／＼I对f＼P1p2p3图5．2迭代法重建过程中射线投影图Fig5．2projectionofradialintheprocessofART 5ART与SIRT算法可见正交优化层次的迭代比一般的顺序迭代只是多了一个选取线性相关度最小的超平面最为下次投影的超平面。在实现中，我们超平面的选取转化为投影分度的选取，即选择两个投影分度使两者的空间夹角尽可能的接近900。5．3．2并行迭代算法在数据采集与处理的章节中，我们已经明确地指出窄角扇束扫描方式下实际上可以看作多探测器的平行束扫描方式。下面我们以图5．2来说明顺序迭代重建图像的过程，并找出重建过程内在的并行性“1。从图中我们可以看出，在某个投影方向做修正时，比如垂直方向，这里共有三根射线(P。，P：，P。)，在顺序迭代时我们总是先将一根射线处理完了后，再处理第二、第三根射线。其实这三根射线的修正互不影响，完全可以同时进行，也即存在并行性。因而我们可以厩时修正这三根射线以加快迭代过程。 6计算机实现与结果讨论6．1系统综述本课题在微机上完成了用代数重建法(ART)和同时迭代法(SIRT)实现图像重建的整个过程，并对二者进行了系统的比较研究，尤其是在不完全数据和数据噪声污染特别严重时重建图像的比较。本系统在W'mdows平台下，借助C++Builder开发完成的。C++Builder””是Inprise公司开发的一款优秀编程开发工具。它是一个面向对象的可视化编程平台，应用于32位Windows应用程序的快速开发，具有良好的集成开发环境(mE)，它提供了所有你需要的用于设计、开发、试验、调试和发布应用程序的工具。C++Builder：[酸大程度体现了C++语言的灵活性和适用性，大大提高了程序的运行效率，可以实现用最小的代码开销编写出高效率的Windows应用程序。本系统提供了良好的人机交互界面，完成了仿真数据采集，噪声模拟，算法选择，图像重建，结果显示与分析的全过程。整个图像重建系统的基本框架如图6．1所示：图6．1图像重建系统基本框图Fig6．1fundamentalSketchmapoflmageRecons咖cdonSystem系统的主要组成部分包括：(1)仿真数据采集(DataAcquisition)：本课题由于实验条件所限，无法取得窄角扇束CT扫描的实测数据，因此投影数据是由采用计算机仿真的方式对测试模型生成的。我们选取256X256的图像作为数据源，窄角扇束的张角为120，探测器数目为128。(2)数据预处理回atapre-Processing)：该部分对模拟采集到的数据进行预处理，在这里我们检验算法的性能，不进行数据校正，而模拟产生噪声源，人为的加入噪声(包含加性噪声和乘性噪声)和进行数据缺损(包含连续性数据缺损和间隔性数据缺损)处理。(3)图像重建(ImageReconstruction)：这部分通过人机对话的方式选择重建参数和希望采用的重建算法，利用采集到的数据和经过处理的数据进行图像重建，以检验算法的性能。(4)重建结果分析(AnalysisofReconstructionResult)：这部分包括了重建图像的显示，查看图像内任意一条线段上的灰度PROFILE曲线，以比较重建图像与源图像之 重庆大学硕士学位论文问的差异，显示重建质量指标。6．2仿真数据采集∞ataAcquisition)6．2．1数据采集过程度仿真实现为了验证我们所提出的算法可行性及计算效率，计算机模拟是一种非常有效的检验手段。根据数据投影的定义，在对计算机测试模型进行数据采集时，首先需要确定方式源、被检物以及各探测器的位置。放射源和探测器两点之间的连线就是穿过测试模型的实现轨迹，找出射线轨迹所经过的各点坐标，并将这些点的灰度值积分，就获得了投影数据n”。在我们的重建实验中，假设在窄角扇束扫描系统中的探测器阵列由M个探测器，并按等角间距方式排列。取被检物每一分度的旋转角度等于射线柬的窄角n，则接个扫描过程物体的旋转分度数为Ns=石／口。在每个旋转分度下源与探测器阵列同步平移K步，这样整个扫描过程中采集的中投影量为MXPXK。一般在实际窄角扇束扫描工业CT中，源—探测器只作平移运动，被测工件作选择运动。为了计算机实现的方便，我们将其转化为源—探测器既作平移运动又作旋转运动，被测工件不进行任何运动的形式。设扫描覆盖区域的直径为D，我们取射线源与探测器间的距离为L=(3／2)D。如图6．2所示，在旋转角度为6时，放射源沿直线平行移动，该直线的方程为(法线式方程)：xsin6+ycosJ—D=0(6．1)为计算机实现方便起见，将上式改写为参数t的方程：Ⅲ一等瞄,inq虬匀,@z，其中p=三一善=D为已知量。如前所述，在任一旋转分度下，射线源与探测器阵列沿直线同步平移需走足步，着平移轨迹的中间点为坐标原点，射线源—探测器阵列的平移步距为d，则显然中间点两侧的相邻平移位置的坐标(矗力应满足：(z—psinJ)2+(Y—pcosJ)2=d2(6．3)将参数方程代入此式可得：r2xct928+t2=d2=>，2(1+c培2J)=d2jr2=d2sin2占(6．4)j，=±dsinJ 6计算机实现与结果讨论∥．．．·--’JLY-·--·。X分一≮义x一图6．2数据采集过程的仿真F培6．2imitationofdataacquisition[习=[一‘挚d。蠹d][：：耋]cs．s，当，=：t2dsin占时阱一d。当d][捌@s，以此类推，我们可以数学归纳法证明，当ra按整数变化时，有下列等式成立：[；]=[一‘苫’d。±三d][：：；]cs．，，上式中Ⅲ表示平移位置的序号，d为探测器阵列的平移步进距离，D为扫描区域的直径。很显然，只需将历从一芝}到兰；逐次递增，就可确定反射源在平移步进过程中，所处的K个位置。，为了便于计算，按角度递增顺序依次将旋转分度编号为j=D，』，2,3⋯，N-1。则式(5．7)中的旋转角度为占=亨a。由此可得，在第_，个旋转分度下，在第m个平移位鼍r。放射源的举标为：47 重庆大学硕士学位论文黜H一≯冀黝]@s，在放射源位置确定后，我们就能确定每个对应位置上探测器的位置。设M个探测器按图6．3等角间距排列，并编号为j=D，』，2,3⋯，34-1，则按第四章中的式(4．6)计算出在第，号旋转分度下，放射源与第i号探测器间连线与坐标轴的夹角西“。经过坐标变换后，可得第j号旋转分度下，对应于射线源的第m次平移，各探测器的坐标为：院帐u’-sin04j,iYT(J)cosL]@。，∽1(别l，pj。7(M／2一玲t凇2})∞D图6．3探测器等角间距排列示意图Fig6．3sketchmapofequal-angledetectorsarray确定了放射源坐标和对应的各探测器的坐标之后，就可以列出射线的直线方程。此外我们还利用放射源坐标与相应的射线与坐标轴的交角求出实现方程。使用该射线方程就可以计算出该射线通过的各像素的线性积分值，就得到了投影矩阵尸，同时可以求出系数矩阵R。6．2．2图像格式在图像重建实验中，用于重建的测试模型都是以位图的形式提供的，重建后的结果也以位图文件的形式进行显示和保存。位图(BMP)‘351是W'mdows下的一种通用图像格式，可以在基于W'mdows的其它有关图像的程序中自由使用，因此使用位图格式开始系统具有良好的通用性；同时由于位图格式是未经任何压缩处理的图像，使用这种格式更方便获得模拟扫描数据。此外，W'mdows提供了一套支持位图图像处理的SDK函数‘36’孔381，从而为软件开发提供了便利。 6计算机实现与结果讨论BMP文件大体上分为四个部分：位图头文件(BrrMAPF几EHEADER)，位图信息头(BITMAPn师OHEADER)，调色板(Palette)，实际的位图数据(ImageData)。第一部分为位图文件头，是一个结构，含有BMP文件的类型、大小和打印格式等信息。这个结构长度是固定的即为14个字节，其定义如下：typedefstructtagBITMAPFILEHEADER{WORDbfType；／／位图文件的类型，必须为0x424D，即字符串“BM”DWORDbfSize；／／位图文件的大小，包含这14个字节wORDbfReservedl；／／位图文件保留字，必须为0WORDbfReserved2；／／位图文件保留字，必须为0DWORDbfOfibits；／／为从文件头到实际的位图数据的偏移字节数，即位图头／／文件，位图信息头和调色板三个部分的长度之和。)BITMAPFILEHEADER；第二部分为位图信息头(BITMAPD盯OH队DER)，它含有位图文件的尺寸和颜色等信息。这个结构长度也是固定的，长度为40个字节，其定义如下：typedefstructtagBITMAPINFOHEADER{DWORDbiSize；／／这个结构的长度，为40LONGbiWidth：／／位图的宽度，以象素为单位LONGbiHeight；／／位图的高度，以象素为单位WORDbiPlanes；／／目标设备的级别，必须为lWORDbiBitCount：／／每个象素所须的位数，必须为l(单色)、／／4(16色)、8(256色)、24(真彩色)之一DWORDbiCompression；／／位图压缩类型，必须是0(不压缩)、／／I(BI压缩类型)或2(BI压缩类型)之一_RLE8RLE8DWORDbiSizelmage；／／位图的大小，以字节为单位LONGbiXPelsPerMeter；∥目标设备水平分辩率，单位是每米象素的个数LONGbiYPelsPerMeter；／，目标设备垂直分辩率，单位是每米象紊的个数DwORDbiClrUsed：／／位图实际使用的颜色数，如果该值为0，则用到的颜／，色数为2大hiBitCount次方DWORDbiClrImportant；脆图显示时被认为重要的颜色数，如果该值为0，则，／认为所以的颜色都重要}BITMAPINFOHEADER：第三部分为调色板(Palette)，当然这里是对那些需要调色板的位图而言的。有 重庆大学硕士学位论文些位图，如真彩色，是不需要调色板的，BITM心矾FOH勘山ER后直接是位图数据。调色板实际上是一个数组，共有biClrUsed个元素(如果该值为0，则用到的颜色数为2的biBitCount次方个元素)。数组中每个元素的类型是一个RGBQUAD结构，占4个字节，RGBQUAD的定义如下：typedefmgRGBQUAD{BYTErgbBlue{／／蓝色的亮度值BYTErgbGreen；／／绿色的亮度值BYTErgbRed；／／红色的亮度值BYTErgbRcsm-wd；／／乇畏留，必须为0}RGBQUAD第四部分就是实际的位图数据(1mageData)。对于一幅灰度图像来说，rgbBlue、rgbGreen和rgbRed均相等。在生成位图文件时，Windows从位图的左下角开始(从左到右从下到上)逐行扫描位图，最后得到的是最上面一行的最右一个元素。，这些象素值的字节组成了位图阵列。在重建实验中我们将采用的是24位位图和256色位图，其中24位位图没有调色板，每个像素对应3个字节，分别代表红、绿、蓝的亮度值。在灰度图像中，每个像素的红、绿、蓝三分量的亮度值相等。256色位图有调色板，它的实际位酎数据就是像素的灰度值。计算投影数据就是计算被各条射线所穿过的像素的灰度和累加，即像素的线性积分。这和射线工业CT中的数据采集有所不同，在工业cT中我们采集到的数据是各条射线穿过工件后剩余的能量值，在图像重建时只需将此值和射线出射能量作差运算即可得到沿射线工件各点的线性积分。投影数据的计算流程图如图6．4。6-3数据预处理∞atapre-Processing)为了模拟真实的采集环境，本课题将对采集到数据进行噪声污染和缺损处理。而在真正的工业CT中，图像重建之前需要对采集的数据进行校正，本文为了检验算法的先进性和可靠性，重建前不对数据进行预先校正，而是对数据进行噪声污染和缺损处理。噪声污染处理主要包括加入乘性噪声和加性噪声，其中乘性噪声主要用来模拟实际数据采集中由高能粒子直线加速器和其他电器设备所产生的电磁干扰噪声；而加性噪声是用来模拟实际数据采集中的高斯白噪声。50 6计算机实现与结果讨论图6．4投影数据的计算流程图F瑭6．4flowch矾ofcalculatingpmjectiondata 重庆大学硕士学位论文数据缺损处理主要包括连续性视角数据缺损和间隔性视角数据缺损。连续性视角数据缺损是指将空间相邻的探测器(比如⋯M．1，M，M+I⋯)的数据丢失；间隔性视角数据缺损是指将空间相隔的探测器(比如⋯M-K，M，M+K⋯，K>I)的数据丢失。进行数据缺损处理主要是检验算法处理不完全数据时的性能。6．4图像重建(ImageReconstruction)6．4．1ARTI型算法图6．5ARTI的计算流程图Fi96．5thecalculatingflowchartofARTI 6计算机实现与结果讨论在此算法中，我们选取图像的初值为T，T是重建图像平均灰度的一个大致估计，这样作可以算法较快的收敛，这也就是我们前面所述的“利用先验知识，可以加快收敛过程”的一个应用。算法中关于‰的计算和计算射线的投影数据是一样的，只不过计算投影数据时，我们计算的是各个像素的灰度和，而计算系数矩阵时我们是将射线所穿过的像素在乍中标记为1，因而计算的复杂度是一样的。从算法的迭代式我们可以看出，算法的图像的修正实际上是按照射线的顺序来进行的，即每次将修正值平均涂抹到被该射线所穿过的各点上，判定条件lIx，一Ⅳ旧’112Times，即判断是否超过指定的迭代次数。ARTI算法的流程图如图6．5。6．4．2ARTII型算法在此算法中，我们选取图像的初值也为T，o<五汁’电，算法中。的计算和ARTI相同，算法的流程图如图6．6。与ARTI型相比，ARTⅡ算法在迭代过程中多了一个收敛参数x和一个迭代判断条件‘j+X耻)<-Pik，只有判断条件成立时，才对图像进行修正，因而ARTlI要比ART1的收敛速度会(当爿”；1时)；当ARTlI收敛参数小于l时，越盯II的收敛速度将小于ARTI，但重建图像的质量要比ARTI好。通过实验我们可咀发现当取0．1<∥”O．8)的图像平滑的多。从Profile曲线对比分析我们还可以看出，当^取值较大时(^>O．5)，重建图像与原始模型的差别主要在于少数点的较大误差，因而图像会出现黑点或白点，也就是“盐和胡椒”现象；而^较小时差别主要是多数点的小误差，因而图像要比^取值较大时平滑的多，但迭代次数较少时，图像比较模糊，且图像的亮度要差一些，我们通过调整灰度等级来提高亮度“““。由此我们可以预见当^取值较小时，经过多次迭代会得到一个较为满意的图像，实验结果也证明了这一点，当0．1<^