生物医学工程毕业论文-高精度X-CT投影成像研究.doc

生物生物医学工程毕业论文医学工程毕业论文题目：高精度X-CT投影成像研究专业:生物医学工程高精度高精度X-CT投影成像研究投影成像研究摘要自从1971年第一台CT设备问世以来，计算机断层成像技术（CT）不断取得巨大进步。CT从理论上讲是一个从投影重建图像的反问题。卷积反投影（CBP）、直接傅立叶重建（DF）以及代数重建算法（ART）同为典型的CT重建算法。其中，DF算法在原理上简单，重建速度较快，但是由于缺乏较好的从极坐标系到直角坐标系的插值方法，使得在大多数情况下DF算法重建质量不如CBP以及ART算法重建质量。而由于CBP算法能够重建出质量足够好的图像，同时其耗费的时间也在可以接受的范围内，因此现代CT中的重建算法几乎都是用CBP算法。论文从DF算法入手，试图在投影数据足够的条件下寻找好的网格化方法，改善DF重建算法的质量。在我们的论文中，研究了DF法中应用三种常见插值器方法时，图像重建效果对比。同时，论文也指出了采用一种插值器能使DF法与CBP法等效。最后我们提出了一种新的DF重建算法。在这种算法中采用了对投影数据进行线性调频z变换近似其频谱数据的方法，使得频谱数据更密，同时在DF法的最后一步中采取线性调频z变换求重建图像，减小了插值误差。论文中把这种方法与CBP重建以及ART、SIRT进行了对比，仿真结果显示，本文中采用的方法重建图像的质量至少不比CBP重建以及ART、SIRT重建图像质量差。关键词：计算机断层成像，直接傅立叶变换重建，卷积反投影，代数重建算法，同时迭代重建法RESEARCHONHIGH-PRECISIONX-CTIMAGERECONSTRUCTIONFROMPROJECTIONSABSTRACTSince1971whenthefirstCTequipmentwasmadeCThasbeenmakinggreatprogressallthetime.IntheoryCTistheinverseproblemofreconstructinganimagefromitsprojectiondata.Theconvolutionbackprojects(CBP)directFourierreconstructs(DF)aswellasalgebrareconstructionalgorithm(ART)arethetypicalCTreconstructionalgorithm.Amongthem，theDFisverysimpleinprinciplebutbecauseoflackingagoodinterpolationinterpolatingdatafromthepolarcoordinatesystemtorectangularcoordinatesystemwhichcausesthattheimagereconstructedbyDFisinferiortotheimagereconstructedbytheCBPandARTinmostcases.WhiletheCBPontheothersidecanyieldtogoodqualityreconstructionsandonlytakeashorttimesotheCBPiswidelyusedinmodernCTequipment.InthispaperwestartwiththeDFalgorithmandtrytoimprovethereconstructionqualityundertheconditionofenoughprojectiondata.WeinvestigatethreetypesofinterpolatorswidelyusedintheDFandconstructtheimagesreconstructedbyusingtheseinterpolators.WealsostudyaninterpolatorwhichcanmaketheDFalgorithmbeequivalenttotheCBPalgorithm.AtlastweprovideanewwaytocompletetheDF.InthisweusetheCZToftheprojectiondatatoapproximateitsfrequencyspectrumdatainwhichwaycanmakethefrequencyspectrumdatabedenser.InthesametimeweusetheChirpztranstoreconstructtheimageinordertoreducetheinterpolationerror.WecomparethiswiththeCBPARTandSIRT.ThesimulationresultsshowsthattheimagereconstructedbythisisatleastnotworsethantheimagesreconstructedbytheCBPARTandSIRT.Keywords:ComputerizedTomographyDirect-FourierConvolutionBackProjectionAlgebraicReconstructionTechniquesSimultaneousIterativeReconstructionTechnique.目录目录第一章绪论.11.1断层成像技术发展简介.11.2CT图像重建.21.3研究目的及论文结构.2第二章从投影重建图像算法概述.32.1投影定理.32.2卷积反投影法2.32.3直接傅里叶重建法.52.4代数重建算法.62.5线性调频z变换（CZT）.8第三章高精度直接傅里叶重建算法研究.93.1常见插值方法.93.2与卷积反投影算法等效的直接傅里叶重建插值方法.123.3一种新的高精度DF重建算法.153.3.1投影数据的获得.153.3.2投影数据一维傅里叶频谱数据的计算.163.3.3频谱数据的网格化.173.3.4由频谱数据重建图像.173.3.5重建图像对比.17第四章仿真结果与算法比较.194.1图像仿真结果.194.2重建后图像与原始图像差异评价.214.3灰度值比较.224.4计算复杂度比较.23第五章结论.25参考文献.26谢辞.27译文及原文.28第1页第一章绪论CT（ComputerizedTomography），即计算机断层成像，是用来获取观测目标断层图像的一门技术。它广泛应用于诊断医学、射电天文学、电子显微和雷达等多科学领域。CT是由在多个观测角度获得的有关目标的一系列投影数据，通过图像重建技术来得到目标的断层图像的。近年来，X-CT无论在基本技术方面，还是在新的临床应用方面都取得了巨大的发展。在CT的各个主要组成部分，如光管、探测器、滑环、数据获取系统和算法等方面，都取得了很大进步。自从螺旋CT和多层CT问世，出现了许多新的临床应用。CT经过三十多年的发展以后，再次成为医学图像领域中最令人兴奋地诊断方法之一。1.1断层成像技术发展简介自从1895年，伦琴在试验阴极射线管时发现了X射线后，X射线便发展为广泛使用的检测工具。用传统的照相方法，三维的人体沿X射线方向被压缩成了两维的图像。体内所有骨骼结构和组织都重叠在一起，使得感兴趣对象的清晰程度大为下降。由于传统照相的这种限制，导致了传统断层成像技术的出现。Bocage是传统断层成像技术的先驱之一。早在1921年，他描述了这样的设备，能使感兴趣平面以外上下的结构模糊得看不清楚。虽然传统断层成像术有了一定的发展，而且这些技术在生成感兴趣平面的图像方面取得某种程度的成功，却都有着最基本的限制，即它们并没有增加物体的反差，也不能从根本上去除焦平面以外的其他结构。加上病人须接受很大的X射线剂量，传统断层成像现在很少应用到临床上了。而从投影重建图像的努力也早在1940年已经开始，再1940年颁布的专利中，GabrielFrank描述了现代断层成像技术的基本思想。1963年，Kuhl和Edward应用放射性同位素提出了横向断层成像方法，该方法被后来进一步发展和改进为今天的发射计算机断层成像（ECT）。从多个投影数据重建图像的数学公式可以追溯到澳大利亚的数学家Radon，他在1917年用数学证明从无限多投影数据可以重建出原来的图像。在医学上，AllanMCormack报道了可能是第一台实际建成的CT扫描机的研究结果。英国EMI的中心研究实验室的GodfreyNHounsfield于1967年开始第一个临床CT扫描机的研制。再进一步改进了数据获取和重建技术后，第一台可供临床应用的CT机安装在Atkinson-Morley医院。虽然第一代扫描机的临床检查结果还可以就接受，但是在长达4分半钟的数据获取过程中病人的运动引起的图像质量问题是严重的。第二代扫描设备仍然是一台平移旋转扫描机，然而由于利用了多个笔形束使得所需旋转的步数减少了。最流行的扫描机类型是第三代扫描机。这种结构将大量探测器布置在以X射线源为中心的圆弧上，探测器的尺寸足够大使得整个检测对象始终落在探测器的视场内。X射线源和探测器在整个设备围绕病人旋转时保持相对静止。直线移动的取消显著减少了数据获取时间。第四代扫描机中的探测器形成一个闭合的圆环，在整个扫描过程中X光管围绕病人旋转，而探测器保持静止。第五代扫描机也称作电子束扫描机，于1980年到1984年建造，用于心脏检查。为了减少心脏跳动形成的伪迹，要求采集一套完整的投影数据必须在2050ms内完成。因此在电子束扫描机中，射线源的旋转由电子束扫描运动代替了传统的X光管的机械运动来完成。自从第一台临床扫描机采用以来，CT技术取得了巨大的进展。CT机已经从第一代发展至第五代，扫描时间也从5分钟缩短至几十毫秒，空间分辨率也达到1mm以下1。第2页CT起源于医学成像，随着CT技术的发展，CT技术的应用也深入到重建图像的科学和技术的其他各个领域。因此，从广义上讲，CT是这样一个过程，它通过某种辐射源（如X射线，放射性核素、超声波、磁场等）对观察目标进行作用并检测投影数据，然后利用计算机完成图像重建、数据处理并和图像显示技术结合，得到观测目标内部的（通常具有某种物理性质的）断层图像或三维立体图像。CT的成像大致过程可用图1-1所示方框图来表示12。图图1-11-1CTCT成像过程流程图成像过程流程图图1-1显示的CT图像成像流程图中，为简便起见，我们把“预处理”模块放到“负对数运算”模块前面，然而，事实上有些预处理步骤在对数步骤之后完成。在后处理模块中，主要包括伪像校正、图像增强、3D以及其他各种重构方法。1.2CT图像重建图像重建技术是CT技术中的一个重要问题。目前已经发现许多种从投影重建图像的方法，这些方法中，往往不同的方法有不同的适应范围，可以在不同的场合下考虑使用不同的重建方法，以获得最好的图像重建质量或者缩短计算时间。所有的图像重建方法一般分为两种类型：变换方法和级数展开法。变换法也称解析法，其基本特点是目标函数与投影函数之间的关系可以用某种解析公式来表达，一旦这种解析公式建立，便寻找一种合适的离散实现这一解析公式的具体算法。变换法中两类最具代表的算法是滤波反投影重建算法和傅立叶重建算法。由于变换法实现比较简单，速度快，这类方法（尤其是卷积反投影法）是商用CT机中普遍采用的方法。而级数法的特点是，从一开始就把连续图像离散化，从而建立起离散重建图像和离散投影之间的代数方程组，然后在某种最优准则的指导下解这个线性方程组。由于这个方程组一般都非常大，很难采取直接的方法求解，因此这种求解过程一般是一个迭代过程。1.3研究目的及论文结构论文中主要研究高精度的直接傅立叶变换重建算法（DF）。在论文中先概述了CBP法、ART法以及DF法，并表明现在的CT机中广泛使用的算法是CBP法。在讨论DF算法中，我们先研究了比较常见的一些插值方法（最近点插值、双线性插值以及混合二阶样条线性插值法），通过仿真表明双线性插值以及混合二阶样条线性插值法要比最近点插值法中取得更好的效果。接着论文证明了可以采用一种插值方法使得DF法与CBP法等效，从而表明有可能找到更好的插值方法使得DF法重建图像质量优于CBP法。最后我们介绍了我们重建图像的方案。经过大量图像仿真的试验，我们采用了增加频谱数据密度的方法，减小网格化误差，取得了较好的实验结果。论文的章节内容具体安排如下：第一章，简单介绍计算机断层成像地发展；第二章，综述了四种常见CT重建算法（CBP、DF、ART、SIRT）以及从投影重建图像的基础知识；第三章，详细介绍了DF算法特别是我们对DF的数字实现方法；第四章，对第三章的算法进行了仿真试验，并与CBP以及ART、SIRT进行了对比；数据采集预处理或标定负对数运算断层图像重建后处理窗口显示第3页第五章，在前面几章的基础上，总结了本文主要内容以及结论。第二章从投影重建图像算法概述从投影重建图像的算法很多。常见的有：反投影重建算法；滤波反投影重建算法；直接傅立叶变换重建算法；迭代重建算法等。在本章，我们将在第一小节介绍投影定理，投影定理是卷积反投影重建算法以及直接傅立叶重建法的基础。然后，为了与我们采用的算法进行对比，我们简单介绍其中四种最常见的重建算法：卷积反投影法（CBP）、直接傅立叶重建法（DF）、代数重建法（ART）、同时迭代重建算法（SIRT）。由于在DF法中，我们需要使用线性调频z变换（CZT），因此在本章最后一小节，我将简单介绍CZT。2.1投影定理投影定理或中心切片定理是滤波反投影重建算法与直接傅立叶变换重建算法的基础。在没有衍射源的情况下，其内容是：某图像f（x，y）在视角为时投影的一维傅立叶变换给出f（x，y）的二维傅()rPx立叶变换的一个切片。切片与轴相交成角，且通过坐标原点2。12()()FF即（2-:121()()()()()rAAFpxPP1）根据投影定理，投影重建图像的方法原则上可按如下流程求解：（1）采集不同时脚下的投影，理论上应为1800范围的连续无穷多投影；（2）求出个投影的一位傅立叶变换，这就是图像二维傅立叶变换过原点的各切片，理论上应为连续的无穷多片；（3）将（2）求出的各切片组成图像的二维傅立叶变换；（4）用（3）得出的图像傅立叶频谱数据求二维傅立叶反变换得到重建图像。在实际的重建过程中，由于反变换时的数学处理不同，又可将变换法分为滤波反投影重建算法和直接傅立叶反变换算法。其中滤波反投影重建算法又可分为卷积反投影法和Radon反变换法。Radon反变换法由于对误差很敏感，因此在商用CT机上未得到应用。从数学表达式看，Radon反变换法与卷积反投影法等效，因此我们没有仿真Radon反变换重建方法，而只使用卷积反投影法重建图像。在下文的讨论中，我们也只介绍卷积反投影法，而不再介绍Radon反变换法。2.2卷积反投影法2由投影定理，可以通过在不同视角下的投影的一维傅立叶变:()A()axy()rPx换求得，因此待见图像：第4页:12()1212121222cos()2cos()12002cos()2co1()()()()4=()()=()()ixyiriririraraxyFAAeddAeddPeddPedddPes()00.(2-2)d从公式（2-2）的第二个积分，即（2-2cos()ired3）该式可以改写成空域变量xr的傅里叶反变换式：2cos()2cos()cos()()()()cos()(2-4)rririrxrrrxredPedhxpxgr式中（2-()()()rrrgxpxhx5）而，。将式（2-4）代入式（2-2）后，得1()()rhxF11()()rpxFP（2-0()cos()argrd6）式（2-6）的物理意义是：经过给定点的所有滤波后的射线投影在范围内累()r0:加反投影重建，得出点的像素值。()r由式（2-6）离散化可以得出滤波反投影重建算法数字实现的具体步骤：（1）对采集到的离散数据进行离散滤波；:()()()()()()()()NNrrrrrnNnNpxmpndmhxndpndmxndhxpndmhx（2-7）式中为滤波函数。为了更清楚的对的选择进行讨论，我们转到频域去讨论。()rhx()rhx设为的傅里叶变换，则理想的滤波函数，但在种滤波函数是不可()H()rhx()H能实现的。在实际中，的高频分量幅度很小，因此近似为带限信号。（事实上，()rpx()P在物体尺寸有限的情况下，投影数据分布在有限的范围内，因此严格的说，其频带是无限的。然而如果物体的密度在空间范围内变化平稳，则高频分量幅度确实不大，并且检测器在接收X射线时，有平均作用，相当于低通滤波，可以认为高频分量很小。）利用这()P第5页一事实，我们可以把强制为带限，即可以对理想的滤波函数加一个窗函()H()H数。常用的窗函数有矩形窗、汉明窗等。在我们的仿真试验中采用简单的矩形窗，也就是使用R-L滤波函数。R-L滤波函数的系统函数为：()RLH（2-()(2)RLHrectB8）式中112(2)0BdrectB；其他。相应的冲激响应函数为：()RLrhx（2-2222()2sinc(2)sinc()rBixRLrrrBhxedBxBBxB9）同时，由于投影数据的采集是离散的，因此，滤波函数也应相应的取其离散形式。对进行采样，其中d为采样间隔对应的最高不失真空间频率为，以()RLrhx1(2)Bd代入式（2-9），将得到的离散形式如（2-10）所示：rxnd()RLrhx（2-2222104()01RLndhndnnnd为偶数为奇数10）如果对进行线性内插，则得到另一种连续的空域函数，可以把()RLhnd()RLrhx看作是的一次近似。也可以对进行精确的sinc函数内插，得()RLrhx()RLrhx()RLhnd到，但是由于运算复杂，在实践中一般不使用，在我们的仿真试验中也采取对()RLrhx进行线性内插重建图像。另外，由于R-L滤波器使用理想的矩形窗，导致重建图()RLhnd像有Gibbs现象，表现为明显的震荡响应。但是由于该滤波器形式简单，因此在我们的仿真试验中，利用CBP重建图像都是采取R-L滤波器。（2）式中，一般此值不一定为的整数倍，即（cossinrxxmymdrxnd为整数），。因此n()rxndnnd()()rhxndhnnd。这样如何求任意的时的:()()rPxmPndm()()pndmhndrx便成为问题的关键。最常见的方法是用内插去求解。插值的过程可以在滤波以:()rPxm前，利用离散的，应用内插得到在轴上连续取值的，就可使用式（2-()RLhndrx()rhx7）求解。也可以在滤波后求得后，在通过内插的方法求得:()rPxm:()Pndm。两种方法的效果一样。最近点内插和线性内插是最常见的插值方法，其具体:()rPxm过程可以参见3.1节。在我们的仿真实验中采取的插值方法是线性插值法，虽然线性插值法是不精确的内插方法，然而它的运算方法颇简单。由步骤（1）所介绍精确的内插方法是函数，运算要比线性插值法复杂得多，因此我们没有采用。sinc（3）设重建图像坐标为，对每个待重建点进行反投影迭加有()iixy（2-:10()cos()sin()NiiijnfxypxnynN第6页11）2.3直接傅里叶重建法直接傅里叶重建算法思想更为直接，即根据投影定理有投影数据求出原图像二维傅里叶变换的每一个切片，将所有切片组合起来即可得到原图像的二维傅里叶变换，在对其进行二维傅里叶反变换就可以重建出原图像。由上面描述可以得出直接傅里叶重建算法的数学描述如下3：（1）假设投影数据为其对t的傅里叶变换为，则可以按式()pt()Pr()Pr（2-12）求得：（2-2()()irtPrptedt12）（2）假设原图为，其二维傅里叶变换为，根据投影定理可以按式（2-()fxy()Fuv13）求得:()Fuv(2-13)()()FuvPr（3）对进行二维傅里叶反变换求得()Fuv()fxy（2-2()2()()()()iuxvyiuxvyfxyFuvedudvPredudv14）式（2-14）就是直接傅里叶重建法的计算公式。但在实际的实现中，由于采集数据的天然离散性，需要对式（2-14）进行离散实现，得到直接傅里叶重建算法的数字实现步骤如下：（1）对投影数据进行傅里叶变换，对式（2-12）离散化得到：()pntm（2-121212()()MintndnMPndmtpntme15）（2）式（2-15）给出了原图呈离散极坐标状态的原图二维傅里叶频谱数据，假设原图在直角坐标系下的二维傅里叶频谱数据为，则可以通过对()Fsurv()Fsurv的插值求得。()PNdm（3）对式（2-14）进行离散化，则可以得到重建图像的离散公式如下：（2-2()()()KKjskuxrlvysKrKfkxlyuvFsurve16）2.4代数重建算法卷积反投影算法与直接傅立叶算法都是变换法重建图像。变化法的突出优点是实现简第7页单，速度快，对足够精度的投影数据能获得很好的重建质量。除了变换法以外，另一类完全不同的重建方法是迭代重建算法5。迭代重建算法是基于模型的重建算法。它能和特定的成像设备及数据采集物理过程结合起来，并能利用某些先验信息，因此即使对于比较粗糙的投影数据，这类方法也能给出较好的重建质量。代数重建算法是迭代重建算法的一种。迭代算法已经被广泛研究了很多年。这种方法假设图像有一组未知的阵列组成，然后得到一组关于此阵列的等式。虽然这种方法在概念上比前面的变换法要简单，但是在医学应用上，这些方法往往不能重建出质量足够好的图像，并且耗费的时间太多。这些算法的实施方法如下：用一个二维矢量表示一个物体，用p表示它的投影。两个变量被一个系统矩阵A和一个误差矢量e连接起来：p=A+e(2-17)重建过程是根据一个测量矢量p估计。估计是通过要求和e满足测定的最优准则来进行。估计的过程是一个迭代的步骤，也就是说，我们产生这样一个矢量序列，(0)(1)(n)以使序列收敛于这里是基于p的最优估计。对于每次迭代j，我们计算p(j):p(j)=A(j)+e（2-18）根据计算的投影p(j)与测量投影p之差，修正估计p(j)这样两个之差减小。通常，在估计过程中在上加一定的限制条件。最常使用的一个限制条件是非负性1。文献中有多种形式的ART（AlgebraicReconstructionTechniques），由于各种ART算法本质都相似，我在算法比较中只采取最常见的一种ART加型ART（ARTII）。（乘型ART与加型ART类似，只是调整方法有所不同。）ARTII满足线性不等式pA（2-19）其中p为投影矢量，A为系统矩阵，由图像形状以及投影位置共同决定。设(0)为迭代初始值，则迭代步骤按式（2-20）进行。（2-(1)2kkkkiiipAAA20）上述校正过程是逐条射线依次进行的，故又称为逐线校正。文献表明，ART在MN时，按照式（2-20）进行迭代并不能保证收敛。其中M为射线条数，N为图像网格数目。当MN时，按照式（2-20）进行迭代将收敛于最小方差准则估计的值5。在式（2-19）中，迭代初始值(0)在理论上可以选择任意值，但是如果我们根据经验，选择某些初始值，可以加速式（2-20）的收敛。加权系数一般选择小于1。而且根据经验表明，当选择较小（0.0250.25）时，根据式（2-20）进行重建，经过3I8I（I为射线数目）次迭代后，结果相当令人满意。加型ART解的形式为：。而乘型ART解的形式为(1)kkk，它的迭代步骤如式（2-21）所示2：(1)kkk（2-(1)iAkkkkipA21）其中初始值一般选取全1矢量（尽量不要取0，否则在迭代时出现除0的情况）。0第8页可以证明按式（2-21）迭代求解，得到的序列收敛于p=A的最大熵解。在我们的仿真实验中，只采用加型ART与DF法对比，而不再采用乘型ART法与DF法对比。SIRT（SimultaneousIterativeReconstructionTechnique）同时迭代重建法是与ART并行的另一种迭代重建算法。这种方法与ART不同的是，它试图通过减慢迭代收敛的速度，来取得比ART法更好的重建图像。SIRT与ART相比，最大的优势在于对测量误差不敏感。虽然在SIRT中，也通过计算其投影误差值，但是与ART不同，在SIRT中，不是每一根射线误差都会调整射线通过网格点的像素值，而是把所有的射线投影误差值计算之后，再去修改迭代网格的像素值。正因为它不是逐个射线调整，而是逐点调整，使得SIRT对单个射线误差不敏感2。SIRT迭代算法的过程与使结果满足最小二乘方准则有关。即SIRT努力使解满足满足最小二乘方准则，也就是说使得式（2-22）有最小值。（2-1()()()TpApA22）一般采用迭代算法去求解，具体的迭代方法如式（2-23）所示：（2-(0)(1)()()()TkkTkApApA23）式（2-23）的意义是：取测量值p作为初始图像，在求第k+1次迭代求(k+1)时，利用第k次的估计值(k)加上校正图像。校正图像正比于第k次估计的误差矢量反投影。因而每个像素的校正值实在是通过该像素的所有射线和的误差值的累加，()()TkApA而不是只与一条射线有关。由于每一像素的校正值是所有通过该像素射线的共同贡献，一些随机误差被平均掉了。为了便于迭代计算，一般将式（2-23）改写为式（2-24）的形式：（2-(0)(1)()()1()TIkkTkiiiiAppaa24）其中I为射线条数，在式（2-24）中很容易看出表示经过j号像素()1()ITkiiijipaa的所有射线的投影误差对该像素值的校正。正因为如此，SIRT的校正过程被称为逐点校正。其中初始值选择，是因为这样选择能够加快SIRT的收敛速度。其中(0)TAp为图像的反投影。反投影重建图像会引起星状误差，在CBP算法中，图像的星状伪TAp迹通过滤波除去；在SIRT算法中，图像的星状伪迹由迭代调整消除。由于在CBP算法中调整方法更加直接，因此其重建速度也更快。2.5线性调频z变换（CZT）离散傅里叶变换（DFT）可以看作是点的序列的变换在单位圆上等间隔取样N()xnz点。同时DFT也可以看作是序列的傅里叶变换的等间隔取样点。但是在实际中N()xnN也许：不需要计算整个单位圆上变换的取样，如对于窄带信号，只需要对信号所在的一z段频带进行分析，这时希望频谱的采样集中在这一频带内，以获得较高的分辨率，而频带以外的部分可不考虑；第9页或者对其它围线上的变换取样感兴趣，例如语音信号处理中，需要知道变换的zz极点所在频率，如极点位置离单位圆较远，则其单位圆上的频谱就很平滑，这时很难从中识别出极点所在的频率，如果采样不是沿单位圆而是沿一条接近这些极点的弧线进行，则在极点所在频率上的频谱将出现明显的尖峰，由此可较准确地测定极点频率。CZT在一定程度上解决了上述问题。序列的变换公式如式（2-25）所示6：()xnz（2-10()()NnkknXzxnz01.1kL25）取，则式（2-25）可以化为：0000()iikkzreRe01.1kL（2-001000()()()()NiinnkknXzxnreRe01.1kL26）公式（2-26）就是线性调频变换的公式。z由公式（2-26）可以看出CZT的一些性质：当时，CZT变换路径为单位圆上一段圆弧，但是变换之后的点数不一定001rRL等于原序列点数；N当时，CZT变换起始点在单位圆外，当时，CZT变换起始点在单位圆内；01r01r当时，CZT变换向外螺旋扩展，当时，CZT变换向内螺旋扩展。01R01R在我们的试验中，希望得到的是信号的频谱分析，故应在单位圆上去实现CZT，并且由于我们希望得到信号傅立叶频谱比较密集的采样，以减小网格化的误差，因此在我们实现的DF重建算法中，取，且，在本文的后面几章将看到采取CZT近001rRLN:似信号傅立叶变换重建图像取得了较好的效果。第三章高精度直接傅里叶重建算法研究按照直接傅里叶重建算法的思想，必须对频谱数据进行从极坐标空间到直角坐标空间的插值，这主要是因为两个方面的原因：第一，为了应用FFT求，因为二维FFT()fxy是基于矩形网格的，目前还没有有效的基于极网格下的二维FFT算法存在；第二，图像的显示大多是基于直角坐标系的矩形网格描述。因此在我们的仿真实验中按照公式（2-16）对图像进行重建，关键是对的求解，由公式（2-15）得到的傅里叶频谱数据()Fsurv呈极坐标形式，因而直接傅里叶重建算法需要对对频谱数据从极坐标向直()PNdm角坐标的插值。如图3-1所示，显示的是极坐标与直角坐标采样的不同特点。二维频域中的插值不像真实空间中的插值一样直接。在真实空间中，一个插值误差局限于像素所在的小区域。然而，对于频域插值，这个特性并不正确，因为二维傅里叶空间中每个采样表示某一个空间频率。于是，在傅里叶空间中一个单独采样点的误差会影响整个图像（经过傅里叶反变换后）的外貌，这也是使得大多数DF重建方法很难得到与CBP法相同重建质量第10页的主要原因。因此，利用直接傅里叶重建算法要产生高精度的重建图像，需要在傅里叶空间内非常精确的插值，插值技术也成为了DF重建算法的关键技术。（a）极坐标采样示意图（b）直角坐标采样示意图图图3-13-1极坐标到直角坐标插值示意图极坐标到直角坐标插值示意图3.1常见插值方法在重建算法中最为常见的插值方法有：最近点插值法、双线性插值法以及三次样条插值法。其中最近点插值法最简单，但是误差也最大，三次样条插值法重建得到的图像质量最好。这几种插值方法均属于所谓的常态插值法，即插值后的函数在取样点的取值保值不变。另一种类型的插值法被称作网格化过程，它是基于采样定理导出的，并且不具备常态插值的特点。这种网格化过程的关键是选择卷积函数。总的来讲，已经发现了一些精心设计的插值函数，可以使得才有DF法重建图像质量不亚于采用CBP法重建图像质量，但是这些插值函数往往较复杂，在本文中不予讨论。在我们的仿真试验中，采用的是最简单的最近的插值法，实现高精度重建图像。为了比较应用最近的插值法、双线性插值法以及三次样条插值法的直接傅立叶变换法重建图像质量，我们在这里将简要介绍这三种插值方法，并且给出了仿真图像的对比。图图3-23-2极坐标空间网格化示意图极坐标空间网格化示意图如图3-2所示，当采用最近点插值法时，在直角坐标系的P点的值，将采用离P点最近的极坐标系中点的值替代。也就是采用使取值的那个点作minmindPAPBPCPD为P点值的估计。这种方法虽然最简单，但是误差较大，只有在A，B，C，D点的距离足第11页够小时，才能实现比较精确的插值。在我们的试验中，正好满足这一条件，所以在后文可以看到我们采用最近点插值方法实现了较高精度的图像重建。在双线性插值法中，先需要估计E、F点的值，再求P点的值。其中E，F点可采用式（3-1）估计。（3-BAEEAFEBFFABCDFFDFFCFFCD1）然后采用式（3-2）计算的值作为估计点P的值。（3-:EFPPFFPEFFEF2）在大多数情况下，利用双线性插值法估计得到的值比利用最近值插值法估计的值要更准确，但是其算法也更复杂。在我们的试验中由于数据足够密集，因此没有采用双线性插值法，而是采用最近点插值法，也能取得很好的结果。在这三种插值法中，重建图像质量最好的插值方法是三次样条插值法，但是三次样条插值法计算量较大。由于重建图像的质量对与径向方向的插值精确性敏感性远大于对与相角方向的插值精确性。因此这里介绍的是混合二阶样条线性插值法7。混合二阶样条线性插值法的计算步骤为：首先沿径向方向进行三次样条曲线插值，其次沿相角方向进行线性插值即可。如图3-2所示，设、分别为A、B、C、D沿径向方向的偏导ApBpCpDp数，则沿径向方向的点E、F估计值为：23(3()2)+(2()(3-3)EAABABABABAFFpEAFFppEAFFppEA23(3()2)+(2()(3-4)FCCDCDCDCDCFFpEAFFppEAFFppEA然后利用式（3-2）沿相角方向进行线性插值得到P点的值。选取Shepp-Logan头部模型作为仿真原模型，投影数据为367360，重建图像大小为256256。选用以上三种插值算法的DF重建图像以及局部放大图像请参见图3-3。其中图（a）、（b）、（c）分别是应用最近点插值、双线性插值以及混合二阶样条线性插值法重建的图像，图像（d）、（e）、（f）分别是对应的三种插值法重建图像的局部放大图。有图3-3的对比可以看出，由最近点插值得到的图像，质量最差，边缘的震荡最明显以及图像的混淆误差都明显比用双线性插值和混合二阶样条线性插值法重建图像误差大。而应用双线性插值法和混合二阶样条线性插值法重建的图像虽然在视觉上难以区分其图像质量，然而应用一些客观评价标准能够判别图像质量好坏。（理论上，三次样条插值法在多数情况下都是一种比较好的插值方法，尤其对平滑的曲线插值8，因此，可以假定由混合二阶样条插值的DF法重建图像的质量应优于用线性插值法重建图像的质量，然而在这里，由于频谱数据较为充足，因此没能显示其重建图像的质量优势。）在这里我们采用第4.3节中介绍的几种评价图像质量的客观标准作为评判这三种方法重建图像质量的标准。表格3-1显示了三种插值法重建图像的这几种判据大小。对于归一化平均绝对距离判据，双线性插值法在三种方法中表现得最出色。在我们的其他试验仿真中表明，当投影射线条数较少时，二阶样条线性插值法在三种插值法中表现得最好。由于我们的论文主要考虑在投影数据充足时的重建效果，因此，稀疏投影数据的仿真图没有收到本论文中。第12页表表3-13-1重建图像质量评价参数对比重建图像质量评价参数对比插值方法最近点插值双线性插值二阶样条线性插值法归一化均方距离0.27460.23910.2398归一化平均绝对距离0.0520.03910.0395最坏情况距离1477（a）最近点插值（b）双线性插值（c）混合二阶样条插值（d）最近点插值局部放大第13页（e）双线性插值局部放大（f）混合二阶样条插值局部放大图图3-33-3三种插值法重建图像以及其局部放大图三种插值法重建图像以及其局部放大图3.2与卷积反投影算法等效的直接傅里叶重建插值方法由于卷积反投影法不但能重建出较好的图像质量，并且它的计算复杂度也在能够接受的范围内，因此几乎所有现代的CT机中使用的重建算法都是卷积反投影法。虽然DF重建算法能够提供较好的计算效率，但是DF算法重建图像质量在大多数情况下都比CBP差。就像3.1节中讨论的几种插值方法，其重建图像都比CBP差。因此，可能有人会怀疑CBP算法会比所有的DF重建算法要好，但是在HyeokhoChoi和DavidC.MunsonJr.的论文中，提供了一种插值方法使得DF算法与CBP等效，因此可以把CBP看作DF算法的一种有效实现。在这里，我将简单讨论对比这两种重建方法9。在平行束CT中，图像的投影是从各种不同方向采集到的。如果角度采集足够密集，就能忽略图像混叠。因此在CBP法中，可以看作是极坐标下的二维傅里叶反变换，如（3-5）式所示：（3-5）200()()exp(cossin)PfxyFjxydd我们假设f(