（信号与信息处理专业论文）三维锥束ct中滤波反投影算法的研究.pdf

中北大学硕士学位论文 三维锥束c t 中滤波反投影算法的研究 摘要 近年来，三维锥束c t 在医学和工业上得到越来越广泛的关注，同时锥束投影图像重 建算法也在迅速发展。在各种三维锥束重建算法中，滤波反投影算法仍然是一类最常用 的重建算法，在医学和工业检测上有着重要的意义。 本论文针对三维锥束c t 中两种滤波反投影算法，以提高重建质量，降低重建图像 的噪声为目的，分别对其中的f d k 算法和k a t s e v i c h 算法进行了一些研究和探索。 在本文中首先分析了f d k 算法，研究了其性质，并进行了仿真验证。然后，在滤 波函数方面，采取了一种介于r - l 和s l 之间的混合滤波函数，由实际的重建结果可知， 该滤波函数的重建结果无论在横向还是在纵向上的分辨率都比较高，并且在空间分辨率 和时间分辨率上都有较好的效果。接着，在插值函数方面，针对f d k 算法的反投影过程 中，各个体素在探测器上投影分布的特点，提出了一种新的插值方法。实际实验结果表 明，在该算法中这个新的插值方法能更好地抑制噪声。 最后，针对螺旋扫描轨迹的精确重建算法叫a t s e v i c h 算法的公式，对其步骤( 投 影数据求导，滤波过程和反投影重建物体三个步骤) 进行了一定的分析和研究，并做了 仿真实验，为以后相关工作打下了基础。 关键词：f d k 算法，性质，混合滤波函数，插值，k a t s e v i c h 算法 中北大学硕士学位论文 r e s e a r c ho nf i l t e r e db a c kp r o je c t i o na l g o r i t h m i n3 dc o n e - b e a mc t a b s t r a c t i nr e c e n ty e a r s ，t h et h r e e - d i m e n s i o n a lc o n e - b e a mc th a v em o r ea n dm o r ew i d e s p r e a d a p p l i c a t i o ni nt h em e d i c i n ea n dt h ei n d u s t r y , s i m u l t a n e o u s l yt h ec o n e - b e a mr e c o n s t r u c t i o n a l g o r i t h ma l s oh a v et h er a p i d l ye x p a n d a m o n ga l lk i n d so fr e c o n s t r u c t i o na l g o r i t h m sf o r3 d c o n e - b e a m c t ，f i l t e r e db a c kp r o j e c t i o na l g o r i t h mi sat y p em o s tc o m m o n l yu s e d r e c o n s t r u c t i o na l g o r i t h m ，s oi ti sv e r yi m p o r t a n ti nm e d i c a la n dn o n d e s t r u c t i o n b a s eo nt h es o m ef i l t e r e db a c kp r o j e c t i o na l g o r i t h mi nt h r e e - d i m e n s i o n a lc o n e - b e a mc t ，i n o r d e rt oi m p r o v et h ei m a g er e c o n s t r u c t e da n dr e s t r a i ni t sn o i s e ，w ed i ds o m es t u d i e sa n d e x p l o r a t i o n sa b o u tt h ef d ka l g o r i t h ma n dt h ek a t s e v i c ha l g o r i t h mi nt h i sp a p e r f i r s to fa l l ，w ec o n d u c t e ds o m er e s e a r c hi n t h ef i l t e rf u n c t i o na n dt h ei n t e r p o l a t i n gf u n c t i o n f i r s to fa l l ，w ed r a w no u tt h ef d ka l g o r i t h mn a t u r ea n dc o n f i r m e dt h e s en a t u r eb yt h e s i m u l a t i o nr e s u l ta n dm a t h e m a t i c si n f e r e n t i a lr e a s o n i n g t h e n , i nt h ef i l t e rf u n c t i o na s p e c t ，w e p u tf o r w a r dam i x e df i l t e rf u n c t i o n , w h i c ha r ei n t e g r a t i o no ft h er - la n ds lf i l e rf u n c t i o n f r o mt h er e s u l t so ft h ea c t u a lr e c o n s t r u c t i o n , w ek n o wt h a tt h er e s u l t sf r o mt h i sk i n df i l t e r f u n ，c t i o nh a v eah i g h e rs p a c er e s o l u t i o ni nb o t hh o r i z o n t a la n dv e r t i c a l ，a n di nt h es p a t i a la n d t e m p o r a lr e s o l u t i o n ，i th a sb e t t e rr e s u l t s t h e n , i ni n t e r p o l a t i n gf u n c t i o na s p e c t ，i nt h e b a c k - p r o j e c t i o np r o c e s so ft h ef d ka l g o r i t h m ，e v e r yv o x e l sp r o j e c td i s t r i b u t i o ni nt h ed e t e c t o r h a v ei t s e l fc h a r a c t e r i s t i c s ，w ep r e s e n t san e wi n t e r p o l a t i o nm e t h o d ，a n da c t u a le x p e r i m e n t a l r e s u l t ss h o wt h a tt h ef d ka l g o r i t h mu s et h i sn e wi n t e r p o l a t i o nm e t h o dc a nr e s t r a i nn o i s e e f f e c t a tl a s t ，t h i sp a p e ra l s oh a sc a r r i e do nc a 恤a n a l y s i sa n dt h er e s e a r c ht ot h eh e l i c a ls c a n p a t h sp r e c i s er e c o n s t r u c t i o na l g o r i t h m - 一k a t s e v i c ha l g o r i t h m sf o r m u l aa n dt h es t e p ( i n c l u d i n g d e r i v a t i v e so f c o n e - b e a md a t a ，f i l t r a t i o no f d e r i v a t i v ed a t aa n db a c k p r o j e c t i o no f l a s td a t a ) ，a n dh a sd o n et h e s i m u l a t i o ne x p e r i m e n t i th a sb u i l tt h ef o u n d a t i o nf o rt h el a t e rr e l a t e dw o r k k e yw o r d ：f d ka l g o r i t h m , n a t u r e , m i x e df i l t e rf u n c t i o n ，i n t e r p o l a t i o n , k a t s e v i c ha l g o r i t h m 原创性声明 本人郑重声明：所呈交的学位论文，是本人在指导教师的指导下， 独立进行研究所取得的成果。除文中已经注明引用的内容外，本论文 不包含其他个人或集体已经发表或撰写过的科研成果。对本文的研究 作出重要贡献的个人和集体，均已在文中以明确方式标明。本声明的 法律责任由本人承担。 论文作者签名：斗日期：上_ 盘查且灶 关于学位论文使用权的说明 本人完全了解中北大学有关保管、使用学位论文的规定，其中包 括：学校有权保管、并向有关部门送交学位论文的原件与复印件： 学校可以采用影印、缩印或其它复制手段复制并保存学位论文； 学校可允许学位论文被查阅或借阅；学校可以学术交流为目的，复 制赠送和交换学位论文：学校可以公布学位论文的全部或部分内容 ( 保密学位论文在解密后遵守此规定) 。 签名： 导师签名： 日期：磁堡竺旦sll i l 日期：3 霜j ( 笺生国三! 国 中北大学学位论文 1 1 课题研究背景 第一章引言 该课题来源于国家自然基金项目“三维c t 的信息获取与处理理论和方法”。x 射线 断层成像技术( c o m p u t e r i z e dt o m o g r a p h y 称c t ) 的物理原理是：x 射线穿过物体时会 发生强度衰减，衰减的程度和物体的密度、原子序数以及入射射线的能量和初始强度有 关。利用这个原理，让x 射线从不同的角度穿过被检测物体并被探测器接收，再利用 各种不同的算法重建出被检测物体各处的射线衰减系数并得到其分布图像，以此反映被 测物体的密度分布。c t 的断层图像是由计算机对物体的投影数据进行计算而获得的， 因此c t 的成像过程可以分为两大步，一是进行投影数据的测量，获得充分的投影数据； 二是用计算机对投影数据进行计算来获得断层图像。x 光的投影测量结构可以分为三 种，分别为平行束投影、扇束投影和锥束投影。平行束投影在工程上是不能直接实现的。 和扇束投影相比，锥束投影一次性获得的信息量更大，因此锥束投影结构有利于提高扫 描速度和图像质量。三维锥束c t 拥有数据采集时间短，射线利用效率高，一次扫描就 能对物体进行三维重建并且可以保证在各个方向获得同样的分辨率等优点，另一个显著 优点就是可以很容易地得到被测对象的任意角度的切片图像 1 】【2 1 。因此，三维锥束c t 已成为目前的研究热点。 现有的重建算法可分为两大类，即解析法和迭代法。而在所有的重建算法中，最实 用的还是滤波反投影算法【3 - 7 。三维锥束中的滤波反投影算法又分为近似算法和精确算 法，它们在各自的领域有着重要的应用价值。 1 2 国内外研究现状 早在1 9 6 0 年，数学家k i r i l l o v 已在文献8 1 中提出了从复数值的锥束投影数据重建 n - 维复数值函数的理论，从而为真正意义上的三维重建技术奠定了理论基础。后来， s m i t h 将k i r i l l o v 的公式改写到实数空间中并且提出了用于无限长扫描轨迹的反变换公 1 中北大学学位论文 式【9 】，在几乎每一个通过物体的超平面都包含一个源的情况下，t u y 推导了一个用于两 个垂直扫描轨迹的实函数重建公式 1 们。锥形束三维重建算法可分为近似重建和精确重建 两种【1 1 1 。其中，以f d k 算法【1 】【1 2 】为代表的近似重建算法，由于具有源点轨迹简单( 不 要求满足重建的完全性条件) 、可以方便地进r o i ( r e g i o no fi n t e r e s t ) 重建以及利于 并行处理等优点，而因此在医学成像和无损探伤领域有着广泛的应用前景。精确重建算 法，由于源点轨迹要求满足重建的完全性条件，其源点轨迹必须为垂直双圆或圆加直线 等非平面轨迹，因此在处理长物体的重建时，物理实现较为麻烦。 1 2 i 近似重建算法 近似锥束重建在实际中有着重要的作用，它的主要优点为：可使用不完全的扫描轨 迹、可允许射束锥角覆盖很少一部分的重建对象、只需要维滤波而不是二维滤波操作。 h 锄姐【1 4 1 、l e g i t 和m c k a y t l 5 1 曾直接利用扇束算法进行锥束扫描的三维重建，当锥束 角较小时，重建结果很好。g r a n g e a t 1 6 】【1 7 1 算法利用二维投影数据估计三维逆雷当变换， 也得到了较好的重建结果。在众多的近似算法中，其中f e l d k a m p ，d a v i sa n dk r e s s 在1 9 8 4 年提出的f e l d k a m p 1 8 】算法最为成功，该类型的重建算法也一直是实际应用中的主流， 在小锥束角的情况下，该算法重建结果可以与二维扇束重建图像媲美。但在大锥角情况 下远离中心层的重建图像误差很大【1 9 】【2 0 1 ，会产生纵向模糊。经典的f e l d c a m p 算法要求 扫描轨迹为圆形，g u l l b e r g t 2 1 】【2 2 】、王革【2 3 1 等推广了f e l d c a m p 算法，使之可应用于各种 扫描轨迹，可对球状、棒状的样品进行重建并得到良好的重建图像质量。如今，f d k 算法已经被广泛采用，而且得到了很大的发展，大家在此基础上提出了很多的改进算法 【2 铊7 】。 2 1 2 精确重建算法 近似算法有着不可逾越的精度问题，所谓的“长物体问题( 1 0 n go b j e c tp r o b l e m ) t 2 s l 也严重地限制了该算法的发展。与此同时，精确重构的理论研究硕果累累。t u y 建立了 实现精确重构的锥形束投影轨道的充要条件2 9 1 ，g m n g e a t 推导了锥形束投影与三维 2 中北大学学位论文 r a d o n 变换的关系【3 们，建立起锥形束精确重构坚实的数学理论基础。该领域的开拓者之 一b d s m i t h 对重构条件做出了几何解析，并描述了实现重构的滤波反投影算法系统【3 l 】。 然而锥形束重构不仅是一个计算机断层成像的工程问题，也是一个较为复杂的弱病态求 逆的数学问题吲。s m i t h t 3 3 1 ，g r a n g e a t 3 4 】； 年l t u y t 3 5 1 在1 9 8 4 年前后分别提出t 3 种类型的 锥束精确重建算法。至此，奠定了锥束重建算法的理论基础。后来，研究人员又基于这3 种算法发展了多种算法。在此基础上，2 0 0 2 年k a t e s v i c h 提出了基于螺旋轨道的移不变的 滤波反投影算法【3 6 1 ，并进一步将其推广到普适轨道的通用系统【3 7 1 。锥形束重构研究由此 进入新阶段。但是，由于精确算法的数学复杂及计算量大，很难得到实际应用。在这些 诸多算法中，美国芝加哥大学以潘晓川教授为首的课题组提出了一种重要的b p f 精确重 建算法【3 8 - 4 0 l 该算法是继k a t s e v i c h 算法【3 4 】1 4 1 1 以后的另一个重要的锥束c t 精确重建算法。 b p f 算法具有鲜明的数学和物理上的特点，该算法的数学证明和推导自成体系，同时具 有直观、明了的物理概念支持。另外，该算法虽然是针对螺旋锥束c t 提出的，但由于 b p f 算法本身的特点，使得该算法能够方便、直观地扩展到其他形式的锥束扫描轨道的 c t 精确重建中【4 2 】【4 3 1 。作为锥束重建的一种特殊情况，b p f 算法也被应用到扇束、平行束 c t 精确重建中m 】。b p f 算法的最大优点在于它第一次在数学上解决了横向截断的投影数 据的c t 精确重建问题，能够利用理论上最少的投影数据精确重建出物体图像。 1 3 论文的主要内容 出于对工程实现难度和算法计算量的考虑，采用的主要是数据不完备的基于圆轨迹 扫描的f d 幻匠似三维重建算法。本论文主要是针对该算法中的滤波部分的滤波函数和插 值部分的插值方法进行的研究，同时也对精确重建中的k a t s e v i c h 算法进行了一定的探索。 由于这一目的，论文的主要工作由以下几个部分组成： l 、简单地介绍了标准f d k 算法以及其性质，并由仿真的结果和数学上的推导来验 证或证明了该算法的性质。 2 、在f d k 算法中，首先对真实的投影数据分别用r - l 和s l 滤波重建，然后采取 了一种介于r - l 和s l 之间的混合滤波函数，由实际的重建结果可知，该滤波函数的重 建结果无论在横向还是在纵向上的空间分辨率都比较高。最后分析了这种滤波函数的优 3 中北大学学位论文 越性。 3 、提出了f d k 算法中的一种新的插值方法。该方法根据体素投影的特点，采用在 重建过程中，考虑在不同扫描角度下待重建体素落在探测器上一定的投影面积，并且在 反投影重建这个体素点的时候，考虑那些探元上的信息量。实际实验结果表明，在f d k 算法中这个新的插值方法比传统的插值方法( 如：最近邻插值，双线性插值) 重建出来 的图像边缘清晰，而且能更好地抑制噪声。 4 、对精确算法中的k a t s e v i c h 算法进行分析，主要是从算法的公式以及运行步骤上 研究和总结归纳，并得到了一些有用的结论，接着对模拟数据进行了仿真试验。 5 、分析了在k a t s e v i c h 算法中的几个问题，如p i 线的求解方法等，最后比较了这 两种算法( f d k 算法和k a t s e v i c h 算法) 的内在联系和区别。 4 中北大学学位论文 2 1f i ) k 算法的公式 第二章f d k 算法简介 c t 锥束重建算法发展至今，近似算法由于数学上形式简单，并且实现起来容易， 所以是实际c t 系统中的主流算法。f d k 算法是圆轨迹锥束重建中用处最广泛也最实用 的重建算法，就是现在的f d k 算法，扫描结构如图2 1 所示。该算法实际上是二维扇 束反投影算法的三维推广，它首先对偏离中心平面的投影数据进行加权，然后利用类似 与滤波反投影的算法对物体进行重建。 笺 一羔7 图2 1 锥宋圆周扫猫结构不葸图 探测源的扫描轨迹是： c = y r 3 ；y 。= r c o s ( s ) ，y ：= r s i n ( s ) ，y ，= o ，j 【0 ，2 万) ) 其具体过程可以分为以下两步： 凇马功2 b 志砌刊宰矿国 z ) = r ”瓦8 习2 。歹以咖凇 5 ( 2 1 1 ) ( 2 1 2 ) ( 2 1 - 3 ) 中北大学学位论文 其中p ( 红a ，易) 是采集到的投影数据，g p ( 口) 是一维滤波器，r 是轨道半径，是扫描角 度，a 、b 是在探测器上的位置。 步骤的具体描述如下： ( 1 ) 首先利用类似于余弦的函数对投影数据进行加权，适当地修正体素到源点的距 离和角度差； ( 2 ) 然后对不同投影角度的投影数据进行水平方向上的一维滤波； ( 3 ) 最后沿x 射线方向进行三维反投影。重建的体素值是通过该体素的所有投影角 度的射线的贡献之和。 2 2 算法的性质 作为实际应用较为广泛的f d k 算法，有其自身的特点。在这一节里，引出了标准 f d k 算法的几个重要性质并说明了各个性质在实际中的应用价值。 性质1 在中心一层，即z - - 0 时，是精确重建。 t u y 、s m i t h 、g r a n g e a t 等基于各自的工作得出了一致的结论：如果与物体相交的每个平 面上至少存在一个锥束源点，则物体可被精确重建【州。标准的f e l d k a m p 方法就是基于 单圆的。然而，除扫描轨迹所在平面( 就是z = 0 面) 外的那些与物体相交的平面都不包 含锥束源点，因此，f e l d k a m p 方法是一种近似的重建方法。但对于中心平面的重建来说， 可以演化成为二维扇束重建，也就是精确重建。离中心平面越远的平面，锥角越大，重 建误差也越大。该算法的重建误差和轨道半径的关系见参考文献【4 5 】【蛔。 这个性质是对于中心平面( 即扫描轨迹的平面) 的一个说明。若在实际中想要得到 长物体的多层“精确 信息，可围绕重建物体的不同横断位置多扫描采集几次。诚然， 这样会使采集数据的时间和重建时间增加。 性质2 对z 方向一致的物体，即厂k y ，z ) = 厂g ，y ) ，f d k 是精确重建。 这里，主要是由模型【4 7 】仿真的结果验证的。其中，z 轴均匀物体的横截面的模型参 数见下表2 1 所示： 6 中北大学学位论文 表2 1z 方向均匀物体的仿真参数 设计虚拟探测器，将探测器镜像到物体中心的位置。令射线源到物体中心的距离和到探 测器的距离各是8 0 0 。 其中一个扫描角度下的投影值如图2 2 图2 2 某一角度下的投影值 对于2 5 6 2 5 6 2 5 6 的三维图象重建结果如下： 可见：中间层和非中间层的任一层重建结果是一致的。由性质1 直接可得f d k 是 精确重建。 这个性质是对特殊物体，当然，该物体只是理想的一个假设，在实际中很少存在。 若通过先验知识，已知了被重建物体的各个横切面都相似，重建一层即可得到整个物体 的图像。 7 中北大学学位论文 里 图2 3 ( a ) 均匀物体中间层重建图2 3 ( b ) 距离中间1 0 0 层的重建 置 ( a ) 过图2 3 ( a ) 中心的直线灰度图( b ) 过图2 3 ( a ) 中心的直线灰度图 图2 4 当x = 1 2 8 时过图2 3 所示切面的中心直线的灰度图 性质3 ：z 方向的积分少k y ，z 皿是保持不变。 设： p 慨口，6 ) 是在投影方向= ( c o s , 8 , s i n , a , o ) 0 0 点, c o = ( 蠢，y 。，z 0 ) 处的投影数据。 p = 胎忆，韶仁蒜南b ( 2 2 1 ) 由f d k 方法从二维投影p 饵口，6 ) 重建出的物体是： 图 中北大学学位论文 缸伍) = f 7 西考扣p ( 口一口b 协口，6 b 志妇妒 ( 2 z 2 ， 舯一尚 慨南4 - x r 8 即要证明肚k y ，z 皿= p g ，y ，z 皿 把( 2 2 1 ) 、( 2 2 2 ) 代入( 2 2 5 ) 中得 肚伍弛= 厅南g p ) 肌阢) 需 南a aa p a z - i - b = = = = = = = = = = = = = = = = r 2 + 口2 2 。 d ( a - r 熹南 r + x 。ji s ( x o ) 南巾 对z 的积分 f 纠6 一 i z o r r 七x o p z o r r + x n 芦 卜h z 南一 尺2 z o r r + x q p 妣赢妒 z o 尺 r + x a 8 9 卜砜 ( 2 2 3 ) ( 2 2 4 ) ( 2 2 5 ) ( 2 2 6 ) ( 2 2 7 ) 1 d z ：墨垄：丝 ( 2 2 8 ) jr ( 2 2 9 ) 一 罗一 鬻 一尺 弦 万若 得 比 勰 啦 榭 舡 。口 村 肠 一 岁一 等舞 一伍 为0 变 弘 皿 7 x z “ 幺 肋 汜 肠 式 中北大学学位论文 最后对投影角度积分 肛伍k = 少阢) m j c 0 ) 础一y 。) 咖。d z 。= p 伍协 得证。 这个性质从数学上证明了：无论扫描轨迹的半径如何改变，用f d k 算法重建出来的 物体在z 方向上的积分和真实物体的在z 方向上的积分值是一致的。这就说明，重建结 果在整体上保持了被重建物体的信息，f d k 算法只是对物体内部信息的重新分配。 2 3 本章小结 通过上面对基于圆轨迹扫描方式的f d k 算法性质的深刻研究，进一步了解了f d k 算法。这种方法具有简单，有效，快速的特点，但它只能近似的重建物体。之所以说 f d k 算法是一种近似的算法，是因为无论测量时的分辨率如何，重建结果和实际情况或多 或少都会有偏离。在锥角比较大的情况下，f d k 算法的伪影就变得很严重了，即重建 结果与实际情况的偏差很大。 1 0 中北大学学位论文 第三章混合滤波函数在f d k 算法中的应用 滤波反投影重建算法的实现，理论上要求滤波器的系统函数日= 纠。这是一个 频带无限的滤波函数。因用日纠2 咖= 月刊2 咖j 按佩利一维纳准则4 3 1 ，这一理想 滤波器是不能实现的。 3 1 滤波函数的选择 滤波函数的选取若一味的追求理想函数则不能实现。但若结合具体成像过程考虑， 则滤波函数既可以实现，又可有足够的精度。一般成像过程的具有如下特点： 1 投影数据的高频( 空间域频率) 分量幅度很小； 2 投影数据是离散采集的； 3 存在噪音。 所以滤波器可以通过对理想的斜坡函数日如) = i 爿，加上一个频率范围在h b 的窗函 数来实现。则滤波器的频域表示形式为：日= h 形式中形窗函数。所以 使用不同的窗函数，就可以得到不同的滤波器1 4 9 由此看来，滤波器的选择实际上是窗函数的选取，为了得到较好的重建图像分辨率， 窗函数的傅立叶逆变换应当是具有一个高而窄的中央突出部分，而侧面突出部分引起的 吉布斯现象，被看作是由重建图像相邻部分的沾染性歪曲，通常，有高而窄的中央突出 ( 即主瓣) 的窗函数，其旁瓣也突出，也就是说，若一味地要求分辨率高会引起较严重 的吉布斯现象。因此，窗函数的选择取决于实际构件的内部结构和重建目标。 一般为了考虑其性能的优越，通常要求窗函数的幅频具有以下特征： ( 1 ) 主瓣宽度窄，以获得较陡的过渡带； ( 2 ) 最大旁瓣相对主瓣尽可能小，以改善通带的平稳度和增大阻带的衰减。 工业c t 图像质量都是通过密度分辨率与空间分辨率来评价的，因此，两者都是影响 成像的主要因素【5 0 】【5 l 】。而滤波函数对重建误差的影响主要有以下三个因素：主瓣，临近 1 1 中北大学学位论文 旁瓣，及远处旁刺5 2 】【5 3 】。主瓣越高越窄，空间分辨率越好；旁瓣越小，数值精度越高， 通常密度分辨率就越好。近年来，不少学者提出了许多滤波函数，并对一些滤波函数做 了分析和改进【蚓。 在通常情况下，我们用到的滤波函数有r l 滤波函数嘲、s - l 滤波函数5 6 1 。近年来， 不少相关人员提出了许多滤波函数如h a r m i n g 滤波函数、h a m m i n g 滤波函数等等。 3 2 滤波函数 在确定了实际构件的内部结构和重建目标时，常用到下面两种滤波函数。 3 2 1r - l ( r a m a c h a n d r a n - l a k s h m i n a r a y a n a ) 滤波函数 该函数的采样序列为： h r - l ( 刀) = 三 刀：0 3 刀为偶数 ( 3 2 1 ) 一1 l 一刀为奇数 玎万 用这种滤波函数计算简单、实用，避免了大量的正弦、余弦计算，得到的采样序列是分 段性的，并没有明显地降低图像的质量，所重建图像轮廓清楚，空间分辨率高，但有 g i b b s 现象，表现为明显的振荡响应。 3 2 2s - l ( s h e p p - l o g a n ) 波函数 该函数的采样序列为： 五“0 ) = 硐- 2 ( 3 2 2 ) 该滤波函数比较平滑，相对于r l 来说能更快地完全趋于零，因此该滤波器对投影的中 高频成分具有抑制作用，从而重建图像的振荡响应减小。但是其分辨率不高。 1 2 中北大学学位论文 3 2 3 混合滤波函数 鉴于以上两种滤波函数有各自的优缺点，为了根据具体问题，提高空间分辨率的同 时并对噪声有很好的抑制，可将上述提到的两种常用的滤波函数同时使用，也就是根据 具体频率响应函数可以考虑他们的加权和，将r - l 和s l 滤波函数各取其一部分，可能 会发挥各自的优点，或许可以使重建图像提高分辨率的同时减少图像的振荡现象。我们 称之为混合滤波函数【5 刀，其函数的采样序列时是： j i l g ) = k l h “( n ) + k 2 h “g ) ( 3 2 3 ) 也就是： j i l g ) = 其中：k l + 后2 = 1 竺堡堕 4 8 2以= 0 确：裟 2 4 ) 一万一刁网 显然可见：当岛= o 时，该函数是s - l 滤波函数；当哎= 0 时，该函数是r - l 滤波函数。 该滤波函数相当于对r - l 和s l 滤波函数的一个折中，也就是与对低频和高频信息的一 个折中。 3 3 实验结果 实验采用2 1 5 k v ，2 5 m a 的x 射线源。单个探测器大小为0 1 2 7 m m 。采用的探测 器为p a x s c a n 2 5 2 0 。工作模式：数字视频。数据类型：u n s i g n e ds h o r t 。a d 1 2 b i t 。 射线源检测工件探测板距离：8 5 0 m m 一1 5 0 m m ，探测器大小2 5 6 * 2 5 6 。旋转一周采样 间隔为1 度，投影数3 6 0 个，被测工件的成分为铝合金。图3 1 是是被测工件在某个角 度下的投影。大小为2 5 6 2 5 6 。下面采用了上节中所提到的三种滤波函数进行滤波， 并对其结果进行了比较。 1 3 中北大学学位论文 图3 1 某一角度下的实际投影 ( 1 ) 在对上述投影进行重建时采用r l 滤波函数结果如图3 2 ： ( a ) 图3 2 ( a ) 和( b ) 分别是用r - l 滤波函数三维重建后z = 3 8 和x = 1 2 8 切面图像 从该重建后的图像中可以看出，横切面( a ) 与纵切面( b ) 的质量都不够满意，有 明显的振荡响应，表现为g i b b s 现象。这是因为在某一点处是理想的矩形，而理想的矩 形窗是产生g i b b s 现象的根源所在。 ( 2 ) 改用s 。l 滤波函数对上述投影进行重建，结果如图3 3 ： 从该重建后的图像中可以看出在横切面( a ) 和纵切面( b ) 上比较平滑，振荡响应 减小，对含有噪声的投影数据有很好的抑制，但是该图像的分辨率不高。 ( 3 ) 以上两种的重建结果都达不到理想的效果，下面改用文中提到混合滤波函数，对 于混合滤波函数中的两个参数k 和也的确定，采用比较得到，当k 。= 0 8 5 ，k 2 = o 1 5 时 重建结果的效果最好，结果如图3 4 ： 1 4 中北大学学位论文 ( a ) ( b ) 图3 3( a ) 和( b ) 分别是s l 滤波函数三维重建后z = 3 8 和x = 1 2 8 切面图像 图3 4 ( a ) 和( b ) 分别是混合滤波函数三维重建后z = 3 8 和x = 1 2 8 切面图像 从该重建结果中，可以看到图像在噪声和分辨率上都有所改善。 在混合滤波函数的实际应用中，k l ，k ：的值该如何确定，为什么在k 。= 0 8 5 ， k ：= o 1 5 时，图像的重建效果较好，可根据当k l , k ：分别取不同值时的滤波函数在各自 空域的主瓣宽和幅值以及第一旁瓣的幅值例大小( 见表3 1 ) 来说明。 1 5 中北大学学位论文 表3 1 主瓣宽和幅值以及第一旁瓣的幅值 这里，我们从判别重建图像质量的空间分辨率和密度分辨率两个指标上来分析。空 间分辨率是指在c t 图像中能够辨别的最小物体的能力。密度分辨率是指利用图像灰度 取分辨被检物质的基础方法。在投影数据一定的情况下，二者是相互矛盾的指标。通常， 主瓣高且窄，图像的空间分辨率就高，如r l 滤波函数的重建结果；旁瓣衰减迅速，可 提高其密度分辨率，如s l 滤波函数的重建结果。而混合滤波函数的效果介于两者之 间，对该物体的重建，在k l = 0 8 5 ，k 2 = 0 1 5 时，滤波函数得到了较好的折中。 由上面的三种滤波函数的比较表明混合滤波函数在图像重建过程中运用，对图像的 噪声有着很好抑制，同时对图像的分辨率有很大的改善 3 4 本章小结 综上所述，关于r l 和s l 相互结合的混合滤波函数对重建图像的噪声有很好的抑 制作用，而且重建图像的x y 平面或是在与z 轴平行方向在空间分辨率和密度分辨率上 有较大的提高。并通过实验证明了该方法的有效性及可行性。但是关于系数毛和如的 确定要根据具体情况来确定。 1 6 中北大学学位论文 第四章f d k 算法中一种新的插值方法 什么是插值? 文献【5 9 】倾向于下面的定义，在己知的坐标范围内，一种基于模型的从 离散数据估计连续数据的方法。 插值有三个重要假设： ( 1 ) 潜在的数据是连续定义的； ( 2 ) 给定数据采样值，就可以计算连续函数在任何坐标点的值； ( 3 ) 连续函数在采样点的值等于给定的数据采样值。 4 1f o k 算法中的插值 在f d k 算法的反投影过程中，由于数据的离散性，会出现象素的投影地址“对不 准 现象，一般需要插值运算。这里的插值是指在己知的坐标范围内，一种基于模型的 从离散数据估计连续数据的方法。经典的线性插值技术包括最邻近插值( r e p l i c a t i o n ) 6 0 】， 双线性插值( b i l i n e a r ) 【6 1 1 ，双三次( b i e u b i c ) 插值删t 6 2 - 6 5 等。特别是对于图像的边缘。最 邻近插值会产生块状和锯齿，而双线性插值和双三次插值则会使图像的边缘变得模糊。 于是，研究者引入了许多非线性的插值方法，如局部适应插值方法【蛔【6 7 1 ，b a y e s i a n 估计 1 6 8 笺 -jo 4 2 常用插值函数 一般认为，像素的反投影值是由点的插值取得。在本节中，我们考虑上诉投影图形 的加权值。任取某一待重建体素上的顶点k y ，z ) ，为扫描角度，l s o 为探源到物体中 心距离，l s 0 2 为探源到探测器距离，那么它在探测器上的落点p 的坐标【6 9 1 ： x = ( 1 s s s + y s i n ( x s i n 一y c o s ( 4 2 1 ) y = ( 1 s 0 2 ( 1 s o x c o s c 功+ y xs i n ( f 1 ) ) ) xz；(422) 1 7 中北大学学位论文 、 因为得到的g ，y ) 不一定是整数，所以需要通过插值得到其对应的反投影值。 在这一节中将分别介绍其中的最近邻插值和双线性插值这两种方法。 4 2 1 最近邻插值 这是一种最简单的插值方法，又可以称为零阶插值。即将得到的g ，j ，) 用离它最近 的整数g ，夕) 处的投影值厂g ，少) 来代替。若令厂g ，y 7 ) 为两个变量的函数，其在单位正 方形定点的值厂伉多) ，厂 + 1 ，夕+ 1 ) ，厂g + 1 ，夕) ，厂仁夕+ 1 ) 己知。假设我们希望通过 插值得到正方形内任意点b 7 ，少) 的值厂g ，少) ，就需要双线性插值。 4 2 2 双线性插值 双线性插值又称为一阶插值，它的效果要好于最邻近插值算法。只是程序相对复杂 一些，运行时间稍长一些。由于通过四个点确定一个平面是一个过约束问题，所以在一 个矩形栅格上进行的一阶插值就需要用到双线性函数。 我们可由双线性方程： g ，y ) = 戤+ 砂+ c r y + d ( 4 2 3 ) 来定义的一个双曲抛物面与四个已知点拟合。 从a 到d 这四个系数需由己知的四个顶点的厂k y ) 值来选定。有一个简单的算 法可以生成一个双线性插值函数，并使之与四个顶点的厂kj ，) 值拟合。首先，我 们对上端的两个顶点进行线性插值可得 g ，o ) = s ( o ，o ) + x l s o ，o ) 一s ( o ，o ) 1 ( 4 2 4 ) 类似地，对于底端两个顶点进行线性插值有 厂g ，1 ) - s ( o ，1 ) + x k o ，1 ) 一s ( o ，1 ) 】 ( 4 2 5 ) 最后，我们做垂直方向的线性插值，以确定 厂g ，y ) = 厂g ，o ) + j ，纱k 1 ) 一g ，o ) 1 ( 4 2 6 ) 1 8 中北大学学位论文 将以上三个等式展开并合并同类项得： 厂g ，y ) = i t 0 ，o ) 一s ( o ，o ) b + l r ( o ，0 - r ( o ，o ) k + 陟( 1 ，1 ) + 厂( o ，o ) 一厂( o 1 ) 一厂( 1 ，o ) b y + f ( o ，0 ) ( 4 2 7 ) 很明显，该式是双线性的。 4 3 新的插值方法 在本节中考虑了三维重建图像中的每个像素在不同扫描角度时在探测器上投影的 近似图形，( 若像素的大小为一个单位，则此近似图形的面积大于等于1 ) ，也就是考虑 重建像素的大小和探测器的探元的大小，根据穿过各个像素的射线落在探测器上的面积 大小，将在不同探测器上的投影值的加权和作为该像素在此扫描角度下的反投影值。 4 3 1 像素投影图形的取法 通常情状下，我们认为物体的像素在探测器上的投影是个点，但在实际中，在某 一个扫描角度下，经过每个体素的射线在探测器上形成一个几何图形。在重建过程中， 这个几何面的形状不规则而且其面积很不易求出，因此，我们可近似考虑每个像素的一 部分点在探测器上所形成的图形。具体步骤如下：如图2 1 所示，在锥束圆轨迹扫描结 构中，定义锥束的投影角为夕，扇角为，锥角为k 。设, m 3 c d e f g h 是要重建图像的 某一个体素( f ， 七) ，如图4 1 所示，a ( i - 互1 ，+ j 1 ，后+ 三) ，b ：( f 一圭，j 一互1 ，七+ 三1 ) c ： 19 f 、，j1 2 一七 12) + l 一2 一 ， 一r 囊 一 l 一2 l 一 ： l 一2 e + ，、 、j 日 小到 h + ” m “列 1 2 一 十 l z 1 2 1 2 一 + ，、 一2 + ： 、 d ：g 1 2 1 2 + 一 后 七 一2 1 2 一 一 。j ，一2 l 一2 + 一 rp一 -；一 中北大学学位论文 y 图4 1 某一像素各项点在坐标系中的位置 当射线源介于x 正半轴到y 正半轴之间这去万弧度( 即0 。夕 9 0 。) 时，考虑a c g e z 1 这个对角面在探测器上的投影图形。当射线源介于y 正半轴到x 负半轴之间这去万弧度 二 ( 即9 0 。 1 8 0 。) 时，考虑b d h f 这个对角面在探测器上的投影图形。当射线源介 1 于x 负半轴到y 负半轴之间这去万弧度( 即1 8 0 。 2 7 0 。) 时，考虑a c g e 这个对角 二 1 面在探测器上的投影图形。当射线源介于y 负半轴到x 正半轴之间这去刀弧度( 即 二 2 7 0 。 o 物体所在范围： u = 备r 3 ：# + ，2 ，o ， 足 其中p ：螺距； r ：探测源的旋转半径； 2 6 ( 5 1 1 ) ( 5 1 2 ) 中北大学学位论文 0 ，v ) ：探测器阵列的坐标； 在扫描点y 处、方向为夕的投影值是： d ，b ，国= 譬f b + t 舢，芦s 2 风一= 酬v l s 膳u 夕眯c 防一j 5 2k a t s e v i c h 算法描述 考虑圆柱体的感兴趣区域定义为 r o i = ( x ，y ，z ) lx 2 + y 2 r 2 , z m i n z z m a x ) 图5 2 对感兴趣区域的局部重建示意图 k a t s e v i c h 算法公式： s ( d = 一专i ，网1 矿g ，风，x 渺 g ，力= r 厅毒( g ) ，c o s 。夕+ s i n g - e g ，z ) b 面1 砂 “g ，力= ”g ，s ：g ，夕) ) p g ，夕) = 夕“g ，穗s 2 2 7 ( 5 1 3 ) ( 5 1 4 ) ( 5 2 1 ) ( 5 2 2 ) ( 5 2 3 ) ( 5 2 4 ) ( 5 2 5 ) 中北大学学位论文 5 3k a t s e v i c h 算法的执行步骤 具体的步骤可分为三步。分别是： ( 一) 根据扫描的参数对锥束投影数据求导； ( 二) 一维线性斜滤波； ( 三) 反投影重建三维物体。 5 3 1 投影数据求导 令4 g 鸬v ) = 忑dg g 鸬y )