CT图像金属伪影校正算法研究

CT图像金属伪影校正算法研究

姓名：陈豫申请学位级别：硕士专业：模式识别与智能系统

指导教师：周荷琴

20090501

摘要摘要计算机断层扫描成像（ＣｏｍｐｕｔｅｒｉｚｅｄＴｏｍｏｇｒａｐｈｙ，ＣＴ）是一ｌ＇－Ｊ结合了机械、自动化、计算机科学和Ｘ射线等多门学科的新成像技术，由于ＣＴ重建的图像具有空间和密度分辨率高等特点，使它在医疗无损害诊断和工业无损检测领域得到了广泛的应用。由于ＣＴ的被扫描物体中的金属物质造成了投影数据的不准确，使通过其重建的图像中产生了明暗相间的条状伪影，或称为金属伪影。这些伪影严重损害图像质量，影响了图像的后续处理和分析。因此，金属伪影问题是ＣＴ研究中的重要问题，也是自Ｘ射线ＣＴ问世以来一直寻求解决的热点问题。本文首先介绍了ＣＴ成像和重建的理论基础，分析了金属伪影的成因，然后对国内外目前的多种关于金属伪影的校正方法进行了归纳综合和仿真实验，并对其中一些主流算法的实现方法进行了优化，分析了各类方法的利与弊，最后提出了基于自适应衰减和滤波的ＣＴ去金属伪影混合法。基于自适应衰减和滤波的ＣＴ去金属伪影混合法首先选取合适的阈值分割出金属区域，对金属投影数据进行自适应衰减和滤波，然后通过滤波反投影（ＦＢＰ）重建图像。再利用原始投影数据对金属区域进行最大期望值（ＥＭ）局部重建，最后根据迭代重建的结果，对衰减和滤波后重建的图像中的金属区域进行补偿。根据Ｍａｔｌａｂ数值仿真实验结果表明，该算法能有效去除金属伪影，完整地呈现金属物体的结构，并能保证金属物体周围的成像不失真，降低了阈值选取精确度对算法重建结果的影响。而且本算法的计算复杂度较小，即使被扫描物体对包含多个金属物体，图像的重建效果也很好。此外，本文最后还提出了进一步研究和完善该算法方向，通过这些研究将有可能解决ＣＴ新技术和新应用中的金属伪影问题。关键词：金属伪影投影数据滤波反投影混合法自适应衰减和滤波补偿ＡｂｓｔｒａｃｔＡｂｓｔｒａｃｔＣｏｍｐｕｔｅｄＴｏｍｏｇｒａｐｈｙ（ｃａ３ｉｓａｎｅｗｉｍａｇｉｎｇｔｅｃｈｎｏｌｏｇｙｗｈｉｃｈｃｏｍｂｉｎｅｓａｎｍｅｃｈａｎｉｃａｌ，ａｕｔｏｍａｔｉｏｎ，ｃｏｍｐｕｔｅｒｓｃｉｅｎｃｅａｎｄＸ－ｒａｙａｎｄｏｔｈｅｒｓｕｂｊｅｃｔｓ．ＩｔｈａｓｅｘｔｅｎｓｉｖｅａｐｐｌｉｃａｔｉｏｎｏｎｎｏｎｈａｒｍｍｅｄｉｃａｌｄｉａｇｎｏｓｔｉｃｓａｎｄｎｏｎｄｅｓｔｒｕｃｔｉｖｅｔｅｓｔｉｎｇｂｅｃａｕｓｅｏｆｔｈｅｈｉｇｈｓｐａｔｉａｌａｎｄｄｅｎｓｉｔｙｒｅｓｏｌｕｔｉｏｎｏｆｔｈｅＣＴｉｍａｇｅｓ．ＡｓｔｈｅｍｅｔａｌｓｉｎＣＴｔｈｅｓｃａｎｎｅｄＳＯｏｂｊｅｃｔｓｃａｕｓｅｔｈｅｉｎａｃｃｕｒａｃｙｏｆｔｈｅｐｒｏｊｅｃｔｄａｔａｔｈａｔｔｈｅｒｅｋｎｏｗｎａｒｅｌｉｇｈｔａｎｄｗｈｉｔｅｓｔｒｉｐｅａｒｔｉｆａｃｔｓｉｎｔｈｅＣＴｒｅｃｏｎｓｔｒｕｃｔｉｏｎｉｍａｇｅｓ，ｗｈｉｃｈａｒｅａｓｔｈｅｍｅｔａｌａｒｔｉｆａｃｔｓ．Ｔｈｅｒｅａｒｔｉｆａｃｔｓｓｅｒｉｏｕｓｌｙｄａｍａｇｅｉｍａｇｅｑｕａｌｉｔｙａｎｄｉｍｐａｃｔｏｆｔｈｅｆｏｌｌｏｗ—ｕｐｉｍａｇｅｐｒｏｃｅｓｓｉｎｇａｎｄａｎａｌｙｓｉｓ．Ｔｈｅｒｅｆｏｒｅ，ｔｈｅｓｔｕｄｙｏｆｍｅｔａｌａｒｔｉｆａｃｔｓｉｎＣＴｉｍａｇｅｉｓａｎｉｍｐｏｒｔａｎｔｉｓｓｕｅ，ｂｕｔａｌｓｏｈａｓｂｅｅｎａｈｏｔｉｓｓｕｅｔｏｆｉｎｄａｓｏｌｕｔｉｏｎｓｉｎｃｅｔｈｅａｄｖｅｎｔｏｆＸ—ｒａｙＣＴ．ＴｈｉｓａｒｔｉｃｌｅｆｉｒｓｔｉｎｔｒｏｄｕｃｅｓｔｈｅｔｈｅｏｒｅｔｉｃａｌｂａｓｉｓｏｆＣＴｉｍａｇｉｎｇａｎｄｉｍａｇｅｒｅｃｏｎｓｔｒｕｃｔｉｏｎ，ａｎａｌｙｓｉｓｏｆｔｈｅｃａｕｓｅｓｏｆｍｅｔａｌａｒｔｉｆａｃｔｓ，ａｎｄｔｈｅｎｓｙｎｔｈｅｓｉｚｅｓａｎｄｓｉｍｕｌａｔｅｓｔｈｅｍａｉｎｈｏｍｅａｎｄａｂｒｏａｄｍｅｔａｌａｒｔｉｆａｃｔｃｏｒｒｅｃｔｉｏｎｍｅｔｈｏｄ．Ａｌｓｏｏｐｔｉｍｉｚｅｓｓｏｍｅｏｆｔｈｅｍａｉｎｓｔｒｅａｍｏｆｔｈｅａｌｇｏｒｉｔｈｍｓａｎｄａｎａｌｙｚｅｓｔｈｅｐｒｏｓａｌｇｏｒｉｔｈｍｓ．ＦｉｎａｌｌｙｔｈｅａｒｔｉｆａｃｔｓａａｎｄｃｏｎｓｏｆｔｈｅｓｅｎｅｗａｄａｐｔｉｖｅｌｙｈｙｂｒｉｄｆｉｌｔｅｒｉｎｇａｌｇｏｒｉｔｈｍｉｓａｎｐｒｏｐｏｓｅｄｔｏｒｅｄｕｃｅｃａｕｓｅｄｂｙｍｅｔａｌｉｎＣＴｉｍａｇｅ．ａｄａｐｔｉｖｅｓｃａｌｉｎｇＩｎｔｈｅｎｅｗａｄａｐｔｉｖｅｓｃａｌｉｎｇｈｙｂｒｉｄａｌｇｏｒｉｔｈｍ，ｆｉｒｓｔｌｙｆｉｌｔｅｒｉｎｇｍｅｔｈｏｄｉｓｕｓｅｄｔｏｐｒｅｐｒｏｃｅｓｓｔｈｅｒｅｄｕｃｅｄｉｍａｇｅｉｓａｎｄｐｒｏｊｅｃｔｉｏｎｄａｔａｏｆｍｅｔａｌａｎｄｔｈｅａｒｔｉｆａｃｔｓｔｈｅｒｅｃｏｎｓｔｒｕｃｔｅｄｂｙｆｉｌｔｅｒｅｄｂａｃｋｐｒｏｊｅｃｔｉｏｎ（ＦＢＰ）ｍｅｔｈｏｄ．Ｔｈｅｎｅｘｐｅｃｔａｔｉｏｎｍａｘｉｍｉｚａｔｉｏｎａｌｇｏｒｉｔｈｍ（ＥＭ）ｉｓｐｅｒｆｏｒｍｅｄｏｎｔｈｅｉｔｅｒａｔｉｖｅｏｒｉｇｉｎａｌｍｅｔａｌｐｒｏｊｅｃｔｄａｔａ．Ｆｉｎａｌｌｙ，ａｃｏｍｐｅｎｓａｔｉｎｇｐｒｏｃｅｄｕｒｅｉｓａｐｐｌｉｅｄｔｏｔｈｅｒｅｃｏｎｓｔｒｕｃｔｅｄｍｅｔａｌｒｅｇｉｏｎ．Ｔｈｅｓｉｍｕｌａｔｉｏｎｒｅｓｕｌｔｓｈａｖｅｄｅｍｏｎｓｔｒａｔｅｄｔｈａｔｔｈｅｐｒｏｐｏｓｅｄａｌｇｏｒｉｔｈｍｃａｎｒｅｍｏｖｅｔｈｅｍｅｔａｌａｒｔｉｆａｃｔｓａｎｄｋｅｅｐｔｈｅｓｔｒｕｃｔｕｒｅｅｆｆｅｃｔｉｖｅｌｙ．Ｉｔｅｎｓｕｒｅｓｉｎｆｏｒｍａｔｉｏｎｏｆｍｅｔａｌｏｂｊｅｃｔｔｈａｔｔｈｅｔｉｓｓｕｅｓａｒｏｕｎｄｔｈｅｍｅｔａｌｗｉｌｌｎｏｔｂｅｄｉｓｔｏｒｔｅｄａｎｄｒｅｄｕｃｅｓｔｈｅｔｈｒｅｓｈｏｌｄｏｆａｃｃｕｒａｃｙｏｆｔｈｅｉｍｐａｃｔ．ＴｈｅｍｅｔｈｏｄｉｓａｌｓｏｃｏｍｐｕｔａｔｉｏｎａｌｅｆｆｉｃｉｅｎｔａｎｄｅｆｆｅｃｔｉｖｅｆｏｒｔｈｅＣＴｉｍａｇｅｓｗｈｉｃｈｃｏｎｔａｉｎｓｓｅｖｅｒａｌｍｅｔａｌｏｂｊｅｃｔｓ．Ｉｎｏｆａｄｄｉｔｉｏｎ，ａｔｅｎｄｏｆｔｈｅａｒｔｉｃｌｅｔｈｅｄｉｒｅｃｔｉｏｎｓｏｆｔｈｅｆｕｒｔｈｅｒｓｔｕｄｙｔｈｅａｌｇｏｒｉｔｈｍａｒｅａｎｄｉｍｐｒｏｖｅｍｅｎｔｐｒｏｐｏｓｅｄ，ｔｈｒｏｕｇｈｔｈｅｓｅｓｔｕｄｉｅｓｗｉｌｌｂｅｐｏｓｓｉｂｌｅｔｏｒｅｓｏｌｖｅｔｈｅｍｅｔａｌａｒｔｉｆａｃｔｓｐｒｏｂｌｅｍｓｉｎｔｈｅＣＴｎｅｗｔｅｃｈｎｏｌｏｇｙａｎｄｎｅｗａｐｐｌｉｃａｔｉｏｎｓ．ＫｅｙＷｏｒｄ：ＭｅｔａｌＡｒｔｉｆａｃｔｓ，ＰｒｏｊｅｃｔＤａｔａ，ＦｉｌｔｅｒｅｄＢａｃｋＰｒｏｊｅｃｔｉｏｎ，ＨｙｂｒｉｄＡｄａｐｔｉｖｅＳｃａｌｉｎｇａｎｄＦｉｌｔｅｒｉｎｇ，ＣｏｍｐｅｎｓａｔｉｎｇＡｌｇｏｒｉｔｈｍ，Ⅱ中国科学技术大学学位论文原创性声明本人声明所呈交的学位论文，是本人在导师指导下进行研究工作所取得的成果。除已特别加以标注和致谢的地方外，论文中不包含任何他人已经发表或撰写过的研究成果。与我一同工作的同志对本研究所做的贡献均已在论文中作了明确的说明。作者签名：签字日期：中国科学技术大学学位论文授权使用声明作为申请学位的条件之一，学位论文著作权拥有者授权中国科学技术大学拥有学位论文的部分使用权，即：学校有权按有关规定向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅，可以将学位论文编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。本人提交的电子文档的内容和纸质论文的内容相一致。保密的学位论文在解密后也遵守此规定。口公开作者签名：口保密（——年）隘【拯作者签名：盔逸导师签名：Ｊ蛰雄一．签字日期：—二翌鹭上』！垒签字隰缉：三：竺．第１章绪论第１章绪论１．１ＣＴ概述计算机层析成像技术（ＣｏｍｐｕｔｅｄＴｏｍｏｇｒａｐｈｙ，ＣＴ）是２０世纪人类在放射诊断领域取得的重大突破。ＣＴ结合了机械、自动化、计算机科学和ｘ射线等多门学科，为医疗无损害诊断和工业无损检测提供了一种新技术。其重建的断层图像具有无影像重叠、密度分辨率和空间分辨率高、可直接进行数字化处理等优点。ＣＴ技术还广泛应用于生物学、材料学等学科领域，而且ＣＴ的理论与方法也可应用于物探、地球物理、射电天文学等学科领域。ＣＴ的原理基于不同的物质对射线具有不同的衰减性质。ＣＴ检测时，用放射线从各方向照射被测物体，并测量穿过物体的射线强度，再通过一定的重建算法计算出物体内部各点处的物质关于射线的线性衰减系数，从而得到被测物体的断层图像。ＣＴ系统由硬件和软件两部分组成。硬件部分主要包括系统的射线源、检测器、电子学部分以及机械转动装置。软件部分主要包括数据预处理、图像重建以及图像显示和应用分析部分。Ｌ１．１ＣＴ的发展历史ＣＴ理论基础是由奥地利科学家Ｊ．Ｒａｄｏｎ提出的Ｒａｄｏｎ变换，利用Ｒａｄｏｎ变换可以通过二维或三维物体各个方向的投影，采用数学方法重建物体图像。这一理论最早应用在无线电天文学的图像重建中。１９３８年，德国的ＧａｂｒｉａｌＦｒａｎｋ第一次在Ｘ射线诊断中应用直接反投影重建法，但当时未获得较好的图像质量。１９６１年，无线电天文学家ｗ．Ｈ．Ｏｌｄｅｎｄｏｆｆ用“旋转一迁移”法实现了最早的图像重建。在医学中应用Ｒａｄｏｎ变换是由Ａ．Ｍ．ｃｋＣｏｒｍａ在１９６４年提出。１９６３年，Ａ．Ｍ．Ｃｏｒｍａｃｋ进一步发展了从Ｘ射线投影重建图像的解析数学方法。１９７２年在Ｇ．Ｎ．Ｈｏｕｎｓｆｉｅｌｄ的直接贡献下，诞生了世界上第一台可以用于临床诊断的ＥＭＩ扫描机，基于Ｇ．Ｎ．Ｈｏｕｎｓｆｉｅｌｄ与Ａ．Ｍ．Ｃｏｒｍａｃｋ在ＣＴ研制中作出的开创性的工作，他们荣获了１９７９年的诺贝尔医学奖和生理学奖。１９７４年Ｌｅｄｌｅｙ研制成功了全人体扫描ＣＴ，并安装在美国乔治镇大学医疗中心。此后，在医学方面，西方发达国家先后研制出具有高分辨率的螺旋ＣＴ、可超高速成像的电子束ＣＴ等设备。工业方面，７０年代末起美国利用研制的透射式工业ＣＴ设备对军工产品的关键部件做无损检测，美制造，开展逆向工程学的研究。国科学测量系统公司的工业咖还应用于ＣＡＤ／洲方面，进行加工元件的仿型第１章绪论自上世纪７０年代ＣＴ问世以来，各国的ＣＴ生产商不断研制新技术，新产品如雨后春笋，ｘ射线ＣＴ已从第一代发展到第四代，扫描时间也重最初的十几分钟大大缩短Ｎ５秒左右，重建算法也日趋成熟，图像质量也得到了很大的改善。国内２０世纪９０年代初自行研制成功了第一台透射式Ｙ射线工业ＣＴ样机，清华大学等单位研制成功ｘ射线微焦点ＩＣＴ机，东大阿尔派先后研制成功普通医用ＣＴ机和螺旋ＣＴ机。由于医用和工业ＣＴ的巨大市场需求，使得国产ＣＴ的研发受到高度重视，进入了一个良好的发展时期。１．１．２ＣＴ成像原理简介Ｘ射线ＣＴ通过控制Ｘ光球管围绕被测物体扫描，Ｘ射线沿各方向的直线穿过物体的横切面，ｉ贝９得ｘ射线的衰减数据，在计算机上利用图像重建原理计算得到物体横切面上的各点对Ｘ射线的吸收系数，并在计算机或是胶片上显示出来。ＣＴ扫描射线如图１．１（ａ）所示，光源强度为Ｉｏ的ｘ射线穿透一部分均匀的物体，得到衰减后的光强为Ｉ，若射线穿过物体的长度为Ｙ，则根据Ｂｅｅｒ定理有：ＴＴｆｕｄｌ＝ｌｎ（专ｑ－）或ｕｙ＝ｌｎ－孥－ｔ１１（１．１）其中，ｕ为物体对ｘ射线的吸收系数。如果物体内部分段均匀且各段的吸收系数分别为Ｕｌ，Ｕ２，Ｕ３，…Ｕｎ，相应的长度分别为ｙ１，ｙ２，Ｙ３，…Ｙｎ，则公式（１．１）改为：ｕｌｙｌ＋ｕ２ｙ２＋ｕ３ｙ３＋…＋Ｕｎｙｎ＝ｌｎ（Ｉｏ／１）。更一般地，如果物体在（ｘ，”平面内分布不均匀，吸收系数ｕ＝Ｕ（ｘ，ｙ），那么沿着某路经Ｌ的衰减为：Ｔｔｆｕｒｌｌ＝ｌｎ（、＇．Ｏ）ｌ（１．２）我们将式（１．２．）中的线积分称为ＣＴ投影定理，若不指定路线，则称为投影，它是某一角度下射线引起的射线投影的集合。ＣＴ设备克服了线积分的弊病，得到了反映人体组织或是物体结构分布的图像。总的来讲，ＣＴ的基本原理是采用准直后的ｘ射线束对人体的某一层面从不同的角度进行照射，在射线穿过的另一端利用探测器接收多组原始数据，经过计算机重建后得到用于显示的二维数据矩阵，通过显示设备显示图像。图１．１（ｂ）简单说明了ＣＴ的扫描方式。图中扫描器发出一组ｘ射线，射线穿过被测区域，在对面通过探测器接收衰减后的Ｘ射线得的射线强度，即投影值。然后再旋转一个小角度进行扫描，如果以此方式作Ｍ次扫描将得到Ｍ组投影数据。２第１章绪论器Ｉ（ａ）射线透过物体图１．１ＣＴ扫描示意图在ＣＴ成像技术中，有两个很重要的概念，就是ＣＴ值与ＣＴ图像和窗宽与窗位。ＣＴ图像是通过数学方法对ＣＴ原始数据进行重建，得到图像矩阵。计算机把重建图像矩阵中的各个像素转变为不同明暗的相应光点，通过显示设备显示出来。ＣＴ值是ＣＴ采用的标准是根据各种组织和物质对ｘ射线的吸收系数来决定的。Ｈｏｕｎｓｆｉｅｌｄ将线性衰减系数划分为２０００个单位，医学上称为ＣＴ值。由于在物理过程中，物质的密度是由物质对ｘ射线的衰减系数来体现的，在研究ＣＴ图像时，能提供诊断信息的是组织之间的密度差异，而不是绝对密度，所以一般定义ＣＴ值如下：ＣＴ值．１０００×—／＇／ｘ－—／＇／ｈ２０纥２Ｄ（１．３）即特定物质的ＣＴ值等于该物质的衰减系数与水的吸收系数之差再与水的吸收系数相比，然后乘以１０００。物质的ＣＴ值越高，表明其密度越大，所以金属物体的ＣＴ数远远大于骨骼和其它组织的ＣＴ数。另外，ＣＴ值的大小还与ｘ射线的能量有关，Ｘ射线的能量越低则相同物质的ＣＴ值就越大。目前，大多数ＣＴ的Ｃ－Ｉ＂数变化范围是．１０００到＋３０００，而计算机显示器一般只能显示２５６个灰度级，人眼只能分辨３０个灰度级，因此，为了能用有限的灰度级较好地表示ＣＴ图像，人们引入了窗口技术，即通过正确地选择窗宽和窗位来显示符合人眼视觉要求的图像，尤其是在临床医学上，能更好的显示医生感兴趣的区域。如图１．２所示，窗宽是ＣＴ数的显示范围，窗位是显示灰度级的中间位置所对应的ＣＴ数。例如，如果窗宽时１００，窗位是１００，那么ＣＴ数的范围是５０至１Ｊ１５０，其３第１章绪论中ＣＴ数１００对应的显示器灰度级是１２８，ＣＴ数５０对应的灰度级为０，而ＣＴ数１５０对应的灰度级为２５５。小于５０的ＣＴ数的灰度级全是０，大于１５０的ＣＴ数灰度级全是２５５。如果增大窗宽可以使一个灰度级对应多个ＣＴ数，这样会使ＣＴ数接近的组织对比度降低，反之，如果减小窗宽会使多个灰度级对应一个ＣＴ数，导致对比度增加，使被扫描物体结构更清晰，但也会丢失部分信息，所以在运用时，要反复调节窗宽和窗位，来获得最佳视觉效果。Ｏ窗一窗唑～图１．２窗宽与窗位ＣＴ数把如图１．２这种将显示的灰度范围和窗宽内的ＣＴ数按线性关系一一对应的窗口称为线形窗，线性窗实现简单，可手动动态改变窗宽和窗位，是商用ＣＴ控制软件中常见的功能。１．１３ＣＴ设备的主要技术指标ＣＴ的主要技术指标包括：ＣＴ的工作时间、ＣＴ系统部件和ＣＴ图像质量指标，这些指标与ＣＴ性能息息相关。（１）ＣＴＴ作时间ＣＴ工作时间包括：扫描时间、图像重建时间、周期时间。扫描时间短，一方面可以缩短病人暴露在ｘ射线中的时间；另一方面可以减少由于病人运动而引起的运动伪影，确保图像质量。ＣＴ扫描时间在医学中采用的指标，是指在能够获得设备所提供的图像质量的前提下，所需要的最短扫描时间；图像重建时间是指根据原始数据进行计算机处理，从而获得ＣＴ图像所需要的时间；重建时间短有利于及时观察扫描图像，迅速调整扫描策略。一般地，重建时间与重建矩阵４第１章绪论的大小密切相关，显然重建矩阵为５１２ｘ５１２时所需要的时间比重建矩阵为２５６ｘ２５６时要长。对病人的某一断层从扫描开始，经重建、显示到拍片完成整个过程所需要的时间称为周期时间。一般，周期时间＝扫描时间＋重建时间＋拍片时间。（２）ＣＴ系统部件ＣＴ设备的系统部件包括：扫描机架和病床、Ｘ射线系统、数据采集系统、计算机及显示系统、照相系统、数据存储系统等。ｏ）ｃｒ图像质量指标空间分辨率：指ＣＴ对物体空间大小及几何尺寸的分辨能力。通常有两种表示方式：一种是用单位厘米内的线数来表示；另一种是用可辨别最小物体的直径用毫米来表示。密度分辨率：表示ＣＴ设备对密度差别的分辨能力，通常以百分数表示。比如０．４５％表示两种物质的密度差大于０．４５％时，ＣＴ可以将它们区分出来。伪影：在被扫描取样的物体中并不存在的，但在图像中却出现的各种不同类型的影像。一般地，按伪影产生的方式分为ＣＴ设备伪影和病人运动伪影两种。１．１．４ＣＴ图像的特点介绍了ＣＴ的成像原理及主要技术指标后，在此基础上来分析ＣＴ断层扫描形成的图像的特点。ＣＴ图像是通过计算机计算出来的ｘ射线衰减值的二维分布图，是由一定数目的像素按矩阵排列而成的二维断层图像，这些像素反映了相应单位容积的ｘ射线吸收系数。目前常用的ＣＴ装置的像素大小为２５６ｘ２５６、３２０ｘ３２０和５１２ｘ５１２。ＣＴ图像一般是８位或１６位的灰度图像。与其它很多成像系统一样，ＣＴ成像也受到噪声的干扰。深入了解ＣＴ成像过程中受到的各种干扰以及各种干扰在ＣＴ图像上的表现，有助于在ＣＴ图像的处理过程中消除伪影的影响，提高图像处理的质量。ＣＴ伪影和噪声是评价ＣＴ成像质量的重要指标，下面就对伪影和噪声对ＣＴ图像的影响进行分析，以最大程度地减少两者对ＣＴ图像处理的影响。在ＣＴ图像上非真实的阴影或干扰称为伪影。它降低图像的质量，易造成误判或不可诊断。它可分为由被扫描物人或物引起的伪影与ＣＴ设备本身所造成的伪影两大类。由被扫描物人或物引起的伪影中，最常见的是在扫描取样时，病人可能体内有金属物体引起的放射状金属伪影，如图１．３（ａ）所示。病人在扫描时自主或不自主的运动引起的运动伪影。均匀物体的ＣＴ影像中各像素的ＣＴ值参差不齐，它使图像呈颗粒性，直接影响其密度分辨率，尤其以低密度的可见度为甚。把这种现象用统计学上的标准偏５第１章绪论差方式表示出来，即为ＣＴ的噪声，可分为随机噪声和统计噪声。它们的产生机理和对图像质量的影响各不相同，严重的随机噪声往往形成了伪影，如图１．３（ｂ）所示。在ｃｒＹ２像中，ｘ射线光予数是一个随机起伏的统计量，即使电压非常稳定．其它一切条件都不变，表示ｘ射线强度的光子数也是随着时间统计而起伏。此外检测器和放大器中也有随机起伏的噪声（张伟，２００６）。所有的这些构成了随机测量误差，反映在投影值上，它使投影值曲线上产生单个向上或向下的小尖峰。由于在ＣＴ图像的重建中，每个测量值都被反投影到整个图像上，所以这个单尖峰就变成重建图像上的一条条纹。当有许多随机噪声的小尖蜂时，重建国像上就出现许多杂乱的细小条纹和由许多条纹交叉形成的颗粒斑点，成为图像随机噪声。统计噪声又称像素噪声，指服从泊松分布的噪声，主要由ｘ射线的强度和数量的起伏，以及ｘ射线穿透人体时要产生散射和吸收所形成的ｆＹａＶＵＺＭｅｔａｌ，１９９９）。糜“）含金属伪影的ｃＴ豳像佃）含严重噪声的ｃＴ图像图１．３古伪影和噪声ＣＴ图像１．２研究意义和国内外研究现状随着ｃ１技术的发展，ＣＴ已经成为医疗诊断和工业检测等行业的一种非常重要的技术。但是在ｃｒ普及的过程中出现很多有待解决的新问题，ｃｒ图像中的金属伪影就是其中之一。理论上伪影可以定义为图像中被重建数值和物体真实的衰减系数之间的差异（Ｋａｌｅｎｄｅｒｄａｌ，１９８７）。尽管定义足够宽泛，但它几乎包括所有的非理想图像，所以我们需要对特定的伪影表现进行说明。条状伪影经常表现为横穿图像的明显直线ｆ－傲不是平行的）’它们可以是亮的也可以是暗的。由于重建算法中滤波器的特性，明暗条纹是成对出现的（ＧｌｏｖｅｒＧＨａｎｄＰｅｌｔＮＪ，１９８１）。由于这种伪影般由金属物体引起，所以又常称为金属伪影。无论在医疗诊断还是工业检测都经常遇到被检测物体中含金属物质，或是其它吸收系数很高的物质，导致投影数据出现了跃变，重建后的图像就会包含明暗６第１章绪论相问的金属伪影，这些伪影一般数量众多、幅值又高，严重降低图像质量（Ｂａｉｓｓａｌｏｖ的热点。ｅｔａｌ，２０００），无论是在医疗或是工业检测都影响了图像的后续处理和分析，如何消除这些条状伪影以提高图像质量，已经成为当前ＣＴ应用研究领域中用于消除ＣＴ图像中金属伪影的方法来自学术机构或是工业界，在不同的杂志或是会议论文集上可以找到许多研究论文和文章。近几年一些比较有效的ＣＴ去金属伪影算法的得到了快速发展，但都存在一定的缺陷。国外宾夕法尼亚大学的ＬｅｗｉｔｔＲＭ和ＢａｔｅｓＲＨ（１９８３）提出了用一种特殊的函数进行插值的方法，后来又提出了Ｃｈｅｂｙｓｈｅｖ多项式的差值法；Ｃｒａｗ－ｆｏｒｄ（１９８６）在线性插值的基础上添加了一些辅助处理；华盛顿大学ＷＡＮＧＧ等人（１９９６）提出的迭代代数重建法。芝加哥大学的ＸＩＡＤａｎ（２００５）Ｊ丕提出了一种ＥＭ（ＥｘｐｅｃｔａｔｉｏｎＭａｘｉｍｉｚａｔｉｏｎ）局部迭代的混合算法。还有一些人提出了马尔可夫随机场方法，算法基于传播最大似然的思想，假设测量用的探测器读数服从泊松分布，通过最优似然来重建图像。国内这方面的研究还比较少，主要是清华大学的谷建伟提出的快速线性插值法（谷建伟等，２００５）和基于伪影权重的非线性插值法（谷建伟等，２００６ｂ），以及林宙辰等（２００１）提出的采用多项式插值。现在主流去金属伪影算法分为三类：插值重建算法、迭代重建算法和混合重建算法。插值算法简单快速，但是只适用于金属结构简单的情况，而且会使金属周围的区域失真。迭代法也是一种常用的ＣＴ图像重建算法，能有效去除金属伪影和抑制噪声，而且能很好的呈现金属物体的结构，但其运算量非常大，速度很慢，难以实用化。混合法是近几年提出的一种新算法，ＥＭ混合法（ＸｎＤｅｔａｌ，２００５）是其中代表，该算法速度较快，能有效反映金属物体结构，不足之处是金属周围的区域失真，使其在医用咖用中存在很大的局限性。了详细的研究，并设计了一种新的去金属伪影算法。为了克服以上这些算法的不足，本文对金属伪影成因和去金属伪影算法进行１．３本文的主要工作及内容安排论文通过研究ＣＴ的成像和重建原理、金属伪影成因和现在主流的去金属伪影重建算法，并在这些研究的基础设计一种新的金属伪影消除算法，各章安排如下：第１章介绍了ＣＴ的发展历史、成像原理、主要技术指标和课题的研究背景、意义及国内外研究现状，针对目前现在主流去金属伪影算法的不足，确定了本文的研究方向。７第１章绪论第２章介绍了ＣＴ重建的理论和数学基础，并说明了如何通过变换法和数级展开法进行ＣＴ图像重建。变换法中主要介绍了滤波反投影重建算法，数级展开法（迭代法），重点介绍了代数迭代重建算法（ＡＲＴ）、联合代数迭代重建算法（ＳＡＲ＇ｒ）、期望最大值迭代算法（ＥＭ）和有序子集期望最大值迭代法（ＯＳＥＭ）。第３章重点研究了ＣＴ中金属伪影的成因和现在主流的去金属伪影算法：插值重建法、迭代重建法和ＥＭ混合重建法。通过Ｍａｔｌａｂ构建数值仿真平台，在仿真平台的基础上对这些算法进行了去金属伪影实验并对它们的实现方法进行了优化，最后通过实验结果总结了这些算法的优点和不足，并提出了去金属伪影算法的改进思路。第４章针对第三章中提到的主流去金属伪影算法的不足，设计出一种去ＣＴ图像中金属伪影自适应衰减滤波混合法。对自适应衰减滤波混合法每一个步骤进行了详细的介绍，并通过数值模拟对比实验证明，该方法可以有效的去除金属伪影，并很好的保留了金属及其周围组织的信息，其效果优于插值法和ＥＭ混合法。特别是在多金属物体的情况下，算法的效果更是优于其它算法。第５章对本文研究的主要内容、创新点进行了总结，并对未来的研究工作的方向进行了展望。８第２章ＣＴ成像基本原理与图像重建算法第２章ＣＴ成像基本原理与图像重建算法２．１经典的断层成像技术为了克服传统的Ｘ射线摄影技术中的图象重叠问题，在ＣＴ发明以前已有经典的断层成像技术出现。这种成像技术沿用至今已有６０余年，并发展成许多模式。下面介绍其中的一种——纵轴直线断层成像技术。纵轴直线断层成像技术是这样一种成像技术，它对被选中的感兴趣断层清晰的成像，而是这一层的上下层模糊可见或是基本看不见。传统的断层成像技术如图２．１，它表示经典的纵轴直线断层成像技术的基本原理（黄庆华，２００２）。被扫描物体或病人位于点状ｘ射线源和检测器之间。ｘ射线源的运动被固定在某一平面上，这一平面被称为源平面。检测器所在平面也是被固定在某一平面，一般这个平面就是感光底板Ｐ。图２．１中，源平面和感光底板相互平行，ｘ射线源Ｘ与检测器ＡＢ作反方向运动，它们的速度应满足ｘ射线源Ｘ到检测器上的某点Ａ所作的直线始终通过断层中的点Ａ。同理Ｘ射线源Ｘ到检测器上的某点Ｂ所作的直线始终通过断层中的点Ｂ，则ＡＢ就是感兴趣断层。图中Ｘ１和Ｘ２分别是射线源Ｘ在ｔ１、ｔ２时刻的位置，Ａ１、Ａ２和Ｂ１、Ｂ２分别是检测器平面上两点Ａ、Ｂ在ｔ１、ｔ２时刻的位置。由它们运动速度的约束可知，直线Ｘ１Ａ１和Ｘ２Ａ２必交于点Ａ，直线Ｘ１８１和Ｘ２８２必交于点Ｂ。这样曝光期内，可以使断层ＡＢ平面的图像在感光板Ｐ上不断加强而清晰成像。Ｘ１Ｘ２射线源Ｘ移动方向Ａ２ＡｌＢ２Ｂ１图２．１传统的断层成像技术９第２章ＣＴ成像基本原理与图像重建算法２．２ＣＴ图像重建的问题描述图２．２表示图像重建中的投影坐标系，设一个具有某种物理性质（吸收系数）的二维分部函数ｆ（ｘ，Ｙ）定义在平面域０内，ｐ（ｘ’，口）表示ｆ（ｘ，Ｙ）沿着由Ｘｗ和０定义的一条直线的线积分，穿过厂＠，ｙ）的所有线积分的集合为：ｐ（ｘ’，日）２Ｊ＝。ｆ（ｘ，ｙ）ｄｙ’。Ｊ＝。ｆ（ｘ’ｃｏｓＯ—Ｙｓｉｎ０，∥ｓｉｎ０＋ｙ’ｃｏｓ０）ｄｙ’（２．１）‘＝ｒｃｏｓ（Ｏ一驴）ｘ’Ｅ（－ｏｏ，＋∞）ｐ∈【０，万）对于一个固定的角度０，ｐＯ’，口）是ｆ（ｘ，ｙ）沿口方向在ｚ’轴上的一个一维平行投影。如果目标函数厂＠，Ｙ）未知，而它的所有投影（基线积分）却能通过某种方式获得，则ｆ（ｘ，Ｙ）便可由ｐＯ’，０）唯一的确定，这个问题是Ｒａｄｏｎ在１９１７提出的，因此，（２．１）式中联系的目标函数和它的投影函数之间的积分变换关系被称为Ｒａｄｏｎ变换。从数学上看，图像重建问题就是求解Ｒａｄｏｎ变换的反变换，即给定投影数据ｐ（ｘ’，０），求出ｆ（ｘ，Ｙ）。而这个问题在工程上称为从投影重建图像，或ＣＴ图像重建。Ｐ（ｘ、，Ｏ）ｙ，／，，Ｘ彦雩、、、、、、、、）７，，，，７７／文■～／＼、、、、、、ｏ／，Ｋ』、≮７，，，，，Ｊ。‘Ｉ／／一、、、、女、、图２．２投影的坐标系２．３ＣＴ重建算法从投影重建图像的算法很多，如：反投影重建算法、滤波反投影重建、直接傅立叶变换重建算法和迭代重建算法。但是反投影重建算法和直接傅立叶变换重建算法都存在明显缺陷，很少使用。反投影重建算法图像中存在密度为零的点，所以重建后不一定为零，从而引入星状伪影。由于傅立叶空间中产生的采样模式１０第２章ＣＴ成像基本原理与图像重建算法不是笛卡尔坐标，为了执行二维傅立叶反变换这些采样不得不重新插值或重新栅格化到一个笛卡尔坐标系中，在傅立叶空间中一个采样点的误差会影响到整个图像，所以直接傅立叶变换重建算法重建图像效果不好。因此将本节重点介绍最常见的滤波反投影重建和迭代重建算法。重建算法在具体实现时，会因为ｃＴ扫描方式的不同略有不同，但整体思想还是不变的。现行ＣＴ主要扫描方式有平行束扫描和扇形束扫描。迭代重建算法不会过分依赖于扫描方式，而滤波反投影重建算法相对较依赖于扫描方式。２．３．１滤波反投影算法２．３．１．１Ｒａｄｏｎ变换和中心切片定理Ｒａｄｏｎ变换及其逆变换和中心切片定理是滤波反投影重建（Ｆｉｌｔｅｒｅｄ础。本节将详细介绍Ｒａｄｏｎ变换和中心切片定理。ＢａｃｋＰｒｏｊｅｃｔｉｏｎ，ＦＢＰ）的理论基础，而且Ｒａｄｏｎ变换及其逆变换还是迭代法的数学基由１．１．２节可知ＣＴ投影数据是衰减系数的线形积分，图像重建问题其实就是是否可以通过投影数据估计物体内部的衰减系数（ＨｅｒｍａｎＧＴ，１９９５），Ｊ．Ｒａｄｏｎ于１９７１年给出了肯定的答案，就是Ｒａｄｏｎ变换及其逆变换。Ｒａｄｏｎ变换可在任意维空间定义，且定义存在多种形式。在ＣＴ扫描中只考虑二维情况。那么二维函数ｆ（ｘ，Ｙ）的Ｒａｄｏｎ变换定义为：Ｒ｛ｆ（ｘ，ｙ））＝ｒｆ／（ｘ，ｙ）６｛ｌ－（ｘｃｏｓＯ＋ｙｓｉｎＯ）｝ｄｘｄｙ一ＲｏＵ）＝Ｒ（１，口）（２．２）式中ｚ代表坐标原点到直线Ｌ的距离，表达式Ｚ－ｘｃｏｓＯ＋ｙｓｉｎ０代表直线Ｌ，沿着一系列平行线（投影线）的积分就组成了投影Ｒｏ（ｆ），所有投影组成的集合｛也Ｕ），目∈【０＇幼））就是Ｒａｄｏｎ变换或者称为投影。其中积分域为整个ｘｙ平面，０为距离Ｚ与ｘ轴的夹角，６为脉冲函数。６ｃｚ，＝ｆ三二二吕函数，如图２．３所示。（２．３，当夹角０为固定值时，那么函数ｆ（ｘ，ｙ）在夹角０下的投影值是距离Ｚ的一维第２章ＣＴ成像基本原理与图像重建算法ｙｌ／／／图２．３Ｒａｄｏｎ变换示意图一ＸＲａｄｏｎ变换可以认为是Ｚ一０空ｆｎ－ｊｑ丁的投影，Ｚ一０空间中的每一点对应一条直线，而Ｒａｄｏｎ变换就是图像像素点在每一条直线上的积分。因此，图像中灰度级高的像素点在Ｚ一０空间中形成亮点（ＣＩ’图像中金属物体是亮点），而低灰度级的像素点在Ｚ一０空间中为暗点。根据Ｒａｄｏｎ变换的定义在Ｚ一０空间中有：ＸＩＤｙ暑瞄ｐ螂口口一ＳＳｎ＋ＳＳｎ口口（２．４）则Ｒａｄｏｎ变换的另一种表达形式为：Ｒ｛ｆ（ｘ，ｙ）】－２ｆｉ（ｐｃｏｓｏ—ｓｓｉｎ０，ｐｃｏｓＯ＋ｐｓｌｎ口灿其示意图如图２．４所示。（２·５）忒／ｓ＼雾弋一＼／＼心图２．４Ｒａｄｏｎ变换的另一种表达示意图１２第２章ＣＴ成像基本原理与图像重建算法为了求解ｆ（ｘ，ｙ），Ｒａｄｏｎ在１９１７年对公式（２．２）进行求逆运算，这就是著名的Ｒａｄｏｎ反变换，求解的，＠，ｙ）为：舷加昙ｒｄ啦南学讲论还不能设计滤波反投影算法，还需要中心切片定理得辅助。（２．６）式中Ｘ＝伍，Ｙ），宇＝（ｃｏｓＯ，ｓｉｎＯ）且歹ｐ，亭）一Ｃ，（ｚ亭＋｝亭上）ｉｔ。只有Ｒａｄｏｎ变换理该定理在非衍射情况下，中心切片定理示意图如图２．５所示，某图像ｆ（ｘ，ｙ）在视角为妒时的投影Ｐ西（Ｘｒ）的一维傅立叶变换给出厂Ｏ，Ｙ）的二维傅立叶变换Ｆ２（ｗ１，ｗ２）－Ｅｐ，妒）的一个切片，切片与Ｗｌ轴相交成驴角，且通过坐标原点。即互【所“）】＝最ｐ，妒）ｋ定（２．７）图２．５中心切片定理示意图式中，Ｆ·【·】表示以为傅立叶变换。下面我们简单的证明一下该定理：先证当７＝Ｏ时，此时Ｘｔ－－Ｘ，所以投影儿以）５ｐｏ（Ｘ）。ｆｆ（ｘ，ｙ）ｄｙ（２·８）ｇ［ｐ。Ｇ）】２ｆＰ。（ｘ）ｅ－ｉＷ，Ｘｄｘ２ｆ［ｆｆ（ｘ，ｙ）ｄｙ］ｅ一叩ｄｘ２ｆｆｆ（ｘ，ｙ）ｅ’呻ｄｘｄｙ（２．９）而ｆ（ｘ，ｙ）的二维傅立叶变换为Ｆ（Ｍ，％）＝Ｅ【厂Ｇ，．），）】＝ｒｆｒＣｘ，ｙ）ｅ一““ｗ’ｄｘｄｙ（２．１０）第２章ＣＴ成像基本原理与图像重建算法比较公式（２．８）和（２．９）可以得到：Ｅ０。＠））＝Ｅ（嵋，％）１ｗ，．。＝Ｅ（Ｍ，０）（２．１１）也就是说ｙ（ｘ，ｙ）沿Ｙ轴的～维傅立叶变换给出了厂Ｏ，ｙ）的二维傅立叶变换Ｖ２（Ｗｌ，ｗ２）在ｗ２－Ｏ处的一个切片。若投影的视角为妒，利用旋转坐标一一Ｙ，，使积分路径总是沿着ｙｒ。此时ｘ，轴与Ｘ轴之间的夹角为≯，如图２．５所示，并有关系［量］。［一ｃｓｏ；ｓｎ＃矽∞ｓｉｎｓ＃≯．｜］。［ｙｘ］以ｙｒ作为投影轴的投影为：ｃ２．１２，功“）。ｒ，＠，ｙ）ａｙ，＝ｒｆ“，Ｙ，）方，它的一位傅立叶变换为：（２．１３）壬式ｐ◆ｂｒⅥ５ｆｐ÷ｑｒｙ却ｒｄｘｒ；｛ｌｆｒｑｒ，Ｙｒ）ｄｙｒ·ｄ２弼ｘｒ缸ｒ。ｆｆｆ（ｚ，ｙ）∥挪ｏ。∞妒郴如妒’ｄｘｄｙ（２．１４）在推导式（２．１３）的过程中用到了关系式ａＥａＥａｘ红方，＝砂砂，砂砂，ａｘ蚴２瞄ｃｏｓ＃妒耋纠蚴；螂ｃ２胚，（２．１６）以及ｆ“，Ｙ，）＝，＠，ｙ）。由式（２．１４）可知，令Ｍ一２死ｐｃｏｓ＃，％一２ｚｐｓｉｎ＃，代入式（２．１４）有可Ｄ“）】＝ｐ≯（ｐ）一Ｆ２（ｗ，，％）＝Ｆ２（ｐ，驴）其中驴＝ａｒｃｔｇ（Ｗ２ＪＷｌ）式（２．１６）说明：Ｗｌ，ｗ２并不是独立的，而是受到直线≯一ａｒｃｔｇ（ｗ２／ｗ１）的约束。这样就证明了投影定理。２．３．１．２反投影算法的数学描述我们把“取投影”一“反投影重建”呻“重建后图像’’这些环节看作一个以原像为输入重建后图像为输出的成像系统，如图２．６所示。为了仔细剖析这个系统，先求该系统的点扩展函数。１４第２章ＣＴ成像基本原理与图像重建算法图像ｉ（ｘ，Ｙ）取投影反投影重建垂Ｌ图２．６反投影重建等效成像系统设一位于坐标原点ｘ＝０，ｙ＝０的点源ｄ（ｘ，ｙ）为ｘ．ｙ断面中的唯一像点。并设扫描方式为平行／旋转。即射线先平行移动，从物体的一侧移动到另一侧；让后旋转一个角度△妒，再加上作一个平行移动。之后又旋转一个角度△妒。如此继续，直到累计的转角达到１８０．△≯度为止。我们还要设置一个旋转坐标系统Ｘｒ—Ｙ，，见图２．３，Ｘ—Ｙ为固定坐标，Ｘｒ—Ｙ，为旋转坐标，（ｒ＇？）为极坐标。它绕原点转动使投影线总时沿着ｙｒ方向，Ｘ—Ｙ坐标的原点与砟一只坐标的原点重合。两者的夹角为≯。不同的≯代表不同的投影角度。射线位置可由坐标化，爹）决定。空间任一点的坐标可用（ｘ，ｙ），“，只）或（ｒ，？）表示，设驴为离散取值，如取≯一唬，则相应的投影为：＾“）‘ｐ“，晚）－。ｆＬ（ｘ，，Ｙｒ）ｄＹｒ。』６似，Ｙ，）ｄｙ，。ｆａ似）６（ｙｒ）ｄｙ，一６（ｊｒ，）ｂ．霸一ｔ５（ｒｃｏｓ（０一唬”（２．１７）根据反投影的定义原像中某点的像素值等于通过该点的所有射线的平均投影值。所以点（Ｘｒ，Ｙ，）的图像在所述坐标系统中表示为：Ⅳ．１Ⅳ．，（，，口）＝ｆ“以）。袁善ｐ呜Ｌ，·ｒ。ｏｓ（¨）。妻善ｐ呜【，ｃｏｓｐ一谚）】△驴（２．１８）式中，△≯＝ｐ／Ｎ：，Ｍ为投影数。若在有限区间内将射线增至不相重的无限条，即连续取投影，则式但．４）过渡到更一般的连续情况下的反投影表达式为：１耳ｆ（ｒ，口）一二ｐ≯【，ｃｏｓ（Ｏ·妒）ｐ妒刀％（２．１９）式中，积分区间为（０，ｐ），因为在忽略射束硬化的情况下，妒在Ｑ，２ｐ）区ｌ＇ｆｉ］内的投影值等于（０，ｐ）区间内的投影值。２．３．１．３正弦图正弦图是反投影重建的一个重要概念，上一节中提到（ｒ，？）表示断层中某一点的坐标，而“，驴）表示投影线的位置，它对应于一个射线投影值。算法实现时，１５第２章ＣＴ成像基本原理与图像重建算法这些投影值存储在计算机的存储区中。反投影算法要求对任一点（ｒ＇？）找出经过该点的所有射线，从而求得对该点的所有射线投影之均值。给出定点（ｒ’目）的线坐标（Ｘｒ，妒）轨迹如图２．７所示，过给定点（ｒ，？）的射线阮，≯）均以ｒ为直径，以（ｒ／２，？）为圆心的圆周上。换句话说，过（ｒ，？）的射线＠，，≯）满足方程：‘＝ｒｃｏｓ（Ｏ一驴）（２．２０）图２．７过给定点（ｒ，口）的线坐标（‘，≯）轨迹如以≯为纵坐标，ｔ为横坐标，即在Ｏ，，≯）空间中给定（ｒ，？）的情况下，由式（２．２０）画出的一条正（余）弦曲线，见图２．８。曲线经过之点就是存储器中的某一地址。实际使用时，妒以？为增量，ｘ，以为增量，则（肌△，ｎＯ）对应于存储器中的地址（ｍ，ｎ均为整数），若？，ｄ足够小，把这些地址中的内容取出相加再取平均，就可足够精确地认为是该点的重建像素值。图象中的每一点，对应于＠，，妒）空间中的一条正弦曲线。全部图像对应于一簇正弦曲线。即是正弦图（ｓｉｎｏｇｒａｍ）。正弦图是一个十分重要的概念，是从投影重建图像的基础（庄天戈，１９９２）。Ｊｌ由７ｃｎ△～．／△／７０，ｄ‘。‘。’—ｎｄＸｒ／／／７ｉ，，Ｊｌ、、、～～图２．８空间点（ｒ，？）＝（１，ｐ／４）及相应的正弦图１６第２章ＣＴ成像基本原理与图像重建算法２．３．１．４滤波反投影重建在Ｒａｄｏｎ变换和正弦图的基础上直接反投影重建的图像，会引入星状伪影，即原来图像中密度为零的点，重建后不一定为零，导致图像失真。但是我们可以通过对重建后图像进行滤波，去除星状伪影。图２．９（ａ）表示成像系统框图，但是即使使用近似，二维滤波器的实现还是很困难的。如果能如图Ｚ９（ｂ）先滤波再反投影，在反投影重建以前对投影数据进行滤波，再把滤波后的数据进行反投影重建，那么不仅消除了星状伪影，而且滤波器变成了一维，简化了运算。这个假设是可以实现的，因为反投影重建属于积分运算，而对投影数据进行滤波属于卷积运算，两者都是线性运算，运算次序是可以互换的。先对投影数据滤波，再进行反投影重建就是常说的滤波反投影重建算法。（ａ）先“反投影”再滤波系统框图（ｂ）先滤波再“反投影”系统框图图２．９滤波反投影重建方案滤波反投影重建的关键是滤波运算和反投影运算的顺序互换，它的理论基础就是２．３．１．１节中进行了介绍投影定理或称中心切片定理。根据上面的投影定理，滤波反投影算法重建可以按以下流程进行：（１）根据采集不同角度下的投影数据，若探测器数量较多可只采０到１８０度。（２）对每个投影角度的投影数据进行傅立叶变换。（３）将每个角度的傅立叶变换后的投影数据组成原图像的二维傅立叶变换。（４）通过二维傅立叶反变换得到原图像。具体实现算法时，关键是选择哪种数学方法来处理反变换，滤波反投影重建算法一般通过变换法求解反变换。如图２．１０中所示，待建图像为ａ（ｘ，ｙ），它的二维傅立叶变换表示为Ａ：（ｗ。，Ｗ：）－－ｈ：ｐ，≯），根据中心切片定理，Ａ：∽，≯）可通过ａ（ｘ，ｙ）在不同的视角≯－－Ｆ的投影ｐ。纯）的一维傅立叶变化求得即：Ａ：（ｗ。，Ｗ２）＝Ａ：（ＪＤ，妒）＝Ｆｐ≯（‘）］－－ｐ≯（ｐ）＝ｐ（ｐ，妒）（２．２１）１７第２章ＣＴ成像基本原理与图像重建算法Ｘ图２．１０滤波反投影所用坐标系这时待建图像可表示为：枷）－口㈨Ｍ愀懒）】－专玢（恍妙唧，ａｗｄｗ：；ｆｆＡ（Ｐ，驴∥砌咧疗一’’ＩＰＩａＰ，枷５ｆｆｐ（ｐ，≯矽枷州８一声’ＩＰｌｄｐｄ≯２户驴．厂ｐ（ｐ，∥砌鲫－’Ｉｐｌｄｐ此有：（２·２）式（２．２）推导中用到了在ＣＴ重建问题数学描述中提到的等式‘＝ｒｃｏｓ（Ｏ一驴），因ｗｌｘ＋ｗ２ｙ＝２ｘｐ（ｘｃｏｓ≯＋ｙｓｉｎ）一２：ｒｐｘ，＝２死ｐｒｃｏｓ（Ｏ一妒）而且ｄｗ。ｄｗ２＝ｌＪＩｄｐｄｑＯａｗｌ】；Ｏｐ帆ａ驴Ｏｗ２Ｏ＝ａ％Ｏｐ妨纫Ｓ｛弓ｍａ砂对硪２．１７），我们先研究第二个积分嚣盆嚣矽）－Ｅ—ｐ∽力）】ｆＰ（Ｐ，≯矽抑７酬８—９’Ｉｐｌｄｐ；ｆｌｐｌｐ（ｐ，≯∥挪Ｘ，ｄｐ＝ＪＩｌ似）水ｐ（ｘ，≯）Ｌ．，ｃｏｓ（州）１８第２章ＣＴ成像基本原理与图像重建算法＝ｇ（ｘｒ，≯），，．，ｃｏｓ（日．ｐ）式中：（２．２３）ｇ以，≯）＝Ｐ（毫，妒）母＾（‘）（２．２４）而Ｊ｝ｌ以）一互“【ｐ】，ｐ以，≯）一互“ｐ（ｐ，妒）】，式（２．１８）的物理意义就是投影数据ｐ“，≯）经过传递函数为ＩＰＪ＝日＾佴）】的滤波器后的修正投影值，且满足‘；ｒｃｏｓ（Ｏ一妒）时的值。而‘一ｒｃｏｓ（Ｏ一妒）恰好是通过给定点（，．，口）的射线方程。将式（２．１８）代入式（２．１７）中得玎石（，，日）＝ｆｇｃｏｓ（０一≯），≯Ｆ≯ｒ（２．２４）气该式的物理意义是：对于给定图像点（，．，口）的像素值等于通过该点的所有射线滤波后的投影值的累加。这也是滤波反投影重建算法的核心思想。在计算机中实现的各个步骤可以分解出来用图２．１１表示，具体地说分如下三步：（１）选定滤波函数（一般选用Ｓ．Ｌ滤波函数，其去噪效果好，重建质量高），对每个固定视角谚的投影数据ｐ“，谚）进行滤波，得到滤波后的投影数据ｇ“，谚）。（２）对于每个谚，如果点＠，谚）满足‘一ｒｃｏｓ（Ｏ一谚），则把滤波后的投影值ｇ阮，谚）累加到该点的像素值中。（３）将步骤（２）中得到的累加像素值就是重建后的图像。图２．１１滤波反投影重建算法过程框图２．４．２迭代重建算法由上一节公式（２．６），Ｒａｄｏｎ逆变换是由关于投影数据的导数算子、Ｈｉｌｂｅｒｔ算子以及反投影算子复合而成，其中Ｈｉｌｂｅｒｔ算子为无限区域上的奇异积分算子，１９第２章ＣＴ成像基本原理与图像重建算法难以直接计算，导数算子对数据噪声有放大作用（ＬａＲｉｖｉｅｒｅＰＪａｎｄＰａｎＸ，２０００），反投影算子会引起伪影。此外实际层析成像计算中采集的数据，有时还是不完全的。因此Ｒａｄｏｎ逆变换不便用直接数值求解，而采用离散重建算法（ＤｅｍｐｓｔｅｒＡＰ，１９７７）。ＣＴ投影数据实际上是沿一组平行束或扇形束进行采集的，所能得到的只是有限个离散角度，有限位置上的投影数据集（本论文都用平行束数据）。由Ｒａｄｏｎ变换知道可以通过有限个投影值【对】０，口），计算，的估计值或者说是重建（ＳａｕｅｒＫａｎｄＢｏｕｍａｎＣ，１９９３）。可以通过两类方法进行重建，一类是变换方法，另一类是级数展开法。上一节中介绍的滤波反投影重建属于前者，本节将要介绍迭代发属于展开法。迭代法优于滤波反投影算法之处就是减少了金属伪影，更好的处理了断层投影、有限角度断层成像，并提高了噪声抑制。ｌ生胄匕（ＬｅａｈｙＲＭａｎｄＢｙｒｎｅＣＬ，２０００）。需要说明的是对ＣＴ迭代法的研究还处于初期阶段，其存在很大的潜在应用可能。一直阻碍着迭代法应用于ＣＴ的原因，就是计算的复杂性。不过随着计算机硬件性能和软件效率的提高，使迭代法应用于ＣＴ成为了可能。本章将介绍四种最具代表性的迭代重建方法：代数迭代重建算法（ＡｌｇｅｂｒａｉｃＲｅｃｏｎｓｔｒｕｃｔｉｏｎＲｅｃｏｎｓｔｒｕｃｔｉｏｎＴｅｃｈｎｉｑｕｅｓ，Ａｍ３、联合代数迭代重建算法（ＳｉｍｕｌｔａｎｅｏｕｓＡｌｇｅｂｒａｉｃＴｅｃｈｎｉｑｕｅｓ，ＳＡＲＴ）、期望最大值迭代算法（ＥｘｐｅｃｔａｔｉｏｎＳｕｂｓｅｔｓＥｘｐｅｃｔａｔｉｏｎＭａｘｉｍｉｚａｔｉｏｎ，ＥＭ）和有序子集期望最大值迭代法（ＯｒｄｅｒｅｄＭａｘｉｍｉｚａｔｉｏｎ，ＯＳＥＭ）。下面对迭代法的一些基本概念进行介绍。迭代重建算法的概念与滤波反投影重建最大的区别在于前者把一幅连续的图像离散化，把一幅图像划分为Ｊ＝ｎｘｎ个像素点，并以夕＠，ｙ）＝罗ｘｆｂｆＯ，ｙ），，Ｏ，ｙ）表示，＠，ｙ）的离散值，且在每个像素内部ｉ（ｒ，臼）为常数。现假定一组基图像饿，也，…，６，】－，它的线性组合能够逼近我们重建的任何图像并定义６，ｏ，ｙ）一｛三如果伍＇ｙ’妻君个像素内则重建图像的离散值可以表示为（２．２６）／（ｘ，ｙ）＝罗ｘ—ｂＯ，ｙ）耳（２．２７）其中ｘ，表示第ｊ个像素点的像素值（或为灰度，或为密度），式（２．２２）通常简记为，ｏ，ｙ）＝∑ｘ，％被向量ｆ：１，２，３，…，ｊ＝１，２，３，…Ｊ２ｎ（２．２８）对给定者组瞽本胃德柳，如，…，６Ｊ，如果它是线性无关的，那么数字化图像，就唯一确定，这一列向量称为图像向量。于是估计图像，的第２章ＣＴ成像基本原理与图像重建算法问题就变成了寻找一个图像向量ｘ－－－｛ｘ，，ｘ：，…，工，），使其尽可能的接近，。定义Ｒａｄｏｎ变换，可得倪，厂；巩；歹＝∑ｘ，既ｂ，（２．２９）ｂ是用户所定义的函数，通常矾；ｂ，能用解析方法进行计算，一般为了简化运算把射线的宽度定义为零，假设每个像素的宽度为ｆ，如图２．１２所示射线穿过像素ｊ的长度为ｄ，则该像素对射线的投影值贡献为＿一∑勺ｐｆ。令投影数据相量为Ｐ＝｛ｐｌ，Ｐ：，…，Ｐ，｝，因此投影向量的第ｊ个分量的值另：ＰＪ一∑勺ｘ，（２．３０）式中ｉ－１，２，３，…，，ｊ『一１’２，３，…Ｊ其中Ｉ为总射线数，Ｊ为总像素点数。我们可以将式（２．２５）表示为矩阵形式Ｐ—ＲＸ（２．３１）其中Ｐ为Ｉ】【ｌ阶投影矩阵，ｘ为待求的Ｊｘｌ阶图像矩阵，Ｒ为ＩｘＪ阶系数矩阵。＼、＼＼ｌ‘＼一。。．．．Ｉ≮ｊ＼．，．＼图２．１２迭代重建法加权系数的定义现在我们的任务就是根据传感器测量得到的投影数据和已知的系数矩阵Ｒ，只要规定了像素的布置和几何结构，Ｒ就可以求出来。求解图像向量，可以通过两种方法求解图像向量：方法之一是求解系数矩阵Ｒ的逆尺～，从而得到ｘ—Ｒ‘１Ｐ（２．３２）这种方法很直观，但是一般系数矩阵很大，即使使用稀疏矩阵求其逆的运算量也是非常巨大的。第二种方法是通过迭代法把经过ｊ号像素的所有射线和累加，得到ｊ号像素值，公式如下：２１第２章Ｃ－ｑ＂成像基本原理与图像重建算法ｐｉｘｆ＝＞：‘ｆ符（２．３３）该公式直接迭代重建的图像伪影很严重，这显然不是我们所希望看到的，但是它为我们深刻理解迭代重建算法很有帮助。在实际情况下传感器的测量误差是不可避免的，噪声影响也是不可避免的。所以等式（２．２６）总是使方程无解。这时必须考虑误差并估计一组解，使它在某个最优准则下达到最佳（ＤｅＭａｎＢｅｔａｌ，２００１）。于是（２．２６）改为：（２．３４）Ｐ＝ＲＸ＋Ｐ这里ｅ是误差矢量，实际使用迭代算法估计图像矢量ｘ时，总是尽量逼近真实值，使某种最优准则下的误差最ｄｘ（ＢｙｎｅＣＬ，１９９６）。迭代算法的主要缺点就是计算量大，随着计算机技术的快速发展，迭代法将越来越受重视。迭代算法根据用观测数据更新中间图像的方法不同来分类，一般的标准是最ｄ＇－－乘函数或者Ｉ．发散，这分别等于在高斯和泊松分布下的最大似然。ＡＲＴ简单、有效，因此成为ＣＴ中的基本重建算法。在最近十几年中又出现了几种基于ＡＲＴ和ＥＭ的改进算法，其中一个比较有意义的是ＯＳＥＭ算法，本章将会介绍这些算法的原理。２．４．２．１ＡＲＴ和ＳＡＲＴ迭代重建算法ＡＲＴ代数重建算法是迭代算法的一种型式（ＧｏｒｄｏｎＲａｌ，１９７０），是最早用在ＣＴ上的迭代算法。ＡＲＴ的基本过程是给定重建区域的一个图像初值，再将所得ｅｔ的投影值沿其射线方向均匀的反投射回去，不断的对图像进行修正，直到满足所需要求，然后结束迭代过程。其迭代公式为：ｘ，（ｎ＋１）－－－－Ｘ抄九寿（Ａ坤∽）式中ｉ＝ｎｍ。ｄ（Ｉ）＋１，九为松弛因子，ＩＩｇ，１２（２．３５）２善芎是Ｒ的第ｉ行的欧几里得范数。在ＡＲＴ方法中，每条射线都要对Ｘ；的值修正一次，也就是说，第ｉ条射线对各ｘ；修正完后，再用第（ｉ＋１）条射线进行修正，所有射线修正完之后，即完成一轮迭代。这是通过最优准则判别是否达到收敛要求，如果否，以上次迭代后的图像矩阵ｘ为初值，则进行下一轮迭代。由于ＡＲＴ每次修正只考虑一条射线，如果某个探测器工作不正常，其采样的数据将对重建结果有很大的影响，而且其重建结果还与射线的迭代顺序有关（ＧｕａｎＨＱａｎｄＲｉｃｈａｒｄＧ，１９９４）。ＳＡＲＴＪ，失＿代算（ＡｍｄｅｒｓｅｎａｎｄＫａｋＡＣ，１９８４）法克服第２章ＣＴ成像基本原理与图像重建算法这些缺点，它是利用通过一个像素点的所有射线的修正值来确定该像素的平均修正值。这样取平均值就可以抑制干扰，而且其迭代结果与射线的迭代顺序无关，其迭代公式为：ｘ卜ｚ；町＋九寿荟景慨却ｎ’）（２．３６）式中Ｒ．，表示通过像素点ｊ的所有投影系数和。Ｒ，＋表示射线ｆ通过的所有的像素点的投影系数之和。其迭代步骤与刽ＲＴ一样。对于ＡＲＴ和ＳＡＲＴ算法一般选取收敛准则都是为Ｉｘ（肘ｎ—ｘ（“’０＜亭２．４．２．２ＥＭ迭代重建算法最大似然估计是最为常用和有效的估计方法之一，最早该迭代算法被用于ＰＥＴ（ｉＥ甘巳子断层成像）和ＳＰＥＣＴ（单光子发射∞，近年来它们在ＣＴ应用上７１起了越来越多的注意。最大似然估计的基本是思想是：利用已知的若干观测值，对没有任何先验知识的参数进行估计（Ｂａｒｒｅｔｔ的基大似然估计可表示为：ｅｔａｌ，１９９４）。因此，在使用最大似然估计方法时，被估计参数常被假定为常数，但是未知，而观测值被认为是随机变量，则未知参数ｘ露，＝ａｒｇｍａｘＰ（Ｙｌ，ｙ２，…，ＹⅣＩｚ）一ａｒｇｍａｘＰ（ｙｌｚ）知参数ｘ决定（李群等，２００５），记为Ｐ（ＹＩｚ）。（２．３７）式中Ｙｌ，），：，…，ＹⅣ表示随即变量Ｙ的Ｎ个观测值，随机变量Ｙ的概率密度函数由未图像重建的最大似然估计模型和ＥＭ算法是ｆｌ了ＳｈｅｐｐＬＡ（１９８２）和ＬａｎｇｅＫ（１９８４）在８０年代独立发展起来的。该模型假定ｘ射线的光子数是服从Ｐｏｓｓｉｏｎ分布，将图像分为Ｊ个像素，像素ｊ的估计参数为Ｘ，。最表示沿着第ｉ条射线所探测到光子数。，如果口｛ｆ（０ｓ％ｓ１）表示第ｊ个探测器探测到的第ｉ条射线的光子数，且满足等式∑％；１，那么只是Ｐｏｓｓｉｏｎ分布的一个样本（ＤｅＭａｎＢｅｔａｌ，２０００）。∑口玎ｚ，ｔ只基于Ｐｏｉｓｓｏｎ分布的似然函数为（２．３８）ｍ）＝驴等＾；∑口｛『ｚＪ（２．３９）其中凡是第ｉ条射线投影值的预测参数值。此外投影Ｐｏｉｓｓｏｎ分布模型满足：（２．４０）第２章ＣＴ成像基本原理与图像重建算法把式（２．３５）代入式（２．３４）有：砌啼扣－掣（２．３６）最大，（２．３６）的对数似然函数ｌ（ｘ）＝１０９（“ｘ））如下：仁４１，ＣＴ图像重建的极大似然方法就是求出图像矢量ｘ＝怯，Ｘ２，…，ｘ，＞，使似然函数７＠）２荟ｚ，善口。＋著见１。ｇ（善口｛『ｚ，）一善１。ｇ（Ｐｉ！）（２·４２）最大似然估计的具体实现过程中，基本是以对数似然函数作为目标函数，如式（２．３７），连续凹函数，０）存在全局极大值（ＯａｎＪＫ，２００１）。因此求ｘ的极大似然解可以通过求对数似然函数的偏导数，由于翌盟是非线性函数，方程里盟。ｏ不存在解析解，只能通过迭代的方法就去近似解。ＳｈｅｐｐＬＡ和Ｌａｎｇｅ各自提出下面的迭代公式：善白箭乏白ｘＪ＠）掷灿爹妻茹仁４３，式中∑勺为经过像素点ｊ的所有射线在ｊ点的加权系数之和，乏白ｘ，＠）为射线ｆ通过的所有像素点的像素值，与该点对射线ｉ的权值的乘积之和。该迭代公式有两个明显特点（张书明等，２００７）：（１）每次迭代之后，对数似然函数ＺＯ）都会增大，即有ＺＯ‘”１’）＞Ｚ＠‘“’）。（２）给定非负的初始估计ｘ（们，以后迭代的ｘ（１），工（劲，…，ｘ（”）也是非负。ＥＭ算法一般选取收敛准则也是为Ｉｘ（”ｎ—ｘ（“’ＩＩ＜亭。２．４．２．３ＯＳＥＭ迭代重建算法ＥＭ重建算法具有很好的抗噪声能力，尤其是在处理统计性差的数据时，更能显示出它相对于解析法的优越性，但是这种方法在每一次迭代过程中，使用所有的投影数据对重建图像每一个像素点的值进行校正，重建图像只被替换一次，所以迭代法的运算量大，运算时间长，难以满足快速重建要求。ＯＳＥＭ算法（ＫａｍｐｈｕｉｓＣＫａｎｄＢｅｃｋｍａｎＦＪ，１９９８）是在ＥＭ算法的基础上引入有序子集来加速收敛的一种有效的实用算法。ＯＳＥＭ在迭代过程中将投影数据分成：Ｔ个子集，每一个子集对重建图像各像素点值校正以后，重建图像便被更新一次，所有的子集运算一遍，称为一次迭代过程。在ＥＭ算法一次迭代过程中，重建图像被更新一第２章ＣＴ成像基本原理与图像重建算法次，而在ＯＳＥＭ算法的一次迭代过程中重建图像被更新Ｔ次，所以ＯＳＥＭ算法可以加快收敛的速度。在０ＳＥＭ中投影数据（ｎ，Ｐ２，…，所）被分为Ｔ个有顺序的子集（墨，Ｓ：，…Ｓｒ），子集的划分很重要的，不能随意的划分，子集太多时就会出现不收敛的情况，太少的话又不能快速收敛，因此我们需要合理的划分子集和选取松弛系数（刘力等，２００２，吴朝霞等，２００４，李健等，２００７）。一般来说划分子集时，先要将投影数据按某个规律划分，这些子集必须有一定的对称性，即要以一种平衡的方法选择子集（印胤等，２００３），选取子集的一般规律如下：（１）包含在每个子集内的的重建图像信息和统计信息应该充分。（２）子集的选择应该是平衡的，每个子集对投影数据的贡献应该是相等的。（３）连续子集将的断层图像距离应该最大。（４）子集的个数应该是在收敛范围内的最小值（张书文等，２００３）。子集划分后，ＯＳＥＭ的迭代公式为：工，＠＋１）。型∑芦（２．㈣三勺怠荟勺ｚＪ＠）式中薹～为第ｉ个子集中经过像素点ｊ的所有射线在ｊ点的加权系数和，Ｚ勺Ｘｉ＠）为射线ｆ通过的所有像素点的像素值，与该点对射线ｉ的权值的乘积之和。２．４．２ＣＴ重建算法的一些补充说明迄今为止我们的讨论几乎都是局限于平行束扫描的情况，可是现在的ＣＴ大多数都是扇形束扫描，为了解决重建问题，需要进一步观测扇形束的几何关系。扇形束ＣＴ又包括两种设计：等角扇形束和等距扇形束。其区别主要是由不同的探测器设计导致的，见图２．１３。每次采样之间的角度间隔相等蕊≤＼＼ｊ１．ｊ汐一≮潮（ａ）等角扇形束图２．１３扇形束扫描图解２５蕊＼＼■ｊｒ、采样角度．ｊ问距不等距（ｂ）等距扇形束第２章ＣＴ成像基本原理与图像重建算法由于射束不再平行，故其重建算法需要重新讨论。一种方法是重排算法，即把一个视图中采集到的扇形数据重新组合成平行数投影数据，在用上面介绍的重建算法进行重建。这种方法对扇形束的投影的视角与各扇形束间的增角有一定的约束条件。第二种方法是扇形束投影直接重建算法，这种算法不必把数据重排，只需要适当的加权即可运用于上一章类似的算法重建图像。最后一种方法是直接使用迭代法，与平行束的区别仅在于系数矩阵的计算略有不同。第３章金属伪影成因及其校正方法研究第３章金属伪影成因及其校正方法研究ＣＴ的应用越来越广泛，随之而来的是一系列的新问题，重建图像中的金属伪影就是其中之一。金属伪影在图像中表现为明暗相间的条状伪影。这些伪影降低图像质量，严重影响图像的后续处理和分析。如何消除这些条状伪影以提高图像质量，是当前ＣＴ应用研究领域中的热点。３．１金属伪影的成因金属伪影的成因相当复杂，随着金属物体的形状和密度不同，伪影的形状可能略有变化（谷建伟等，２００４），但其主要原因是射束硬化。射束硬化是由多色的Ｘ射线能谱与能量相关的衰减系数造成的，低能的Ｘ射线比高能的ｘ射线更容易被大多数材料吸收（Ｓｅｉｔｅ定律可以按下面的形式重写：Ｐｅｔａｌ，１９８５），这主要是光电吸收的原理。为了强调衰减的能量相关特性，对于单能量ｘ射线，Ｌａｍｂｅｒｔ．ＢｅｅｒｓＩｏ正。Ｉｌ’Ｅ泸以ｂｑ其中ＬＪ和易Ｊ表示能量为Ｅ的火射和透射ｘ射线强度，如Ｊ是物体在能量Ｅ下的衰减系数。一个给定物体的衰减系数定义为透射与入射强度之比的负对数。实际上，ＣＴ实际产生的ｘ射线覆盖了一个很宽的频谱。图３．１显示了一个工作在１２０ｋＶｐ的Ｘ光管输出ｘ射线的能谱的例子。能谱显示工作在１２０ｋＶｐ的射线管输出，Ｘ光子的能量从２０ｋｅＶ到１２０ｋｅＶ。归化输出Ｘ射线能量（ｋｅＶ）图３．１Ｘ射线光管的输出能谱２７第３章金属伪影成因及其校正方法研究前面提到过，大多数材料的线性衰减系数是和能量相关的。一般的材料对于低能光子具有很高的吸收系数，对于高能光子具有较低的吸收系数（ＨｅｒｍａｎＧＴ，１９７９）。例如ＣＴ应用于医疗时，当具有多能谱的ｘ射线射向一软组织时，能量较低的光子比能量较高的光子被吸收得更多。能量在１０ｋｅＶ以下的光子基本全被皮肤吸收，由于低能光子被大量吸收，导致高能光子的比率变大，这就是人们称这种现象为“射束硬化＂的原因（张泽宏，２００５）。图３．２显示了骨头和水的线性衰减系数和Ｘ射线能量间的关系。，Ｌ一一一姗ｍｍ。衰减系数１／ｅ１ｌｍ¨ＯＺＵ４０６０８ＵｌＵＵ１２０１４０Ｘ射线能量（ｋｅＹ）图３．２骨和水的线性衰减系数和Ｘ射线能谱的关系图通过上面射束硬化的介绍，将有助于我们理解金属伪影的成因。实际ＣＴ应用中，为了减小射束硬化的影响，尽量使ＣＴ的ｘ射线逼近单一能量（ＤｅＭａｎＢｅｔａ１．１９９９），但是当ｘ射线在透射物体的过程中总是会由于衰减呈现多能量，特别是射线通过金属物质这种高吸收系数的物质后，能量大大降低，这时射束硬化效果会越来越明显，最后的透射强度位于数据电子设备动态范围的底部，导致投影数出现跳跃变化。由ＣＴ成像原理可知，如果断层射线扫描系数的变化出现非连续性的跳跃变化，如被扫描物体中包含金属物体，其衰减系数远远大于其它物体的衰减系数，就会导致投影数据的一阶导数在某一段表现为弱连续性，经过滤波处理后，这种弱连续性迸一步扩大，就在图像中形成明暗相间的“条状”伪影（ＰａｔｒｉｃｋＪ，２００５），即我们常说的金属伪影。３．２金属伪影校正算法随着凹技术的发展，人们提出许多方法来校正金属伪影。这些方法基本上可以分为三类：投影插值方法、迭代重建方法和混合法。投影插值方法理论简单、第３章金属伪影成因及其校正方法研究计算快速和容易算法实现等优点而具有较高的应用价值，但是它只能处理简单的金属物体。迭代法能有效去除金属伪影和抑制噪声，而且能很好的呈现金属物体的结构，但其运算量非常大，速度很慢，难以实用化。但是随着计算机硬件的发展和人们对算法的不断优化，使得迭代法用于实际的校正方案的可行性越来越高。混合法一般是将插值法和迭代法混合使用，该方法吸收了插值法和迭代法的优点，具有算法速度较快，能有效反映金属物体结构，但是金属周围的区域失真，使其在医用ＣＴ应用中存在很大的局限性。下面具体介绍这三种常用的去金属伪影算法。３．２．１插值重建法插值重建算法是工业ＣＴ中最常用的去金属伪影算法，该算法理论简单、计算快速和容易算法实现。插值算法有很多种，如谷建伟等提出了快速线性插值法和基于伪影权重的非线性插值法，林宙辰等采用多项式插值，ＬｅｖｉｔｔＲＭ和ＢａｔｅｓＲＨ（１９８３）提出了用一种特殊的函数进行插值的方法，后来又提出ＴＣｈｅｂｙｓｈｅｖ多项式的差值法。２０００年ＺＨＡＯＳｈｉｙｉｎｇ等人还提出了对投影数据的小波系数线性插值的方法。各种插值算法的步骤基本都是一样的（谷建伟等，２００６ａ），唯一区别就是通过何种插值方法对金属投影数据区域进行插值。３．２．１．１插值算法步骤插值算法一般分为如下几步：（１）确定金属区域通过插值重建算法先对原始数据通过ＦＢＰ直接重建，再通过选取的阈值从原始图像中分标记出金属区域Ｑ４。由第二章的介绍可知从探测器获取投影数据Ｉ和１０，得到，－ｔｏｅ—Ｊ肚，投影数据Ｐ生成原理的公式为：Ｐ＝ｌｏｇＩｏ—ｌｏｇＩ＝阻出（３．２）对投影数据直接使用滤波反投影重建，得到一幅含金属伪影的ＣＴ图像，Ｑ＠，Ｙ）。选取阈值ｓ（一般Ｓ的选取需要作一些试验，才能使分割的金属区域最佳），分割金属区域：弘ｙ）伴力’笔∥净ｓ式中‘＠，Ｙ）表示金属区域的像素值。（３．３）第３章金属伪影成冈及其校正方法研究（２）对图像厂Ｑ。＠，ｙ）使用旋转坐标系的方法进行重投影重投影是插值算法中很重要的步骤，即对通过阈值分割出来的金属区域Ｑ“进行再投影，确定金属投影数据的区域，这样才能对该范围内的投影数据进行插值。在计算金属区域的投影坐标时需要知道，旋转射线对金属区域进行投影和射线固定不变旋转金属区域是等价［孔Ｊ（ＪｏｓｅｐｈＰＭ，１９８２）。我们首先记下每个金属像素点的坐标（Ｍｘ，Ｍｙ），用如下公式计算旋转后的金属区域的横坐标：Ｔｘ＝Ｍ。ｃｏｓ（Ｃ）－Ｍ，ｓｉｎ（Ｓ）（３）对金属区域投影数据进行插值校正这里主要是针对平行束投影数据，如果是扇形投影数据，我们假设已经进行（３．４）式中妒为旋转角度。通过对这些横坐标进行标记计算出金属投影数据区域｛ｇ）。过重排。插值校正可以使用很多插值法如：线性插值、多项式插值以及三次样条插值等等。线性插值可以消除由结构简单的金属物体引起的绝大多数金属伪影，并且不会引入新的伪影，使图像整体质量得到提高。样条插值校正随着噪声的提高，会在重建图像中引入彗星状伪影，使图像模糊。多项式插值的效果介于线性插值和样条插值之间。综合考虑去金属伪影的效果和对噪声抑制等因素，在实际的ＣＴ图像金属伪影校正所使用的这几种插值方法中，线性插值是一种最为简单、最常用而且校正效果较好的金属伪影的校正方法。下面将如何使用线性插值。上一步中通过旋转坐标系对金属区域进行重投影确定的投影数据区域｛酚，｛酚由金属投影ｇＩｎ和相应扫描角度下的非金属投影ｇｎ两部分组成，投影区域｛ｇ）线性插值后的结果用ｇ，表示，理想中的结果是ｇ，一ｇ。，假定投影角度妒下的金属投影数据范ｍ为（ｉｘ，ｉ２），假定ｉＥ［ｉｌ，ｉ２］，则插值公式为：ｇ，Ｏ）一÷互—÷ｇ（‘）＋÷—寻ｇ（ｉ２）０．５）‘２一王ｌ‘２一ｈ则插值后投影角度驴下的投影数据｛ｇ’）表示为：ｇ’ｏ，＝ｆ；善了；主：耋：乏；（４）重建图像并恢复图像金属信息ｃ３．６，对校正后的投影数据｛ｇ’】．直接使用ＦＢＰ莺＿建得到图像尼０，Ｙ），并恢复图像部分金属信息。在图像尼Ｏ，Ｙ）中，原来的金属区域如图３．４（ｃ）所示变成了一个暗斑，这是由于线性插值引起的，这样虽然基本消除了金属伪影，但是却丢失了金属信息。所以需要对图像金属信息进行恢复，只需要将原图像，Ｑ＠，ｙ）中的金属第３章金属伪影成因及其校正方法研究区域取代尼＠，ｙ）金属区域即可，’。，ｙ，２｛童￡富ｏ，ｙ其）它∈ｆｌ”整个算法的流程图如图３．３所示。原始投影数据ＦＢＰ重建原始图像（３．７）分割出金属区域图像标记每个金属区的重投影获取金属投影像素点对金属投影数据进行线性插值计算ＦＢＰ重建处理后图像对金属区域进行补偿图３．３算法流程图３．２．１．２插值算法仿真实验结果及分析为了考察线性插值算法的效果，建立仿真实验模型参数见表３．１，１到４组数据代表组织，第５组参数代表一个金属区域。实验采用平行光束解析法模拟投影数据，投影角度为ｏ．１８０度，步长为２度，共９０＋投影度。探测器数目为９０个，重建图像大小为１８０Ｘ１８０。表３．１数值仿真模型参数在仿真实验中除了加入高斯噪声。通过函数Ｈ（ｘ）来模拟射束硬化（ＺｈａｏＳＹ，２０００），使射线载经过金属物体后被进一步吸收，更好的模拟金属伪影：３１第３章金属伪影成冈及其校正方法研究Ｈ（ｘ）－ｘ－塑丝笔丕曼匝咿㈣）果如图３４所示。（３８）式中ａ为硬化闽值，ｃ为硬化幅度，仿真实验ｑ诹ａ＝８，ｃ＝０．５。此外电脑配置为单核（主频１＿４劬ＨＺ），５１２眺Ｊ存（所有仿真试验均在同一机器上完成）。仿真实验的结（ｃ）线性插值后ＦＢＰ重建图３４线性插值仿真结果其中图３．４（曲是需要扫描的理想模型原图，白色的是金属物体；图３．４Ｃｏ）给出了直接ＦＢＰ重建后的结果，由图３．４可知其金属伪影明显，图像质量非常差，而且金属结构不清晰；图３．４化１是图像对金属投影数据区域进行线性插值校正后用ＦＢＰ重建后的结果，其基本消除了金属伪影，当是缺失了所有的金属信息，而且原金属区域变为一个暗斑，这是由线性插值的不光滑性引起的。国３．４（ｄ）是图像３．４（ｃ）经过金属区域补偿后得到的结果，虽然其恢复了太部分的金属信息，但是边缘仍存在一些失真，这主要是由于判别金属区域的阈值选取没有达到最优，导致部风金属像素点或是非金属像素点被误判。不过总体来说．对于结构简单的金属引起的伪影，线性插值重建法的效果还是比较令人满意的。ｉ万一ｆｄｌ插值重建后金属补偿第３章金属伪影成因及其校正方法研究３．２．２迭代重建法由于投影数据本来就易受噪声的影响，再加上金属这种高吸收系数物质，使采集的投影数据出现了偏差，因此直接通过解析法重建总是不能很好的消除金属伪影。在第二章中我详细了几种迭代重建算法，迭代法优于滤波反投影算法之处只要在投影数据充分并且迭代次数足够多的情况下就能有效减少了金属伪影，更好的处理了断层投影、有限角度断层成像（ＦｅｓｓｌｅｒＪＡａｎｄＨｅｒｏＡＯ，１９９５），并对各种噪声有良好的抑制性能。虽说的ＣＴ迭代法的研究还处于初期阶段，但其存在很大的潜在应用可能，而且一直阻碍着迭代法应用于ＣＴ的计算量大的原因，随着计算机硬件性能和软件效率的提高，也得到了一定程度的解决。因此近年来迭代法越来越受人们重视，特别是在去金属伪影这个问题上的优越性，使它更受学者和工程师的亲睐。３．２．２．１系数矩阵的定义和计算迭代法已经被广泛的研究了很多年，这些算法的实施方法如下。为了简化讨论，让我们把讨论限制到一个二维物体以及它的投影测量上。由２．２．２节可知迭代法都是通过等式Ｐ；ＲＸ＋ｅ（ｅ是误差矢量，实际使用迭代算法估计图像矢量Ｘ）来估计图像像素，使其尽量逼近真实值。由等式可知，要使用迭代必须先求出投影矩阵Ｒ，该系数只与像素的布置和几何结构有关。我们使用的是平行束射线，且如果是扇形扫描认为已经对射线进行过重排。要求投影矩阵Ｒ必须先定义射线和定义投影权重系数％。第一种表示射线和％的方法是把射线看作宽为ｆ的粗线（常数了等于像素宽度Ｐ），这样射线经过某个像素点时会覆盖像素点的部分面积ｓ∥此时％定义为：％ｔ￥ｉｔ／Ｐ。（３．９）第二种表示射线和％的方法是把射线宽度认为是为零，当射线经过像素点时，其与像素点相交的长度为厶，此时％定义为：ｒｏ；２：ｆ『／（４２ｐ）第三种表示方法如Ｆ白（３．１０）Ｌ墨ｉ号射线经粥号像素其它Ｑ（３．１１）不难看出，射线宽度ｆ和％的定义会影响投影矩阵Ｒ的值，从而影响到重建后图像的精度，一般情况我们选择第二种定义。迭代法的系数矩阵都是提前计算好，并存放于内存中，在迭代时直接从内存中读取。３３第３章金属伪影成因及其校正方法研究３．２．２．２迭代算法实现时的一些优化方法使用迭代算法去金属伪影，其实就是通过迭代法重建图像，步骤就是先计算系数矩阵，并根据需要以相应的形式进行存储，并读入内存，再是根据迭代公式进行迭代计算。迭代法实现相对简单，但是计算量很大、运算时间长。我们可以在实现算法时尽可能对算法进行优化，达到提速目的。在对几种迭代法进行仿真实验时，都根据迭代法不同的特点，在实现上作了一些优化。可以通过一下的方法提高算法运算速度：（１）不要直接进行矢量运算由迭代公式（２．３５）可知，ＡＲＴ的迭代过程如下，即第ｉ条射线都要对所有的ｘ，修正一次，第ｉ条射线对各ｘ，修正完后，再用第（ｉ＋１）条射线进行修正，所有射线修正完之后，即完成一轮迭代，即每次需要进行矢量计算Ｒ工（…。Ｒｚ（“）表示第ｉ条射线的估计投影值，Ｒ是第ｉ条射线的系数矩阵（Ｊｘｌ），当射线ｉ没有经过像素点ｊ时Ｒｉ等于０，由此可知Ｒ是一个稀疏矩阵，当多数值为零，只有射线ｉ经过的像素点所对应的系数不为零，而每条射线经过的像素点最多２Ｎ个（图像分辨率为ＮｘＮＯＰＪ＝ＮｘＮ），为了避免最少Ｎ术（Ｎ．２）次无意义的计算，我们需要记下每条射线当然这样会使程序便得复杂点，因为需要跟据不同的迭代公式计算其它的附经过的点，这将大大减小运算量，而且随着Ｎ的增加效果会更明显。加信息，如因为ＡＲＴ的迭代过程是一根射线一根的逐线校正，所以对于ＡＲＴ迭代就需要计算二维矩阵Ｌｉｎｅ［ｉ］［ｊｌｉｌ己录第ｉ条直线经过的第ｊ个点的系数权值。而ＳＡＲＴ、ＥＭ和ＯＳＥＭ迭代法都是利用通过该像素点内的所有射线的估计投影值Ｒｚ（“’来确定一个像素的平均修正值，因此还需要在计算并记录下经过每个像素点的所有射线序号，这些数据只要像素的布置和几何结构确定，则都是不变的，可以和Ｒ矩阵一起计算并存入内存。（２）更具迭代法的修正顺序，对射线估计值足ｚｏ）的计算进行调整估计投影值足ｚ（“）计算是很有讲究的。ＡＲＴ迭代是需要把第ｉ条射线校正后的像素点更新，才能进行第（ｉ＋１）条射线的校正，因ｌｋ蚣ｄｌＴ迭代法需要在每次校正时动态计算Ｒ工ｏ’（鼍动态更新的），其余的迭代法可以在每轮迭代开始时计算好每条射线的估计投影值，以后每次校正直接使用结果即可。这大大提高了程序的运算效率和速度，使ＳＡＲＴ、ＥＭ和ＯＳＥＭ每次迭代速度远快于ＡＲＴ迭代。（３）通过迭代初值得选取提高收敛速度迭代初值对收敛速度有一定程度影响，如果初值越接近真实值则收敛越快，除了ＡＲＴ迭代算法，其余算法初值取投影数据平均值或是直接ＦＢＰ重建后像素值都比初值取零收敛更快。第３章金属伪影成因及其校正方法研究３．２．２．３仿真实验结果及分析为了考察迭代算法的去金属伪影效果，建立仿真实验模型参数见表３．２，１Ｎ４组数据代表组织，第５组参数代表一个金属区域。实验采用平行光束解析法模拟投影数据（如果是扇形束扫描，迭代算法不变，只是系数计算略有不同），投影角度为０－１８０度，步长为３度，共６０个投影度。探测器数目为１２８个，重建图像大小为１２８ｘ１２８．此外通过公式（３．８）模拟射束硬化。表３．２数值仿真模型参数迭代法去金属伪影仿真实验结果如图３．５所示０）理想模犁（ｂ）ＦＢＰ虱Ｉ（０ＡＲ幢建ｆｄ）ＳＡＲＴ重建（ｅ）ＥＭ重建（００ＳＥＭ重建图３．５选代法去金属伪影仿真实验结果图３５（ａ１是该仿真试验的理想模型图，下方白色的亮点是一金属物体。图３．５Ｃｏ）给出了直接ＦＢＰ重建后的结果，不仅有明显的会属伪影，中间还有由于射束硬化引起的暗带。图３５（ｃ）是通过Ａ】嘴代１０次后的重建结果，由图可知ＡＲＴ算法大大削弱了金属伪影，但是由金属物体引起的射束硬化喑带还是很明显，而且金属物体周围出现了环状伪影，这主要是由于ＡＩ盯迭代是一根射线一根射线的逐线校正，使其易受噪声和一些随机因素干扰。ＡＩＩＴ算法的收敛速度还受迭代射第３章金属伪影成因及其校正方法研究线顺序的影响，不过由于ＡＲＴ收敛快，这点影响是可以忽略的。图３．５（ｄ）是ＳＡＲＴ算法迭代三轮的结果，其基本消除了金属伪影，很好的还原了金属信息，不过由于ＳＡＲＴ是使用通过该像素的所有射线的估计投影值的平均值来校正像素，这样虽然能有效抑制干扰，但是也会使非金属区域的一些细节变模糊（图像背景的暗纹是由于迭代次数过少，没有达到最佳收敛）。图３．５（ｅ）是通过ＥＭ迭代三次后的结果，由图可知ＥＭ算法能很好的消除金属伪影，还原金属信息。而且相同迭代次数下，ＥＭ算法的重建效果更好，说明了ＥＭ收敛Ｉ：匕ＳＡＲＴ要快。ＥＭ同样在迭代次数不充分的情况下为丢失一些非金属区域的细节。图３．５（ｄ）是通过ＯＳＥＭ迭代两次后的结果，在仿真试验中把投影数据分为三个有顺序的子集，每个子集隔３×△度取一投影射线（△位投影角度旋转步长，本实验中为３度），例如在投影角度驴下，射线集｛Ｏ，６，１２，…，８６）为第一个子集，射线集｛２，８，１４，…，８８｝为第二个子集。从重建结果中可知，ＯＳＥＭ基本消除了金属伪影和射束硬化引起的暗带，而且能基本还原金属信息，但是同样会模糊细节。而且由图我们可知ＯＳＥＭ虽然迭代两次，但是其数据更新了六次，所以其收敛速度较快。３．２．３ＥＭ混合重建法插值法虽然对简单金属物体引入的金属伪影效果明显，但是它不能清晰的呈现金属物体，而且对于对金属物体效果不好。迭代法虽然能有效去除金属伪影，但是要得到理想的重建结果，需要迭代很多次，运算量大，重建时间较长。芝加哥大学的ＸＩＡＤａｎ（２００５）结合了插值法和迭代法的优点，提出了ＥＭ局部迭代法。３．２．３．１ＥＭ混合法主要步骤ＣＴ成像包括两个步骤：数据采集和图像重建。由于金属的吸收系数一般远远大于非金属区域，这样每个投影角度下，总有一段投影数据实际没有检测到任何射线，导致这部分投影数据过大。通过ＦＢＰ重建图像时，就会出现金属伪影，其实金属伪影主要是投影数据不连续的反应映。混合法主要分两步，首先和插值重建法一样，首先确定投影数据中的金属投影数据，然后分别重建两部分数据，非金属区域直接ＦＢＰ重建，金属图像区域通过ＥＭ迭代法重建。具体说明如下：（１）确定金属和非金属投影数据既然金属伪影主要是投影数据不连续的反应映，我们可以先确定每个投影角度下的金属投影数据。确定金属投影数据的方法和插值重建法一样，首先通过ＦＢＰ算法重建原始投影数据，选取适当的灰度级阈值确定金属图像区域。然后利第３章金属伪影成因及其校正方法研究用旋转坐标系确定每个投影角度下的金属投影数据范围。ｆ２１重建非会属区域图像金属投影数据区域确定后，先将这部分投影数据移除，即使是金属投影数据也认为由金属部分和非金属部分（金属部分比重大），因此我们还需要估计非金属部分对投影值得贡献，线性插值重建是常用方案之一。然后通过ＦＢＰ重建非金属区域，因为去除会属投影数据并插值后的原始投影数据较为平滑，通过ＦＢＰ重建不仅效果好而且速度快（在商用ＣＴ中ＦＢＰ的部分算法可用硬件实现，使其重建速度更快）。但是重建图像中的金属区域为一暗斑（如图３．６（ｃ）所示），因此需要下一步对金属信息进行恢复。ｆ３）ＥＭ迭代重建金属区域图像该算法之所ｌ三ｃ被称为混合法就是因为使用两种方法重建图像。因为ＥＭ迭代法不易受数据不连续性的影响，所以我们使用ＥＭ迭代法重建金属物体，而且当金属物体结构较复杂时，迭代法能更好的呈现金属结构。虽然迭代法运算量巨大，但是由于一般ＣＴ图像中的金属物体都比较小，所以该算法的速度略大于插值重建法，远快于整幅图像使用ＥＭ迭代法重建。最后将通过ＥＭ迭代得到的属区域取代步骤（２）中重建的非金属区域图像的暗斑。３．２．３．２仿真实验结果及分析为了验证ＥＭ混合法的去金属伪影效果，建立仿真实验模型参见３２．２３节表２，第５组参数代表一个金属区域。实验采用平行光束解析法模拟投影数据（如果是扇形束扫描，迭代算法不变，只是系数计算略有不同ｋ投影角度为０－１８０度，步长为２度，共９０个投影度。探铡幡数日为９０个，重建图像大小为１８０ｘ１８０，通过公式（３．８１模拟射束硬化。图３．６是仿真实验结果。ｎ）理想模刑酗Ｏ）直接ＦＢＰ重建３７第３章金属伪影成固及其枝正方法研究（ｃ邱Ｐ重建１＃金属区域棚）ＥＭ混台洼重建结果图３．６ＥＭ混合法仿真实验结果图３．６（ａ１给出了该仿真试验的理想模型囤，下方白色的是表示金属物体。图３．６ｆｂｌ是ＦＢＰ重建原始投影数据后的结果，图中含有明显的金属伪影，严重影响了图像的清晰度，而且中间还有由于射束硬化引起的暗带。图３．６（ｃ）通过ＦＢＰ重建去除金属投影数据并插值后的投影数据，基本消除了金属伪影，但是由于线性插值的托连续性，导致金属区域周围有失真。图３．６（ｄ）为完整的ＥＭ混合法重建后的效果，捎除了金属伪影。但是金属周围略有失真，这主要是阈值的选取没有达到最佳和线性插值的不光滑性引起的。ＥＭ混合法图像已经很接近理想模型了，但是它最大的不足金属周围区域发生形变，特别是在金属物体较大时很明显，见图４３（ｃ）。３．３本章小结现在主流去金属伪影算法分为三类：插值重建算法、迭代重建算法和混合重建算法。通过仿真实验可知它们各有各优点和缺点。插值算法首先判断金属投影区域．然后对金属投影进行插值，晟后重建图像。该方法简单快速，但是由于这种算法是通过移出金属区域的投影数据，从而消除投影数据的跃变，去除盒属伪影，显然使用插值的数据代替金属覆盖区域的投影数据显然是不精确的，因此金属区域周围会发生形变，而且金属物体越大，失真越明显。因此插值算法只适用于金属结构简单的情况，而且金属周围的区域结构相对简单。迭代法也是一种常用的ＣＴ图像重建算法，ＥＭ迭代法和ＯＳＥＭ迭代法都是典型代表。迭代法能有效去除金属伪影和抑制噪声，而且能很好的呈现金属物体的结构，但是迭代法巨大的运算量，导致其运算速度很慢，难以实用化，而且选代法会模糊一些局部细节。ＥＭ混合法是近几年提出的一种新算法，该算法仅对金属区域进行ＥＭ迭代第３章金属伪影成因及其校正方法研究重建，对非金属区域直接ＦＢＰ重建，该算法速度较快，能有效反映金属物体结构，但是由于该算法也是通过移出金属区域的投影数据，从而消除投影数据的跃变，并通过线性插值估计金属覆盖区域的投影数据，因此和插值法一样，通过该算法重建的图像金属周围的区域失真，使其在医用ＣＴ应用中存在很大的局限性。通过实验发现，当前已经实用的金属去伪影算法主要存在以下几个缺点：（１）由于使用插值的数据代替金属覆盖区域的投影数据，导致的金属周围区域失真，而且金属物体越大、结构越复杂，失真越明显。（２）不能很好呈现多金属物的结构。０）阈值的选取是否最优，直接影响算法的处理结果。所以我们需要设计一种算法能保证金属周围的区域不失真的情况下，有效去除伪影，呈现金属物体的构造细节，而且该算法尽量降低阈值选取的精确度对算法去伪影效果的影响。３９第４章基于自适应衰减滤波的ＣＴ图像去金属伪影混合法第４章基于自适应衰减滤波的ＣＴ图像去金属伪影混合法上一章中介绍了现在主流的区金属伪影重建算法，插值法是一种常用的去金属伪影算法。无论是谷建伟等人提出了快速线性插值法，还是林宙辰等采用多项式插值都是首先判断金属投影区域，然后对金属投影进行插值，最后重建图像。这类方法简单快速，但是只适用于金属结构简单的情况，而且会使金属周围的区域失真。迭代法是一种常用的ＣＴ图像重建算法，如ＡＲＴ、ＳＡＲＴ、ＥＭ和ＯＳＥＭ迭代法，它们迭代算法都是典型代表。迭代法能有效去除金属伪影和抑制噪声，而且能很好的呈现金属物体的结构，但其运算量非常大，速度很慢，难以实用化。第三章最后实现了一种ＥＭ混合法，仅对金属区域进行ＥＭ迭代重建，对非金属区域直接ＦＢＰ重建，该算法速度较快，能有效反映金属物体结构，不足之处是金属周围的区域失真，使其在医用啦用中存在很大的局限性。本章提出一种去ＣＴ图像中金属伪影自适应衰减滤波混合法，该方法首先判断金属投影区域，对该局域进行自适应衰减和滤波，然后通过ＦＢＰ重建整幅图像，再利用原始投影数据对金属区域进行ＥＭ迭代重建，并对自适应衰减滤波重建后的金属区域进行校正。数值模拟的ＣＴ实验证明，该方法可以有效的去除金属伪影，并很好的保留了金属及其周围组织的信息，其效果优于插值法和ＥＭ混合法，而且减弱了阈值选取精确度对算法性能的影响，特别是在多金属物体的情况下，算法的效果更是优于其它算法。而且该算法计算复杂度小，有很高的实用价值。４．１算法的基本思想投影数据认为由两部分组成：一部分是非金属区域投影数据。另一部分是射线透射金属的投影数据，这部分投影数据一般被认为是“被破坏的”并应该对其进行校正，这样通过ＦＰＢ重建时金属伪影被大大削弱（Ｃｈｅｒｔ崛２００２）。插值法重建法和ＥＭ混合法都是移出金属区域的投影数据，简单的使用邻域插值的数据代替金属区域的组织或背景，但是这样并不能真实反映被金属覆盖区域的衰减系数（ＤａｖｉｓＧＲ．１９９４），因此通过这些算法重建的图像，会出现金属区域附近细节失真。本章的算法将保留金属投影数据的部分特征数据，这样重建后能真实反映金属区域附近细节。由第３章我们知道射线通过高吸收系数的金属物质后，能量大大降低，这时射束硬化效果会越来越明显，最后的透射强度位于数据电子设备动态范围的底部，导致投影数出现跳跃变化，由ＣＴ成像原理可知，这种跃变就会导致投影数第４章基于自适应衰减滤波的ＣＴ图像去金属伪影混合法据的一阶导数在某一段表现为弱连续性，经过滤波处理后，这种弱连续性进一步扩大（ＰａｎＸ，１９９９），就在图像中形成金属伪影。为了削弱金属投影数造成的跃变并保留这部分数据的特征，可以先对金属投影区域进行了自适应衰减（ＨｓｉｅｈＪ，１９９８），消除投影数据的跃变，再进行滤波，进一步改善金属投影区域的平滑性，而且还能降低阈值选取的精确度对算法去伪影效果的影响。但是通过自适应衰减和滤波后的投影数据重建的图像，金属结构一定会模糊。获取金属区域信息最好的方法当然是迭代法，迭代法能有效去除金属伪影和抑制噪声，并能完整的恢复金属信息。迭代法的最大缺点就是计算时间长，但是ＣＴ扫描对象中金属物体一般都较小，所以仅对金属区域进行迭代是运算量是可以接受的。我们可以通过局部迭代法获得较精确的金属区域信息，只要合理的利用这些信息对金属区域进行补偿，就能获得理想的图像结果。算法流程如图４．１所示。原始投影数据ＦＢＰ重建原始图像分割出金属区域图像ＥＭ迭代金属区域，重投影获取金属投影估计金属区域像素估计金属投影值并对其进行自适应衰减和滤波ＦＢＰ重建处理后图像对金属区域加权补偿图４．１算法流程图４．２算法与实现算法的主要步骤是先确定金属投影区域，对金属投影数据进行自适应衰减和滤波，然后通过ＦＢＰ重建图像，最后再根据局部迭代法计算出的金属区域像素点估计值对金属区域进行补偿运算。下面将详细的介绍算法的实现过程。４１第４章基于自适应衰减滤波的ＣＴ图像去金属伪影混合法４．２．１局部迭代算法的选取在算法实现前，必需先选择局部迭代算法。迭代法有很多种（具体参见３．２．２节１，如ＡＲＴ、ＳＡＲＴ、ＥＭ和ＯＳＥＭ，我们应该选择哪一种迭法，这就需要对这些主流的迭代法进行一下比较。比较主要是基于两个性能：时间和重建效果。时间又包括一次迭代需要时间，收敛速度和局部迭代算法复杂度。重建效果主要从去伪影效果和金属结构的清晰度来判断。为了对迭代法的时间性能进行比较，需要建立仿真实验模型（参见３．２．２．３节表３．２），实验采用平行光束解析法模拟投影数据（如果是扇形束扫描，迭代算法不变，只是系数计算略有不同），投影角度为０．１８０度，步长为３度，共６０个投影度。探测器数目为１２８个，重建图像大小为１２８ｘ１２８。（１）时间性能比较算法迭代一次耗时如表４．１所示，从表可知ＡＲＴ迭代一次的时间远大于其它迭代算法，主要是由于ＡＲＴ迭代是需要把第ｉ条射线校正后的像素点更新，才能进行第（ｉ＋１）条射线的校正，即需要动态的计算射线ｉ的估计投影值Ｒ工ｏ），而ＳＡＲＴ、ＥＭ和ＯＳＥＭ都是进行以轮或是一个子集的迭代才更新估计图像矢量Ｘ，因此可以在每轮迭代开始时计算好每条射线的估计投影值，以后每次校正时直接使用结果即可，大大减少了运算量。而且由于ＡＲＴ是逐线迭代，即使是局部迭代也需对每条射线进行校正，所以ＡＲＴ算法局部迭代运算量也最大。因为ＯＳＥＭ需要合理划分子集，所以它的局部迭代运算量也略大于ＳＡＲＴ和ＥＭ算法。表４．１迭代一次时间和局部迭代运算量表图４．２是ＡＲＴ、ＳＡＲＴ、ＥＭ和ＯＳＥＭ的均方差曲线图，从图中我们可以看出，ＡＲＴ收敛慢，而且随着迭代次数的增加，均方差呈发散趋势。ＥＭ和ＯＳＥＭ收敛第一次和第二次迭代时收敛’［＇曼－Ｔ＇ＳＡＲＴ（ＥＭ骨＊１１０ＳＥＭ第二次均方差分别为０．０１３６和０．０１２３），但随着迭代次数的增加，ＥＭ和ＯｓＥＭ的收敛速度远远快于Ｓ灿盯。ＥＭ收敛速度Ｌ匕ＯＳＥＭ慢一些，主要体现在均方差的数值上。４２第４章基于自适应衰减滤波的ＣＴ图像去金属伪影混合法（ａ）ＡＲＴ迭代均方差曲线（ｂ）ＳＡＲＴ迭代均方差曲线（ｄ）ＯＳＥＭ迭代均方差曲线图４．２迭代法均方差曲线对比图（２）重建结果比较由３．２．２．３节的图３．４可知，ＡＲＴ迭代算法能大大削弱了金属伪影，但是对金属物体引起的射束硬化暗带明显不明显，而且由于ＡＲＴ的逐线校正方式易受噪声干扰，导致金属物体周围出现了环状伪影。ＳＡＲＴ能基本消除了金属伪影，而且比ＡＲＴ更好的还原了金属信息，不过由于ＳＡＲＴ是使用通过该像素的所有射线的估计投影值的平均值来校正像素，这样虽然能有效抑制干扰，但是也会使非金属区域的一些细节变模糊。ＥＭ算法能很好的消除金属伪影，还原金属信息，效果要好于ＳＡＲＴ迭＿代法，而．ＫＥＭ收敛Ｉ：ＢＳＡＲＴ要快（迭代次数大于２次），但是ＥＭ同样在迭代次数不充分的情况下为丢失一些非金属区域的细节。ＯＳＥＭ同样能基本消除了金属伪影和射束硬化引起的暗带，而且能基本还原金属信息，但是同样会模糊细节，而且其收敛最快。综合以上算法的局部重建时间和重建效果，选则ＥＭ和ＯＳＥＭ是最好的，ＥＭ虽然收敛慢于ＯＳＥＭ，但是其算法局部迭代运算量小，而且ＯＳＥＭ还受子集分割的影响。因此选择了ＥＭ算法作为局部迭代算法。４３第４章基于自适应衰减滤波的ＣＴ图像去金属伪影混合法４．２．２算法的具体实现步骤４．２．２．１确定金属投影区域本算法和插值重建法与ＥＭ混合法一样，第一步都是确定金属投影区域，假设使用的是平行射线投影数据Ｐ（ｍ，ｎ）（如果是扇形束投影数据，我们认为是已经进行过重排），ｍ＝ｌ，．．．Ｍ，ｎ＝ｌ，．．．Ｎ。Ｍ是探测器总数，Ｎ是投影角度总数。使用ＦＢＰ重建图像，获得原始图像，选择适当的阈值提取金属区域，然后对金属区域利用旋转像素坐标系的方法，在每个投影角度下进行重投影，确定金属区域在投影数据中的范围。Ｕｎ）＝｛【ｓ（ｉ，ｎ），ｅ（ｉ，ｎ）】），ｉ＝１，．．．Ｉ（ｎ）（４．１）Ｉ（ｎ）是金属区域在角度ｎ下总的投影点数，ｓ（ｉ，ｎ），ｅ（ｉ，ｎ）是金属投影数据在每个角度下的起始和结束位置。虽然只有一个金属物体时可以直接从投影数据中分割出金属投影区域，但是重投影方法才是最稳定的分割金属投影数据的方法，它适用于多金属情况。４．２．２．２自适应衰减调节和滤波自适应衰减和滤波是算法中很重要的一步，它不仅减小了金属投影数据对周边区域的影响，使投影数据变得较为平滑（ＫａｃｈｅｒｉｅｓｓＭｅｔａｌ，２００１），降低了阈值选取精确度对算法的影响，大大削弱了金属伪影，而且还保留了被金属覆盖区域的细节特征。通过对４．２．２．１节中求出的金属投影数据区域【ｓ（ｉ，ｎ），ｅ（ｉ，ｎ）】线性插值，得到对金属区域的背景区域的投影值ｂ（ｍ，ｎ）ｍＥ“ｎ），ｂ（ｍ，ｎ）是对金属背景区域的粗略估计，求出金属的投影值Ｍｍ：Ｍｍ（ｍ，ｎ）＝【Ｐ（ｍ，ｎ）．ｂ（ｍ，ｎ）】，ｍＥＵｎ）对取出的金属投影值进行自适应衰减调节：（４．２）（４．３）Ｍｌ（ｍ，１３）＝ＸｃＸｂ（ｍ，ｎ）Ｍｍ（ｍ，ｎ），ｒｅＥＬ（ｎ）Ｘｂ（ｍ，ｎ）是衰减因子，根据投影数据Ｍｍ（ｍ，ｎ）的振幅计算，用于修正由金属物体引起的射束硬化。Ｘｃ是一个常量缩放因子，其值小于１，用于减少金属投影数据的振幅，金属区域投影数据的振幅越小，在ＦＢＰ重建时产生的金属伪影越小，Ｘｃ可取０．１．０．５，Ｘｃ越小，金属伪影也越小，但是会导致金属背景区域的细节失真。在很多实际情况中，ｘ射线的光谱和金属信息是不可知的，所以我们可以简单的令Ｘｂ（ｍ，ｎ）＝１，这样自适应衰减调节其实就蜕变为线性衰减，在仿真程序中也只提供一个缩放因子参数Ｘｃ＝０．１。第４章基于自适应衰减滤波的ｃｒ图像去金属伪影混合法由于在确定金属投影区域时选取的闽值不可能精确的划分出所有金属像素点，所以对金属投影区域进行自适应衰减调节后，金属投影区域附近仍然存在数据波动，需要通过滤波进一步减小数据的波动，我们对衰减后的金属投影数据ＭＩ（ｍ，ｎ）使用自适应滤波，进一步降低噪声，得到较平滑的金属投影数据Ｍ２（ｍ，ｎ１，这个滤波过程可以表示为：Ｍ２（ｍ，ｎ）＝Ｍｌ（ｍ，ｎ）＠Ｋ（ｍＸｍ∈Ｕｎ）（４．４）Ｋｆｍ）是一个滤波函数，可以是均值滤波，中值滤波和裁减中值滤波，根据经验，中值滤波也能很好的去除金属伪影和光子饥渴噪声。中值滤波函数窗宽的选取是由探测器总数决定的。仿真实验中对于１８０个探测器的ＣＴ，窗宽可以取５，随着探测器数的增大，应适当加大窗宽，但也不能过大，否则会使金属物体边缘形变。图４．２为衰减和滤波前后投影数据对比图，其中图４．３（ａ）为原始投影数据，中间突出部分就是金属区域投影数据。图４．３（ｂ）是经过自适因衰减和滤波后的投影数据（ｘｃ－－ｏ．１）。由图可知投影数据的跃变明显削弱，通过ＦＢＰ重建后的图像基本消除金属伪影，但是由于对投影数据进行衰减导致金属细节模糊和边界失真，因此还要进行金属区域补偿。４．２．２．３对金属区域进行补偿烹参一曲原始投影数据（ｂ）衰减和滤波后投影数据闰４３衰减和滤波前后的投影数据把４２．２２修正过的金属投影数据加入背景数据中Ｍ３（ｍ，ｎ）＝Ｍ２（ｍ，ｎ）＋ｂ（ｍ，ｎ），ｍ∈Ｌ（ｎ）（４．５）金属区域外的投影数据保持不变，对Ｐ３使用ＦＢＰ重建图像Ｆｎ。。Ｏ，ｊ），ｉｊ＝１，２，３…Ｍ，Ｍ是探测器总数。由于自适应衰减会模糊金属的细节，所以需要对图像Ｒ—ｉ，ｊ）的金属局域进行补偿．达到细节增强的目的。如何把金属信息恢复，就是这－４，节的主要工作。第４章基于自适应衰减滤波的ＣＴ图像去金属伪影混合法通过ＥＭ局部迭代，得到金属区域的估计像素值之后，如果像ＥＭ混合法一样，简单地将估计的金属区域取代ＦＢＰ重建后的金属区域，会使金属区域边缘失真，而且还会丢失部分由白适应衰减保留下来的金属覆盖区域的细节。因此可以通过金属区域的估计像素值，对衰减滤波后重建的图像进行加权补偿。如对原始投影数据Ｐ的金属区域使用ＥＭ迭代，获取金属区域的迭代像素值Ｆｅｒｎ（ｉ，ｊ），ｉ，ｊＥＳ，Ｓ表示金属区域。则对金属区域的加权补偿公式为：Ｆｎ。“ｉ，ｊ）＝（Ｆｎ附（ｉ，ｊ）＋Ｗａ木Ｆｅｍ（ｉ，ｊ））／Ｗｄｉ，ｊＥＳｗａ为常量权重因子，表示迭代的金属像素值的补偿度，由于ＥＭ迭代法是基于统计理论的迭代法，受噪声影响小，当迭代像素点较多时其计算的像素值可信度高，取Ｗａ大于等于１。但是当金属区域很小时，由于用于统计的像素点很少，导致可信度降低，此时Ｗａ的取值应小于１。Ｗｄ是一个调节显示窗宽的参数，取值一般为１０至ｌＪ２０。４．３仿真实验结果及分析为了考察本文算法的效果，通过仿真实验将其与商用ＣＴ主流去金属伪影算法（线性插值法、ＥＭ混合法）进行了比较。模型参数见表１，序号５、６两组参数代表一个有重合的金属区域。实验采用平行光束解析法模拟投影数据，投影角度为０．１８０度，步长为２度，共９０个投影度。探测器数目为１８０个，重建图像大小为１８０ｘ１８０。在仿真实验中除了加入高斯噪声，也通过公式（３．８）模拟了射束硬化。表１数值仿真模型参数本文算法的仿真结果如图４．４所示，图４．“ａ）是直接ＦＢＰ重建后的结果，可以看到图像金属周围出现条纹状伪影，而且中间出现明显暗带。图４．４（ｂ）是经过自适应衰减和滤波后重建的图像，可以看出有效的去除了金属伪影和硬化引起的暗带，但是金属物体对比度下降，而且边界失真。图４．４（ｃ）是对金属区域进行补偿后的结果，金属结构清晰、对比度明显，而且有效修正了金属物体边缘的失真。第４章基丁白适应衰减滤波的ｃＴ图像去金属伪影混合法■■■■囵囵囹第４章基于自适应衰减滤波的ＣＴ图像去金属伪影混合法ＥＭ混合法相当的情况下明显提高了成像质量。４．４本章小结本章设计了一种基于自适应衰减和滤波的ＣＴ去金属伪影混合法，该算法通过选取适当的衰减系数和滤波函数对金属投影数据进行自适应衰减和滤波，然后使用ＦＢＰ算法对处理后投影数据重建。并通过对原始投影数据进行ＥＭ局部迭代估计金属区域像素值，根据估计值对重建图像中的金属区域进行加权补偿。该算法能有效去除金属伪影，抑制噪声和减小射束硬化的影响，降低了算法处理效果对阈值选取精确度的依赖程度，并弥补了插值法和ＥＭ混合金属区域失真的不足，而且重建的金属的结构更清晰。该算法计算复杂度与ＥＭ混合法相当，能比较容易的移植到各种结构的ＣＴ中，因此具有很高的实用价值。此外，该算法可以进一步优化，当金属物体较大时，金属区域重建估计需要较长时间，可以根据图像的先验知识分布结合ＥＭ迭代或者使用分块迭代，都可以缩短重建时间。第５章总结和展望第５章总结和展望５．１总结ＣＴ是现在应用最为广泛的医学和工业检测成像设备之一，但是金属伪影的存在严重降低了图像质量，无论是在医疗或是工业检测都影响了图像的后续处理和分析。本文首先介绍ＣＴ的各种重建算法，然后详细介绍了三种主流的去金属伪影算法。并对这些算法进行了数值模拟仿真实验，通过研究算法的理论和进行仿真实验，对这些算法的优缺点进行了总结，在此基础上提出了一种新的去金属伪影算法。本文的主要研究工作有：（１）为了实现和验证ＣＴ去金属伪影重建算法，使用Ｍａｔｌａｂ开发工具，实现了一个简单的ＣＴ仿真平台，为进一步深入开展这方面研究打下了基础。（２）实现了３种主流的去金属伪影算法：插值重建法、迭代重建法和ＥＭ混合法，并总结了它们去金属伪影的优缺点。（３）对各种迭代法的实现提出了一些优化措施，并详细比较了这些迭代法的收敛速度和复杂度。（４）设计了一种基于自适应衰减滤波的ＣＴ去金属伪影混合法。该算法能有效去除金属伪影，抑制噪声和减小射束硬化的影响，弥补了插值法和ＥＭ混合金属区域失真的不足，而且使金属的结构更清晰；此外，该算法计算复杂度与ＥＭ混合法相当，移植到不同设计的ＣＴ中也比较容易。５．２展望本文的研究主要针对平行束和扇形束ＣＴ展开的，研究的范围及实现的功能还很有限。近些年ＣＴ的新技术发展迅速，应用领域也越来越广泛。对于Ｃ，ｒ中金属伪影的研究需要进一步深入，因为对于新的ＣＴ技术和应用，始终会遇到金属伪影问题，如何在这些新领域中去除金属伪影，提高ＣＴ成像的准确性，将会影响ＣＴ技术和应用的进一步发展。本文设计的自适应衰减滤波混合法，如果能一步完善，将有可能解决ＣＴ新技术和应用中的金属伪影问题，本人认为下面几个方面需要经一步研究：（１）ＣＴ的新技术主要是螺旋ＣＴ和多层ＣＴ，特别是在医疗领域，螺旋ＣＴ和多层ＣＴ已经被广泛的应用。由于螺旋Ｃ１’和多层ＣＴ和扫描方式与传统的平行束和扇形束ＣＴ不同，如何确定螺旋ＣＴ和多层ＣＴ投影数据中的金属投影区域，将会成为研究重４９第５章总结和展望点（ＣＴ的扫描方式对局部迭代法影响较小）。（２）确定金属区域时，为了达到最优划分对阈值的精确定有很高要求，现在还没有一种全自动的分割方法，适当人工干预是不可避免的。金属区域的确定，对金属结构信息、和成像方式的知识均是必不可少的。如何确定一些模型，在这些模型的基础上可以不需要人工干预而确定最佳阈值。因此基于模型的金属区域划分是今后去金属伪影算法的研究重点。（３）ＣＴ三维重建随着计算机性能的提高越来越受到重视，特别是在临床应用中通过ＣＴ三维重建，得到的结果可以充分的挖掘图像中所包含的有用信息，提高医疗诊断的准确性与科学性。因此将自适应衰减滤波混合法应用于ＣＴ－－－维重建，有利于扩大ＣＴ在医疗诊断和治疗中的应用范围。参考文献参考文献谷建伟，张丽．２００４．圆轨迹锥束ＣＴ的伪影成因和校正方法综述【Ｊ】．北京：ＣＴ幂ｎ－－维成像学术年会议文集．１０（２）：１２—１６谷建伟，张丽，陈志强．２００５．ＣＴ图像中金属伪影的快速校ｉＥ［Ｊ１．中国体视学与图像分析．１０（２）：１０８—１１１．谷建伟，张丽，陈志强．２００６．常见的几种插值方法在ＣＴ系统金属伪影校正中的性能比较【Ｊ】．核电子科学与探测技术．砸６）：９０５－９０７谷建伟，张丽，陈志强．２００６．基于投影的ＣＴ图像金属伪影非权重校ｉＥ［Ｊ】．清华大学学报．４６（６）：８２５－８２８黄庆华．２００２．ＣＴ图像重建与医学影像图像增强【Ｄ】．【硕士】．合肥：中国科学技术大学．１４李建，吴志芳．２００５．最大后验概率重建算法在发射ＣＴ中的应用【Ｊ】．核电子学与探测技术．２７（１）：４５·４９李建，吴志芳，安洪振，等．２００７．松弛系数．有序子集．最大似然期望算法在∞Ｃｏ集装箱硎佥测中的应用【Ｊ】．核电子学与探测技术．２７（１）：４５－４９吴朝霞，刘力，柴新禹．２００４．可变约束ＯＳＥＭ图像重建算法仿真研究［Ｊ］．中国生物医学工程学报．２３（１）：１－３林宙辰，石青云．２００１．用四次多项式插值消除医用ｘ射线ＣＴ中的金属伪影【Ｊ】．中国图象图形学报．６（Ａ）（动：１４２－１４７刘力，吴朝霞，赵书俊．２００２．具有松弛因子的ＯＳＥＭ重建算法［Ｊ］．中国图像图形学报．７（ＡＸｌ）－８１４－８１７张泽宏．２００５．多能ｘ射线硬化校正方法的研究【Ｄ】：【硕士１．重庆：重庆大学．１８－２５印胤，刘力．２００３．扇束图像中ＯＳＥＭ算法及子集划分的研究［Ｊ］．ＣＴ理论语应用研究．７（３）：１－８庄天戈主编．１９９２．ＣＴ原理与算法【Ｍ】．上海交通大学出版社．１６－２２张书文，陈英茂，田嘉禾，等．２００３．ＯＳＥＭ子集及迭代速度对ＰＥＴ巨］像质量的影响【Ｊ】．中国医学影像杂志．１１（３）：１９９－２０１张伟．２００６．ＣＴ图像质量评估及图像伪影模拟校正方法【Ｄ】．【硕士】．上海：上海交通大学．７．１０张书明，胡宇光，施元丁，等．２００７．同步辐射光成像中的最大似然迭代法【Ｊ】．ＣＴＮ论与应用研究．１６（２）：８—１３ＡｍｄｅｒｓｅｎＡＨ，ＫａｋＡＣ．１９８４．Ｓｉｍｕｌｔａｎｅｏｕｓａｌｇｅｂｒａｉｃｒｅｃｏｎｓｔｒｕｃｔｉｏｎｔｅｃｌｍｉｑｕｅ（ＳＡＲＴ）：ａｓｕｐｅｒｉｏｒｉｍｐｌｅｍｅｎｔａｔｉｏｎｏｆｔｈｅＡＲＴａｌｇｏｒｉｔｈｍ［Ｊ］．Ｕｌｔｒａｓｏｎｉｃ５１Ｉｍａｇｉｎｇ．６：８１－９４参考文献ＢａｒｒｅｔｔＨＨ，ＷｉｌｓｏｎＤＷ，ＴｓｕｉＢＭ．１９９４．ＮｏｉｓｅｐｒｏｐｅｒｔｉｅｓｏｆＥＭａｌｇｏｒｉｔｈｍ．ＩｍａｇｅＴｈｅｏｒｙＰｈｙｓｉｃａｌＭｅｄｉｃａｌＢｉｏｌｏｇｙ．ＢｙｒｎｅＣＬ１９９６．Ｂｌｏｃｋ－ｉｔｅｒａｔｉｖｅｍｅｔｈｏｄｓｆｏｒｉｍａｇｅｒｅｃｏｎｓｔｒｕｃｔｉｏｎｆｒｏｍＴｒａｎｓａｔｉｏｎｏｎｐｒｏｊｅｃｔｉｏｎｓ［Ｊ］．ＩＥＥＥＩｍａｇｅＰｒｏｃｅｓｓｉｎｇ．５（５）：７９２—７９４ＢａｉｓｓａｌｏｖＲ，ＳａｎｄｉｓｏｎＧＡ，ＤｏｎｎｅｌｌｙＢＪ．２０００．Ｓｕｐｐｒｅｓｓｉｏｎｏｆｈｉｇｈ·ｄｅｎｓｉｔｙａｒｔｅｆａｃｔｓｉｎｘ．ｒａｙＣＴｉｍａｇｅｓｕｓｉｎｇｔｅｍｐｏｒａｌｄｉｇｉｔａｌｓｕｂｔｒａｃｔｉｏｎｗｉｔｈａｐｐｌｉｃａｔｉｏｎｔｏｃｒｙｏｔｈｅｒａｐｙ［Ｊ］．ＰｈｙｓｉｃｓｉｎＭｅｄｉｃａｌａｎｄＢｉｏｌｏｇｙ．４５：Ｎ５３·Ｎ５９．ＣＲ．１９８６．ＲｅｐｒｏｊｅｃｔｉｏｎｕｓｉｎｇａＣｒａｗｆｏｒｄ４８０．４８３．ｐａｒａｌｌｅｌｂａｃｋｐｒｏｊｅｃｔｏｒ［Ｊ］．ＭｅｄｉｃａｌＰｈｙｓｉｃｓ。１３（４）：ＣｈｅｎＬＧＬｉａｎｇａｔｔｅｎｕａｔｉｏｎＹＧｅｏｒｇｅＡＳａｎｄｉｓｏｎ，ｅｔａ１．２００２．ＡｎｏｖｅｌｍｅｔｈｏｄｆｏｒｒｅｄｕｃｉｎｇｈｉｇｈｏｂｊｅｃｔａｒｔｉｆａｃｔｓｉｎＣＴｒｅｃｏｎｓｔｒｕｃｔｉｏｎｓ［Ｃ］．ＰｒｏｃｅｅｄｉｎｇｓｏｆＳＰＩＥ，４６８４：ＰＰ．８４１—８５０．ＤｅｍｐｓｔｅｒＡＰ，ＬａｉｒｄＮＭ，ＲｕｂｉｎＤＢ．１９７７．ＭａｘｉｍｕｍｌｉｋｅｈｏｏｄｆｒｏｍｉｎｃｏｍｐｌｅｔｅｄａｔａｖｉａｔｈｅＥＭａｌｇｏｒｉｔｈｍ［Ｊ］．ＪｏｕｒｎａｌｏｆｔｈｅＲａｙａｌＳｔａｔｉｓｔｉｃａｌＳｏｃｉｅｔｙＳｅｒｉｅｓＢ．３９（１）：１－１５．ｐｒｏｊｅｃｔｉｏｎｓｏｎＤａｖｉｓＧＲ．１９９４．ＥｆｆｅｃｔｏｆｌｉｎｅａｒｉｎｔｅｒｐｏｌａｔｉｏｎｏｆｔｈｅｆｉｌｔｅｒｅｄＸ－Ｒａｙｃｏｍｐｕｔｅｄｔｏｍｏｇｒａｐｈｅ．ＪｏｕｒｎａｌｏｆＸ－ｒａｙＳｃｉｅｎｃｅａｎｄＤｅｉｍａｇｅｎｏｉｓｅｉｎＴｅｃｈｎｏｌｏｇｙ［Ｊ］．４１：１９２－１９９．ＭａｎＢ，ＮｕｙｔｓＪ，ＤｕｐｏｎｔＰ，ｅｔａ１．１９９９．ＭｅｔａｌｓｔｒｅａｋａｒｔｉｆａｃｔｓｉｎＸ—ｒａｙｃｏｍｐｕｔｅｄＴｒａｎｓａｃｔｉｏｎｓｏｎｔｏｍｏｇｒａｐｈｙＡｓｉｍｕｌａｔｉｏｎＤｅｓｔｕｄｙ［Ｊ］．ＩＥＥＥＮｕｃｌｅａｒＳｃｉｅｎｃｅ４６（３）：６９１—６９６．ＭａｎＢ，ＮｕｙｔｓＪ，ＤｕｐｏｎｔａＰ，ｅｔａ１．２０００．Ｒｅｄｕｃｔｉｏｎｏｆｍｅｔａｌｓｔｒｅａｋａｒｔｉｆａｃｔｓｉｎｘ—ｒａｙｃｏｍｐｕｔｅｄｔｏｍｏｇｒａｐｈｙｕｓｉｎｇＴｒａｎｓａｃｔｉｏｎｓｏｎｔｒａｎｓｍｉｓｓｉｏｎｍａｘｉｍｕｍａｐｏｓｔｅｒｉｏｒｉａｌｇｏｒｉｔｈｍ【Ｊ】．ＩＥＥＥＮｕｃｌｅａｒＳｃｉｅｎｃｅ．４７（３）：９７７－９８１．ＤａｎＪＫ．２００１ＳｔａｔｉｓｔｉｃａｌｌｙｒｅｇｕｌａｔｅｄａｎｄａｄａｐｔｉｖｅＥＭｒｅｃｏｎｓｔｒｕｃｔｉｏｎｆｏｒｅｍｉｓｓｉｏｎｃｏｍｐｕｔｅｄｔｏｍｏｇｒａｐｈｙ［Ｊ］．ＩＥＥＥＤｅＴｒａｎｓａｃｔｉｏｎｏｎＮｕｃｌｅａｒＳｃｉｅｎｃｅ．镐（３）：７９０—７９７ＭａｎＢ，ＮｕｙｔｓＪ，ＤｕｐｏｎｔＰ，ｅｔａ１．２００１．Ａｎｉｔｅｒａｔｉｖｅｍａｘｉｍｕｍ—ｌｉｋｅｌｉｈｏｏｄｐｏｌｙｃｈｒｏｍａｔｉｃｏｎａｌｇｏｒｉｔｈｍｆｏｒｏｒ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓＭｅｄｉｃａｌＩｍａｇｉｎｇ．２０（１０）：９９９－１００８．ＦｅｓｓｌｅｒＪＡ，ＨｅｒｏＡＯ．１９９５．Ｐｅｎａｌｉｚｅｄｍａｘｉｍｕｍ－ｌｉｋｅｌｉｈｏｏｄｉｍａｇｅｒｅｃｏｎｓｔｒｕｃｔｉｏｎｕｓｉｎｇｓｐａｃｅＯｎＩｍａｇｅａｌｔｅｒｎａｔｉｎｇｇｅｎｅｒａｌｉｚｅｄＥＭ１４１７．１４１８ａｌｇｏｒｉｔｈｍ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓＰｒｏｃｅｓｓｉｎｇ．４（１０）：ＧｏｒｄｏｎＲ，ＢｅｎｄｅｒＲ，Ｈｅｒｍａｎｔｈｒｅｅ－ｄｉｍｅｎｓｉｏｎａｌｅｌｅｃｔｉｏｎＧＴ．１９７０．Ａｌｇｅｂｒａｉｃｒｅｃｏｎｓｔｒｕｃｔｉｏｎａｎｄｘ—ｒａｙｔｅｃｈｎｉｑｕｅｓ（ＡＲＴ）ｆｏｒｍｉｃｒｏｓｃｏｐｙｐｈｏｔｏｇｒａｐｈｙ［Ｊ］．ＪｏｕｒｎａｌｏｆＴｈｅｏｒｅｔｉｃａｌＢｉｏｌｏｇｙ．２９（３）：４７１－４７５ＧｌｏｖｅｒＧＨ，ＰｅｌｃＮＪ．１９８１．ＡｎａｌｇｏｒｉｔｈｍｆｏｒｔｈｅｒｅｄｕｃｔｉｏｎｏｆｍｅｔａｌｃｌｉｐａｒｔｉｆａｃｔｓｉｎＣＴｒｅｃｏｎｓｔｒｕｃｔｉｏｎｓ［Ｊ］．ＭｅｄｉｃａｌＰｈｙｓｉｃｓ．８（６）：７９９—８０７，ＧｕａｎＨＱ．ＲｉｃｈａｒｄＧ．１９９４．ＡｐｒｏｊｅｃｔｉｏｎａｃｃｅｓｓｏｒｄｅｒｆｏｒｓｐｅｅｄｙｃｏｎｖｅｒｇｅｎｃｅｏｆＡＲＴｉｎ（ａｌｇｅｂｒａｉｃｒｅｃｏｎｓｔｒｕｃｔｉｏｎｔｅｃｈｎｉｑｕｅ）：ａｍｕｌｔｉｌｅｖｅｌｓｃｈｅｍｅｆｏｒｃｏｍｐｕｔｅｄｔｏｍｏｇｒａｐｈｙ［Ｊ］．Ｐｈｙｓｉｃｓ５２参考文献ＭｅｄｉｃａｌａｎｄＢｉｏｌｏｇｙ．３９：２００５－２０１０ＨｅｒｍａｎＧＴ．１９７９．Ｃｏｒｒｅｃｔｉｏｎｆｏｒｂｅａｍｈａｒｄｅｎｉｎｇｉｎｃｏｍｐｕｔｅｄｔｏｍｏｇｒａｐｈｙ［Ｊ］．ＰｈｙｓｉｃｓｉｎＭｅｄｉｃａｌａｎｄＢｉｏｌｏｇｙ．２４（１）：８１－１０６．ＨｅｒｍａｎＧＴ．１９９５．Ｉｍａｇｅｒｅｃｏｎｓｔｒｕｃｔｉｏｎｆｒｏｍｐｒｏｊｅｃｔｉｏｎｓ［Ｊ］．Ｒｅａｌ－ＴｉｍｅｃｏｍｐｕｔｅｄＩｍａｇｉｎｇ．１：３—１８ＨｓｉｅｈＪ．１９９８．Ａｄａｐｔｉｖｅｓｔｒｅａｋａｒｔｉｆａｃｔｒｅｄｕｃｔｉｏｎｉｎｔｏｍｏｇｒａｐｈｙｒｅｓｕｌｔｉｎｇｆｒｏｍｅｘｃｅｓｓｉｖｅｘ－ｒａｙｐｈｏｔｏｎＰｎｏｉｓｅ［Ｊ］．ＭｅｄｉｃａｌＰｈｙｓｉｃｓ．２５（１１）：２１３９—２１４７．ｉｍｐｒｏｖｅｄａｌｇｏｒｉｔｈｍＪｏｓｅｐｈＭ．１９８２．Ａｎｏｎｆｏｒｒｅｐｒｏｊｅｃｔｉｎｇｒａｙｓｔｈｒｏｕｇｈｐｉｘｅｌｉｍａｇｅｓ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓＫａｌｅｎｄｅｒＭｅｄｉｃａｌＩｍａｇｉｎｇ．１（３）：１９２—１９６Ｊ．１９８７．ＲｅｄｕｃｔｉｏｎｏｆＣＴａｒｔｉｆａｃｔｓｃａｕｓｅｄｂｙｍｅｔａｌｌｉｃＷＡＨｅｂｅｌＲ，Ｅｂｅｒｓｂｅｒｇｅｒｉｍｐｌａｎｔｓ［Ｊ］．Ｒａｄｉｏｌｏｇｙ．１６４（２）：５７６－５７７．ＫａｍｐｈｕｉｓａｌｌＣＫＢｅｃｋｍａｎＦＪ．１９９８．ＡｃｃｅｌｅｒａｔｅｉｔｅｒａｔｉｖｅｔｒａｎｓｍｉｓｓｉｏｎＣＴｒｅｃｏｎｓｔｒｕｃｔｉｏｎｕｓｉｎｇｏｎｏｒｄｅｒｅｄｓｕｂｓｅｔＫａｃｈｅｌｒｉｅｓｓｃｏｎｖｅｘａｌｇｏｒｉｔｈｍ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓＯ，Ｋａｌｅｎｄｅｒｍｅｄｉｃａｌｉｍａｇｉｎｇ．１７（８）：１１０１－１１０５．Ｍ，ＷａｔｚｋｅＷＡ．２００１．Ｇｅｎｅｒａｌｉｚｅｄｍｕｌｔｉ－ｄｉｍｅｎｓｉｏｎａｌａｄａｐｔｉｖｅｆｉｌｔｅｒｉｎｇｆｏｒＰｈｙｓｉｃｓｃｏｎｖｅｎｔｉｏｎａｌａｎｄｓｐｉｒａｌｓｉｎｇｌｅ－ｓｌｉｃｅ，ｍｕｌｔｉ－ｓｌｉｃｅ，ａｎｄｃｏｎｅ－ｂｅａｍＣＴ［Ｊ］．Ｍｅｄｉｃａｌ２８（４）：４７５－４９０．ＲＭ，ＢａｔｅｓＲＨ．１９８３．Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎａｌｇｏｒｉｔｈｍｓ：ｔｒａｎｓｆｏｒｍｍｅｔｈｏｄｓ［Ｊ］．ＰｒｏｃｅｅｄｉｎｇｏｆＬｅｗｉｔｔｔｈｅＩＥＥＥ．７１（３）：３９０－４０８．ＬａｎｇｅＫＣａｒｓｏｎＲ．１９８４．ＥＭｒｅｃｏｎｓｔｒｕｃｔｉｏｎｏｆＣｏｍｐｕｔｅｒＡｓｓｉｓｔｅｄａｌｇｏｒｉｔｈｍｓｆｏｒｅｍｉｓｓｉｏｎａｎｄｔｒａｎｓｍｉｓｓｉｏｎｔｏｍｏｇｒａｐｈｙ［Ｊ］．ＪｏｕｒｎａｌＬｅａｈｙａｎｄＴｏｍｏｇｒａｐｈｙ．８（２）：３０６－３１６ＲＭ，ＢｙｒｎｅＣＬ．２０００．ＲｅｓｅｎｔｄｅｖｅｌｏｐｍｅｎｔｓｉｎｉｔｅｒａｔｉｖｅｉｍａｇｅｒｅｃｏｎｓｔｒｕｃｔｉｏｎｆｏｒＰＥＴｏｎＳＰＥＣＴ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓＬａＲｉｖｉ６ｒｅｍｅｄｉｃａｌｉｍａｇｉｎｇ．１９（４）：２５７－２６０．ｒｅｇｒｅｓｓｉｏｎＰＪ，ＰａｎＸ．２０００．Ｎｏｎｐａｒａｍｅｔｒｉｃｌｉｋｅｌｉｈｏｏｄｓｉｎｏｇｒａｍｓｍｏｏｔｈｉｎｇｏｎｕｓｉｎｇａｒｏｕｇｈｎｅｓｓ－ｐｅｎａｌｉｚｅｄｐｏｉｓｓｏｎｏｂｊｅｃｔｉｖｅＦｕｎｃｔｉｏｎ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｍｅｄｉｃａｌｉｍａｇｉｎｇ．１９（８）：７７３－７８６．ＰａｎＸ．１９９９．Ｏｐｔｉｍａｌｔｏｍｏｇｒａｐｈｙｎｏｉｓｅｃｏｎｔｒｏｌｉｎａｎｄｆａｓｔｒｅｃｏｎｓｔｒｕｃｔｉｏｎｏｆｆａｎ·ｂｅａｍｃｏｍｐｕｔｅｄｉｍａｇｅ［Ｊ］．ＭｅｄｉｃａｌＰｈｙｓｉｃｓ．２６（５）：６８９－６９７．Ｎｏｉｓｅ－ＩｎｄｕｃｅｄＳｔｒｅａｋＡｒｔｉｆａｃｔｓｉｎＰａｔｒｉｃｋＪ，Ｌ丑Ｒ，ＤａｖｉｄＭＢ，ｅｔａ１．２００５．ＲｅｄｕｃｔｉｏｎｏｆＸ－ＲａｙＣｏｍｐｕｔｅｄＴｏｍｏｇｒａｐｈｙＴｈｒｏｕｇｈＳｐｌｉｎｅ·－ＢａｓｅｄＰｅｎａｌｉｚｅｄ··ＬｉｋｅｌｉｈｏｏｄＳｉｎｏｇｒａｍＳｍｏｏｔｈｉｎｇ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｏｎｍｅｄｉｃａｌｉｍａｇｉｎｇ．２“１）：１０５—１１１ＳｅｉｔｚＥＲｕｅｇｓｅｇｇｅｒＰ．１９８５．ＣＴｂｏｎｅｄｅｎｓｉｔｏｍｅｔｒｙｏｆｔｈｅａｎｃｈｏｒａｇｅｏｆａｒｔｉｆｉｃｉａｌｋｎｅｅｊｏｉｎｔｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔｅｒＡｓｓｉｓｔｅｄＴｏｍｏｇｒａｐｈｙ．９：６１２－６２２．ＲｅｃｏｎｓｔｒｕｃｔｉｏｎｉｎＥｍｉｓｓｉｏｎＳｈｅｐｐＬＡ，ＶａｒｄｉＹ．１９８２．ＭａｘｉｍｕｍＬｉｋｅｌｉｈｏｏｄＩＥＥＥＴｒａｎｓａｃｔｉｏｎＴｏｍｏｇｒａｐｈｙ［Ｊ］．ＭｅｄｉｃａｌＩｍａｇｉｎｇ．１（２）：１１３·１２２．５３参考文献ＳａｎｅｒＫＢｏｕｍａｎＣ１９９３．ＡＴｒａｎｓａｃｔｉｏｎｓｏｎｌｏｃａｌｕｐａａｔｅｓｔｒａｔｅｇｙｆｏｒｉｔｅｒａｔｉｖｅｒｅｃｏｎｓｔｒｕｃｔｉｏｎｆｒｏｍｐｒｏｊｅｃｔｉｏｎｓ［Ｊ］．ＩＥＥＥＳｉｇｎａｌＰｒｏｃｅｓｓｉｎｇ．４１（２）：５４０－５４８ＷａｎｇＧＳｎｙｄｅｒＤＬ．，Ｏ’ＳｕｌｌｉｖａｎＪＡ，ｅｔａ１．１９９６．ＩｔｅｒａｔｉｖｅＤｅｂｌｕｒｒｉｎｇｆｏｒＣＴＭｅｔａｌＡｒｔｉｆａｃｔＲｅｄｕｃｔｉｏｎ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓＸＩＡＤ，ＲｏｅｓｋｅＪｏｈｎｏｎｍｅｄｉｃａｌｉｍａｇｉｎｇ．１５（５）：６５７—６６４．ｒｅｄｕｃｉｎｇＣ，ＹＵ崛ｅｔａ１．２００５．ＡｈｙｂｒｉｄａｐｐｒｏａｃｈｔｏｃｏｍｐｕｔｅｄｔｏｍｏｇｒａｐｈｙｍｅｔａｌａｒｔｉｆａｃｔｓｉｎｉｎｔｒａｃａｖｉｔａｒｙＹａｖｕｚｂｒａｃｈｙｔｈｅｒａｐｙ【Ｊ】．Ｂｒａｃｈｙｔｈｅｒａｐｙ．４（１）：１８－２３．ａｎｄｎｏｉｓｅａｎａｌｙｓｉｓＭ，ＦｅｓｓｌｅｒＪＡ．１９９９．Ｐｅｎａｌｉｚｅｄ·ｌｉｋｅｌｉｈｏｏｄｅｓｔｉｍａｔｏｒｓＰＥＴｆｏｒｒａｎｄｏｍｓ．ｐｒｅｃｏｒｒｅｃｔｅｄ６６５—６７４．ｔｒａｎｓｍｉｓｓｉｏｎｓｃａｎｓ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｏｎｍｅｄｉｃａｌｉｍａｇｉｎｇ．１８（８）：ＺｈａｏＳＹ，ＲｏｂｅｒｔｓｏｎＤＤ，ＷａｎｇＧ，ｅｔａ１．２０００．Ｘ－ｒａｙＣＴｍｅｔａｌａｒｔｉｆａｃｔｒｅｄｕｃｔｉｏｎｕｓｉｎｇｗａｖｅｌｅｔｓ：ａｎａｐｐｌｉｃａｔｉｏｎｆｏｒｉｍａｇｉｎｇｔｏｔａｌｈｉｐｐｒｏｓｔｈｅｓｅｓ［Ｊ］．ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｏｎｍｅｄｉｃａｌｉｍａｇｉｎｇ．１９（１２）：１２３８－１２４７．攻读硕士学位期间发表的学术论文攻读硕士学位期间发表的学术论文已发表论文：［１］陈豫，罗海，周荷琴．基于自适应滤波的ＣＴ图像去金属伪影混合算法．中国医疗器械杂志．２００９，３３（２）：８７－９０５５致谢致谢在本论文完成之际，谨向所有给予我指导和帮助的老师，同学，朋友和家人致以最诚挚的感谢。首先要特别感谢导师周荷琴教授，在信息工程研究实验室四年难忘的时光里，正是在她耐心的指导和悉心的关怀下，使我充实和愉快地完成了本科、研究生论文和一系列的工程项目。此外，周老师渊博的知识，开阔的视野，前瞻的眼光以及敏锐的思维给了我很深的印象，并在项目选题和研究过程中给了我很多指导和意见，不仅开拓了我的思维，更培养了我的学习能力和实际解决问题的能力。在为人方面，周老师严谨的治学精神，一丝不苟的工作态度和积极的生活态度深深的影响了我，让我养成了严肃认真的科研态度及乐观的生活态度，这将对我今后的工作和生活都大有益处。其次要感谢信息工程研究实验室的谭裴师兄，罗海及实验室的其他成员，大家的交流和讨论，极大的促进了我的学习和研究。还要感谢王雷老师和０６研的所有同学，感谢你们对我生活上的帮助，三年的早夕相处给我留下了很多美好的回忆。最后感谢父母和亲人，你们的支持是我最大的动力。２００９年５月５６CT图像金属伪影校正算法研究

作者：

学位授予单位：

陈豫

中国科学技术大学

本文链接：http://d.g.wanfangdata.com.cn/Thesis_Y1497509.aspx

授权使用：中北大学(zbdxtsg)，授权号：505296ac-4552-416f-a8f5-9e99016592f7

下载时间：2011年3月1日