生物医学工程毕业论文

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

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

2、够好的图像，同时其耗费的时间也在可以接受的范围内，因此现代 CT 中的重建算法几乎都是用 CBP 算法。论文从 DF 算法入手，试图在投影数据足够的条件下寻找好的网格化方法，改善 DF 重建算法的质量。在我们的论文中，研究了 DF 法中应用三种常见插值器方法时，图像重建效果对比。同时，论文也指出了采用一种插值器能使 DF 法与 CBP 法等效。最后我们提出了一种新的 DF 重建算法。在这种算法中采用了对投影数据进行线性调频 z 变换近似其频谱数据的方法，使得频谱数据更密，同时在 DF 法的最后一步中采取线性调频 z 变换求重建图像，减小了插值误差。论文中把这种方法与 CBP 重建以及ART、S

3、IRT 进行了对比，仿真结果显示，本文中采用的方法重建图像的质量至少不比 CBP重建以及 ART、SIRT 重建图像质量差。关键词：计算机断层成像，直接傅立叶变换重建，卷积反投影，代数重建算法 ，同时迭代重建法 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

4、 inverse problem of reconstructing an image from its projection data. 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

5、 of lacking a good interpolation method interpolating data from the polar 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 c

6、an yield to good quality reconstructions and only take a short time , 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 inte

7、rpolators widely used in the DF method and construct the images reconstructed 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

8、 projection data to approximate its frequency spectrum data, in which 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 simulat

9、ion results shows that the image reconstructed by this method is at least 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 Reconstructio

10、n Technique. 目录目录第一章 绪论.11.1 断层成像技术发展简介.11.2 CT 图像重建.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.

11、4 由频谱数据重建图像.173.3.5 重建图像对比.17第四章 仿真结果与算法比较.194.1 图像仿真结果.194.2 重建后图像与原始图像差异评价.214.3 灰度值比较.224.4 计算复杂度比较.23第五章 结论.25参考文献.26谢辞.27译文及原文.28 第 1 页第一章 绪论CT（Computerized Tomography），即计算机断层成像，是用来获取观测目标断层图像的一门技术。它广泛应用于诊断医学、射电天文学、电子显微和雷达等多科学领域。CT是由在多个观测角度获得的有关目标的一系列投影数据，通过图像重建技术来得到目标的断层图像的。近年来，X- CT无论在基本技术方面，还

12、是在新的临床应用方面都取得了巨大的发展。在CT的各个主要组成部分，如光管、探测器、滑环、数据获取系统和算法等方面，都取得了很大进步。自从螺旋CT和多层CT问世，出现了许多新的临床应用。CT经过三十多年的发展以后，再次成为医学图像领域中最令人兴奋地诊断方法之一。1.1 断层成像技术发展简介自从 1895 年，伦琴在试验阴极射线管时发现了 X 射线后，X 射线便发展为广泛使用的检测工具。用传统的照相方法，三维的人体沿 X 射线方向被压缩成了两维的图像。体内所有骨骼结构和组织都重叠在一起，使得感兴趣对象的清晰程度大为下降。由于传统照相的这种限制，导致了传统断层成像技术的出现。Bocage 是传统断层

13、成像技术的先驱之一。早在 1921 年，他描述了这样的设备，能使感兴趣平面以外上下的结构模糊得看不清楚。虽然传统断层成像术有了一定的发展，而且这些技术在生成感兴趣平面的图像方面取得某种程度的成功，却都有着最基本的限制，即它们并没有增加物体的反差，也不能从根本上去除焦平面以外的其他结构。加上病人须接受很大的 X 射线剂量，传统断层成像现在很少应用到临床上了。而从投影重建图像的努力也早在 1940 年已经开始，再 1940 年颁布的专利中，Gabriel Frank 描述了现代断层成像技术的基本思想。1963 年，Kuhl 和 Edward 应用放射性同位素提出了横向断层成像方法，该方法被后来进一

14、步发展和改进为今天的发射计算机断层成像（ECT）。从多个投影数据重建图像的数学公式可以追溯到澳大利亚的数学家Radon，他在 1917 年用数学证明从无限多投影数据可以重建出原来的图像。在医学上，Allan MCormack 报道了可能是第一台实际建成的 CT 扫描机的研究结果。英国 EMI 的中心研究实验室的 Godfrey NHounsfield 于 1967 年开始第一个临床 CT 扫描机的研制。再进一步改进了数据获取和重建技术后，第一台可供临床应用的 CT 机安装在Atkinson-Morley 医院。虽然第一代扫描机的临床检查结果还可以就接受，但是在长达 4 分半钟的数据获取过程中病

15、人的运动引起的图像质量问题是严重的。第二代扫描设备仍然是一台平移旋转扫描机，然而由于利用了多个笔形束使得所需旋转的步数减少了。最流行的扫描机类型是第三代扫描机。这种结构将大量探测器布置在以 X 射线源为中心的圆弧上，探测器的尺寸足够大使得整个检测对象始终落在探测器的视场内。X 射线源和探测器在整个设备围绕病人旋转时保持相对静止。直线移动的取消显著减少了数据获取时间。第四代扫描机中的探测器形成一个闭合的圆环，在整个扫描过程中 X 光管围绕病人旋转，而探测器保持静止。第五代扫描机也称作电子束扫描机，于 1980 年到 1984 年建造，用于心脏检查。为了减少心脏跳动形成的伪迹，要求采集一套完整的投

16、影数据必须在 2050ms内完成。因此在电子束扫描机中，射线源的旋转由电子束扫描运动代替了传统的 X 光管的机械运动来完成。自从第一台临床扫描机采用以来，CT 技术取得了巨大的进展。CT 机已经从第一代发展至第五代，扫描时间也从 5 分钟缩短至几十毫秒，空间分辨率也达到 1mm以下1。 第 2 页CT 起源于医学成像，随着 CT 技术的发展，CT 技术的应用也深入到重建图像的科学和技术的其他各个领域。因此，从广义上讲，CT 是这样一个过程，它通过某种辐射源（如 X 射线，放射性核素、超声波、磁场等）对观察目标进行作用并检测投影数据，然后利用计算机完成图像重建、数据处理并和图像显示技术结合，得到

17、观测目标内部的（通常具有某种物理性质的）断层图像或三维立体图像。CT 的成像大致过程可用图 1-1 所示方框图来表示1,2。图图 1-11-1 CTCT 成像过程流程图成像过程流程图图 1-1 显示的 CT 图像成像流程图中，为简便起见，我们把“预处理”模块放到“负对数运算”模块前面，然而，事实上有些预处理步骤在对数步骤之后完成。在后处理模块中，主要包括伪像校正、图像增强、3D 以及其他各种重构方法。1.2 CT 图像重建图像重建技术是 CT 技术中的一个重要问题。目前已经发现许多种从投影重建图像的方法，这些方法中，往往不同的方法有不同的适应范围，可以在不同的场合下考虑使用不同的重建方法，以获

18、得最好的图像重建质量或者缩短计算时间。所有的图像重建方法一般分为两种类型：变换方法和级数展开法。变换法也称解析法，其基本特点是目标函数与投影函数之间的关系可以用某种解析公式来表达，一旦这种解析公式建立，便寻找一种合适的离散实现这一解析公式的具体算法。变换法中两类最具代表的算法是滤波反投影重建算法和傅立叶重建算法。由于变换法实现比较简单，速度快，这类方法（尤其是卷积反投影法）是商用 CT 机中普遍采用的方法。而级数法的特点是，从一开始就把连续图像离散化，从而建立起离散重建图像和离散投影之间的代数方程组，然后在某种最优准则的指导下解这个线性方程组。由于这个方程组一般都非常大，很难采取直接的方法求解

19、，因此这种求解过程一般是一个迭代过程。1.3 研究目的及论文结构论文中主要研究高精度的直接傅立叶变换重建算法（DF）。在论文中先概述了CBP法、ART法以及DF法，并表明现在的CT机中广泛使用的算法是CBP法。在讨论DF算法中，我们先研究了比较常见的一些插值方法（最近点插值、双线性插值以及混合二阶样条线性插值法），通过仿真表明双线性插值以及混合二阶样条线性插值法要比最近点插值法中取得更好的效果。接着论文证明了可以采用一种插值方法使得DF法与CBP法等效，从而表明有可能找到更好的插值方法使得DF法重建图像质量优于CBP法。最后我们介绍了我们重建图像的方案。经过大量图像仿真的试验，我们采用了增加频

20、谱数据密度的方法，减小网格化误差，取得了较好的实验结果。论文的章节内容具体安排如下：第一章，简单介绍计算机断层成像地发展；第二章，综述了四种常见CT重建算法（CBP、DF、ART、SIRT）以及从投影重建图像的基础知识；第三章，详细介绍了DF算法,特别是我们对DF的数字实现方法；第四章，对第三章的算法进行了仿真试验，并与CBP以及ART、SIRT进行了对比；数据采集预处理或标定负对数运算断层图像重建后处理窗口显示 第 3 页第五章，在前面几章的基础上，总结了本文主要内容以及结论。第二章 从投影重建图像算法概述从投影重建图像的算法很多。常见的有：反投影重建算法；滤波反投影重建算法；直接傅立叶变换

21、重建算法；迭代重建算法等。在本章，我们将在第一小节介绍投影定理，投影定理是卷积反投影重建算法以及直接傅立叶重建法的基础。然后，为了与我们采用的算法进行对比，我们简单介绍其中四种最常见的重建算法：卷积反投影法（CBP）、直接傅立叶重建法（DF）、代数重建法（ART）、同时迭代重建算法（SIRT）。由于在DF法中，我们需要使用线性调频z变换（CZT），因此在本章最后一小节，我将简单介绍CZT。2.1 投影定理投影定理或中心切片定理是滤波反投影重建算法与直接傅立叶变换重建算法的基础。在没有衍射源的情况下，其内容是：某图像 f（x，y）在视角为时投影的一维傅立叶变换给出 f（x，y）的二维傅()rP

22、x立叶变换的一个切片。切片与 轴相交成角，且通过坐标原点2。12(,)( , )FF 即 （2-121(,)( , )()( )( , ),rAAF pxPP 1）根据投影定理，投影重建图像的方法原则上可按如下流程求解：（1）采集不同时脚下的投影，理论上应为 1800范围的连续无穷多投影；（2）求出个投影的一位傅立叶变换，这就是图像二维傅立叶变换过原点的各切片，理论上应为连续的无穷多片；（3）将（2）求出的各切片组成图像的二维傅立叶变换；（4）用（3）得出的图像傅立叶频谱数据求二维傅立叶反变换得到重建图像。在实际的重建过程中，由于反变换时的数学处理不同，又可将变换法分为滤波反投影重建算法和直接

23、傅立叶反变换算法。其中滤波反投影重建算法又可分为卷积反投影法和Radon 反变换法。Radon 反变换法由于对误差很敏感，因此在商用 CT 机上未得到应用。从数学表达式看，Radon 反变换法与卷积反投影法等效，因此我们没有仿真 Radon 反变换重建方法，而只使用卷积反投影法重建图像。在下文的讨论中，我们也只介绍卷积反投影法，而不再介绍 Radon 反变换法。2.2 卷积反投影法2由投影定理，可以通过在不同视角下的投影的一维傅立叶变( , )A ( , )a x y()rP x换求得，因此待见图像： 第 4 页12()1212121222cos()2cos()12002cos()2co1(

24、, )( , )( (,)(,)4 =(,)( , ) =( , )( , )ixyiriririra ra x yFAAeddAed dPed dPed ddPe s()00. (2-2)d 从公式（2-2）的第二个积分，即 （2-2cos()ired 3）该式可以改写成空域变量 xr的傅里叶反变换式：2cos()2cos()cos()( , ) ()(, ) cos(), , (2-4)rririrxrrrxredPedh xp xg r 式中 （2-(, )(, )* ()rrrg xp xh x5）而，。将式（2-4）代入式（2-2）后，得1()( )rh xF11(, ) ( , )

25、rp xFP （2-0( , ) cos(), a rg rd 6）式（2-6）的物理意义是：经过给定点的所有滤波后的射线投影在范围内累( , )r0加反投影重建，得出点的像素值。( , )r由式（2-6）离散化可以得出滤波反投影重建算法数字实现的具体步骤：（1）对采集到的离散数据进行离散滤波； (,)(,) ()(,) ()* ()(,)* ()NNrrrrrnNnNp x mp nd mh xndp nd mxndh xp nd mh x （2-7）式中为滤波函数。为了更清楚的对的选择进行讨论，我们转到频域去讨论。()rh x()rh x设为的傅里叶变换，则理想的滤波函数，但在种滤波函数是

26、不可( )H()rh x( )H能实现的。在实际中，的高频分量幅度很小，因此近似为带限信号。 （事实上，()rp x( )P在物体尺寸有限的情况下，投影数据分布在有限的范围内，因此严格的说，其频带是无限的。然而如果物体的密度在空间范围内变化平稳，则高频分量幅度确实不大，并且检测器在接收 X 射线时，有平均作用，相当于低通滤波，可以认为高频分量很小。 ）利用这( )P 第 5 页一事实，我们可以把强制为带限，即可以对理想的滤波函数加一个窗函( )H( )H数。常用的窗函数有矩形窗、汉明窗等。在我们的仿真试验中采用简单的矩形窗，也就是使用 R-L 滤波函数。R-L 滤波函数的系统函数为：( )R

27、LH （2-( )(/ 2 )R LHrectB8）式中1, 1/ 2(/ 2 )0, BdrectB；其他。相应的冲激响应函数为：()R Lrhx （2-2222()2sinc(2)sinc ()rBixR LrrrBhxedBx BBx B9）同时，由于投影数据的采集是离散的，因此，滤波函数也应相应的取其离散形式。对进行采样，其中 d 为采样间隔对应的最高不失真空间频率为，以()R Lrhx1/(2 )Bd代入式（2-9） ，将得到的离散形式如（2-10）所示：rxnd()R Lrhx （2-22221 04()0 1 R Lndhndnnnd为偶数为奇数10）如果对进行线性内插，则得到另

28、一种连续的空域函数，可以把()R Lhnd()R Lrhx看作是的一次近似。也可以对进行精确的 sinc 函数内插，得()R Lrhx()R Lrhx()R Lhnd到，但是由于运算复杂，在实践中一般不使用，在我们的仿真试验中也采取对()R Lrhx进行线性内插重建图像。另外，由于 R-L 滤波器使用理想的矩形窗，导致重建图()R Lhnd像有 Gibbs 现象，表现为明显的震荡响应。但是由于该滤波器形式简单，因此在我们的仿真试验中，利用 CBP 重建图像都是采取 R-L 滤波器。（2）式中，一般此值不一定为的整数倍，即（cossinrxxmymdrxn d 为整数） ，。因此n()rxndn

29、n d()() rh xndh nn d。这样如何求任意的时的(,)(,)rP x mP n d m (,) ()p n d mh n drx便成为问题的关键。最常见的方法是用内插去求解。插值的过程可以在滤波以(,)rP x m前，利用离散的，应用内插得到在轴上连续取值的，就可使用式（2-()R Lhndrx()rh x7）求解。也可以在滤波后求得后，在通过内插的方法求得(,)rP x m(,)P nd m。两种方法的效果一样。最近点内插和线性内插是最常见的插值方法，其具体(,)rP x m过程可以参见 3.1 节。在我们的仿真实验中采取的插值方法是线性插值法，虽然线性插值法是不精确的内插方法

30、，然而它的运算方法颇简单。由步骤（1）所介绍精确的内插方法是函数，运算要比线性插值法复杂得多，因此我们没有采用。sinc（3）设重建图像坐标为，对每个待重建点进行反投影迭加有( ,)iix y （2-10( ,)cos()sin()Niiijnf x yp xnynN 第 6 页11）2.3 直接傅里叶重建法直接傅里叶重建算法思想更为直接，即根据投影定理有投影数据求出原图像二维傅里叶变换的每一个切片，将所有切片组合起来即可得到原图像的二维傅里叶变换，在对其进行二维傅里叶反变换就可以重建出原图像。由上面描述可以得出直接傅里叶重建算法的数学描述如下3：（1）假设投影数据为,其对 t 的傅里叶变换为

31、，则可以按式( , )p t( , )P r( , )P r（2-12）求得： （2-2( , )( , )irtP rp tedt12）（2）假设原图为，其二维傅里叶变换为，根据投影定理可以按式（2-( , )f x y( , )F u v13）求得:( , )F u v (2-13)( , )( , )F u vP r（3）对进行二维傅里叶反变换求得( , )F u v( , )f x y （2-2 ()2 ()( , )( , )( , )iux vyiux vyf x yF u v edudvP redudv 14）式（2-14）就是直接傅里叶重建法的计算公式。但在实际的实现中，由于采

32、集数据的天然离散性，需要对式（2-14）进行离散实现，得到直接傅里叶重建算法的数字实现步骤如下：（1）对投影数据进行傅里叶变换，对式（2-12）离散化得到：(,)p n t m （2-1/2 121/2(,)(,)Min tndnMP nd mtp n t me 15）（2）式（2-15）给出了原图呈离散极坐标状态的原图二维傅里叶频谱数据，假设原图在直角坐标系下的二维傅里叶频谱数据为，则可以通过对(,)F s u r v(,)F s u r v的插值求得。(,)P N d m（3）对式（2-14）进行离散化，则可以得到重建图像的离散公式如下： （2-2 ()(,)(,)KKjsk u x rl

33、 v ysK rKf k x l yu vF s u r v e 16）2.4 代数重建算法卷积反投影算法与直接傅立叶算法都是变换法重建图像。变化法的突出优点是实现简 第 7 页单，速度快，对足够精度的投影数据能获得很好的重建质量。除了变换法以外，另一类完全不同的重建方法是迭代重建算法5。迭代重建算法是基于模型的重建算法。它能和特定的成像设备及数据采集物理过程结合起来，并能利用某些先验信息，因此即使对于比较粗糙的投影数据，这类方法也能给出较好的重建质量。代数重建算法是迭代重建算法的一种。迭代算法已经被广泛研究了很多年。这种方法假设图像有一组未知的阵列组成，然后得到一组关于此阵列的等式。虽然这种

34、方法在概念上比前面的变换法要简单，但是在医学应用上，这些方法往往不能重建出质量足够好的图像，并且耗费的时间太多。这些算法的实施方法如下：用一个二维矢量 表示一个物体，用 p 表示它的投影。两个变量被一个系统矩阵 A 和一个误差矢量 e 连接起来： p = A + e (2-17)重建过程是根据一个测量矢量 p 估计 。估计是通过要求 和 e 满足测定的最优准则来进行。估计的过程是一个迭代的步骤，也就是说，我们产生这样一个矢量序列，(0), (1), (n),以使序列收敛于 *,这里 *是基于 p 的最优估计。对于每次迭代 j，我们计算 p(j): p(j) = A(j) + e （2-18）根

35、据计算的投影 p(j)与测量投影 p 之差，修正估计 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

36、）进行。 （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

37、解的形式为：。而乘型 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（Simultaneous Iterative Reconstruction Technique）同时迭代重建法是与 ART并行的另一种迭代重建算法。这种方法与 ART 不

38、同的是，它试图通过减慢迭代收敛的速度，来取得比 ART 法更好的重建图像。SIRT 与 ART 相比，最大的优势在于对测量误差不敏感。虽然在 SIRT 中，也通过计算其投影误差值，但是与 ART 不同，在 SIRT 中，不是每一根射线误差都会调整射线通过网格点的像素值，而是把所有的射线投影误差值计算之后，再去修改迭代网格的像素值。正因为它不是逐个射线调整，而是逐点调整，使得 SIRT 对单个射线误差不敏感2。SIRT 迭代算法的过程与使结果满足最小二乘方准则有关。即 SIRT 努力使解 满足满足最小二乘方准则，也就是说使得式（2-22）有最小值。 （2-1( )() ()TpApA 22）一般

39、采用迭代算法去求解 ，具体的迭代方法如式（2-23）所示： （2-(0)(1)( )( )()TkkTkA pApA23）式（2-23）的意义是：取测量值 p 作为初始图像，在求第 k+1 次迭代求 (k+1)时，利用第 k 次的估计值 (k)加上校正图像。校正图像正比于第 k 次估计的误差矢量反投影。因而每个像素的校正值实在是通过该像素的所有射线和的误差值的累加，( )()TkApA而不是只与一条射线有关。由于每一像素的校正值是所有通过该像素射线的共同贡献，一些随机误差被平均掉了。为了便于迭代计算，一般将式（2-23）改写为式（2-24）的形式： （2-(0)(1)( )( )1()TIkk

40、TkiiiiA ppaa24）其中 I 为射线条数，在式（2-24）中很容易看出表示经过 j 号像素( )1()ITkiiijipaa的所有射线的投影误差对该像素值的校正。正因为如此，SIRT 的校正过程被称为逐点校正。其中初始值选择，是因为这样选择能够加快 SIRT 的收敛速度。其中(0)TA p为图像的反投影。反投影重建图像会引起星状误差，在 CBP 算法中，图像的星状伪TA p迹通过滤波除去；在 SIRT 算法中，图像的星状伪迹由迭代调整消除。由于在 CBP 算法中调整方法更加直接，因此其重建速度也更快。2.5 线性调频 z 变换（CZT）离散傅里叶变换（DFT）可以看作是点的序列的变换

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

42、极点频率。CZT 在一定程度上解决了上述问题。序列的变换公式如式（2-25）所示6：( )x nz （2-10()( )NnkknX zx n z0,1,.,1kL25）取 ，则式（2-25）可以化为：0000()iikkzr eR e0,1,.,1kL （2-001000()( )()()NiinnkknX zx n r eR e0,1,.,1kL26）公式（2-26）就是线性调频变换的公式。z由公式（2-26）可以看出 CZT 的一些性质：当时，CZT 变换路径为单位圆上一段圆弧，但是变换之后的点数不一定001rRL等于原序列点数；N当 时，CZT 变换起始点在单位圆外，当时，CZT 变换

43、起始点在单位圆内；01r 01r 当 时，CZT 变换向外螺旋扩展，当时，CZT 变换向内螺旋扩展。01R 01R 在我们的试验中，希望得到的是信号的频谱分析，故应在单位圆上去实现 CZT，并且由于我们希望得到信号傅立叶频谱比较密集的采样，以减小网格化的误差，因此在我们实现的 DF 重建算法中，取，且，在本文的后面几章将看到采取 CZT 近001rRLN似信号傅立叶变换重建图像取得了较好的效果。第三章 高精度直接傅里叶重建算法研究按照直接傅里叶重建算法的思想，必须对频谱数据进行从极坐标空间到直角坐标空间的插值，这主要是因为两个方面的原因：第一，为了应用 FFT 求，因为二维 FFT( , )f

44、 x y是基于矩形网格的，目前还没有有效的基于极网格下的二维 FFT 算法存在；第二，图像的显示大多是基于直角坐标系的矩形网格描述。因此在我们的仿真实验中按照公式（2-16）对图像进行重建，关键是对的求解，由公式（2-15）得到的傅里叶频谱数据(,)F s u r v呈极坐标形式，因而直接傅里叶重建算法需要对对频谱数据从极坐标向直(,)P N d m角坐标的插值。如图 3-1 所示，显示的是极坐标与直角坐标采样的不同特点。二维频域中的插值不像真实空间中的插值一样直接。在真实空间中，一个插值误差局限于像素所在的小区域。然而，对于频域插值，这个特性并不正确，因为二维傅里叶空间中每个采样表示某一个空

45、间频率。于是，在傅里叶空间中一个单独采样点的误差会影响整个图像（经过傅里叶反变换后）的外貌，这也是使得大多数 DF 重建方法很难得到与 CBP 法相同重建质量 第 10 页的主要原因。因此，利用直接傅里叶重建算法要产生高精度的重建图像，需要在傅里叶空间内非常精确的插值，插值技术也成为了 DF 重建算法的关键技术。 （a）极坐标采样示意图 （b）直角坐标采样示意图图图 3-13-1 极坐标到直角坐标插值示意图极坐标到直角坐标插值示意图3.1 常见插值方法在重建算法中最为常见的插值方法有：最近点插值法、双线性插值法以及三次样条插值法。其中最近点插值法最简单，但是误差也最大，三次样条插值法重建得到的

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

47、-23-2 极坐标空间网格化示意图极坐标空间网格化示意图如图 3-2 所示，当采用最近点插值法时，在直角坐标系的 P 点的值，将采用离 P 点最近的极坐标系中点的值替代。也就是采用使取值的那个点作minmin,dPA PB PC PD为 P 点值的估计。这种方法虽然最简单，但是误差较大，只有在 A，B，C，D 点的距离足 第 11 页够小时，才能实现比较精确的插值。在我们的试验中，正好满足这一条件，所以在后文可以看到我们采用最近点插值方法实现了较高精度的图像重建。在双线性插值法中，先需要估计 E、F 点的值，再求 P 点的值。其中 E，F 点可采用式（3-1）估计。 （3-*BAEEA FEB

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

49、行三次样条曲线插值，其次沿相角方向进行线性插值即可。如图 3-2 所示，设、分别为 A、B、C、D 沿径向方向的偏导ApBpCpDp数，则沿径向方向的点 E、F 估计值为： 23(3()2) +( 2() (3-3)EAABABABABAFFp EAFFppEAFFppEA23(3()2) +( 2() (3-4)FCCDCDCDCDCFFp EAFFppEAFFppEA然后利用式（3-2）沿相角方向进行线性插值得到 P 点的值。选取 Shepp-Logan 头部模型作为仿真原模型，投影数据为 367360，重建图像大小为256256。选用以上三种插值算法的 DF 重建图像以及局部放大图像请参

50、见图 3-3。其中图（a） 、 （b） 、 （c）分别是应用最近点插值、双线性插值以及混合二阶样条线性插值法重建的图像，图像（d） 、 （e） 、 （f）分别是对应的三种插值法重建图像的局部放大图。有图 3-3 的对比可以看出，由最近点插值得到的图像，质量最差，边缘的震荡最明显以及图像的混淆误差都明显比用双线性插值和混合二阶样条线性插值法重建图像误差大。而应用双线性插值法和混合二阶样条线性插值法重建的图像虽然在视觉上难以区分其图像质量，然而应用一些客观评价标准能够判别图像质量好坏。 （理论上，三次样条插值法在多数情况下都是一种比较好的插值方法，尤其对平滑的曲线插值8，因此，可以假定由混合二阶样

51、条插值的 DF 法重建图像的质量应优于用线性插值法重建图像的质量，然而在这里，由于频谱数据较为充足，因此没能显示其重建图像的质量优势。 ）在这里我们采用第 4.3 节中介绍的几种评价图像质量的客观标准作为评判这三种方法重建图像质量的标准。表格 3-1 显示了三种插值法重建图像的这几种判据大小。对于归一化平均绝对距离判据，双线性插值法在三种方法中表现得最出色。在我们的其他试验仿真中表明，当投影射线条数较少时，二阶样条线性插值法在三种插值法中表现得最好。由于我们的论文主要考虑在投影数据充足时的重建效果，因此，稀疏投影数据的仿真图没有收到本论文中。 第 12 页表表3-13-1 重建图像质量评价参数

52、对比重建图像质量评价参数对比插值方法最近点插值双线性插值二阶样条线性插值法归一化均方距离0.27460.23910.2398归一化平均绝对距离0.0520.03910.0395最坏情况距离1477 （a）最近点插值 （b）双线性插值 （c）混合二阶样条插值 （d）最近点插值局部放大 第 13 页 （e）双线性插值局部放大 （f）混合二阶样条插值局部放大图图 3-33-3 三种插值法重建图像以及其局部放大图三种插值法重建图像以及其局部放大图3.2 与卷积反投影算法等效的直接傅里叶重建插值方法由于卷积反投影法不但能重建出较好的图像质量，并且它的计算复杂度也在能够接受的范围内，因此几乎所有现代的 C

53、T 机中使用的重建算法都是卷积反投影法。虽然 DF 重建算法能够提供较好的计算效率，但是 DF 算法重建图像质量在大多数情况下都比 CBP 差。就像 3.1 节中讨论的几种插值方法，其重建图像都比 CBP 差。因此，可能有人会怀疑 CBP算法会比所有的 DF 重建算法要好，但是在 Hyeokho Choi 和 David C. Munson, Jr.的论文中，提供了一种插值方法使得 DF 算法与 CBP 等效，因此可以把 CBP 看作 DF 算法的一种有效实现。在这里，我将简单讨论对比这两种重建方法9。在平行束 CT 中，图像的投影是从各种不同方向采集到的。如果角度采集足够密集，就能忽略图像混

54、叠。因此在 CBP 法中，可以看作是极坐标下的二维傅里叶反变换，如（3-5）式所示： （3-5）20 0( , )( , )exp( cossin )Pf x yFjxydd 我们假设 f(x,y)在时域上限制在一个方形区域中，特别地： （3-( , )0 if or . f x yxayb6）因此可以把公式（3-5）改写为 （3-0( , )( , )exp(),pf x yFj t dd7）其中。在现实中，由于只能在离散点取得，不妨假设cossintxy( , )pF能在，并且，其中，( , )pF,.pp ,0,1,.,1qqQ /Q那么我们把式（3-7）离散化可以得到 （3-10( ,

55、 )(,)exp(),Qpqpf x yFpqpjpt 8）定义 fn,m = f(nx,my)为 f(x,y)的采样，并且 n = -N,0,N; m = -N,0,N. 可以 第 14 页求得 fn,m的近似值，计算式如下：10 , ,(,)exp(cos()sin().QPpxyqpPf n mf n mFpqjp nqmq （3-9）由公式（2-2）可知 CBP 法是实现（3-5）式重建图像的一种方法，而式（3-9）是根据式（3-5）代数变换得到的，因此可以看出由式（3-9）重建图像与 CBP 法重建图像等效。因为在 DF 重建算法中，最后一个步骤是对频谱数据进行二维离散傅里叶反变换。

56、由于式（3-10）式插值得到的是式（3-9）中的二维 DFT。因此可以得出用式（3- , F s r ,f n m10）进行插值得到的频谱数据再进行二维离散傅里叶反变换重建得到的图像与 CBP 法重建图像等效。1022 , ,exp(21212 =(,)exp(21(21) cos()exp2(21)2 (sin()212NNnN m NQPNpqpPnNNxmNyF s rf n mjsnNNnFpqpjNNspqNmjrpqN . (3-10) 由式（3-10）我们可以得到按以下两个步骤重建图像与 CBP 算法重建图像等效9：步骤一：利用式（3-10）对傅里叶频谱数据进行极坐标到直角坐标系

57、的插值；步骤二：计算插值好的频谱数据的二维离散傅里叶反变换。显而易见，经过式（3-10）插值的 DF 重建算法与 CBP 算法中采取矩形窗的滤波函数等效。同时，也很容易按照上面的推导方法找到与 CBP 算法中采取其他窗函数的滤波函数等效的插值方法，只是插值函数的加权系数不会是雅克比加权系数。当然，两种方法也有微小的不同。在 CBP 算法中，在反投影的步骤中需要对滤波后投影进行插值，而应用式（3-10）进行插值的 DF 重建算法不需要进行插值。因此，确切地说，应用式（3-10）进行频域插值的 DF 重建算法与在反投影步骤中对滤波后投影进行精确插值的 CBP 算法等效。按照式（3-10）重建图像的

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

59、10）中影响较大的均为靠近的点。在文献9中表明，对式（3-10）利用 , F s r ,f n m , F s r166 的窗函数时，重建效果与式（3-10）插值效果类似。虽然对插值器进行加窗处理后计算量大大降低，同时对图像质量影响不大，但是与 CBP 法相比，即使是加窗的插值器也没有比 CBP 法节省计算量。因此式（3-10）的最大意义在于证实存在一种非最优的插值器使得 DF 重建法与 CBP 法等效，也就表明了如果利用最优的插值器有可能使得 DF 法重建图像质量要优于用 CBP 法重建图像质量。文献7中同时也表明利用这种插值法的 DF 法重建图像之所以复杂度比 CBP 法要高，很可能是由于

60、这种方法没有考虑到极坐标中频谱数据的网格特性，而 CBP 法利用了极坐标的网格特性。 第 15 页 （a）CBP 法重建图像 （b）等效 DF 法重建图像 （c）CBP 法重建图像局部放大图 （d）等效 DF 法重建图像局部放大图图3-43-4 CBPCBP法与等效法与等效DFDF法重建图像以及其局部放大图法重建图像以及其局部放大图表表3-23-2 CBPCBP法与等效法与等效DFDF法重建图像质量评价参数对比法重建图像质量评价参数对比重建方法CBP等效 DF 法归一化均方距离0.24240.2336归一化平均绝对距离0.03880.0358 续表 3-2重建方法CBP等效 DF 法最坏情况距

61、离983.3 一种新的高精度 DF 重建算法3.3.1 投影数据的获得实际的投影数据是由检测器测得的，在计算机模拟时，则取自仿真模型，在我们的试验中，采取 Shepp-Logan 模型。由于 Radon 在他的经典论文中描述的问题由投影重建图像，对科学的众多领域产生了重大影响。因此，人们把联系目标函数空间到投影函数空间的积分变换称为 Radon 变换10。Radon 变换的公式如式（3-11）所示： 第 16 页 （3-( , )( , ) ( cossin)p xf x yxyx dxdy11）其中，为 Dirac delta 函数。叫做(,)x 0, )(,)( , )p x沿方向的投影。

62、如果注意到坐标与它的旋转坐标的关系，即式（3-( , )f x yxoyx oy12） （3-cossinsincosxxyy12）那么式（3-11）的 Radon 变换式可以写作式（3-13）： （3-( , )( , )(cossin ,sincos )p xf x y dyf xyxydy13）其中，。式中，和是定义平面上一条直线的位置的两个(,)x 0, )x变量，该直线由方程决定，对一个固定的和，这个直线方程决cossinxxy x定了一条与轴平行的直线，为沿该直线的线积分。由式（3-12）可以y( , )p x( , )f x y发现，空间的一个点对应于空间的一条正弦曲线段11。(

63、 , )x y( , )x在数学上 Radon 变换是一个线积分变换，它可以看做是一个算子。该算子对空( , )x y间的函数作用产生空间的另一个函数。与 CT 问题联系起来，Radon( , )f x y( , )x( , )p x在不同的 CT 系统中具有不同的物理含义。称作被观测目标的表征函数，( , )f x y则是与这一表征函数有关的测量值。在试验中模拟投影数据有一般两种( , )p x( , )p x方法，一种通过解析法计算，另一种通过数值累加法计算。解析法计算实际上就是利用式（3-11）对模拟图像求其 Radon 变换；而通过数值累加法求其 Radon 变换需要做某些方面的近似。

64、为了对 CBP、DF 以及 ART 算法进行对比，在我们的试验中，采取数值累加法求投影数据（因为在 ART 算法中，需要对算法开始即对原图进行离散化处理，故需要采用数值累加法求投影数据） ，这样三种方法都可以采用相同的投影数据，使得重建图像更具有可比性。在我们的仿真实验中，图像在某个方向的投影数据是所有像素点在此方向投影数据之和。因此我们采取的待重建图像为数字图像（在实际的重建过程中，也可以把原图采样为数字图像，因为重建后的图像是天然离散的） 。图图 3-53-5 投影示意图投影示意图如图 3-5 所示，每个像素点的投影都按比例分到离它最近的两个投影点上（也可以更精确的把每个像素再平均分为多个

65、小的像素的，然后再求投影） 。虽然在 CT 机上采集的数据方法与此不同，但是这种方法较精确的模拟了 CT 投影的采集值。按照这种算法思想，我们可以在每个投影角度遍历原图的每一个像素点，得到每个像素点在每个投影角度的投影值，累加之后就可以得到每个投影角度的投影值12。 第 17 页图图 3-63-6 投影数据图像投影数据图像图 3-6 为修改过的 Shepp-Logan 模型的投影数据，其中原图为 256256,投影数据为367180，在本文后面三种方法的对比中，采用的投影数据均采用这种方法获得，并且投影数据都为 367360.3.3.2 投影数据一维傅里叶频谱数据的计算对投影数据频谱值的取得，

66、我们采取式（2-14）计算。观察式（2-14） ，在我们的实验中获取的投影数据，可以令t = 1, d = 1/(128*N),其中N为重建图像大小，在我们的试验中N均取 256。重建图像时，d取值越小，重建图像质量越高，但是图像重建速度也越慢。因此代入数据式（2-14）可以化为式（3-14）的形式 （3-1/22/(128)11/2 , , 64,64)MinnNnMP n mp n m enNN 14）式（3-14）实际上是对序列求解其线性调频 z 变换（CZT） 。由 CZT 的性质， ,p n m我们可以知道当 CZT 变换路径为单位圆上一段圆弧时，能够得到信号的频谱分析13。在这里我们没有采用 DFT 作为投影的频谱分析，因为 DFT 对频段的抽样密度不够，这样导致插值误差较大。而 CZT 能够满足我们对提高频段抽样密度的要求，因此我们选择使用 CZT 作为对数据的频谱分析。对每一个角度的投影数据，按照式（3-14）求得其频谱数据。根据投影定理，每个角度投影数据的傅立叶变换对应原图二维傅立叶变换的一个切片。把所有角度的频谱数据组合起来，也就得到了原图二维傅立叶变换的很多切片，