生物医学工程毕业论文

《生物医学工程毕业论文》由会员分享，可在线阅读，更多相关《生物医学工程毕业论文（31页珍藏版）》请在装配图网上搜索。

1、 生物医学工程毕业论文题目：高精度X-CT投影成像研究专业: 生物医学工程 高精度X-CT投影成像研究摘要自从1971年第一台CT设备问世以来，计算机断层成像技术（CT）不断取得巨大进步。CT从理论上讲是一个从投影重建图像的反问题。卷积反投影（CBP）、直接傅立叶重建（DF）以及代数重建算法（ART）同为典型的CT重建算法。其中，DF算法在原理上简单，重建速度较快，但是由于缺乏较好的从极坐标系到直角坐标系的插值方法，使得在大多数情况下DF算法重建质量不如CBP以及ART算法重建质量。而由于CBP算法能够重建出质量足够好的图像，同时其耗费的时间也在可以接受的范围内，因此现代CT中的重建算法几乎都

2、是用CBP算法。论文从DF算法入手，试图在投影数据足够的条件下寻找好的网格化方法，改善DF重建算法的质量。在我们的论文中，研究了DF法中应用三种常见插值器方法时，图像重建效果对比。同时，论文也指出了采用一种插值器能使DF法与CBP法等效。最后我们提出了一种新的DF重建算法。在这种算法中采用了对投影数据进行线性调频z变换近似其频谱数据的方法，使得频谱数据更密，同时在DF法的最后一步中采取线性调频z变换求重建图像，减小了插值误差。论文中把这种方法与CBP重建以及ART、SIRT进行了对比，仿真结果显示，本文中采用的方法重建图像的质量至少不比CBP重建以及ART、SIRT重建图像质量差。关键词：计算

3、机断层成像，直接傅立叶变换重建，卷积反投影，代数重建算法 ，同时迭代重建法RESEARCH ON HIGH-PRECISION X-CT IMAGE RECONSTRUCTION FROM PROJECTIONSABSTRACTSince 1971 when the first CT equipment was made, CT has been making great progress all the time. In theory, CT is the inverse problem of reconstructing an image from its projection data.

4、 The convolution back projects (CBP), direct Fourier reconstructs (DF) as well as algebra reconstruction algorithm (ART) are the typical CT reconstruction algorithm. Among them，the DF method is very simple in principle, but because of lacking a good interpolation method interpolating data from the p

5、olar coordinate system to rectangular coordinate system which causes that the image reconstructed by DF method is inferior to the image reconstructed by the CBP and ART method in most cases. While the CBP method on the other side can yield to good quality reconstructions and only take a short time ,

6、 so the CBP method is widely used in modern CT equipment. In this paper, we start with the DF algorithm and try to improve the reconstruction quality under the condition of enough projection data. We investigate three types of interpolators widely used in the DF method and construct the images recon

7、structed by using these interpolators. We also study an interpolator which can make the DF algorithm be equivalent to the CBP algorithm. At last, we provide a new way to complete the DF method. In this method, we use the CZT of the projection data to approximate its frequency spectrum data, in which

8、 way can make the frequency spectrum data be denser. In the same time, we use the Chirp z transform to reconstruct the image in order to reduce the interpolation error. We compare this method with the CBP, ART and SIRT. The simulation results shows that the image reconstructed by this method is at l

9、east not worse than the images reconstructed by the CBP, ART and SIRT.Key words: Computerized Tomography, Direct-Fourier, Convolution Back Projection Method, Algebraic Reconstruction Techniques, Simultaneous Iterative Reconstruction Technique.目录第一章 绪论11.1 断层成像技术发展简介11.2 CT图像重建21.3 研究目的及论文结构2第二章 从投影重

10、建图像算法概述32.1 投影定理32.2 卷积反投影法232.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 计算复杂度

11、比较23第五章 结论25参考文献26谢辞27译文及原文28第一章 绪论CT（Computerized Tomography），即计算机断层成像，是用来获取观测目标断层图像的一门技术。它广泛应用于诊断医学、射电天文学、电子显微和雷达等多科学领域。CT是由在多个观测角度获得的有关目标的一系列投影数据，通过图像重建技术来得到目标的断层图像的。近年来，X- CT无论在基本技术方面，还是在新的临床应用方面都取得了巨大的发展。在CT的各个主要组成部分，如光管、探测器、滑环、数据获取系统和算法等方面，都取得了很大进步。自从螺旋CT和多层CT问世，出现了许多新的临床应用。CT经过三十多年的发展以后，再次成为医

12、学图像领域中最令人兴奋地诊断方法之一。1.1 断层成像技术发展简介自从1895年，伦琴在试验阴极射线管时发现了X射线后，X射线便发展为广泛使用的检测工具。用传统的照相方法，三维的人体沿X射线方向被压缩成了两维的图像。体内所有骨骼结构和组织都重叠在一起，使得感兴趣对象的清晰程度大为下降。由于传统照相的这种限制，导致了传统断层成像技术的出现。Bocage是传统断层成像技术的先驱之一。早在1921年，他描述了这样的设备，能使感兴趣平面以外上下的结构模糊得看不清楚。虽然传统断层成像术有了一定的发展，而且这些技术在生成感兴趣平面的图像方面取得某种程度的成功，却都有着最基本的限制，即它们并没有增加物体的反

13、差，也不能从根本上去除焦平面以外的其他结构。加上病人须接受很大的X射线剂量，传统断层成像现在很少应用到临床上了。而从投影重建图像的努力也早在1940年已经开始，再1940年颁布的专利中，Gabriel Frank描述了现代断层成像技术的基本思想。1963年，Kuhl和Edward应用放射性同位素提出了横向断层成像方法，该方法被后来进一步发展和改进为今天的发射计算机断层成像（ECT）。从多个投影数据重建图像的数学公式可以追溯到澳大利亚的数学家Radon，他在1917年用数学证明从无限多投影数据可以重建出原来的图像。在医学上，Allan MCormack报道了可能是第一台实际建成的CT扫描机的研究

14、结果。英国EMI的中心研究实验室的Godfrey NHounsfield于1967年开始第一个临床CT扫描机的研制。再进一步改进了数据获取和重建技术后，第一台可供临床应用的CT机安装在Atkinson-Morley医院。虽然第一代扫描机的临床检查结果还可以就接受，但是在长达4分半钟的数据获取过程中病人的运动引起的图像质量问题是严重的。第二代扫描设备仍然是一台平移旋转扫描机，然而由于利用了多个笔形束使得所需旋转的步数减少了。最流行的扫描机类型是第三代扫描机。这种结构将大量探测器布置在以X射线源为中心的圆弧上，探测器的尺寸足够大使得整个检测对象始终落在探测器的视场内。X射线源和探测器在整个设备围绕

15、病人旋转时保持相对静止。直线移动的取消显著减少了数据获取时间。第四代扫描机中的探测器形成一个闭合的圆环，在整个扫描过程中X光管围绕病人旋转，而探测器保持静止。第五代扫描机也称作电子束扫描机，于1980年到1984年建造，用于心脏检查。为了减少心脏跳动形成的伪迹，要求采集一套完整的投影数据必须在2050ms内完成。因此在电子束扫描机中，射线源的旋转由电子束扫描运动代替了传统的X光管的机械运动来完成。自从第一台临床扫描机采用以来，CT技术取得了巨大的进展。CT机已经从第一代发展至第五代，扫描时间也从5分钟缩短至几十毫秒，空间分辨率也达到1mm以下1。CT起源于医学成像，随着CT技术的发展，CT技术

16、的应用也深入到重建图像的科学和技术的其他各个领域。因此，从广义上讲，CT 是这样一个过程，它通过某种辐射源（如X射线，放射性核素、超声波、磁场等）对观察目标进行作用并检测投影数据，然后利用计算机完成图像重建、数据处理并和图像显示技术结合，得到观测目标内部的（通常具有某种物理性质的）断层图像或三维立体图像。CT的成像大致过程可用图1-1所示方框图来表示1,2。负对数运算预处理或标定数据采集窗口显示后处理断层图像重建图1-1 CT成像过程流程图图1-1显示的CT图像成像流程图中，为简便起见，我们把“预处理”模块放到“负对数运算”模块前面，然而，事实上有些预处理步骤在对数步骤之后完成。在后处理模块中

17、，主要包括伪像校正、图像增强、3D以及其他各种重构方法。1.2 CT图像重建图像重建技术是CT技术中的一个重要问题。目前已经发现许多种从投影重建图像的方法，这些方法中，往往不同的方法有不同的适应范围，可以在不同的场合下考虑使用不同的重建方法，以获得最好的图像重建质量或者缩短计算时间。所有的图像重建方法一般分为两种类型：变换方法和级数展开法。变换法也称解析法，其基本特点是目标函数与投影函数之间的关系可以用某种解析公式来表达，一旦这种解析公式建立，便寻找一种合适的离散实现这一解析公式的具体算法。变换法中两类最具代表的算法是滤波反投影重建算法和傅立叶重建算法。由于变换法实现比较简单，速度快，这类方法

18、（尤其是卷积反投影法）是商用CT机中普遍采用的方法。而级数法的特点是，从一开始就把连续图像离散化，从而建立起离散重建图像和离散投影之间的代数方程组，然后在某种最优准则的指导下解这个线性方程组。由于这个方程组一般都非常大，很难采取直接的方法求解，因此这种求解过程一般是一个迭代过程。1.3 研究目的及论文结构论文中主要研究高精度的直接傅立叶变换重建算法（DF）。在论文中先概述了CBP法、ART法以及DF法，并表明现在的CT机中广泛使用的算法是CBP法。在讨论DF算法中，我们先研究了比较常见的一些插值方法（最近点插值、双线性插值以及混合二阶样条线性插值法），通过仿真表明双线性插值以及混合二阶样条线性

19、插值法要比最近点插值法中取得更好的效果。接着论文证明了可以采用一种插值方法使得DF法与CBP法等效，从而表明有可能找到更好的插值方法使得DF法重建图像质量优于CBP法。最后我们介绍了我们重建图像的方案。经过大量图像仿真的试验，我们采用了增加频谱数据密度的方法，减小网格化误差，取得了较好的实验结果。论文的章节内容具体安排如下：第一章，简单介绍计算机断层成像地发展；第二章，综述了四种常见CT重建算法（CBP、DF、ART、SIRT）以及从投影重建图像的基础知识；第三章，详细介绍了DF算法,特别是我们对DF的数字实现方法；第四章，对第三章的算法进行了仿真试验，并与CBP以及ART、SIRT进行了对比

20、；第五章，在前面几章的基础上，总结了本文主要内容以及结论。第二章 从投影重建图像算法概述从投影重建图像的算法很多。常见的有：反投影重建算法；滤波反投影重建算法；直接傅立叶变换重建算法；迭代重建算法等。在本章，我们将在第一小节介绍投影定理，投影定理是卷积反投影重建算法以及直接傅立叶重建法的基础。然后，为了与我们采用的算法进行对比，我们简单介绍其中四种最常见的重建算法：卷积反投影法（CBP）、直接傅立叶重建法（DF）、代数重建法（ART）、同时迭代重建算法（SIRT）。由于在DF法中，我们需要使用线性调频z变换（CZT），因此在本章最后一小节，我将简单介绍CZT。2.1 投影定理投影定理或中心切片

21、定理是滤波反投影重建算法与直接傅立叶变换重建算法的基础。在没有衍射源的情况下，其内容是：某图像f（x，y）在视角为时投影的一维傅立叶变换给出f（x，y）的二维傅立叶变换的一个切片。切片与轴相交成角，且通过坐标原点2。即 （2-1）根据投影定理，投影重建图像的方法原则上可按如下流程求解：（1） 采集不同时脚下的投影，理论上应为1800范围的连续无穷多投影；（2） 求出个投影的一位傅立叶变换，这就是图像二维傅立叶变换过原点的各切片，理论上应为连续的无穷多片；（3） 将（2）求出的各切片组成图像的二维傅立叶变换；（4） 用（3）得出的图像傅立叶频谱数据求二维傅立叶反变换得到重建图像。在实际的重建过程

22、中，由于反变换时的数学处理不同，又可将变换法分为滤波反投影重建算法和直接傅立叶反变换算法。其中滤波反投影重建算法又可分为卷积反投影法和Radon反变换法。Radon反变换法由于对误差很敏感，因此在商用CT机上未得到应用。从数学表达式看，Radon反变换法与卷积反投影法等效，因此我们没有仿真Radon反变换重建方法，而只使用卷积反投影法重建图像。在下文的讨论中，我们也只介绍卷积反投影法，而不再介绍Radon反变换法。2.2 卷积反投影法2由投影定理，可以通过在不同视角下的投影的一维傅立叶变换求得，因此待见图像：从公式（2-2）的第二个积分，即 （2-3）该式可以改写成空域变量xr的傅里叶反变换式

23、：式中 （2-5）而，。将式（2-4）代入式（2-2）后，得 （2-6）式（2-6）的物理意义是：经过给定点的所有滤波后的射线投影在范围内累加反投影重建，得出点的像素值。由式（2-6）离散化可以得出滤波反投影重建算法数字实现的具体步骤：（1）对采集到的离散数据进行离散滤波； （2-7）式中为滤波函数。为了更清楚的对的选择进行讨论，我们转到频域去讨论。设为的傅里叶变换，则理想的滤波函数，但在种滤波函数是不可能实现的。在实际中，的高频分量幅度很小，因此近似为带限信号。（事实上，在物体尺寸有限的情况下，投影数据分布在有限的范围内，因此严格的说，其频带是无限的。然而如果物体的密度在空间范围内变化平稳，

24、则高频分量幅度确实不大，并且检测器在接收X射线时，有平均作用，相当于低通滤波，可以认为高频分量很小。）利用这一事实，我们可以把强制为带限，即可以对理想的滤波函数加一个窗函数。常用的窗函数有矩形窗、汉明窗等。在我们的仿真试验中采用简单的矩形窗，也就是使用R-L滤波函数。R-L滤波函数的系统函数为： （2-8）式中相应的冲激响应函数为： （2-9）同时，由于投影数据的采集是离散的，因此，滤波函数也应相应的取其离散形式。对进行采样，其中d为采样间隔对应的最高不失真空间频率为，以代入式（2-9），将得到的离散形式如（2-10）所示： （2-10）如果对进行线性内插，则得到另一种连续的空域函数，可以把看

25、作是的一次近似。也可以对进行精确的sinc函数内插，得到，但是由于运算复杂，在实践中一般不使用，在我们的仿真试验中也采取对进行线性内插重建图像。另外，由于R-L滤波器使用理想的矩形窗，导致重建图像有Gibbs现象，表现为明显的震荡响应。但是由于该滤波器形式简单，因此在我们的仿真试验中，利用CBP重建图像都是采取R-L滤波器。（2）式中，一般此值不一定为的整数倍，即（为整数），。因此。这样如何求任意的时的便成为问题的关键。最常见的方法是用内插去求解。插值的过程可以在滤波以前，利用离散的，应用内插得到在轴上连续取值的，就可使用式（2-7）求解。也可以在滤波后求得后，在通过内插的方法求得。两种方法的

26、效果一样。最近点内插和线性内插是最常见的插值方法，其具体过程可以参见3.1节。在我们的仿真实验中采取的插值方法是线性插值法，虽然线性插值法是不精确的内插方法，然而它的运算方法颇简单。由步骤（1）所介绍精确的内插方法是函数，运算要比线性插值法复杂得多，因此我们没有采用。（3）设重建图像坐标为，对每个待重建点进行反投影迭加有 （2-11）2.3 直接傅里叶重建法直接傅里叶重建算法思想更为直接，即根据投影定理有投影数据求出原图像二维傅里叶变换的每一个切片，将所有切片组合起来即可得到原图像的二维傅里叶变换，在对其进行二维傅里叶反变换就可以重建出原图像。由上面描述可以得出直接傅里叶重建算法的数学描述如下

27、3：（1）假设投影数据为,其对t的傅里叶变换为，则可以按式（2-12）求得： （2-12）（2）假设原图为，其二维傅里叶变换为，根据投影定理可以按式（2-13）求得: (2-13)（3）对进行二维傅里叶反变换求得 （2-14）式（2-14）就是直接傅里叶重建法的计算公式。但在实际的实现中，由于采集数据的天然离散性，需要对式（2-14）进行离散实现，得到直接傅里叶重建算法的数字实现步骤如下：（1）对投影数据进行傅里叶变换，对式（2-12）离散化得到： （2-15）（2）式（2-15）给出了原图呈离散极坐标状态的原图二维傅里叶频谱数据，假设原图在直角坐标系下的二维傅里叶频谱数据为，则可以通过对的插

28、值求得。（3）对式（2-14）进行离散化，则可以得到重建图像的离散公式如下： （2-16）2.4 代数重建算法卷积反投影算法与直接傅立叶算法都是变换法重建图像。变化法的突出优点是实现简单，速度快，对足够精度的投影数据能获得很好的重建质量。除了变换法以外，另一类完全不同的重建方法是迭代重建算法5。迭代重建算法是基于模型的重建算法。它能和特定的成像设备及数据采集物理过程结合起来，并能利用某些先验信息，因此即使对于比较粗糙的投影数据，这类方法也能给出较好的重建质量。代数重建算法是迭代重建算法的一种。迭代算法已经被广泛研究了很多年。这种方法假设图像有一组未知的阵列组成，然后得到一组关于此阵列的等式。虽

29、然这种方法在概念上比前面的变换法要简单，但是在医学应用上，这些方法往往不能重建出质量足够好的图像，并且耗费的时间太多。这些算法的实施方法如下：用一个二维矢量表示一个物体，用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之

30、差，修正估计p(j),这样两个之差减小。通常，在估计过程中在上加一定的限制条件。最常使用的一个限制条件是非负性1。文献中有多种形式的ART（Algebraic Reconstruction Techniques），由于各种ART算法本质都相似，我在算法比较中只采取最常见的一种ART加型ART（ART II）。（乘型ART与加型ART类似，只是调整方法有所不同。）ART II满足线性不等式 p A （2-19）其中p为投影矢量，A为系统矩阵，由图像形状以及投影位置共同决定。设(0)为迭代初始值，则迭代步骤按式（2-20）进行。 （2-20）上述校正过程是逐条射线依次进行的，故又称为逐线校正。文献

31、表明，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解的形式为，它的迭代步骤如式（2-21）所示2： （2-21）其中初始值一般选取全1矢量（尽量不要取0，否

32、则在迭代时出现除0的情况）。可以证明,按式（2-21）迭代求解，得到的序列收敛于p = A的最大熵解。在我们的仿真实验中，只采用加型ART与DF法对比，而不再采用乘型ART法与DF法对比。SIRT（Simultaneous Iterative Reconstruction Technique）同时迭代重建法是与ART并行的另一种迭代重建算法。这种方法与ART不同的是，它试图通过减慢迭代收敛的速度，来取得比ART法更好的重建图像。SIRT与ART相比，最大的优势在于对测量误差不敏感。虽然在SIRT中，也通过计算其投影误差值，但是与ART不同，在SIRT中，不是每一根射线误差都会调整射线通过网格点

33、的像素值，而是把所有的射线投影误差值计算之后，再去修改迭代网格的像素值。正因为它不是逐个射线调整，而是逐点调整，使得SIRT对单个射线误差不敏感2。SIRT迭代算法的过程与使结果满足最小二乘方准则有关。即SIRT努力使解满足满足最小二乘方准则，也就是说使得式（2-22）有最小值。 （2-22）一般采用迭代算法去求解，具体的迭代方法如式（2-23）所示： （2-23）式（2-23）的意义是：取测量值 p作为初始图像，在求第k+1次迭代求(k+1)时，利用第k次的估计值(k)加上校正图像。校正图像正比于第k次估计的误差矢量反投影。因而每个像素的校正值实在是通过该像素的所有射线和的误差值的累加，而不

34、是只与一条射线有关。由于每一像素的校正值是所有通过该像素射线的共同贡献，一些随机误差被平均掉了。为了便于迭代计算，一般将式（2-23）改写为式（2-24）的形式： （2-24）其中I为射线条数，在式（2-24）中很容易看出表示经过j号像素的所有射线的投影误差对该像素值的校正。正因为如此，SIRT的校正过程被称为逐点校正。其中初始值选择，是因为这样选择能够加快SIRT的收敛速度。其中为图像的反投影。反投影重建图像会引起星状误差，在CBP算法中，图像的星状伪迹通过滤波除去；在SIRT算法中，图像的星状伪迹由迭代调整消除。由于在CBP算法中调整方法更加直接，因此其重建速度也更快。2.5 线性调频z变

35、换（CZT）离散傅里叶变换（DFT）可以看作是点的序列的变换在单位圆上等间隔取样点。同时DFT也可以看作是序列的傅里叶变换的等间隔取样点。但是在实际中也许：不需要计算整个单位圆上变换的取样，如对于窄带信号，只需要对信号所在的一段频带进行分析，这时希望频谱的采样集中在这一频带内，以获得较高的分辨率，而频带以外的部分可不考虑；或者对其它围线上的变换取样感兴趣，例如语音信号处理中，需要知道变换的极点所在频率，如极点位置离单位圆较远，则其单位圆上的频谱就很平滑，这时很难从中识别出极点所在的频率，如果采样不是沿单位圆而是沿一条接近这些极点的弧线进行，则在极点所在频率上的频谱将出现明显的尖峰，由此可较准确

36、地测定极点频率。CZT在一定程度上解决了上述问题。序列的变换公式如式（2-25）所示6： （2-25）取 ，则式（2-25）可以化为： （2-26）公式（2-26）就是线性调频变换的公式。由公式（2-26）可以看出CZT的一些性质：当时，CZT变换路径为单位圆上一段圆弧，但是变换之后的点数不一定等于原序列点数；当 时，CZT变换起始点在单位圆外，当时，CZT变换起始点在单位圆内；当 时，CZT变换向外螺旋扩展，当时，CZT变换向内螺旋扩展。在我们的试验中，希望得到的是信号的频谱分析，故应在单位圆上去实现CZT，并且由于我们希望得到信号傅立叶频谱比较密集的采样，以减小网格化的误差，因此在我们实现

37、的DF重建算法中，取，且，在本文的后面几章将看到采取CZT近似信号傅立叶变换重建图像取得了较好的效果。第三章 高精度直接傅里叶重建算法研究按照直接傅里叶重建算法的思想，必须对频谱数据进行从极坐标空间到直角坐标空间的插值，这主要是因为两个方面的原因：第一，为了应用FFT求，因为二维FFT是基于矩形网格的，目前还没有有效的基于极网格下的二维FFT算法存在；第二，图像的显示大多是基于直角坐标系的矩形网格描述。因此在我们的仿真实验中按照公式（2-16）对图像进行重建，关键是对的求解，由公式（2-15）得到的傅里叶频谱数据呈极坐标形式，因而直接傅里叶重建算法需要对对频谱数据从极坐标向直角坐标的插值。如图

38、3-1所示，显示的是极坐标与直角坐标采样的不同特点。二维频域中的插值不像真实空间中的插值一样直接。在真实空间中，一个插值误差局限于像素所在的小区域。然而，对于频域插值，这个特性并不正确，因为二维傅里叶空间中每个采样表示某一个空间频率。于是，在傅里叶空间中一个单独采样点的误差会影响整个图像（经过傅里叶反变换后）的外貌，这也是使得大多数DF重建方法很难得到与CBP法相同重建质量的主要原因。因此，利用直接傅里叶重建算法要产生高精度的重建图像，需要在傅里叶空间内非常精确的插值，插值技术也成为了DF重建算法的关键技术。 （a）极坐标采样示意图 （b）直角坐标采样示意图图3-1 极坐标到直角坐标插值示意图

39、3.1 常见插值方法在重建算法中最为常见的插值方法有：最近点插值法、双线性插值法以及三次样条插值法。其中最近点插值法最简单，但是误差也最大，三次样条插值法重建得到的图像质量最好。这几种插值方法均属于所谓的常态插值法，即插值后的函数在取样点的取值保值不变。另一种类型的插值法被称作网格化过程，它是基于采样定理导出的，并且不具备常态插值的特点。这种网格化过程的关键是选择卷积函数。总的来讲，已经发现了一些精心设计的插值函数，可以使得才有DF法重建图像质量不亚于采用CBP法重建图像质量，但是这些插值函数往往较复杂，在本文中不予讨论。在我们的仿真试验中，采用的是最简单的最近的插值法，实现高精度重建图像。为

40、了比较应用最近的插值法、双线性插值法以及三次样条插值法的直接傅立叶变换法重建图像质量，我们在这里将简要介绍这三种插值方法，并且给出了仿真图像的对比。图3-2 极坐标空间网格化示意图如图3-2所示，当采用最近点插值法时，在直角坐标系的P点的值，将采用离P点最近的极坐标系中点的值替代。也就是采用使取值的那个点作为P点值的估计。这种方法虽然最简单，但是误差较大，只有在A，B，C，D点的距离足够小时，才能实现比较精确的插值。在我们的试验中，正好满足这一条件，所以在后文可以看到我们采用最近点插值方法实现了较高精度的图像重建。在双线性插值法中，先需要估计E、F点的值，再求P点的值。其中E，F点可采用式（3

41、-1）估计。 （3-1）然后采用式（3-2）计算的值作为估计点P的值。 （3-2）在大多数情况下，利用双线性插值法估计得到的值比利用最近值插值法估计的值要更准确，但是其算法也更复杂。在我们的试验中由于数据足够密集，因此没有采用双线性插值法，而是采用最近点插值法，也能取得很好的结果。在这三种插值法中，重建图像质量最好的插值方法是三次样条插值法，但是三次样条插值法计算量较大。由于重建图像的质量对与径向方向的插值精确性敏感性远大于对与相角方向的插值精确性。因此这里介绍的是混合二阶样条线性插值法7。混合二阶样条线性插值法的计算步骤为：首先沿径向方向进行三次样条曲线插值，其次沿相角方向进行线性插值即可。

42、如图3-2所示，设、分别为A、B、C、D沿径向方向的偏导数，则沿径向方向的点E、F估计值为： 然后利用式（3-2）沿相角方向进行线性插值得到P点的值。选取Shepp-Logan头部模型作为仿真原模型，投影数据为367360，重建图像大小为256256。选用以上三种插值算法的DF重建图像以及局部放大图像请参见图3-3。其中图（a）、（b）、（c）分别是应用最近点插值、双线性插值以及混合二阶样条线性插值法重建的图像，图像（d）、（e）、（f）分别是对应的三种插值法重建图像的局部放大图。有图3-3的对比可以看出，由最近点插值得到的图像，质量最差，边缘的震荡最明显以及图像的混淆误差都明显比用双线性插值

43、和混合二阶样条线性插值法重建图像误差大。而应用双线性插值法和混合二阶样条线性插值法重建的图像虽然在视觉上难以区分其图像质量，然而应用一些客观评价标准能够判别图像质量好坏。（理论上，三次样条插值法在多数情况下都是一种比较好的插值方法，尤其对平滑的曲线插值8，因此，可以假定由混合二阶样条插值的DF法重建图像的质量应优于用线性插值法重建图像的质量，然而在这里，由于频谱数据较为充足，因此没能显示其重建图像的质量优势。）在这里我们采用第4.3节中介绍的几种评价图像质量的客观标准作为评判这三种方法重建图像质量的标准。表格3-1显示了三种插值法重建图像的这几种判据大小。对于归一化平均绝对距离判据，双线性插值

44、法在三种方法中表现得最出色。在我们的其他试验仿真中表明，当投影射线条数较少时，二阶样条线性插值法在三种插值法中表现得最好。由于我们的论文主要考虑在投影数据充足时的重建效果，因此，稀疏投影数据的仿真图没有收到本论文中。表3-1 重建图像质量评价参数对比插值方法最近点插值双线性插值二阶样条线性插值法归一化均方距离0.27460.23910.2398归一化平均绝对距离0.0520.03910.0395最坏情况距离1477 （a）最近点插值 （b）双线性插值 （c）混合二阶样条插值 （d）最近点插值局部放大 （e）双线性插值局部放大 （f）混合二阶样条插值局部放大图3-3三种插值法重建图像以及其局部放

45、大图3.2 与卷积反投影算法等效的直接傅里叶重建插值方法由于卷积反投影法不但能重建出较好的图像质量，并且它的计算复杂度也在能够接受的范围内，因此几乎所有现代的CT机中使用的重建算法都是卷积反投影法。虽然DF重建算法能够提供较好的计算效率，但是DF算法重建图像质量在大多数情况下都比CBP差。就像3.1节中讨论的几种插值方法，其重建图像都比CBP差。因此，可能有人会怀疑CBP算法会比所有的DF重建算法要好，但是在Hyeokho Choi 和 David C. Munson, Jr.的论文中，提供了一种插值方法使得DF算法与CBP等效，因此可以把CBP看作DF算法的一种有效实现。在这里，我将简单讨论

46、对比这两种重建方法9。在平行束CT中，图像的投影是从各种不同方向采集到的。如果角度采集足够密集，就能忽略图像混叠。因此在CBP法中，可以看作是极坐标下的二维傅里叶反变换，如（3-5）式所示： （3-5）我们假设f(x,y)在时域上限制在一个方形区域中，特别地： （3-6）因此可以把公式（3-5）改写为 （3-7）其中。在现实中，由于只能在离散点取得，不妨假设能在，并且，其中，那么我们把式（3-7）离散化可以得到 （3-8）定义fn,m = f(nx,my)为f(x,y)的采样，并且n = -N,0,N; m = -N,0,N. 可以求得fn,m的近似值，计算式如下： （3-9）由公式（2-2）

47、可知CBP法是实现（3-5）式重建图像的一种方法，而式（3-9）是根据式（3-5）代数变换得到的，因此可以看出由式（3-9）重建图像与CBP法重建图像等效。因为在DF重建算法中，最后一个步骤是对频谱数据进行二维离散傅里叶反变换。由于式（3-10）式插值得到的是式（3-9）中的二维DFT。因此可以得出用式（3-10）进行插值得到的频谱数据再进行二维离散傅里叶反变换重建得到的图像与CBP法重建图像等效。 由式（3-10）我们可以得到按以下两个步骤重建图像与CBP算法重建图像等效9：步骤一：利用式（3-10）对傅里叶频谱数据进行极坐标到直角坐标系的插值；步骤二：计算插值好的频谱数据的二维离散傅里叶反

48、变换。显而易见，经过式（3-10）插值的DF重建算法与CBP算法中采取矩形窗的滤波函数等效。同时，也很容易按照上面的推导方法找到与CBP算法中采取其他窗函数的滤波函数等效的插值方法，只是插值函数的加权系数不会是雅克比加权系数。当然，两种方法也有微小的不同。在CBP算法中，在反投影的步骤中需要对滤波后投影进行插值，而应用式（3-10）进行插值的DF重建算法不需要进行插值。因此，确切地说，应用式（3-10）进行频域插值的DF重建算法与在反投影步骤中对滤波后投影进行精确插值的CBP算法等效。按照式（3-10）重建图像的仿真结果如图3-4所示，其中CBP法在反投影步骤中应用了的插值方法为线性插值法：由

49、图像仿真结果同样可以看出，利用式（3-10）重建图像结果与利用CBP法重建图像结果很类似，以至于直接用视觉观察难以区分二者差别，但是通过第4.3章介绍的评判重建图像的标准，可以看出二者的区别。表3-2显示了这两种方法重建图像的客观评价标准指标。然而式（3-10）的插值法计算复杂度为，这使得利用这种插值法的DF算法其复杂度要高于CBP算法，因此这种算法没有很大的实用性。不过在文献9中提出的一种对式（3-10）的改进算法，即对式（3-10）描述的插值器进行加窗处理，因为对（3-10）中影响较大的均为靠近的点。在文献9中表明，对式（3-10）利用166的窗函数时，重建效果与式（3-10）插值效果类似

50、。虽然对插值器进行加窗处理后计算量大大降低，同时对图像质量影响不大，但是与CBP法相比，即使是加窗的插值器也没有比CBP法节省计算量。因此式（3-10）的最大意义在于证实存在一种非最优的插值器使得DF重建法与CBP法等效，也就表明了如果利用最优的插值器有可能使得DF法重建图像质量要优于用CBP法重建图像质量。文献7中同时也表明利用这种插值法的DF法重建图像之所以复杂度比CBP法要高，很可能是由于这种方法没有考虑到极坐标中频谱数据的网格特性，而CBP法利用了极坐标的网格特性。 （a）CBP法重建图像 （b）等效DF法重建图像 （c）CBP法重建图像局部放大图 （d）等效DF法重建图像局部放大图3

51、-4 CBP法与等效DF法重建图像以及其局部放大图表3-2 CBP法与等效DF法重建图像质量评价参数对比重建方法CBP等效DF法归一化均方距离0.24240.2336归一化平均绝对距离0.03880.0358 续表3-2重建方法CBP等效DF法最坏情况距离983.3 一种新的高精度DF重建算法3.3.1 投影数据的获得实际的投影数据是由检测器测得的，在计算机模拟时，则取自仿真模型，在我们的试验中，采取Shepp-Logan模型。由于Radon在他的经典论文中描述的问题由投影重建图像，对科学的众多领域产生了重大影响。因此，人们把联系目标函数空间到投影函数空间的积分变换称为Radon变换10。Ra

52、don变换的公式如式（3-11）所示： （3-11）其中，为Dirac delta函数。叫做沿方向的投影。如果注意到坐标与它的旋转坐标的关系，即式（3-12） （3-12）那么式（3-11）的Radon变换式可以写作式（3-13）： （3-13）其中，。式中，和是定义平面上一条直线的位置的两个变量，该直线由方程决定，对一个固定的和，这个直线方程决定了一条与轴平行的直线，为沿该直线的线积分。由式（3-12）可以发现，空间的一个点对应于空间的一条正弦曲线段11。在数学上Radon变换是一个线积分变换，它可以看做是一个算子。该算子对空间的函数作用产生空间的另一个函数。与CT问题联系起来，Radon在

53、不同的CT系统中具有不同的物理含义。称作被观测目标的表征函数，则是与这一表征函数有关的测量值。在试验中模拟投影数据有一般两种方法，一种通过解析法计算，另一种通过数值累加法计算。解析法计算实际上就是利用式（3-11）对模拟图像求其Radon变换；而通过数值累加法求其Radon变换需要做某些方面的近似。为了对CBP、DF以及ART算法进行对比，在我们的试验中，采取数值累加法求投影数据（因为在ART算法中，需要对算法开始即对原图进行离散化处理，故需要采用数值累加法求投影数据），这样三种方法都可以采用相同的投影数据，使得重建图像更具有可比性。在我们的仿真实验中，图像在某个方向的投影数据是所有像素点在此

54、方向投影数据之和。因此我们采取的待重建图像为数字图像（在实际的重建过程中，也可以把原图采样为数字图像，因为重建后的图像是天然离散的）。图3-5投影示意图如图3-5所示，每个像素点的投影都按比例分到离它最近的两个投影点上（也可以更精确的把每个像素再平均分为多个小的像素的，然后再求投影）。虽然在CT机上采集的数据方法与此不同，但是这种方法较精确的模拟了CT投影的采集值。按照这种算法思想，我们可以在每个投影角度遍历原图的每一个像素点，得到每个像素点在每个投影角度的投影值，累加之后就可以得到每个投影角度的投影值12。图3-6投影数据图像图3-6为修改过的Shepp-Logan模型的投影数据，其中原图为

55、256256,投影数据为367180，在本文后面三种方法的对比中，采用的投影数据均采用这种方法获得，并且投影数据都为367360.3.3.2 投影数据一维傅里叶频谱数据的计算对投影数据频谱值的取得，我们采取式（2-14）计算。观察式（2-14），在我们的实验中获取的投影数据，可以令t = 1, d = 1/(128*N),其中N为重建图像大小，在我们的试验中N均取256。重建图像时，d取值越小，重建图像质量越高，但是图像重建速度也越慢。因此代入数据式（2-14）可以化为式（3-14）的形式 （3-14）式（3-14）实际上是对序列求解其线性调频z变换（CZT）。由CZT的性质，我们可以知道当C

56、ZT变换路径为单位圆上一段圆弧时，能够得到信号的频谱分析13。在这里我们没有采用DFT作为投影的频谱分析，因为DFT对频段的抽样密度不够，这样导致插值误差较大。而CZT能够满足我们对提高频段抽样密度的要求，因此我们选择使用CZT作为对数据的频谱分析。对每一个角度的投影数据，按照式（3-14）求得其频谱数据。根据投影定理，每个角度投影数据的傅立叶变换对应原图二维傅立叶变换的一个切片。把所有角度的频谱数据组合起来，也就得到了原图二维傅立叶变换的很多切片，如果切片足够多，根据采样定理，就可以精确的得到原图的二维傅立叶频谱数据。根据式（3-14）求频谱数据，其频率间隔，其中为图像大小。根据试验，我们可

57、以发现当，频谱抽样密度已经足够大，能够使网格化误差控制在可接受的范围内。当然如果需要提高精度，可以继续减小，但是这样会增加计算复杂度，延长了计算时间，同时对计算机的存储空间有较大的要求。当时，在计算所需时间和网格化误差大小作了较好的权衡。3.3.3 频谱数据的网格化在3.3.2 节得到的频谱数据是极坐标下的频谱数据，而DF算法的最后一个步骤进行的二维傅立叶反变换需要的频谱数据需要的是直角坐标系下的，因此需要经过频谱数据从极坐标到频谱数据网格化的过程。如在3.1节中讨论的所述，大多数插值方法都不能满足精确重建图像的要求。虽然到目前为止，已经发现了一些插值方法，使得DF法重建图像质量不差于CBP法

58、重建图像质量。然而这些插值方法往往很复杂，实际运算中并不方便。在我们的试验方法中并不对插值方法改进，而是对DF重建算法的最后一步进行改进以减小插值误差。在传统DF算法中，最后一步是对频谱数据进行二维DFT逆变换，也就是说，在网格化过程中，得到的频谱数据是原图像DFT频谱数据的近似。这样，由于频谱数据的间隔较大，加大了插值的难度，从而增大了插值过程的误差。在我们的试验中，DF算法的最后一步是对频谱数据进行Counter-z变换，这样能减小频谱数据之间的间隔，从而减小在插值步骤的误差。网格化方法采取最近点插值方法（具体步骤参见3.1节）。当然在这个步骤的插值过程中，也可以采用其它比最近点插值法更好的插值法。然而由于在我们的仿真试验中，我们使得频谱数据足够密，这样使得其它的高阶插值方法与最近点插值法相比没有很大的优势。而且我们的仿真试验也显示利用双线性插值法以及混合二阶样条线性插值法重建的图像没有比采用最近点插值重建图像质量有明显改善，却耗费了比最近点插值法多得多的计算时间，因此我们选择最近点插值法作为频谱数据网格化的方法。（限于空间的缘故，在这里我们没有把利用双线性插值法和利用混合二阶样条线性插值法重建的图像与利用最近点插值重建的图像列出来进行对比。）3.3.4 由频谱数据重建图像由3.3.3节得到了原图像二维傅立叶频谱