（模式识别与智能系统专业论文）ct图像金属伪影校正算法研究.pdf

摘要 摘要 计算机断层扫描成像( c o m p u t e r i z e dt o m o g r a p h y ，c t ) 是一l - j 结合了机械、自 动化、计算机科学和x 射线等多门学科的新成像技术，由于c t 重建的图像具有 空间和密度分辨率高等特点，使它在医疗无损害诊断和工业无损检测领域得到了 广泛的应用。 由于c t 的被扫描物体中的金属物质造成了投影数据的不准确，使通过其重 建的图像中产生了明暗相间的条状伪影，或称为金属伪影。这些伪影严重损害图 像质量，影响了图像的后续处理和分析。因此，金属伪影问题是c t 研究中的重 要问题，也是自x 射线c t 问世以来一直寻求解决的热点问题。 本文首先介绍了c t 成像和重建的理论基础，分析了金属伪影的成因，然后 对国内外目前的多种关于金属伪影的校正方法进行了归纳综合和仿真实验，并对 其中一些主流算法的实现方法进行了优化，分析了各类方法的利与弊，最后提出 了基于自适应衰减和滤波的c t 去金属伪影混合法。 基于自适应衰减和滤波的c t 去金属伪影混合法首先选取合适的阈值分割出 金属区域，对金属投影数据进行自适应衰减和滤波，然后通过滤波反投影( f b p ) 重建图像。再利用原始投影数据对金属区域进行最大期望值( e m ) 局部重建， 最后根据迭代重建的结果，对衰减和滤波后重建的图像中的金属区域进行补偿。 根据m a t l a b 数值仿真实验结果表明，该算法能有效去除金属伪影，完整地呈现 金属物体的结构，并能保证金属物体周围的成像不失真，降低了阈值选取精确度 对算法重建结果的影响。而且本算法的计算复杂度较小，即使被扫描物体对包含 多个金属物体，图像的重建效果也很好。此外，本文最后还提出了进一步研究和 完善该算法方向，通过这些研究将有可能解决c t 新技术和新应用中的金属伪影 问题。 关键词：金属伪影投影数据滤波反投影混合法自适应衰减和滤波补偿 a b s t r a c t a b s t r a c t c o m p u t e dt o m o g r a p h y ( c a 3i san e wi m a g i n gt e c h n o l o g yw h i c hc o m b i n e s m e c h a n i c a l ，a u t o m a t i o n ，c o m p u t e rs c i e n c ea n dx - r a ya n do t h e rs u b j e c t s i th a sa n e x t e n s i v ea p p l i c a t i o no nn o nh a r mm e d i c a ld i a g n o s t i c sa n dn o nd e s t r u c t i v et e s t i n g b e c a u s eo ft h eh i g hs p a t i a la n dd e n s i t yr e s o l u t i o no ft h ec t i m a g e s a st h em e t a l si nc tt h es c a n n e do b j e c t sc a u s et h ei n a c c u r a c yo ft h ep r o j e c td a t a s ot h a tt h e r ea r el i g h ta n dw h i t e s t r i p ea r t i f a c t si nt h ec t r e c o n s t r u c t i o ni m a g e s ，w h i c h a r ek n o w na st h em e t a la r t i f a c t s t h e r ea r t i f a c t ss e r i o u s l yd a m a g ei m a g eq u a l i t ya n d i m p a c to ft h ef o l l o w u pi m a g ep r o c e s s i n ga n da n a l y s i s t h e r e f o r e ，t h es t u d yo fm e t a l a r t i f a c t si nc ti m a g ei sa ni m p o r t a n ti s s u e ，b u ta l s oh a sb e e nah o ti s s u et of i n da s o l u t i o ns i n c et h ea d v e n t o fx r a yc t t h i sa r t i c l ef i r s ti n t r o d u c e st h et h e o r e t i c a lb a s i so fc ti m a g i n ga n di m a g e r e c o n s t r u c t i o n ，a n a l y s i so ft h ec a u s e so fm e t a la r t i f a c t s ，a n dt h e ns y n t h e s i z e sa n d s i m u l a t e st h em a i nh o m ea n da b r o a dm e t a la r t i f a c tc o r r e c t i o nm e t h o d a l s oo p t i m i z e s s o m eo ft h em a i n s t r e a mo ft h ea l g o r i t h m sa n da n a l y z e st h ep r o sa n dc o n so ft h e s e a l g o r i t h m s f i n a l l yan e wa d a p t i v e l yh y b r i df i l t e r i n ga l g o r i t h mi sp r o p o s e dt or e d u c e t h ea r t i f a c t sc a u s e db ym e t a li nc t i m a g e i nt h en e wa d a p t i v es c a l i n gh y b r i da l g o r i t h m ，f i r s t l ya na d a p t i v es c a l i n ga n d f i l t e r i n gm e t h o di su s e dt op r e p r o c e s st h ep r o j e c t i o nd a t ao fm e t a la n dt h ea r t i f a c t s r e d u c e di m a g ei sr e c o n s t r u c t e db yf i l t e r e db a c kp r o j e c t i o n ( f b p ) m e t h o d t h e nt h e e x p e c t a t i o nm a x i m i z a t i o na l g o r i t h m ( e m ) i sp e r f o r m e do nt h ei t e r a t i v eo r i g i n a lm e t a l p r o j e c td a t a f i n a l l y , ac o m p e n s a t i n gp r o c e d u r ei sa p p l i e dt ot h er e c o n s t r u c t e dm e t a l r e g i o n t h es i m u l a t i o nr e s u l t sh a v ed e m o n s t r a t e dt h a tt h ep r o p o s e da l g o r i t h mc a n r e m o v et h em e t a la r t i f a c t sa n dk e e pt h es t r u c t u r ei n f o r m a t i o no fm e t a lo b j e c t e f f e c t i v e l y i te n s u r e st h a tt h et i s s u e sa r o u n dt h em e t a lw i l ln o tb ed i s t o r t e da n d r e d u c e st h et h r e s h o l do fa c c u r a c yo ft h ei m p a c t t h em e t h o di sa l s oc o m p u t a t i o n a l e f f i c i e n ta n de f f e c t i v ef o rt h ec t i m a g e sw h i c hc o n t a i n ss e v e r a lm e t a lo b j e c t s i n a d d i t i o n ，a te n do ft h ea r t i c l et h ed i r e c t i o n so ft h ef u r t h e rs t u d ya n di m p r o v e m e n to f t h ea l g o r i t h ma r ep r o p o s e d ，t h r o u g ht h e s es t u d i e sw i l lb ep o s s i b l et or e s o l v et h em e t a l a r t i f a c t sp r o b l e m si nt h ec tn e w t e c h n o l o g ya n dn e wa p p l i c a t i o n s k e yw o r d ：m e t a la r t i f a c t s ，p r o j e c td a t a ，f i l t e r e db a c kp r o j e c t i o n ，h y b r i da l g o r i t h m ， a d a p t i v es c a l i n ga n df i l t e r i n g , c o m p e n s a t i n g 中国科学技术大学学位论文原创性声明 本人声明所呈交的学位论文，是本人在导师指导下进行研究工作所取得的成果。除已特 别加以标注和致谢的地方外，论文中不包含任何他人已经发表或撰写过的研究成果。与我一 同工作的同志对本研究所做的贡献均已在论文中作了明确的说明。 作者签名：签字日期： 中国科学技术大学学位论文授权使用声明 作为申请学位的条件之一，学位论文著作权拥有者授权中国科学技术大学拥有学位论 文的部分使用权，即：学校有权按有关规定向国家有关部门或机构送交论文的复印件和电子 版，允许论文被查阅和借阅，可以将学位论文编入有关数据库进行检索，可以采用影印、缩 印或扫描等复制手段保存、汇编学位论文。本人提交的电子文档的内容和纸质论文的内容相 一致。 保密的学位论文在解密后也遵守此规定。 口公开 口保密( 年) 作者签名：盔逸 作者签名：隘【拯 签字隰缉：三：竺 导师签名：j 蛰雄一 签字日期：二翌鹭上! 垒 第1 章绪论 1 1c t 概述 第1 章绪论 计算机层析成像技术( c o m p u t e dt o m o g r a p h y ，c t ) 是2 0 世纪人类在放射诊断 领域取得的重大突破。c t 结合了机械、自动化、计算机科学和x 射线等多门学 科，为医疗无损害诊断和工业无损检测提供了一种新技术。其重建的断层图像具 有无影像重叠、密度分辨率和空间分辨率高、可直接进行数字化处理等优点。 c t 技术还广泛应用于生物学、材料学等学科领域，而且c t 的理论与方法也可 应用于物探、地球物理、射电天文学等学科领域。c t 的原理基于不同的物质对 射线具有不同的衰减性质。c t 检测时，用放射线从各方向照射被测物体，并测 量穿过物体的射线强度，再通过一定的重建算法计算出物体内部各点处的物质关 于射线的线性衰减系数，从而得到被测物体的断层图像。c t 系统由硬件和软件 两部分组成。硬件部分主要包括系统的射线源、检测器、电子学部分以及机械转 动装置。软件部分主要包括数据预处理、图像重建以及图像显示和应用分析部分。 l 1 1c t 的发展历史 c t 理论基础是由奥地利科学家j r a d o n 提出的r a d o n 变换，利用r a d o n 变换可 以通过二维或三维物体各个方向的投影，采用数学方法重建物体图像。这一理论 最早应用在无线电天文学的图像重建中。1 9 3 8 年，德国的g a b r i a lf r a n k 第一次在 x 射线诊断中应用直接反投影重建法，但当时未获得较好的图像质量。1 9 6 1 年， 无线电天文学家w h o l d e n d o f f 用“旋转一迁移”法实现了最早的图像重建。在医 学中应用r a d o n 变换是由a m c k c o r m a 在1 9 6 4 年提出。1 9 6 3 年，a m c o r m a c k 进 一步发展了从x 射线投影重建图像的解析数学方法。1 9 7 2 年在g n h o u n s f i e l d 的 直接贡献下，诞生了世界上第一台可以用于临床诊断的e m i 扫描机，基于 g n h o u n s f i e l d 与a m c o r m a c k 在c t 研制中作出的开创性的工作，他们荣获了 1 9 7 9 年的诺贝尔医学奖和生理学奖。1 9 7 4 年l e d l e y 研制成功了全人体扫描c t ，并 安装在美国乔治镇大学医疗中心。此后，在医学方面，西方发达国家先后研制出 具有高分辨率的螺旋c t 、可超高速成像的电子束c t 等设备。工业方面，7 0 年代 末起美国利用研制的透射式工业c t 设备对军工产品的关键部件做无损检测，美 国科学测量系统公司的工业咖还应用于c a d 洲方面，进行加工元件的仿型 制造，开展逆向工程学的研究。 第1 章绪论 自上世纪7 0 年代c t 问世以来，各国的c t 生产商不断研制新技术，新产品如 雨后春笋，x 射线c t 已从第一代发展到第四代，扫描时间也重最初的十几分钟大 大缩短n 5 秒左右，重建算法也日趋成熟，图像质量也得到了很大的改善。国内 2 0 世纪9 0 年代初自行研制成功了第一台透射式y 射线工业c t 样机，清华大学等单 位研制成功x 射线微焦点i c t 机，东大阿尔派先后研制成功普通医用c t 机和螺旋 c t 机。由于医用和工业c t 的巨大市场需求，使得国产c t 的研发受到高度重视， 进入了一个良好的发展时期。 1 1 2c t 成像原理简介 x 射线c t 通过控制x 光球管围绕被测物体扫描，x 射线沿各方向的直线穿过 物体的横切面，i 贝9 得x 射线的衰减数据，在计算机上利用图像重建原理计算得到 物体横切面上的各点对x 射线的吸收系数，并在计算机或是胶片上显示出来。 c t 扫描射线如图1 1 ( a ) 所示，光源强度为i o 的x 射线穿透一部分均匀的物体， 得到衰减后的光强为i ，若射线穿过物体的长度为y ，则根据b e e r 定理有： tt f u d l = l n ( 专q - ) 或u y = l n - 孥- ( 1 1 ) t 11 其中，u 为物体对x 射线的吸收系数。如果物体内部分段均匀且各段的吸收系数 分别为u l ，u 2 ，u 3 ，u n ，相应的长度分别为y 1 ，y 2 ，y 3 ，y n ，则公式( 1 1 ) 改 为：u l y l + u 2 y 2 + u 3 y 3 + + u n y n = l n ( i o 1 ) 。更一般地，如果物体在( x ，”平面内分 布不均匀，吸收系数u = u ( x ，y ) ，那么沿着某路经l 的衰减为： t f u r l l = l n ( 、 o ) ( 1 2 ) t l 我们将式( 1 2 ) 中的线积分称为c t 投影定理，若不指定路线，则称为投影，它是 某一角度下射线引起的射线投影的集合。 c t 设备克服了线积分的弊病，得到了反映人体组织或是物体结构分布的图 像。总的来讲，c t 的基本原理是采用准直后的x 射线束对人体的某一层面从不 同的角度进行照射，在射线穿过的另一端利用探测器接收多组原始数据，经过计 算机重建后得到用于显示的二维数据矩阵，通过显示设备显示图像。 图1 1 ( b ) 简单说明了c t 的扫描方式。图中扫描器发出一组x 射线，射线穿过 被测区域，在对面通过探测器接收衰减后的x 射线得的射线强度，即投影值。然 后再旋转一个小角度进行扫描，如果以此方式作m 次扫描将得到m 组投影数据。 2 第1 章绪论 ( a ) 射线透过物体 器i 图1 1c t 扫描示意图 在c t 成像技术中，有两个很重要的概念，就是c t 值与c t 图像和窗宽与窗位。 c t 图像是通过数学方法对c t 原始数据进行重建，得到图像矩阵。计算机把 重建图像矩阵中的各个像素转变为不同明暗的相应光点，通过显示设备显示出 来。 c t 值是c t 采用的标准是根据各种组织和物质对x 射线的吸收系数来决定的。 h o u n s f i e l d 将线性衰减系数划分为2 0 0 0 个单位，医学上称为c t 值。由于在物理 过程中，物质的密度是由物质对x 射线的衰减系数来体现的，在研究c t 图像时， 能提供诊断信息的是组织之间的密度差异，而不是绝对密度，所以一般定义c t 值如下： c t 值1 0 0 0 x - h 2 0 ( 1 3 ) 纥2 d 。 即特定物质的c t 值等于该物质的衰减系数与水的吸收系数之差再与水的吸收系 数相比，然后乘以1 0 0 0 。物质的c t 值越高，表明其密度越大，所以金属物体的 c t 数远远大于骨骼和其它组织的c t 数。另外，c t 值的大小还与x 射线的能量有 关，x 射线的能量越低则相同物质的c t 值就越大。 目前，大多数c t 的c - i 数变化范围是1 0 0 0 到+ 3 0 0 0 ，而计算机显示器一般只 能显示2 5 6 个灰度级，人眼只能分辨3 0 个灰度级，因此，为了能用有限的灰度级 较好地表示c t 图像，人们引入了窗口技术，即通过正确地选择窗宽和窗位来显 示符合人眼视觉要求的图像，尤其是在临床医学上，能更好的显示医生感兴趣的 区域。 如图1 2 所示，窗宽是c t 数的显示范围，窗位是显示灰度级的中间位置所对 应的c t 数。例如，如果窗宽时1 0 0 ，窗位是1 0 0 ，那么c t 数的范围是5 0 至1 j 1 5 0 ，其 3 第1 章绪论 中c t 数1 0 0 对应的显示器灰度级是1 2 8 ，c t 数5 0 对应的灰度级为0 ，而c t 数1 5 0 对 应的灰度级为2 5 5 。小于5 0 的c t 数的灰度级全是0 ，大于1 5 0 的c t 数灰度级全是 2 5 5 。如果增大窗宽可以使一个灰度级对应多个c t 数，这样会使c t 数接近的组织 对比度降低，反之，如果减小窗宽会使多个灰度级对应一个c t 数，导致对比度 增加，使被扫描物体结构更清晰，但也会丢失部分信息，所以在运用时，要反复 调节窗宽和窗位，来获得最佳视觉效果。 o c t 数 图1 2 窗宽与窗位 把如图1 2 这种将显示的灰度范围和窗宽内的c t 数按线性关系一一对应的窗 口称为线形窗，线性窗实现简单，可手动动态改变窗宽和窗位，是商用c t 控制 软件中常见的功能。 1 1 3c t 设备的主要技术指标 c t 的主要技术指标包括：c t 的工作时间、c t 系统部件和c t 图像质量指 标，这些指标与c t 性能息息相关。 ( 1 ) c t t 作时间 c t 工作时间包括：扫描时间、图像重建时间、周期时间。扫描时间短，一 方面可以缩短病人暴露在x 射线中的时间；另一方面可以减少由于病人运动而 引起的运动伪影，确保图像质量。c t 扫描时间在医学中采用的指标，是指在能 够获得设备所提供的图像质量的前提下，所需要的最短扫描时间；图像重建时间 是指根据原始数据进行计算机处理，从而获得c t 图像所需要的时间；重建时间 短有利于及时观察扫描图像，迅速调整扫描策略。一般地，重建时间与重建矩阵 4 唑 窗一窗 第1 章绪论 的大小密切相关，显然重建矩阵为5 1 2 x 5 1 2 时所需要的时间比重建矩阵为 2 5 6 x 2 5 6 时要长。对病人的某一断层从扫描开始，经重建、显示到拍片完成整个 过程所需要的时间称为周期时间。一般，周期时间= 扫描时间+ 重建时间+ 拍片 时间。 ( 2 ) c t 系统部件 c t 设备的系统部件包括：扫描机架和病床、x 射线系统、数据采集系统、 计算机及显示系统、照相系统、数据存储系统等。 o ) c r 图像质量指标 空间分辨率：指c t 对物体空间大小及几何尺寸的分辨能力。通常有两种表 示方式：一种是用单位厘米内的线数来表示；另一种是用可辨别最小物体的直径 用毫米来表示。 密度分辨率：表示c t 设备对密度差别的分辨能力，通常以百分数表示。比 如0 4 5 表示两种物质的密度差大于0 4 5 时，c t 可以将它们区分出来。 伪影：在被扫描取样的物体中并不存在的，但在图像中却出现的各种不同类 型的影像。一般地，按伪影产生的方式分为c t 设备伪影和病人运动伪影两种。 1 1 4c t 图像的特点 介绍了c t 的成像原理及主要技术指标后，在此基础上来分析c t 断层扫描 形成的图像的特点。c t 图像是通过计算机计算出来的x 射线衰减值的二维分布 图，是由一定数目的像素按矩阵排列而成的二维断层图像，这些像素反映了相应 单位容积的x 射线吸收系数。目前常用的c t 装置的像素大小为2 5 6 x 2 5 6 、 3 2 0 x 3 2 0 和5 1 2 x 5 1 2 。c t 图像一般是8 位或1 6 位的灰度图像。与其它很多成像 系统一样，c t 成像也受到噪声的干扰。深入了解c t 成像过程中受到的各种干 扰以及各种干扰在c t 图像上的表现，有助于在c t 图像的处理过程中消除伪影的 影响，提高图像处理的质量。c t 伪影和噪声是评价c t 成像质量的重要指标，下 面就对伪影和噪声对c t 图像的影响进行分析，以最大程度地减少两者对c t 图像 处理的影响。 在c t 图像上非真实的阴影或干扰称为伪影。它降低图像的质量，易造成误判 或不可诊断。它可分为由被扫描物人或物引起的伪影与c t 设备本身所造成的伪 影两大类。由被扫描物人或物引起的伪影中，最常见的是在扫描取样时，病人可 能体内有金属物体引起的放射状金属伪影，如图1 3 ( a ) 所示。病人在扫描时自主 或不自主的运动引起的运动伪影。 均匀物体的c t 影像中各像素的c t 值参差不齐，它使图像呈颗粒性，直接影 响其密度分辨率，尤其以低密度的可见度为甚。把这种现象用统计学上的标准偏 5 第1 章绪论 差方式表示出来，即为c t 的噪声，可分为随机噪声和统计噪声。它们的产生机 理和对图像质量的影响各不相同，严重的随机噪声往往形成了伪影，如图1 3 ( b ) 所示。在c r y 2 像中，x 射线光予数是一个随机起伏的统计量，即使电压非常稳定 其它一切条件都不变，表示x 射线强度的光子数也是随着时间统计而起伏。此外 检测器和放大器中也有随机起伏的噪声( 张伟，2 0 0 6 ) 。所有的这些构成了随机测 量误差，反映在投影值上，它使投影值曲线上产生单个向上或向下的小尖峰。由 于在c t 图像的重建中，每个测量值都被反投影到整个图像上，所以这个单尖峰 就变成重建图像上的一条条纹。当有许多随机噪声的小尖蜂时，重建国像上就出 现许多杂乱的细小条纹和由许多条纹交叉形成的颗粒斑点，成为图像随机噪声。 统计噪声又称像素噪声，指服从泊松分布的噪声，主要由x 射线的强度和数量的 起伏，以及x 射线穿透人体时要产生散射和吸收所形成的f y a v u zme ta l ，1 9 9 9 ) 。 糜 “) 含金属伪影的c t 豳像佃) 含严重噪声的c t 图像 图1 3 古伪影和噪声c t 图像 1 2 研究意义和国内外研究现状 随着c 1 技术的发展，c t 已经成为医疗诊断和工业检测等行业的一种非常重 要的技术。但是在c r 普及的过程中出现很多有待解决的新问题，c r 图像中的金 属伪影就是其中之一。理论上伪影可以定义为图像中被重建数值和物体真实的衰 减系数之间的差异( k a l e n d e r da l ，1 9 8 7 ) 。尽管定义足够宽泛，但它几乎包括所有 的非理想图像，所以我们需要对特定的伪影表现进行说明。条状伪影经常表现为 横穿图像的明显直线f - 傲不是平行的) 它们可以是亮的也可以是暗的。由于重 建算法中滤波器的特性，明暗条纹是成对出现的( g l o v e rg h a n dp e l tn j ，1 9 8 1 ) 。 由于这种伪影般由金属物体引起，所以又常称为金属伪影。 无论在医疗诊断还是工业检测都经常遇到被检测物体中含金属物质，或是其 它吸收系数很高的物质，导致投影数据出现了跃变，重建后的图像就会包含明暗 6 第1 章绪论 相问的金属伪影，这些伪影一般数量众多、幅值又高，严重降低图像质量 ( b a i s s a l o ve ta l ，2 0 0 0 ) ，无论是在医疗或是工业检测都影响了图像的后续处理和分 析，如何消除这些条状伪影以提高图像质量，已经成为当前c t 应用研究领域中 的热点。 用于消除c t 图像中金属伪影的方法来自学术机构或是工业界，在不同的杂 志或是会议论文集上可以找到许多研究论文和文章。近几年一些比较有效的c t 去金属伪影算法的得到了快速发展，但都存在一定的缺陷。 国外宾夕法尼亚大学的l e w i t tr m 和b a t e sr h ( 1 9 8 3 ) 提出了用一种特殊的函数 进行插值的方法，后来又提出了c h e b y s h e v 多项式的差值法；c r a w - f o r d ( 1 9 8 6 ) 在 线性插值的基础上添加了一些辅助处理；华盛顿大学w a n gg 等人( 1 9 9 6 ) 提出 的迭代代数重建法。芝加哥大学的x i ad a n ( 2 0 0 5 ) j 丕提出了一种e m ( e x p e c t a t i o n m a x i m i z a t i o n ) 局部迭代的混合算法。还有一些人提出了马尔可夫随机场方法， 算法基于传播最大似然的思想，假设测量用的探测器读数服从泊松分布，通过最 优似然来重建图像。 国内这方面的研究还比较少，主要是清华大学的谷建伟提出的快速线性插值 法( 谷建伟等，2 0 0 5 ) 和基于伪影权重的非线性插值法( 谷建伟等，2 0 0 6 b ) ，以及林 宙辰等( 2 0 0 1 ) 提出的采用多项式插值。 现在主流去金属伪影算法分为三类：插值重建算法、迭代重建算法和混合重 建算法。插值算法简单快速，但是只适用于金属结构简单的情况，而且会使金属 周围的区域失真。迭代法也是一种常用的c t 图像重建算法，能有效去除金属伪 影和抑制噪声，而且能很好的呈现金属物体的结构，但其运算量非常大，速度很 慢，难以实用化。混合法是近几年提出的一种新算法，e m 混合法( x n de ta l ，2 0 0 5 ) 是其中代表，该算法速度较快，能有效反映金属物体结构，不足之处是金属周围 的区域失真，使其在医用咖用中存在很大的局限性。 为了克服以上这些算法的不足，本文对金属伪影成因和去金属伪影算法进行 了详细的研究，并设计了一种新的去金属伪影算法。 1 3 本文的主要工作及内容安排 论文通过研究c t 的成像和重建原理、金属伪影成因和现在主流的去金属伪 影重建算法，并在这些研究的基础设计一种新的金属伪影消除算法，各章安排如 下： 第1 章介绍了c t 的发展历史、成像原理、主要技术指标和课题的研究背 景、意义及国内外研究现状，针对目前现在主流去金属伪影算法的不足，确定了 本文的研究方向。 7 第1 章绪论 第2 章介绍了c t 重建的理论和数学基础，并说明了如何通过变换法和数 级展开法进行c t 图像重建。变换法中主要介绍了滤波反投影重建算法，数级展 开法( 迭代法) ，重点介绍了代数迭代重建算法( a r t ) 、联合代数迭代重建算法 ( s a r r ) 、期望最大值迭代算法( e m ) 和有序子集期望最大值迭代法( o s e m ) 。 第3 章重点研究了c t 中金属伪影的成因和现在主流的去金属伪影算法： 插值重建法、迭代重建法和e m 混合重建法。通过m a t l a b 构建数值仿真平台， 在仿真平台的基础上对这些算法进行了去金属伪影实验并对它们的实现方法进 行了优化，最后通过实验结果总结了这些算法的优点和不足，并提出了去金属伪 影算法的改进思路。 第4 章针对第三章中提到的主流去金属伪影算法的不足，设计出一种去c t 图像中金属伪影自适应衰减滤波混合法。对自适应衰减滤波混合法每一个步骤进 行了详细的介绍，并通过数值模拟对比实验证明，该方法可以有效的去除金属伪 影，并很好的保留了金属及其周围组织的信息，其效果优于插值法和e m 混合法。 特别是在多金属物体的情况下，算法的效果更是优于其它算法。 第5 章对本文研究的主要内容、创新点进行了总结，并对未来的研究工作的 方向进行了展望。 8 第2 章c t 成像基本原理与图像重建算法 第2 章c t 成像基本原理与图像重建算法 2 1 经典的断层成像技术 为了克服传统的x 射线摄影技术中的图象重叠问题，在c t 发明以前已有经 典的断层成像技术出现。这种成像技术沿用至今已有6 0 余年，并发展成许多模 式。下面介绍其中的一种纵轴直线断层成像技术。 纵轴直线断层成像技术是这样一种成像技术，它对被选中的感兴趣断层清晰 的成像，而是这一层的上下层模糊可见或是基本看不见。 传统的断层成像技术如图2 1 ，它表示经典的纵轴直线断层成像技术的基本原 理( 黄庆华，2 0 0 2 ) 。被扫描物体或病人位于点状x 射线源和检测器之间。x 射线 源的运动被固定在某一平面上，这一平面被称为源平面。检测器所在平面也是被 固定在某一平面，一般这个平面就是感光底板p 。图2 1 中，源平面和感光底板 相互平行，x 射线源x 与检测器a b 作反方向运动，它们的速度应满足x 射线 源x 到检测器上的某点a 所作的直线始终通过断层中的点a 。同理x 射线源x 到检测器上的某点b 所作的直线始终通过断层中的点b ，则a b 就是感兴趣断层。 图中x 1 和x 2 分别是射线源x 在t 1 、t 2 时刻的位置，a 1 、a 2 和b 1 、b 2 分别 是检测器平面上两点a 、b 在t 1 、t 2 时刻的位置。由它们运动速度的约束可知， 直线x 1 a 1 和x 2 a 2 必交于点a ，直线x 1 8 1 和x 2 8 2 必交于点b 。这样曝光期 内，可以使断层a b 平面的图像在感光板p 上不断加强而清晰成像。 x 1x 2 射线源x 移动方向 a 2a lb 2b 1 图2 1 传统的断层成像技术 9 第2 章c t 成像基本原理与图像重建算法 2 2c t 图像重建的问题描述 图2 2 表示图像重建中的投影坐标系，设一个具有某种物理性质( 吸收系数) 的二维分部函数f ( x ，y ) 定义在平面域0 内，p ( x ，口) 表示f ( x ，y ) 沿着由xw 和0 定 义的一条直线的线积分，穿过厂 ，y ) 的所有线积分的集合为： p ( x ，日) 2 j = 。f ( x ，y ) d y 。j = 。f ( x c o s o ys i n 0 ，s i n 0 + y c o s 0 ) d y = r c o s ( o 一驴) x e ( - o o ，+ ) p 【0 ，万)( 2 1 ) 对于一个固定的角度0 ，p o ，口) 是f ( x ，y ) 沿口方向在z 轴上的一个一维平行投 影。如果目标函数厂 ，y ) 未知，而它的所有投影( 基线积分) 却能通过某种方式获 得，则f ( x ，y ) 便可由p o ，0 ) 唯一的确定，这个问题是r a d o n 在1 9 1 7 提出的， 因此，( 2 1 ) 式中联系的目标函数和它的投影函数之间的积分变换关系被称为 r a d o n 变换。从数学上看，图像重建问题就是求解r a d o n 变换的反变换，即给定 投影数据p ( x ，0 ) ，求出f ( x ，y ) 。而这个问题在工程上称为从投影重建图像，或 c t 图像重建。 p ( x 、，o ) y 彦雩 ， ，x 文 、) 7 ，7 7 、o 77 ，j 。i 、，7 一 ，k 、女 、 2 3c t 重建算法 图2 2 投影的坐标系 从投影重建图像的算法很多，如：反投影重建算法、滤波反投影重建、直接 傅立叶变换重建算法和迭代重建算法。但是反投影重建算法和直接傅立叶变换重 建算法都存在明显缺陷，很少使用。反投影重建算法图像中存在密度为零的点， 所以重建后不一定为零，从而引入星状伪影。由于傅立叶空间中产生的采样模式 1 0 第2 章c t 成像基本原理与图像重建算法 不是笛卡尔坐标，为了执行二维傅立叶反变换这些采样不得不重新插值或重新栅 格化到一个笛卡尔坐标系中，在傅立叶空间中一个采样点的误差会影响到整个图 像，所以直接傅立叶变换重建算法重建图像效果不好。因此将本节重点介绍最常 见的滤波反投影重建和迭代重建算法。 重建算法在具体实现时，会因为c t 扫描方式的不同略有不同，但整体思想 还是不变的。现行c t 主要扫描方式有平行束扫描和扇形束扫描。迭代重建算法 不会过分依赖于扫描方式，而滤波反投影重建算法相对较依赖于扫描方式。 2 3 1 滤波反投影算法 2 3 1 1r a d o n 变换和中心切片定理 r a d o n 变换及其逆变换和中心切片定理是滤波反投影重建( f i l t e r e db a c k p r o j e c t i o n ，f b p ) 的理论基础，而且r a d o n 变换及其逆变换还是迭代法的数学基 础。本节将详细介绍r a d o n 变换和中心切片定理。 由1 1 2 节可知c t 投影数据是衰减系数的线形积分，图像重建问题其实就是 是否可以通过投影数据估计物体内部的衰减系数( h e r m a ng t , 1 9 9 5 ) ，j r a d o n 于 1 9 7 1 年给出了肯定的答案，就是r a d o n 变换及其逆变换。 r a d o n 变换可在任意维空间定义，且定义存在多种形式。在c t 扫描中只考 虑二维情况。那么二维函数f ( x ，y ) 的r a d o n 变换定义为： r f ( x ，y ) ) = rf ( x ，y ) 6 l - ( x c o s o + y s i n o ) d x d y - 一一 一r o u ) = r ( 1 ，口)( 2 2 ) 式中z 代表坐标原点到直线l 的距离，表达式z - x c o s o + y s i n 0 代表直线l ，沿 着一系列平行线( 投影线) 的积分就组成了投影r o ( f ) ，所有投影组成的集合 也u ) ，目【0 幼) ) 就是r a d o n 变换或者称为投影。其中积分域为整个x y 平面，0 为距离z 与x 轴的夹角，6 为脉冲函数。 6 c z ，= f 三二二吕 ( 2 3 ， 当夹角0 为固定值时，那么函数f ( x ，y ) 在夹角0 下的投影值是距离z 的一维 函数，如图2 3 所示。 第2 章c t 成像基本原理与图像重建算法 j y 一 x 图2 3r a d o n 变换示意图 l r a d o n 变换可以认为是z 一0 空f n - j q 丁的投影，z 一0 空间中的每一点对应一条直 线，而r a d o n 变换就是图像像素点在每一条直线上的积分。因此，图像中灰度 级高的像素点在z 一0 空间中形成亮点( c i 图像中金属物体是亮点) ，而低灰度级 的像素点在z 一0 空间中为暗点。根据r a d o n 变换的定义在z 一0 空间中有： ( 2 4 ) 则r a d o n 变换的另一种表达形式为： r f ( x ，y ) 】- 2 fi ( p c o s o s s i n 0 ，p c o s o + p s l n 口灿 ( 2 5 ) 其示意图如图2 4 所示。 忒 s 雾弋一 心 图2 4r a d o n 变换的另一种表达示意图 1 2 口口 n n s s s s 一 +口口 瞄 螂 d p i 暑 x y rl 第2 章c t 成像基本原理与图像重建算法 为了求解f ( x ，y ) ，r a d o n 在1 9 1 7 年对公式( 2 2 ) 进行求逆运算，这就是著名 的r a d o n 反变换，求解的， ，y ) 为： 舷加昙r d 啦南学讲 ( 2 6 ) 式中x = 伍，y ) ，宇= ( c o s o ，s i n o ) 且歹p ，亭) 一c ，( z 亭+ 亭上) i t 。只有r a d o n 变换理 论还不能设计滤波反投影算法，还需要中心切片定理得辅助。 该定理在非衍射情况下，中心切片定理示意图如图2 5 所示，某图像f ( x ，y ) 在视角为妒时的投影p 西( x r ) 的一维傅立叶变换给出厂o ，y ) 的二维傅立叶变换 f 2 ( w 1 ，w 2 ) - e p ，妒) 的一个切片，切片与w l 轴相交成驴角，且通过坐标原点。即 互【所“) 】= 最p ，妒) k 定 ( 2 7 ) 图2 5 中心切片定理示意图 式中，f 【】表示以为傅立叶变换。下面我们简单的证明一下该定理： 先证当7 = o 时，此时x t - - x ，所以投影 儿以) 5 p o ( x ) 。f f ( x ，y ) d y( 2 8 ) g p 。g ) 】2fp 。( x ) e - i w , x d x2 f ff ( x ，y ) d y e 一叩d x2 f f f ( x ，y ) e 呻d x d y 而f ( x ，y ) 的二维傅立叶变换为 f ( m ，) = e 【厂g ，) ，) 】= rf r c x ，y ) e 一“w d x d y ( 2 9 ) ( 2 1 0 ) 第2 章c t 成像基本原理与图像重建算法 比较公式( 2 8 ) 和( 2 9 ) 可以得到： e 0 。 ) ) = e ( 嵋，) 1 w ，。= e ( m ，0 )( 2 1 1 ) 也就是说y ( x ，y ) 沿y 轴的维傅立叶变换给出了厂o ，y ) 的二维傅立叶变换 v 2 ( w l ，w 2 ) 在w 2 - o 处的一个切片。 若投影的视角为妒，利用旋转坐标一一y ，使积分路径总是沿着y r 。此时x ， 轴与x 轴之间的夹角为，如图2 5 所示，并有关系 量 。 一c s o ；s n # 矽s i n s # | 。 y x c 2 1 2 ， 以y r 作为投影轴的投影为： 功“) 。r ， ，y ) a y ，= rf “，y ，) 方， ( 2 1 3 ) 它的一位傅立叶变换为： 壬式p b r 5f p q r y 却rd x r ； lf r q r ，y r ) d y r d 2 弼x r 缸r 在推导式( 2 1 3 ) 的过程中用到了关系式 红方，= a e a x 砂， a x a e 砂 砂， 砂 。ff f ( z ，y ) 挪o 。妒郴如妒d x d y ( 2 1 4 ) 蚴2 瞄c o s # 妒耋纠蚴；螂c 2 胚， 以及f “，y ，) = ， ，y ) 。由式( 2 1 4 ) 可知，令m 一2 死p c o s # ，一2 z p s i n # ，代入 式( 2 1 4 ) 有 可d “) 】= p ( p ) 一f 2 ( w , ，) = f 2 ( p ，驴) 其中驴= a r c t g ( w 2 j w l ) ( 2 1 6 ) 式( 2 1 6 ) 说明：w l ，w 2 并不是独立的，而是受到直线一a r c t g ( w 2 w 1 ) 的约束。这 样就证明了投影定理。 2 3 1 2 反投影算法的数学描述 我们把“取投影”一“反投影重建”呻“重建后图像这些环节看作一个以原像 为输入重建后图像为输出的成像系统，如图2 6 所示。为了仔细剖析这个系统， 先求该系统的点扩展函数。 1 4 第2 章c t 成像基本原理与图像重建算法 图像 垂 i ( x ，y ) l 取投影 反投影重建 图2 6 反投影重建等效成像系统 设一位于坐标原点x = 0 ，y = 0 的点源d ( x ，y ) 为x y 断面中的唯一像点。并设 扫描方式为平行旋转。即射线先平行移动，从物体的一侧移动到另一侧；让后 旋转一个角度妒，再加上作一个平行移动。之后又旋转一个角度妒。如此继续， 直到累计的转角达到1 8 0 度为止。我们还要设置一个旋转坐标系统x r y ， 见图2 3 ，x y 为固定坐标，x r y ，为旋转坐