当前位置:首页 >> 信息与通信 >>

CT图像环形伪影校正方法研究


CT 图像环形伪影校正方法研究

重庆大学硕士学位论文

学生姓名:黄苏红 指导教师:王 珏 教 授

蔡玉芳 助理研究员 专 业:控制科学与工程

学科门类:工 学

重庆大学自动化学院
二 O 一一年四月

Research on Method

for Removing Ring Artifacts in CT Images

A Thesis Submitted to Chongqing University in Partial Fulfillment of the Requirement for the Degree of Master of Engineering

By Huang Suhong

Supervised by Prof. Wang Jue Assistant Supervised by Cai Yufang Major: Control Science and Engineering

College of Automation of Chongqing University, Chongqing, China April, 2011

中文摘要





CT 技术不仅可以精确地、清晰地显示物体内部细节的结构关系及缺陷状况, 还可以定量地给出细节的辐射吸收数据及物质组成,已在医学的病情诊断和工业 的缺陷识别及属性测量中得到广泛的应用。三代扫描模式具有扫描速度快的优点, 是目前医学 CT 和工业 CT 采用的主要扫描模式。环形伪影作为三代扫描模式的一 个重要问题,不仅影响了图像的质量,还给图像噪声处理以及图像分割等后续处 理造成困扰,降低了图像的识别能力和测量精度。这个问题工业 CT 成像系统检测 中更是突出,因此,环形伪影的去除是 CT 图像预校正的关键一步。 根据重建图像中的环形伪影在投影正弦图中表现为竖直方向直线的特征,从 投影正弦图的角度,提出两种环形伪影校正算法: 第一种是基于投影曲线高通滤波的环形伪影校正法。首先,采用 S-L 滤波函 数对原始投影数据进行滤波,增强伪影信息;其次,对滤波后的投影数据进行投 影积分,并对投影曲线采用微分处理,进一步增大伪影数据与工件轮廓数据的差 异;然后,按照高频采样插值次数对差分后的投影数据进行均值滤波,并根据粗 大误差判断准则选择阈值,确定伪影的位置;最后,分析传统线性插值校正不彻 底的原因,采用线性外推插值来对伪影处投影数据进行校正。 第二种是基于改进的 Canny 算法的环形伪影校正法。在高斯滤波前,首先对 原始投影数据采用 S-L 滤波器进行滤波,增强伪影信息;采用二维高斯函数的一 阶偏导数构造滤波器计算梯度,降低噪声的影响;设定角度阈值对梯度方向进行 限制,排除斜角方向边缘的检测,实现竖直方向边缘点检测;采用梯度阈值和链 长度阈值实现伪影边缘点的检测与连接;分析了各行投影数据插值校正的问题, 提出分段 B 样条拟合法对投影数据进行校正,实现多个连续的环形伪影校正。 通过对实际 CT 图像进行校正实验,结果表明,这两种校正方法都能有效地消 除 CT 图像的环形伪影,又很好保持了图像边缘及分辨率。其中,第一种校正算法 实现简单,对单个分散的环形伪影校正效果好,但是对多个连续的环形伪影校正 受限;第二种实现了多个连续的环形伪影的校正,但是对于最大连接长度小于链 长度阈值的伪影无法检测。它们各有特点,互相补充,但也不能完全互相取代。 关键词:CT 图像,环形伪影,投影正弦图,伪影校正

I

重庆大学硕士学位论文

II

英文摘要

ABSTRACT
CT technology not only is able to show that structure relation and defect status of object interior details accurately and clearly, but also can quantitatively give the radiation density data and material composition. So it is widely applied in medical illness diagnoses and industrial defect recognition and property measurement. Because of the advantages of rapid scanning, the third generations scanning mode has become dominating in MCT and ICT fields. As one of the important problems of the third generations scanning mode, ring artifacts not only effect the image quality, but also hamper subsequent processing, such as image noise processing, image segmentation, and so on. Finally, image recognition ability and measurement accuracy is reduced. This problem is particularly outstanding in ICT imaging system. Therefore, ring artifact removing is the key step to pre-correction of CT images. According to the characteristics that the ring artifacts of CT image show vertical directional lines in the sinogram of CT images, two kinds of ring artifact correction algorithms were proposed based on the sinogram. The first ring artifact correction method was based on projection profile high-pass filter. Firstly, the S-L filter was adopted to filter original projection data and enhance artifact information. Secondly, the filtered projection data was processed by projection integral, and then the projection profile was processed by differential to further increase the difference between artifacts edge data and work piece contour data. Thirdly, the differenced projection data was mean filtered according to the high-frequency sampling interpolation times, and artifacts’ location was automatically searched according to threshold selected by gross error criterion. Finally, through analyzing the reasons why the ring artifacts could not be corrected by traditional linear interpolation, the linear extrapolation interpolation was proposed to correction the projection data of ring artifacts. The second ring artifact correction method was based on the improved Canny algorithm. Firstly, the S-L filter was adopted to filter original projection data and enhance artifact information before using Gaussian filter. Then the gradient was acquired by the filter, which composed by 2d first-order partial derivative to reduce the noise affection. Thirdly, the angle thresholds are established for limiting gradient direction, particularly excluding diagonal direction edge detection and achieve vertical
III

重庆大学硕士学位论文

direction edge detection. Fourthly, a gradient threshold and an angle threshold were used to detect and connect marginal point of the artifacts. Finally, through analyzing the problems that the traditional linear interpolation algorithm is employed to correct the projection data, the subsection B-spline fitting method is proposed to implement multiple continuous ring artifacts correction. Based on the actual CT images calibration experiment, the result indicates that the two proposed correction algorithms is able to remove the ring artifacts effectively and keep image edges and resolution unchanged. On one hand, the first correction algorithm is simple, and effective to remove single dispersed ring artifacts. But it is limited for multiple continuous ring artifacts. On the other hand, the second correction algorithm is able to reduce the multiple continuous ring artifacts, yet not effective to detect the ring artifacts which the maximum connecting length is smaller than chain length threshold. The two proposed methods have different characteristics, complement each other in some aspects,but can not totally replace. Keywords: CT image, Ring artifact, Projection sinogram, Artifact correction

IV









中文摘要 ..........................................................................................................................................I 英文摘要 ....................................................................................................................................... III 1 绪论 .............................................................................................................................................. 1
1.1 课题的研究背景及意义 ........................................................................................................... 1 1.2 CT 图像伪影概述 ...................................................................................................................... 1 1.3 CT 图像环形伪影校正研究现状 .............................................................................................. 3 1.4 论文的研究内容与结构安排 ................................................................................................... 4

2 CT 图像环形伪影的产生及影响 ................................................................................. 5
2.1 CT 系统三代扫描方式及重建算法概述 .................................................................................. 5 2.1.1 投影数据和正弦图 ............................................................................................................ 5 2.1.2 CT 系统三代扫描方式概述 ............................................................................................... 7 2.1.3 扇束滤波反投影重建算法 ................................................................................................ 7 2.2 环形伪影的产生原因 ............................................................................................................... 9 2.3 环形伪影特征及分析 ............................................................................................................. 12 2.4 环形伪影的影响 ..................................................................................................................... 13 2.5 本章小结 ................................................................................................................................. 14

3 环形伪影校正的常用方法 ............................................................................................ 15
3.1 改变扫描方式校正法 ............................................................................................................. 15 3.2 暗场和增益不一致校正法 ..................................................................................................... 16 3.3 投影正弦图的校正 ................................................................................................................. 18 3.3.1 Hough 变换插值校正法 ................................................................................................... 18 3.3.2 小波-傅里叶变换滤波法 ................................................................................................. 21 3.4 重建图像的校正 ..................................................................................................................... 25 3.5 本章小结 ................................................................................................................................. 26

4 基于投影曲线高通滤波的环形伪影校正 ............................................................ 27
4.1 基于投影曲线高通滤波的条状伪影特征识别...................................................................... 27 4.1.1 伪影增强 .......................................................................................................................... 27 4.1.2 投影积分 .......................................................................................................................... 30 4.1.3 一阶微分 .......................................................................................................................... 32 4.1.4 均值滤波 .......................................................................................................................... 33 4.1.5 阈值选择 .......................................................................................................................... 34

V

重庆大学硕士学位论文

4.2 单环伪影投影数据校正 ......................................................................................................... 36 4.2.1 传统的线性插值 .............................................................................................................. 36 4.2.2 存在的问题 ...................................................................................................................... 36 4.2.3 线性外推插值校正 .......................................................................................................... 37 4.3 实验结果与分析 ..................................................................................................................... 38 4.3.1 伪影识别结果 .................................................................................................................. 38 4.3.2 伪影校正结果 .................................................................................................................. 39 4.3.3 分析与总结 ...................................................................................................................... 42 4.4 本章小结 ................................................................................................................................. 43

5 基于改进的 Canny 算法的环形伪影校正 .......................................................... 45
5.1 基于改进 Canny 算法的条状伪影特征识别 ........................................................................ 45 5.1.1 传统 Canny 算法的准则与实现...................................................................................... 45 5.1.2 存在的问题 ...................................................................................................................... 48 5.1.3 改进的 Canny 算法 ......................................................................................................... 48 5.1.4 算法的实现 ...................................................................................................................... 51 5.2 环带伪影投影数据校正 ......................................................................................................... 51 5.2.1 问题分析 .......................................................................................................................... 51 5.2.2 分段 B-样条拟合校正 ..................................................................................................... 53 5.2.3 算法的实现 ...................................................................................................................... 57 5.3 实验结果与分析 ..................................................................................................................... 58 5.3.1 伪影识别结果 .................................................................................................................. 58 5.3.2 伪影校正结果 .................................................................................................................. 59 5.3.3 分析与总结 ...................................................................................................................... 61 5.4 本章小结 ................................................................................................................................. 62

6 结论和展望 ............................................................................................................................ 63
6.1 论文工作总结 ......................................................................................................................... 63 6.2 后续工作展望 ......................................................................................................................... 63

致 谢 ...................................................................................................................................... 65 参考文献 ...................................................................................................................................... 67 附 录 ...................................................................................................................................... 71

VI

1 绪论

1 绪 论
1.1 课题的研究背景及意义
CT(Computed Tomography)技术,即计算机断层扫描技术。它是现代核物理学 和图像技术理论相结合而产生的一门边缘交叉的新兴学科,近年来已发展为无损 检测的一个重要领域[1~2]。CT 技术作为无损检测技术重要的手段之一,与其它无 损检测相比,例如超声检测、射线照相检测、红外检测和微波检测等技术,它具 有检测速度快、图像分辨率高、不受被测工件几何结构限制等优点[3]。不仅可以精 确地、清晰地显示物体内部细节的结构关系及缺陷状况,还可以定量地给出细节 的辐射密度数据及物质组成,最后一点是任何其它无损检测手段所不具有的特点
[1]

。所以,当今国际上已将 CT 技术作为最佳的无损检测技术。 CT 技术首先应用于医学领域, 形成了医学 CT(MCT)技术, 不仅广泛地应用于

人类疾病诊断,还被大量地应用于动物实验去推动生物学和药学发展,被评价为 是医学诊断的革命。由于医用 CT 技术应用的成功,以及其检测速度快、图像质量 高、直观地显示被检测物的形状、内部结构、密度和材质组成等优点,使得 CT 技 术很快地进入了工业领域,形成了工业 CT 技术。随着科学技术的发展,工业 CT 由 X 射线成像技术发展到 γ 射线成像技术, 再发展到新兴的太赫兹(THz)成像技术, 已广泛应用于航空、航天、汽车等工业中精密工件内部结构的测量及其缺陷检测; 兵工工业中对弹药填充密度的检测和武器关键部件的质量检查;钢铁工业中铸件 质量的在线监控和检测;油钻探测中岩芯评估和管道损伤的检测;海关和港口中 对走私、违禁物品的检查等[3~4]。 无论在医学还是工业 CT 系统中,对被检工件进行测量和识别都依赖于 CT 图 像的质量, 只有正确和清晰的 CT 图像才能让检测者做出正确的测量与识别。 但是, CT 技术本身的原理以及 CT 设备的硬件、软件及扫描工艺技术都会导致各种噪声 和伪影。这些伪影不仅影响了图像的质量,同时还给图像噪声处理以及图像分割 等后续处理造成困扰,降低了图像的识别能力和测量精度,这个问题工业 CT 成像 系统检测中更是突出。而环形伪影作为常见伪影的一种,已成为影响 CT 图像质量 的一个重要因素,因此,研究 CT 图像环形伪影校正方法是十分必要的。

1.2 CT 图像伪影概述
伪影是指在 CT 图像上出现了与物体的物理结构不相符的图像特征[1]。这些伪 影产生的原因很多,其中,CT 系统常见的是采样混叠、部分容积效应、射束硬化、 散射、探测器不一致性和金属伪影等。
1

重庆大学硕士学位论文

混叠是对原始投影数据进行采样时,采样密度不够导致谱重叠,从而导致 CT 图像产生混叠伪影;部分容积效应是由被测工件部分地伸入扫描视场所致,随着 切片厚度的增加,CT 系统发生部分容积效应的可能性也增加;射线硬化是由多色 的 X 射线束能谱和与能量相关的衰减系数造成的,高能量 X 射线穿透力强,低能 射线穿透力差,导致图像中圆柱形物体形成环状阴影伪影;散射使得探测器采集 的投影数据偏离了 X 射线测量的真实值,从而导致了重建图像中 CT 数偏移或阴 影伪像;探测器系统失调、增益变化、非线性响应和辐射损伤都会导致探测器响 应不一致,使得 CT 图像中出现环形伪影;金属伪影则是多种产生伪影的因素的综 合结果,往往强度很高,且影响着金属物周边的检测。 从伪影形貌上,伪影主要可分为 4 种类型:条状、阴影、环状及带状伪影等。 条状伪影通常表现为横穿 CT 图像的明显直线,它们可以是亮的,也可以是暗 的。一个条纹通常是由出现在孤立测量中的不一致性引起,这是由于滤波反投影 过程将投影空间中每个数据点映射到图像域中一条直线上。这种不一致性可能是 与数据采集过程相关的内在问题、机械故障或不同观测之间突然变化的结果。单 个探测器单元损坏也会出现明显的直线形伪影。还有一种常见情况就是射线在某 些路径上衰减过大,这时严重的统计噪声也会引起直线状伪影,不过这时可能是 一条或者一簇直线。 阴影伪影经常出现在被测工件的高对比度区域附近,它也是由投影中不一致 性的测量所产生的,但是,与条状伪影不同的是,它是由偏离真实测量结果的一 组通道或投影观测导致的。由于信号中没有明显间断性,这些误差导致图像中的 工件没有清晰的边界。根据错误通道的数量和误差大小,阴影伪影可以局限于局 部,或者影响更大区域。 环状和带状伪影类似叠加在原始图像结构上的环或带,它是第三代CT系统某 个通道在数据采集或测量过程中性能差异所差生的。如果一个通道数据差异,每 一个视角被反投影产生一条直线,在整个采集中持续存在误差产生一组到旋转中 心距离固定的直线,直线尾部被抵消而形成一个同心环状伪影;如果连续若干相 邻的通道存在差异,伪影以带状形式出现在图像中;如果通道数据只有部分视角 下的数据和测量过程中存在差异,环状或带状伪影以部分圆弧形式出现。 完全没有伪影的CT图像是不存在的,在应用过程中,对于这些影响图像质量 的伪影,就需要在设计系统时,根据检测要求考虑如何尽可能降低它们的水平, 或设计专门的处理算法和程序来减轻这些伪影的影响。本文就是要对CT图像中的 环形伪影校正算法进行研究。

2

1 绪论

1.3 CT 图像环形伪影校正研究现状
自第三代 CT 出现以来,探测器不一致性和不稳定性所致的环形伪影,一直都 是面临挑战的技术难题。多年来很多学者都对它的校正进行了深入的研究,得了 到不少有效的校正方法,都取得了一定的效果,但也都有一定的局限性。 对于探测器不一致性的校正,Jiang Hsieh[5]通过对探测器每一通道采用逐一标 定来抑制环状伪影,这种方法对探测通道少的 CT 系统是有效的,但对于探测通道 较多时,实现非常困难。Davis 和 Elliott[6]等人通过改变探测器通道积分时间来减 弱环状伪影,这种方法需要对每个通道进行控制,显然在通道数量多时,实现技 术要求高且造价贵。目前,应用比较成熟的校正方式是暗场和增益不一致校正法
[7~9]

,在探测器暗场偏置和响应灵敏度校正取得非常好的效果。但是,实际工程中, 对于探测器不一致性校正之后残留的环形伪影的校正, 也是 CT 图像伪影校正

由于探测器不稳定性,在完成探测器一致性校正后,图像中还残留环形伪影。 的重要课题之一。这方面也有很多学者进行了研究,提出了各种抑制方法。根据 目前国内外有关 CT 图像环形伪影校正方面的相关文献来看, 这些校正方法大致分 为两大类:第一类是基于 CT 图像的校正方法,它属于重建后图像处理方法。这种 方法主要通过极坐标变换将直角坐标的环形伪影变换成极坐标的线状伪影,然后 再对极坐标图像进行水平方向平滑滤波,最后再反变换到直角坐标中得到校正后 的 CT 图像[10~14], 虽然该方法对环形伪影有不错的校正效果, 但是它只是一种重建 后的处理方法,并不能从根本上解决探测器不稳定问题,而且两次坐标变换都要 用到插值处理,使得图像分辨率下降[15]。 第二类是投影正弦图的校正方法,即图像重建前的处理方法,与第一类方法 相比,校正后的图像在质量和保真度上都要优越一些。 Kowalski[16]通过对原始投 影数据进行高通滤波,从原始信号中剔除响应不一致带来的高频信号,抑制环状 伪影, 但它同时也会过滤掉一些缺陷产生的高频信号, 降低系统成像质量; Haining Sun[17]和 Matani[18]等通过对投影数据采用迭代算法来消除 CT 图像中的环形伪影, 对射束硬化造成非线性伪影也有一定的校正效果,但是迭代算法处理时间长,且 子集及子集之间的处理顺序对校正质量有一定的影响,不易选择;Raven C[19~20] 等学者提出先对投影数据进行小波或傅里叶变换,然后再对水平方向进行低通滤 波的方法来校正,该方法很好地抑制了环形伪影,但同时也给重建后的 CT 图像带 来新的环形伪影;Rainer Raupach[21]等提出了差分法校正环形伪影的构思,即采用 一定的方法提取数据中的环形伪影背景,用原始数据减去伪影背景,从而获得校 正后的投影数据;余晓锷[22~23]和孙少华[24]等利用 Hough 变换或多尺度小波分析来 查找正弦线状伪影的位置,并对该位置上的数据进行线性插值处理进行校正,该 方法伪影查找能力有限,且对复杂的投影正弦图,无法正确查找伪影位置或提取
3

重庆大学硕士学位论文

伪影特征;李保磊[25]对单个独立探测器通道存在突出相应不一致导致环形伪影, 提出了基于特征识别的方法找到投影正弦图的条状伪影位置,该方法的缺点是阈 值不易选择,校正不彻底[15]。 基于上述分析可知,对于探测器一致性校正之后残留的环形伪影的校正,从投 影正弦图的角度来对 CT 图像环形伪影进行校正是合理, 且可避免重建后校正中极 坐标变换所导致图像分辨率下降的问题。总的来说,以上所述方法对 CT 图像环形 伪影校正有一定的抑制效果,但都存在很多不足之处。

1.4 论文的研究内容与结构安排
为了提高 CT 重建图像质量以及后续处理和量化分析精度,本文对 CT 图像环 形伪影的校正方法进行了研究,提出两种去除 CT 图像环形伪影的校正算法。并通 过实验,对本文所提的环形伪影校正算法进行了比较和评价。论文的结构内容安 排如下: 第一章,简单介绍了 CT 技术的概念、发展和应用,以及本文研究的背景及意 义,总结了 CT 系统中常见的伪影类型,并阐述了 CT 图像环形伪影校正算法的研 究现状。 第二章,主要介绍 CT 系统三代扫描方式和扇束滤波反投影重建算法,分析环 形伪影产生的原因,并对环形伪影在重建图像和投影正弦图中表现形式和特征, 以及应用过程所造成的影响进行分析和总结。 第三章,详细介绍一些 CT 图像环形伪影校正的常用算法及相关理论,并分析 这些方法在环形伪影校正方面的优缺点。 第四章,重点介绍本文所提的基于投影曲线高通滤波的环形伪影校正方法, 详细介绍正弦图中伪影的识别过程,以及伪影异常数据的校正,通过实验结果验 证了该算法的有效性,并对它进行分析和总结。 第五章,重点阐述本文所提的基于改进的 Canny 算法的环形伪影校正方法, 对正弦图中条状伪影特征识别进行问题分析,提出改进 Canny 算法,并通过实验 结果验证了该算法的有效性。最后,对实验结果进行分析和总结。 第六章,总结本课题的研究成果,并分析本文所提算法存在的一些问题,对 今后的进一步研究进行展望。

4

2 CT 图像环形伪影的产生及影响

2 CT 图像环形伪影的产生及影响
环形伪影作为 CT 系统三代扫描模式的一个重要问题, 不仅降低了图像的质量, 还给图像噪声处理以及图像分割等后续处理造成困扰。在工业 CT 成像系统检测 中,这个问题更是突出,因此,CT 图像的环形伪影必须加以消除。 为了更好地对 CT 图像的环形伪影消除方法进行研究,下面对 CT 系统三代扫 描方式及重建算法进行了简单的介绍,并分析了环形伪影产生机理及分布特征。

2.1 CT 系统三代扫描方式及重建算法概述
2.1.1 投影数据和正弦图
对最原始的一代扫描模式进行分析,如图 2.1,假定 X 射线被准直成很窄的平 行射束,用一个尺寸很小的探测器测量其强度。这样在一个确定位置就可以测量 到一个投影数据,然后将射线源和射线探测器同步地沿垂直于射线的方向平行移 动,每移动一个小的固定长度,测量一个投影数,直到这些投影数据扫过整个待 测物体,形成一组从同一方向(常称为视角或分度)获得的投影数据。完成这次 直线测量以后,将射线源和探测器相对于待测物体旋转一个小的固定角度,按照 上面的步骤测量第二组投影数据,依次类推,一直旋转 180 o 或 360 o 。这样的一整 套数据通常被称为被测工件的投影数据,一般也可以表示为二维矩阵,其中行代 表探测器通道,列代表旋转角度,强度代表投影采样的幅值[1~2]。
1° 增量

第一次平移

第45次 平移

X射线源 X射线探 测器

图 2.1 最原始的扫描方式 Fig.2.1 The most original scanning mode

如果在平面直角坐标系上,选取样品工作台的旋转中心为坐标原点, x 轴表示
5

重庆大学硕士学位论文

X 射线源和探测器平行移动的位置(简称探测器位置) , y 轴表示 X 射线源和探测 器连线(即 X 射线方向)相对于原始方位转过的角度(简称旋转角度) 。 考虑到反投影重建算法的基础是认为断层图像上任何一点的值等于该平面上 经过该点的全部射线投影之和。如图 2.2 所示,对于断层上任一固定点 A( x, y) ,用 极坐标系表示则为 A(r , θ ) 。假设扫描过程中,检测样品顺时针旋转,用 β 表示旋 转角度。在初始位置 β ? 0 o ,这时射线源与探测器平行移动的位置等于 A 点在 x 轴 上的投影,即 x ? r cos θ 。当检测样品顺时针旋转一个角度 β 时, A( x, y) 点沿圆周 到达了 A?( x?, y?) 点。通过 A 点 X 射线平行移动的位置为 x? ? r cos(θ ? β ) 式中, β 是射线与 x 轴形成的投影角度。 (2.1)

x A(r,θ-β)
旋转

A(r,θ)

θ β

方向

y

O

X射线源

探测器

图 2.2 投影数据的分布图的绘制原理 Fig.2.2 The distribution figure principle of projection data

y

3 π 2

(r , θ ) ? (1, π 4)
r θ x


Δ

π
π 2

d 0
π 2

nd

?



图 2.3 点 (1, π 4) 及相应的正弦图 Fig.2.3 The (1, π 4) point and corresponding sinogram

6

2 CT 图像环形伪影的产生及影响

当被测样品旋转 360 o 时,可以看出在按照前面所述方法绘出的投影数据分布 图上,一个固定点的轨迹是一条余弦曲线。余弦曲线的幅度等于该点在极坐标下 的角度。全部被测工件的图像对应了一簇相互重叠、不同幅度、不同相位的余弦 曲线,因此投影数据分布图可以称为余弦图,但是由于传统的习惯叫法,都将其 称为正弦图[1~2]。对于给定 (r , θ ) ? (1, π 4) 对应的正弦图中的一段如图 2.3 所示。

2.1.2 CT 系统三代扫描方式概述
为了采集全套测试的投影数据,CT 设备中射线源、检测工件和探测器之间所 做的整套相对运动称为扫描。自从 CT 问世以来,早起的发展过程基本上围绕一条 提高扫描速度的主线,其方法不外乎设计不同的运动轨迹、适当增加探测器数量 或改变它们的布置方式等。 根据这些结构的变化, 逐渐形成了 CT 扫描模式中 “代” 的概念,它的发展经历了第一代和第二代的扫描模式,以及现在主流的三代扫描 模式,后来,医学上还发展了四代和五代动态空间扫描模式。 三代扫描模式在二代的基础上发展起来的,为了进一步增加探测器数量,使 射线源和线探测器阵列组成的扇形能够包括整个被检测工件的范围,只要射线源探测器组成的扇形围绕工件旋转 360 o ,射线源-探测器系统不用平行移动就能够获 得重建 CT 图像的完整数据。这样射线的利用效率可以显著提高,省去了二代扫描 中费时的多次平移往复运动,完成一次完整扫描的时间可以显著提高。这就是目 前在医学和工业应用中占有绝对统治地位的三代扫描模式,如图 2.4 所示。

X射线源

探测器

图 2.4 三代扫描模式 Fig.2.4 The third scanning mode

2.1.3 扇束滤波反投影重建算法
CT 图像的环形伪影是三代扫描方式 CT 系统的一个重要问题,目前,三代扫 描的 CT 图像一般都是采用扇束滤波反投影重建算法。 通常三代所用的线阵探测器
7

重庆大学硕士学位论文

从结构上可分为两类:等距射线的直线型探测器和等角射线的弧线型探测器,如 图 2.5 所示。 针对这两种探测器的图像重建算法分别是等距扇形束滤波反投影重建 算法和等角扇形束滤波反投影重建算法。为了对后面环形伪影的分析做准备,本 节以介绍等间距扇形束重建算法为例,对扇束滤波反投影重建算法进行简单介绍。
y S D2 S y D2

a a O x r D1 D1 O r x

(a)等角射线的直线型探测器 (a) Equi-angle ray line detector

(b) 等距射线的弧线型探测器 (b) Equi-spacial ray arc detector

图 2.5 探测器的两种形式 Fig.2.5 Two scan modes of detector

等距扇形束滤波反投影重建算法的示意图如图 2.6 所示。其中, S 0 是 CT 系统 的射线源,K 代表是实际的探测器,S0B 为经过扫描工件 A 点的射线束,它与阵 列探测器的交于 B 点,B 点到探测器阵列中心点 O ? 的距离确定了 S0B 在扇形中的 相对位置, ( s, β ) 则决定了它的绝对位置,其中 s 为该射线与虚拟探测器 L 的交点 到旋转中心 O 点的距离,虚拟探测器 L 平行于实际探测器 K,且通过旋转中心点, 故 L 也称作为中心虚拟探测器。因为 A 是中心虚拟探测器 L 和射线 S0B 的交点, 所以可以采用扇形几何关系由 B 的位置来确定 A 的位置。 等距扇形束滤波反投影算法的基本思想是[26~27]:将射线 S0B 看成平行射束中 的一条射线,即通过加权和修正将 ( s, β ) 确定的射线 S0B 变成平行射线 (0, β ) 的射 线,来对投影数据 p ? ( s) 进行校正,再按照平行束卷积反投影算法来重建图像。它 具体的原理和推导过程见参考文献,在此仅给出重建公式:
f (r , ? ) ? 1 1 ~ P [ s(r ,? )]d? 2 ? 2? U 0
?0 1 Re ? ?e i?s d? , 0 2
2?

(2.2)

其中

~ P? ( s) ?

??

??

?

D D 2 ? s? 2

P? ( s ?) g ( s ? s ?)ds ? , g ( s) ?

8

2 CT 图像环形伪影的产生及影响

s?

D ? r sin( ? ? ? ) y D ? rcos( β - θ ) ,U ? , r ? x 2 ? y 2 , ? ? arctan 。 D x D - rsin( β - θ )
S0

y L K D β θ O A
(r,θ)(x,y)

s

B O’

x

图 2.6 等距扇形束投影示意图 Fig.2.6 Projection schematic diagram of the equi-spacial ray

对公式 2.2 的等距扇形束滤波反投影重建算法进行分解,它的具体步骤如下: ① 加权和修正投影数据 p( ? , s) ,获得修正后的投影数据 p ? ( s, ? )。 D (2.3) p ? ( s, ? ) ? p ( s, ? ) D2 ? s2 ② 对修正后的投影数据用斜坡滤子滤波。 ~ p(s, ? ) ? p ? (s, ? ) * g (s) (2.4) ③ 对滤波后的投影数据加权反投影

f ( x, y) ? f (r ,? ) ? ?

2?

0

1 ~ p ( s, ? )d? U2

(2.5)

2.2 环形伪影的产生原因
环形伪影作为第三代 CT 中常见伪影噪声之一,造成它的来源因素很多[28~30], 主要包括: ① 探测器通道响应不一致性。由于探测器加工工艺的限制,或长期使用的老 化,各个探测器单元的效率和灵敏度不可能存在完全一致,导致在同等强度 X 射 线的照射下不同的阵列探测器单元输出的电流信号也存在差别。 ② 探测器通道故障。探测器某一个通道损坏,会出现严重的单环伪影。 ③ 数据采集系统故障。数据采集系统电源故障、积分板损坏、积分板接触不

9

重庆大学硕士学位论文

良等原因引起,导致采集到的投影数据与真实投影数据不一致,从而使 CT 图像中 产生环形伪影。 ④ X 射线光子不足。这种情况会导致 CT 图像出现浓淡不均的多环伪影。 ⑤ 准直器、补偿器故障及 X 线管窗口有异物等。与探测器通道故障相似,会 导致 CT 图像产生严重的环状伪影。 以上所介绍的因素中,其中数据采集系统(DAS)故障通常是造成 CT 图像环形 伪影最常见的因素来源,也最为复杂。下面从数学上分析三代扫描中环状伪影的 产生机理[31]。 CT 成像的示意图如图 2.7 所示。其中,T 为探测器,F 为 X 射线源。它们彼 此固定,扫面时共同围绕旋转中心 O 旋转。FC 是 X 射线扇束的中心线,FJ 是 FO 的延长线,它们之间的夹角为 α 。若 FC 在 FJ 的反时针方向, α 为正;反之,FC 在 FJ 的顺时针方向, α 为负。坐标系 ( x, y) 固定在被测物体上,原点为 O,即图像 坐标。探测器和球管相对于被测物体旋转了一个角度 Φ ,任取一个探测器 N,它 接收的投影为 FN,经 O 作垂直于 FN 的直线,它与 FN 的交点为 P。
y
J
v C u N
T

D

Φ o θ S1 P α β F S2 x

图 2.7 CT 系统成像示意图 Fig.2.7 CT system imaging schematic diagram

P 点的坐标是

X ? ? f sin( β ? α) cos( β ? α ? Φ) Y ? ? f sin( β ? α) sin( β ? α ? Φ)

(2.6) (2.7)

其中 f 是 FO 的长度, β 是 FN 与 FC 的夹角,如果 FN 在 FC 的反时针方向,

10

2 CT 图像环形伪影的产生及影响

β 为正,反之, β 为负。当 N 选定后,式中除 Φ 外,其他各量均为常数。式 2.6
和 2.7 式是以 Φ 为参变量的圆,其圆心为 O,半径为 OP ? f sin( β ? α) 。即系统扫 描时,P 的轨迹是一个圆。设此圆半径是 r ,得: r ? f sin( β ? α) (2.8)

现在分析 FN 对图像上不同半径的圆环的贡献。在扫描野内作一个圆环,如图 2.8,它由半径为 R 和 ? ΔR ? R ? 的两个同心圆组成,圆心为 O。这里 r<R<S2。当 ΔR 很小时,FN 在此圆环单位面积上经过的像素数目 δ 是 d δ? π R2 ? r 2 (2.9)

N

S2
R O r ΔR

P

图 2.8 FN 对图像上不同半径的圆环的贡献 Fig.2.8 FN contribution of different radius in the image

π

20 15 10 5

rδ d

F

0

1

2

3

4

5

R r

图 2.9 R 圆环上单位面积内有多少像素函数图像 Fig.2.9 The pixel function image number per unit area on the R ring

式中 d 为投影在单位长度上经过的像素数目,当 FN 选定后,d 是常数。上式 表示图像的各个 R 圆环上单位面积内有多少像素接受 FN 的贡献。其变化曲线如

11

重庆大学硕士学位论文

图 2.9 所示。 由此可见:当 R≤r,δ 无意义,意味着 FN 对半径小于 r 的圆环没有贡献,从 物理上讲, FN 永远不经过半径小于 r 的区域, 当然对图像的这一部分不会有贡献; 当 R 接近 r 时,δ 很大,随着 R 上升,δ 急剧下降,这说明 FN 对 P 附近环的贡献 比远离 P 的大得多;扫描时,系统旋转一周,FN 的贡献主要在以 r 为半径,以 O 为圆心的圆上。 通过以上分析,可以看出探测器每个通道对图像的贡献主要是一个圆。如果 其中的一个探测器通道出现了异常, 那么该通道对应的 CT 图像上的圆的像素数据 就会异常,从而产生以重建中心为的圆心的环形伪影。

2.3 环形伪影特征及分析

(a)

(b) 图 2.10 环形伪影在重建图像和投影正弦图中的表现形式

(c)

Fig.2.10 The forms of ring artifacts in reconstruction image and sinogram

重建图像中环形伪影的表现形式与探测器系统异常通道的个数与位置息息相 关。异常探测器通道的个数影响着图像中环形伪影的宽度,它的位置影响着环形 伪影在重建图像中的位置,它们都是一一对应的关系。根据 2.2 节环形伪影产生原

12

2 CT 图像环形伪影的产生及影响

因分析可知,环形伪影在 CT 图像中表现为以重建中心为圆心的同心圆;由于每个 单独探测器通道采集的数据在所有投影角度与中心保持对称,且它们的相对位置 是固定的,所以,环形伪影在投影正弦图中表现为竖直方向上的直线。 如果 CT 系统中有一个探测器通道损坏, 那么对应的重建图像中则存在一个明 显的单环伪影,而正弦图中则出现一根沿着竖直方向的直线,如图 2.10(a);如果 探测器通道响应不一致性会导致全幅 CT 图像产生间隔均匀的环形伪影,如图(b); 如果是数据采集系统积分板电路参数的变化,或者由于射线硬化、强度不一致和 射线通过扫描工件的路径不同,而造成伪影的严重程度也不同,使得图像中伪影 表现为多个不均匀的弧或完整的环,如图(c)。实际应用中,最常见的还是不同亮 度的弧形或环形伪影。 除了以上特征外,通常认为环形伪影还具有以下特性: ① 具有一定的宽度。CT 系统中有一个探测器通道异常,对应着重建图像中 就会出现一个像素宽度的弧或环,那么它的正弦图中也就对应着一个像素的竖线。 而实际工程中,由于工业 CT 中使用的探测器价格昂贵,通常采用高采样频率扫描 技术[32~34]来实现高分辨率 CT 成像。 那么一个异常探测器所致伪影的宽度与高频采 样的插值次数相等,本文将其看作为单个的环形伪影来进行处理。而且如果是多 个连续的探测器异常,那么则在正弦图和重建图像中表现为多个连续的条带和环 带。 ② 同一伪影灰度不均匀。对于数据采集系统积分板电路参数的变化所致的环 形伪影,由于射线通过扫描工件的路径长度不同,伪影的严重程度也不同,从而 导致图像中环形伪影不是均匀的全环,而可能是弧或亮暗不均的环;同理,投影 正弦图中则表现为间断的或者灰度不均匀的竖条。 ③ 可能具有一定的周期性。为了提高探测器通道响应的一致性,通常 CT 系 统将多个探测器组成探测器模块。Klaus J.Engel[35]在对探测器串扰分析时,提出探 测器模块边界与模块内部差异,信号(X 射线和光学信号)在模块边界更容易被 吸收,使得探测器模块边界通道响应不一致,从而导致重建图像上产生了等间距 的单个分散环形伪影,投影正弦图则产生等间距的直线。 从以上特定分析中,可以看出,不同的 CT 图像、甚至是同一幅 CT 图像或同 一个环形伪影,环形伪影的表现都可能不同。在校正时,根据伪影的特征和表现 形式,可以采用不同的校正方法。

2.4 环形伪影的影响
根据伪影的名称和定义可知,只要 CT 图像上出现了伪影,那么它与物体的真 实物理结构就存在差异。 环形伪影作为 CT 图像中最常见伪影之一, 也同样是如此。
13

重庆大学硕士学位论文

它会给图像带来明显的特征改变,从而影响到医学疾病诊断和工业材料检测。 医学上通过 CT 成像获得人体组织和器官的图像,来对病人进行诊断和治疗。 环形伪影的存在,影响了医学 CT 图像的质量,在医学应用中影响非常大。例如环 形伪影可能与人体的某些病变组织图像相似,例如,穿过大动脉血管的黑色弧形 与大动脉破裂非常相似,环形伪影的一个特殊例子是中央模糊,它属于较小范围 的环形伪影,对于这种伪影要格外重视,因为它与病理组织图像表现非常类似, 可能导致误诊[2,23]。 对于工业 CT 领域的应用,环形伪影也同样存在极大地危害。例如,对钢铁铸 件进行检测时,环形伪影就可能影响铸件中裂纹、气泡的识别。对于其它方面应 用,航空航天关键材料结构和属性的检测,港口和海关的旅客行李的检测,军工 器械关键工件的检测和火工品密度测量,这些都会因为图像中环形伪影的干扰, 使得图像对比度降低或工件细节信息被覆盖或模糊,图的可读性变差,影响了这 些工件感兴趣区域的正确识别,以及材料属性和尺寸的精确测量。

2.5 本章小结
本节简单介绍了 CT 系统三代扫描方式和扇束滤波反投影重建算法, 分析了环 形伪影产生的原因,并对环形伪影在重建图像和投影正弦图中表现形式和特征, 以及应用过程所造成的影响进行了分析和总结。

14

3 环形伪影校正的常用方法

3 环形伪影校正的常用方法
自第三代 CT 出现以来,环形伪影一直都是面临挑战的技术难题,很多学者都 对它的消除方法进行了研究。对 CT 图像中环形伪影的校正,按照投影数据获取流 程,总体可分为三大类:数据采集校正、重建图像校正和投影正弦图校正。其中, 数据采集和投影正弦图的校正方法都属于图像重建前的校正方法。数据采集校正 法一般是根据 CT 系统的物理原理来对探测器的不一致性进行校正, 目前应用比较 成熟的是暗场和增益不一致校正法;重建图像校正法一般是将笛卡尔坐标下的环 形伪影变换极坐标下的直线伪影进行滤波,以达到去除环形伪影的目的;投影正 弦图校正法属于数据采集和重建图像校正的中间过程, 根据 CT 图像的环形伪影在 投影正弦图中表现为竖直方向直线的特征,采用 Hough 变换滤波法进行校正,或 通过傅里叶变换及小波变换将投影正弦图到频域,再对高频谱的伪影信息进行滤 波处理。

3.1 改变扫描方式校正法
为了去除CT系统中探测器不一致所导致环形伪影, 傅健[36]从改变CT系统扫面 方式的角度,提出了一类可以抑制探测器响应不一致,从而能消除环状伪影的扫 描采样方式。 文献中提出的校正算法的基本思想是:具有响应不一致性的探测器通道之所 以会影响重构图像,产生环状伪影,关键是其在扫描过程中始终采集同一根射线 受被检物质调制后携带的信号,这样反投影时,该探测器通道具有的响应不一致 性通过这些信号在重构矩阵固定位置形成以这根射线为切线,以旋转中心为中心 的环,所以,如果能在扫描过程当中,采用一定的扫描方式,使具有响应不一致 性的通道能采集不同射线携带的信号,在原始投影采集环节即将探测器系统响应 不一致性平均分布到整个采样序列,避免其集中在投影序列的固定位置,那么环 形伪影可以得到抑制。 作者采用了如图3.1所示扫描方式,并进行了计算机模拟,来验证这几种扫描 方式对环形伪影的减弱效果。通过仿真实验,作者发现图3.1(b)扫描方式抑制环状 伪影没有效果,它只是将探测器响应不一致引起的环状伪影变换了一种形式,这 是由于采用的算法不同产生的,图(c)扫描方式能有效抑制环形伪影。图(a)、(b)所 示的采样扫描方式之所以不能抑制环状伪影,就是因为它们并没有改变探测器通 道与射线的对应关系,在整个扫描采样过程中,每个探测器通道只采集同一根射 线调制的信号,而图(c)所示扫描采样方式正是在扫描采样过程中改变了探测器通

15

重庆大学硕士学位论文

道与射线的对应关系,所以可以抑制环形伪影。同理,图(d)方式亦能实现环形伪 影的抑制。

y Source

y 1 Source

Object 2 Rotation center Detector array x 3

Object x

Rotation center Detector array

(a)物体旋转 (a) Only objects rotate
y Source

(b)物体旋转与平移 (b) Objects rotate and transverse simultaneously
y Source

Object x

Object x

Rotation center Detector array

Rotation center Detector array

(c) 物体旋转,与射线源同步平移 (d) 射线源平移 (c) Objects rotate and transverse with the source simultaneously (d) The source transverses 图 3.1 各种扫描采样方式 Fig.3.1 Various scanning sampling method

根据上述结论,任何一种扫描方式只有在扫描采样过程中改变探测器通道与 射线之间的对应关系,才能抑制图像中的环形伪影。这种方法能相对减弱CT图像 中的环形伪影,但是,实现复杂,且要以增加硬件成本为代价。

3.2 暗场和增益不一致校正法
在探测器制造过程中,受加工工艺限制,不可能制造出完全一致的探测单元。 例如,当探测单元在不同批切割时,它们的尺寸可能略微变化;探测单元之间的 表面粗糙度也会因为加工过程微小改变而有变化;探测单元之间的反射材料不能

16

3 环形伪影校正的常用方法

保证完全一致;链接每个闪烁体单元的光电二极管的谱响应和转换效率可能有所 差别;不同通道之间的数据采集电子器件存在固有的差异,这些因素都导致了探 测器通道响应增益的不一致,从而导致重建图像中产生全环伪影。目前对探测器 通道响应不一致的校正主要采用暗场和增益不一致校正法[7~9]。 在理想情况下,假设每个探测器像元通道的响应度与射束入射到探测器上的 强度成线性关系,它的响应模型为: I 0 ( x) ? 零输入响应B( x) ? 零状态响应R( x) (3.1)

各个探测器像元的 B( x) 和 R( x) 的不同造成了探测器通道对射束的响应不一 致,也导致了探测器通道响应的不一致。所以探测器响应不一致的本质是探测器 响应曲线的斜率与截距不同,对系统来说就是零点与增益的不同,因此必须对探 测器系统进行暗场校正和增益不一致校正。 ① 暗场校正 由于每个探测器通道的 B( x) 都不一样,使得重建图像不能等价于输入的射线 场强,很明显,零输入响应是探测器的一种固有特性。通常也将零输入响应校正 称作为偏置校正。 偏置校正主要是去除探测器系统中暗电流不一致性,它包括放大器相关的DC 分量和暗电流,由于与输入无关,因此,它的校正公式为 I B ( xi ) ? I r ( xi ) ? B( xi ) 束照射条件下,采集的多幅叠加平均的 CT 图像。 ② 增益不一致校正 零状态响应即探测器在工作参数已设定的条件下,当输入响应 D( x) ? 0 时,它 的射线输入响应情况。 虽然探测器输出为 R( x) , 但是零输入响应导致了它产生 B( x) 的偏置,降低了探测器系统的动态范围。所以,由式 3.1 可得 (3.2)

其中, I r ( xi ) 为原始图像; B( xi ) 为暗场基准图像,即在相同曝光时间和无射

R ( x ) ? I 0 ( x ) ? B( x )
测器,在被光敏元件采样之前对输入射线场强的零状态响应为 R( x) ? I i ( x) ? h( x) ? ? I i ( x ? α)h(α)dα
??

(3.3)

假设成像系统的点扩散函数(Point-spread-function, PSF)为 h( x) , 则射线到达探 (3.4)

经模拟和实验验证,此时的 PSF 近似为高斯函数

? ( x ? α) 2 ? h( x ? α) ? A exp ? ln 2? 2 ? ?BWd 2? ?
直径, M 为放大比。

(3.5)

式中 BWd ? d 2 ? [a(M ? 1)]2 为 PSF 的半波宽; d 为像元尺寸,a 为射线焦斑 以幅值为 S ?( x) ,直径为 BWd 的小圆柱,对高斯函数场强等效,则是 3.4 简化
17

重庆大学硕士学位论文

为:

π R( x) ? I i ( x) ? h( x) ? I i ( x) ? ( BWd2 ) S ?( x) ? I i ( x) ? S ( x) 4
其中, S ( x) 为探测器的灵敏度。 根据以上分析,射线输入场强的估计公式为 ? ( x) ? R( x) S ( x) ? ( I ( x) ? B( x)) S ( x) I i 0 经灵敏度校正后的图像可由下式得到 I ( x ) ? B ( xi ) I g ( xi ) ? K ? 0 i G ( xi ) ? B ( xi ) 式中 G( xi ) 为灵敏度基准图像; K 为 G( xi ) 经偏置校正后的归一化因子。

(3.6)

(3.7)

(3.8)

综上分析,采用暗场和增益不一致校正法,对探测器不一致性所导致的环形 伪影有很好的抑制效果。但是,该校正算法的模型是建立在 CT 系统探测器响应在 理想的状态下所提出的线性模型,而实际应用中,由于射线穿透被测工件路径不 同,探测器所采集的射线强度也不同,通常认为探测器对不同能量的射线响应是 非线性的。所以,在对探测器系统所采集的投影数据进行暗场和增益不一致校正 后,通常 CT 图像中还会残留有环形伪影。

3.3 投影正弦图的校正
3.3.1 Hough 变换插值校正法
根据重建图像中环形伪影在投影正弦图中表现为竖直方向直线的特征,作者 依据Hough变换对线检测的优点,提出基于Hough变换插值校正的方法。该算法的 基本思想是[22~23]:利用Sobel检测器对边缘数据进行标记,结合Hough变换,找到 投影正弦图中直线的位置。最后,采用线性插值对检测出的伪影处投影数据进行 校正。它的伪影定位过程如下: ① 对环形伪影投影数据进行标记
-1 0 1 -2 0 2 -1 0 1 -1 -2 -1 0 0 0 1 2 1

Gx ? ( z7 ? 2 z8 ? z9 ) ? ( z1 ? 2 z2 ? z3 )

Gy ? ( z3 ? 2 z6 ? z9 ) ? ( z1 ? 2 z4 ? z7 )

图 3.2 Sobel 掩模 Fig.3.2 Sobel template

18

3 环形伪影校正的常用方法

根据重建图像中的环形伪影在投影正弦图中表现为竖直方向直线的特征,利 用含有Sobel检测器的边缘检测函数edge对边缘数据进行标记。采用如图3.2所示的 Sobel算子来求一阶导数 G x 和 G y ,那么中心邻域梯度计算公式为:
2 2 g ? Gx ? Gy

?

?

12

? [( z 7 ? 2 z8 ? z9 ) ? ( z1 ? 2 z 2 ? z 3 )]2 ? ?( z3 ? 2 z 6 ? z9 ) ? ( z1 ? 2 z 4 ? z 7 )?

?

2 12

?

(3.9) ② 利用Hough函数定位伪影数据 Hough变换时对图像进行某种形式的坐标变换, 现在已成为图像中识别几何形 状的基本方法之一。它的基本思想是将原始图像中给定形状的曲线或直线换成变 换空间中的一个点,及原始图像中给定形状的曲线和直线上的所有点都集中到变 换空间的某个点上形成峰点,使得原始图像中给定形状曲线或直线的检测问题, 变成寻找变换空间中的峰点问题,即把检测整体特性(给定曲线的点集)变成检测 局部特性的问题。由于它是抗干扰性强,对曲线的检测能力基本不受中间断点的 干扰,不用提前对边缘点进行连接,所以常常用于灰度图像中的直线检测。 设在原始图像空间坐标 ( x, y) , 采用极坐标 ( ρ, θ ) 作为变换空间, 其极坐标方程 可写成:

? ? x cos ? ? y sin ?
法线和 x 轴的夹角。

(3.10)

参数 ? 和 ? 可以唯一确定一条直线, ? 表示原点到直线的距离,? 是该直线的

y

θ π l ρ θ O 图像平面 x O Hough 空间 ρ

P

图 3.3 Hough 变换的对偶性原理 Fig.3.3 Dual principle of Hough transformation

对于 ( x, y) 空间中的任一点 ( xi , yi ) 采用极坐标 ( ρ, θ ) 作为变换空间,对应于极 坐标空间中的一条正弦曲线,其初始角和幅值随 x i 和 y i 的值而变。将 ( x, y) 坐标空 间中在同一条直线上的一个点序列,变换到 ( ρ, θ ) 空间,则所有正弦曲线都经过一 点 ( ρ0 , θ0 ) , ? 0 为这条直线到原点的距离,? 0 为法线与 x 轴的夹角,如图 3.3 所示。

19

重庆大学硕士学位论文

通常将以 ( ρ, θ ) 坐标平面称为 Hough 空间,图像平面到 Hough 空间的变换称 为 Hough 变换。图像平面的直线和 Hough 空间的点是一一对应的,如果已知图像 平面上一条直线的参数,那么它对应的 Hough 空间坐标也就可以确定。

(a)Lena 图像 (a) Lena image

(b)Sobel 算子二值图像 (b) Binary image of Sobel detection

(c)Hough变换图像 (c)Hough transform image

(d)检测图像 (d) Detection image

图 3.4 图像的 Hough 变换 Fig.3.4 Hough transform image

图 3.4 表示是 Hough 变换检测过程和结果图像,图(a)为含有条状伪影的 Lena 图像,图(b)为对图(a)采用 Sobel 算子获得的二值图像,然后对(b)进行 Hough 变换 得到图像(c),从检测结果的图像(d)中可以看出,虽然检测了背景中的竖直伪影, 但是检测不完全。这是因为 Hough 变换进行线检测,抗间断能力好,但是存在检 测不精确, 多峰值检测和受噪声干扰严重的问题。 采用 Sobel 算子进行边缘检测时,
20

3 环形伪影校正的常用方法

检测背景的条状伪影边界,但是表现一定的波动情况,最终表现为伪影识别不全, 从而导致重建图像环形伪影校正不彻底,这就是本算法在进行投影正弦图中伪影 检测的缺点所在。

3.3.2 小波-傅里叶变换滤波法
傅里叶变换和小波变换都是从频域对投影数据进行处理, 从而达到去除 CT 图 像环形伪影的目的。它们的基本思想[19~20]都是通过将空间域的投影正弦图,变换 到频域进行平滑滤波,再平滑后图像反变换到空间域,来获得校正后的 CT 图像。 ① 傅里叶变换 傅里叶变换的基本思想是:时域内一般的函数或者信号都可以看作是一种复 杂的波动,而任何一种复杂的波动现象又都是由不同幅度、不同频率的正弦波叠 加而成的[37]。 一个图像尺寸为 M ? N 的函数 f ( x, y) 的离散傅里叶变换由以下等式给出: 1 M ?1 N ?1 (3.11) F (u, v) ? f ( x, y)e ? j 2π (ux M ?vy N ) ?? MN x ?0 y ?0 式中必须对 u 值( u ? 0,1,2,?, M ? 1)和 v 值( v ? 0,1,2,?, N ? 1)计算。同样,给出

F (u, v) ,可以通过傅里叶反变换获得 f ( x, y) ,由以下表达式给出:
f ( x, y ) ? ?? F (u, v)e j 2π (ux M ?vy
u ?0 v ?0 M ?1 N ?1 N)

(3.12)

式中, x ? 0,1,2,?, M ? 1, y ? 0,1,2,?, N ? 1 。式 3.11 和式 3.12 构成了二维离 散傅里叶变换对(DFT)。变量 u 和 v 是变换或频率变量, x 和 y 是空间或图像变量。 由于傅里叶变换后的图像显示的是复数形式的频谱信息,所以,为了方便具 体分析,通常采用幅度谱 F (u, v) 和相位谱 φ(u, v) 来对它们进行分析。它们的计算 公式为
F (u, v) ? F (u, v) exp[ jφ(u, v)] ? R(u, v) ? jI (u, v)
? I (u, v) ? φ(u, v) ? arctan ? ? ? R(u, v) ?

(3.13) (3.14)

F (u, v) ? R 2 (u, v) ? I 2 (u, v)

?

?

1 2

(3.15)

实际应用中,为了更好地分析傅里叶变换后的各频率分量,通常对幅度谱采 用取对数的方法来实现傅里叶变换结果的可视化显示。 如图 3.5 所示, 图(a)是傅里 叶变换前含有伪影的 lena 图像,图(b)为它的傅里叶变换后的对数幅度图像。 傅立叶变换将时间域与频率域联系了起来,在变换后的频谱中,一般信号都 表现为低频信号,而图像的细节和伪影噪声则表现为高频信号,通过滤掉频谱中 的高频成分,来达到去除 CT 图像中环形伪影的目的。但是,傅里叶变换无法对局 部时域信号的频率特征进行分析,所以在滤波的时候可能也滤掉了高频谱的图像

21

重庆大学硕士学位论文

细节,造成图像模糊。

(a) 图 3.5 图像的傅里叶变换 Fig.3.5 Fourier transform image

(b)

② 小波变换 小波变换在二维图像的应用,只考虑尺度函数是可分离的情况,即 Φ( x, y) ? Φ( x)Φ( y) (3.16)

式中, ?( x) 为一维尺度函数,它的基小波是? ( x) 。二维基小波的公式为: ? 1 ( x, y) ? ?( x)? ( y),? 2 ( x, y) ? ?( y)? ( x),? 3 ( x, y) ? ? ( x)? ( y) (3.17) 它构成二维平方可积函数空间 L2 ( R 2 ) 的正交归为一基: ? 1j ,m,n ( x, y) ? 2 j? 1 ( x ? 2 j m, y ? 2 j n) j ? 0, l ? 1,2,3 , j, l , m, n 都为整数 (3.18) 1)正变换 从一幅图像 f ( x, y) 大小为 N ? N ,其中 N 是 2 的幂数。对于 j ? 0 ,尺度,也 就是原始图像的尺度。 j 值的每一次增大都使尺度加倍,而使分辨率减半。在变换 的每一层次,图像都被分解为 4 个四分之一大小的图像,它们都是由原始图像与 一个小波基图像的内积后,再经过行和列方向进行 2 倍的间隔抽样而生成的。对
(j ? 1 ) 于第一个层次 ,可写成 f 20 (m, n) ? f1 ( x, y ), ?( x ? 2m, y ? 2n)
f 21 (m, n) ? f1 ( x, y ),? 1 ( x ? 2m, y ? 2n) f 22 (m, n) ? f1 ( x, y ),? 2 ( x ? 2m, y ? 2n) f 23 (m, n) ? f1 ( x, y ),? 3 ( x ? 2m, y ? 2n)

(3.19)

(j ? 1 ) 后续的层次 ,依次类推,形成如图 3.6 所示的形式。

22

3 环形伪影校正的常用方法

f(x ,y)
(a)

(b)

(c)

(d)

(a)原始图像

(b)第一层

(c)第二层

(d)第三层

(a) Original image (b) First layer (c) Second layer (d) Third layer 图 3.6 二维离散小波变换 Fig.3.6 2-D discrete wavelet transform

如果将内积改写卷积形式,可得 f 20j ?1 (m, n) ? f 20 1 ( x, y ) * ? ( ? x, ? y ) ( 2m, 2m)

f 21j ?1 f 22j ?1 f 23j ?1

? (m, n) ? ? f (m, n) ? ? f (m, n) ? ? f

0 21 0 21 0 21

? ( x, y ) *? (? x,? y )?(2m,2m) ( x, y ) *? (? x,? y )?(2m,2m) ( x, y ) *? (? x,? y )?(2m,2m)

(3.20)

因为尺度函数和小波函数都是可分离的,所以每个卷积都可以分解成行和列 的一维卷积。 如在第一层, 首先用 h0 (? x) 和 h1 (? x) 分别于图像 f ( x, y) 的每行作卷积 并丢弃奇数列;接着这个 ( N ? N ) / 2 矩阵的每列再和 h0 (? x) 、 h1 (? x) 作卷积,丢弃 奇数行;结果就是该层变换所要求的 4 个 ( N / 2) ? ( N / 2) 的数组,如图 3.7 所示。

h0 (? x)

h0 (? x)

f 20j ( x, y ) f 21j ( x, y ) f 22j ( x, y ) f 23j ( x, y )


f

0 2 j ?1

( x, y )

h1 (? x)
h0 (? x)

h1 (? x)
行 列

h1 (? x)


图 3.7 DWT 图像分解步骤 Fig.3.7 DWT image decomposition steps

23

重庆大学硕士学位论文

2)逆变换 与正变换相似,通过在每一层的每一列的左边插入一列零来增频采样前一层 的 4 个矩阵;接着用 h0 ( x) 和 h1 ( x) 来卷积各行,再成对地把这几个 ( N / 2) ? N 的矩 阵加起来;然后通过在每行上面插入一行 0 来将刚才所得的两个矩阵的增频采样 为 N ? N ;再用 h0 ( x) 和 h1 ( x) 与这两个矩阵的每列卷积。这两个矩阵的和就是这一 层重建的结果。逆小波变换图像重建的过程如下图 3.8 所示。
行 列

f 20j ( x, y ) f 21j ( x, y ) f 22j ( x, y ) f 23j ( x, y )

h0 (? x)
+





h0 (? x)
+
X4

h1 (? x)

f 20j?1 ( x, y )

h0 (? x)
+

h1 (? x)

h1 (? x)
图 3.8 DWT 图像分解步骤 Fig.3.8 DWT image decomposition steps

根据傅里叶变换的变换滤波时存在的问题,以及小波分解的优点,文献[19] 中提出小波—傅里叶变换滤波的方法来对 CT 图像的环形伪影, 弥补了傅里叶变换 校正的不足。






(a)

(b) 图 3.9 小波—傅里叶分解示意图 Fig.3.9 Wavelet-FFT decomposition schematic diagram

(c)

24

3 环形伪影校正的常用方法

图 3.9 中的图(b)表示对图(a)(含有伪影的 Lena 图像)进行 3 次小波分解,每 一次小波分解都把原图分解为①垂直细节信息、②水平细节信息、③对角线细节 信息和④平均细节信息等 4 个部分。从分解后的图中可以看出,对投影正弦图进 行小波分解后,竖直方向的伪影信息只分布在垂直伪影信息中,所以只要对垂直 细节部分进行处理,就可以去除图像中的伪影;然后,对小波分解后的垂直细节 信息进行傅里叶变换,变换后伪影的高频信息分布在如图(c)的白色线框所示,并 采用高斯滤波对傅里叶变换后的数据图像进行处理,以去除 CT 图像中的环形伪 影;最后,傅里叶逆变换和小波重构获得校正的 CT 图像。

3.4 重建图像的校正
根据环形伪影在 CT 图像中表现为以重建图像中心为圆心的圆环, 考虑在该圆 环内外选取适当宽度的环形区域,对该圆环数据进行平滑滤波,来消除 CT 图像的 环形伪影。很明显,这种方法不容易确定圆环的位置,且平滑过程的运算量也非 常大。 针对这个问题,很多学者都提出了极坐标变换的思想来对环形伪影进行校正, 他们的基本思路是将 CT 图像从笛卡尔坐标系 f ( x, y) 转换到极坐标系 p(r , θ ) ,基 于环状伪影特点,不难发现,若将直角坐标系中的 CT 图像 f ( x, y) 转换为极坐标 系中图像 p(r , θ ) ,环状伪影将会由同心圆环转换为竖直方向的线状伪影,此时的 特征类似于投影正弦图,因此采用类似上述投影正弦图的处理方法,去掉极坐标 图像 p(r , θ ) 上竖直的线状伪影,得到新的图像 p?(r , θ ) ,然后再将极坐标 p?(r , θ ) 转 换为笛卡尔坐标 f ?( x, y) ,获得校正后的 CT 图像[10~14]。 将笛卡尔坐标下的 f ( x, y) 变换到极坐标 p(r , θ ) ,它的变换公式为: y θ ? arctan( ) , r ? x 2 ? y 2 x 其中, p?(r , θ ) 变换成 f ( x, y) ,所用公式为: x ? r cos Φ , y ? r sin Φ

(3.21)

(3.22)

由于极坐标变换是以重建图像的中心为原点,所以重建图像的环形伪影在极 坐标变换后, 表现为沿着 θ 轴的条形伪影。 图 3.10 为一含有环形伪影 CT 图像的极 坐标变换图像。 对于极坐标下 p(r , θ ) 中伪影的校正,可以采用形如 H (W ) ? (1 W )ones(W ,1) 的 平滑模板,其中 W ? 2,3,4 ?,式中 ones (W ,1) 为一个 1 列 W 行且元素全部为 1 的矩 阵。它的优点是只在沿着一个方向上对图像进行平滑,而图像的另一个方向保持 不变。但是,环形伪影的不均匀性决定了它在极坐标下的条形伪影的强度也是不 均匀的,所以希望在伪影强度大的区域,适当的多进行平滑;而在伪影强度小的

25

重庆大学硕士学位论文

区域,尽量少平滑,甚至不平滑,以尽量保持原图像的细节信息。在这方面,学 者们在文献中都提出了各种方法,且对环形伪影有一定的抑制效果。

(a) 图 3.10 极坐标变换

(b)

Fig.3.10 Polar coordinates transformation

从以上校正算法原理可以看出,由于该算法将笛卡尔坐标下图像变换到极坐 标图像,伪影校正后再反变换的直角坐标,两次用到图像插值处理,必然会导致 图像分辨率下降。而且采用的滤波模板 W ,也局限了它只在对单个分散或少数几 个相连的环形伪影校正有效。

3.5 本章小结
本节对常用的环形伪影校正方法进行介绍,详细介绍了它们的校正原理,并 对它们的校正效果以及存在的问题与不足进行分析。

26

4 基于投影曲线高通滤波的环形伪影校正

4 基于投影曲线高通滤波的环形伪影校正
如果 CT 系统中一个探测器通道异常,那么对应投影数据矩阵中一列数据异 常,投影正弦图对应的列像素出现竖线。所以若想去除 CT 图像中环形伪影,首先 必须识别到投影正弦图中竖线的位置,即数据采集过程中异常探测器的位置,然 后利用周围的投影数据对竖线处的异常数据进行校正。 对一幅图像来讲,由于空间相关性的关系,其相邻位置具有较强的相关性。 相应地,图像沿着列方向(时间轴方向)的投影应具有近似连续的灰度分布,若把 图像的列投影数据序列看成一维信号,则此信号曲线应是近似光滑且连续的。正 因为环形伪影在投影正弦图中沿着竖直方向存在,改变了这些位置列的灰度真实 值,如果沿着列方向进行投影积分,那么积分曲线就会呈现出噪声干扰式的特点。 同理,将二维空间线识别问题转化为一维空间点特征查找,使问题得到简化。

4.1 基于投影曲线高通滤波的条状伪影特征识别
4.1.1 伪影增强
CT 系统在进行成像扫描时, 如果数据采集系统 A/D 转换位数不够或采集的数 据中噪声过大,都会使投影数据的动态范围变小,从而影响了重建图像的对比度, 这个问题在工业 CT 中更是突出。如果投影数据的动态范围太小,很明显在投影正 弦图中表现为图像偏暗,条状伪影在投影正弦图中所表现的特征都不明显和清晰, 所以从投影正弦图对 CT 图像的环形伪影进行校正时, 首先考虑采用图像增强处理 方法来增强正弦图中的伪影信息,以提高伪影与工件弦线的对比度。 目前,数字图像处理中对图像增强的方法有很多,主要分为空间域和频率域 的处理方法,但是这些方法大多针对图像局部或整体轮廓信息的增强,并没有满 足增强本文竖直方向上伪影边界的要求。 根据 CT 图像卷积反投影重建原理以及其 重建滤波器的优点,本文采用重建滤波器来对投影数据进行滤波,以增强正弦图 中的竖直方向上的伪影信息[15]。 卷积反投影重建算法中,理论上要求滤波器的系统函数 H( ρ) ? ρ ,这是一个 频带无限的理想滤波函数。按照佩利-维纳准则,这一理想函数是不能实现的。实 际情况中,频域中的高频分量幅度极小,故认为 H( ρ) 为带限的。为此,令
H ( ρ) ? ρ W ( ρ)

(4.1) (4.2)

式中, W ( ρ) 为带宽为 1 (2d ) 的窗函数, d 为采样间隔,其表达式为
W ( ρ) ? 0, ρ ? 1 (2d )

而 ρ ? 1 (2d ) 时, W ( ρ) 的值是要讨论的内容。最简单的办法是采用矩形窗,
27

重庆大学硕士学位论文

此时
?1, ? W ( ρ) ? ? ?0, ? ρ ? 1 ( 2d ) ρ ? 1 ( 2d )

(4.3)

如不顾投影数据采集时的离散性和噪声的影响,也不计 Gibbs 现象,按式 4.3 处理属合理,而且可以获得很高的分辨率。但实际应用过程中,要考虑离散、噪 声和 Gibbs 现象等影响, H ( ρ) 的选择就需仔细考虑。现在,已经研究出了各种窗 函数的滤波器,常见的滤波器有 R-L 滤波函数和 S-L 滤波函数[26]。 ① R-L 滤波函数 该函数的基本出发点是认为实际的二维图像总有一个上限频率 1 (2d ) , 若以宽 度为 ?? 1 (2d ),1 (2d )? 的矩形窗将函数截断,则滤波函数的频率响应变为 ?ρ ρ ? 1( 2d ) ? H R ? L ( ρ) ? ? ?0 ρ ? 1 ( 2d ) ? 冲激响应为
? sin 2πsB ? 2 ? sin πsB ? hR ? L ( s) ? 2 B 2 ? ??B ? ? ? 2πsB ? ? πsB ?
HR?L(ρ)

(4.4)

(4.5)

hR ? L ( s )
B2

1 4d 2

hR ? L ( s)

1 ? 2d

0

1 2d

ρ

?

2 1 4 3 ? ? ? 2B 2B 2B 2B

1 2 3 2 B 2B 2 B

4 2B
S

?8d ?6d ?4d ?2d

?7d ?5d ?3d ?d

0

d 3d 5d 7d 0

2d 4d 6d 8d

s

(a)频率特性 (a) Frequency characteristic

(b)连续空域特性 (b) Spatial characteristic

(c)离散内插后特性 (c) Discrete characteristic after interpolated

图 4.1 R-L 滤波函数 Fig.4.1 R-L filter

H R ? L ( ρ) 与 hR? L ( s) 的图形见图 4.1(a)和(b)。若探测器的采样间隔为 d,对应的
最高不失真频率 B ? 1 (2d ) ,以 s ? nd 代入式,可得离散化的冲激响应为 ? 1 n?0 ? 4d 2 ? ? hR ? L (nd ) ? ?0 n ? 偶数 ? ? 1 ?? 2 2 2 n ? 奇数 ? nπ d
28

(4.6)

4 基于投影曲线高通滤波的环形伪影校正

R-L 滤波函数形式简单、实用,滤波后图像轮廓清楚。图 4.1(c)为 R-L 滤波器 离散化冲激响应的示意图,从图中可以看出,用它重建图像,轮廓清楚,缺点是 具有 Gibbs 现象,表现为明显的振荡响应。另外若投影数据含有噪声,重建质量也 不够满意,其原因在于当 ? ? 1 /( 2d ) 时, H ( ? ) ? ? 。理想的矩形窗是产生 Gibbs 现象的根源所在。 ② S-L 滤波函数 为了缓解振荡响应,并更好地补偿 ρ ? 1 (2d ) 处的混迭,S-L 滤波函数设法使
ρ ? 1 (2d ) 处的 H ( ρ) 幅值压低。于是,它的冲击响应函数为 ? 2B πρ sin ρ ? 1( 2d ) ? ? π 2B H S? L ( ρ) ? ? ?0 ρ ? 1 ( 2d ) ? ?

(4.7)

冲激响应为
π 1 ? 4 B ? s ? sin( 4 B ? s) 1 ? 4B ? 2 hS ? L ( s) ? ? ? 2? π ? 1 ? (4 B ? s ) 2
2

(4.8)

HS?L(ρ)

hS ? L
1.0 ? 1 ? 4 B ? ? ? 0.8 2 ? π ? 0.6 0.4 0.2

2

hS ? L ( s ) 2 π 2d 2

?2B

?B

0

B

2B

ρ

-10 -8 -6 -4 -2 0

2 4 6

8 10 y ? 4Bs

?4d?3d?2d?d d 2d 3d 4d 0 s

(a)频率特性

(b)连续空域特性

(c)离散内插后特性

(a) Frequency characteristic (b) Spatial characteristic (c) Discrete characteristic after interpolated 图 4.2 S-L 滤波函数 Fig.4.2 S-L filter

H S ? L ( ρ) 与 hS ? L ( s) 的图形见图 4.1(a) 和 (b) 。仍采用采样间隔 d ? 1 (2B) ,对 hS ?L ( s) 进行均匀采样,即以 s ? nd ? n (2B) 代入式中,得 ?2 hS ?L (nd ) ? 2 2 n ? 0,?1,?2, ? π d (4n 2 ? 1)
(4.9)

虽然 R-L 滤波器形式简单、使用方便,但是存在有 Gibbs 现象,会扩大伪影 的边界,给伪影识别带来困扰。S-L 滤波器重建图像,其振荡响应明显减小,抗噪 性能更强,所以本文考虑采用 S-L 滤波器对投影数据进行滤波,以增强投影正弦 图中伪影信息。图 4.3 为某工件投影正弦图采用 S-L 卷积滤波前后的对比图(局部 图像) ,可以看出,滤波前正弦图像(a)对比度低,根本无法看出伪影信息。而滤波
29

重庆大学硕士学位论文

后的(b)图中不仅工件弦线更清晰,而且伪影信息也得到增强,以或亮或暗的竖线 形式明显的表现在正弦图上。

(a) 原始投影图像 (a) Original sinogram

(b) 滤波后投影图像 (b) Filtered sinogram of (a)

图 4.3 某工件投影数据 S-L 滤波前后正弦图 Fig.4.3 Sinogram of one object before and after S-L filtered

4.1.2 投影积分
根据 2.3 节环形伪影的特征分析可知,受工件材料和探测器性能的影响,CT 系统在扫描成像时,探测器的响应通常是非线性且不稳定的,常见的环形伪影在 重建图像中大部分表现为不均匀的(或亮或暗的)弧形或半环伪影,相对应投影正 弦图中的竖直伪影则表示为不均匀的或间断的。所以说同一探测器在不同的时间 分度采集的投影数据,有以下两种情况:正常和异常。正因为不确定某一时间分 度下的投影数据是正常的还是异常的,如果只是按一个时间分度采样的投影数据 值,来识别出导致投影正弦图中条状伪影的位置,那是不可能的。而且一次采样 的投影数据中,高频信息既可能是伪影数据,也可能是工件的细节信息。 所以,需要综合考虑和分析各个探测器在总的扫描时间下,采集的总投影数据 值变化情况。对滤波后的投影数据沿着时间分度方向做一次积分投影,即线积分, 其公式为:
m ~ sum(i) ? ? P (k , i), i ? 1?n k ?1

(4.10)

~ 式中, P (k , i) 为第 i 探测器通道在第 k 时间分度下所采集的投影数据, m 为扫
描旋转一周的时间分度。

30

4 基于投影曲线高通滤波的环形伪影校正

(a)投影正弦图 (a) Original sinogram

(b)投影数据列投影曲线 (b) Projection Profile of integrated projection data in vertical 图 4.4 某工件投影数据正弦图和列投影积分曲线 Fig.4.4 Sinogram and vertical projectin profile of one object

图 4.4(a)为原始投影正弦图,对它的投影数据沿着竖直方向投影积分后,二维 的投影数据矩阵线积分后变成一维的数组,其数据曲线如图 4.4(b)所示,其横坐标 对应着探测器位置,纵坐标表示相应的探测器在整个扫描周期下的总采样幅值。 正常情况下,投影数据的列投影曲线应是近似光滑且连续的,但是由于投影正弦 图中列方向竖直伪影的存在,改变了这些位置的投影数据值,使得投影曲线受到 了较大程度的影响,从而表现出噪声干扰式的特点。根据比尔定律可知,在同等

31

重庆大学硕士学位论文

强度射线源下,探测器采集的投影数据受工件的密度和构造的影响,所以图中数 据恒定区域(平坦段)表示图像均匀灰度区域,曲线中小幅值的波峰和斜坡处的 区域表示的是重建图像的边缘细节信息,而突变的波峰或波谷位置对应着异常探 测器位置。

4.1.3 一阶微分
积分后的投影数据曲线在一定程度上反映了探测器系统在扫描周期内采集数 据的情况,但是如图 4.4(b),如果只是采用低通滤波器来平滑投影曲线波峰的高频 信息,必然会平滑掉小幅值波峰的细节信息,导致图像模糊。而且复杂工件的细 节信息在投影曲线中也表现为一定幅度的波峰,而由探测器异常所导致的突变的 伪影数据并不突出,很难找到合适的阈值来查找伪影的位置。所以,为了更好地 区别伪影和工件细节信息,考虑对投影数据曲线中的突变信号进行锐化,更进一 步增强曲线中伪影的信息,抑制工件细节信息数据的影响。 常用的信号锐化滤波器有:空间域微分算子和频域高通滤波器。本文采用微 分算子来增强投影曲线中伪影数据的突变,同时对图像均匀区域和工件细节区域 进行平滑。 数学函数的微分可以用不同的术语定义,也有各种方法定义这些差别,然而 对于一阶微分的任何定义,都必须保证以下几点:(1)在平坦段(灰度不变的区域) 微分值为零;(2)在灰度阶梯或斜坡的起点处微分值为零;(3)沿着斜坡的微分值为 零[38]。任何二阶微分的定义也同样遵循以上几点。由于一阶微分处理对阶梯具有 较强的响应,所以,本文采用一阶微分对投影曲线进行锐化,其公式为:
diff (i) sum(i) ? sum(i ? 1) ? ?i Δi

(4.11)

因为微分是用来处理数字量,其值是有限的。所以,最大灰度级的变化也是 有限的,变化发生的最短距离是在两相邻像素之间,即 Δi ? 1 。因为只考虑幅值插 值大小,不计方向,为方便判断,式 4.11 简化为
diff (i) ? sum(i) ? sum(i ? 1) ?i

(4.12)

如图 4.5 为对列投影积分曲线数据一阶微分处理后的投影曲线,可以看出,微 分后正弦图均匀区域的投影曲线更加平坦,而曲线中的代表着伪影的高频信息更 加突出,伪影明显得到增强。与图 4.4(b)的列投影曲线相比,经一阶微分锐化后, 曲线中伪影信息更明显,为进一步确定伪影的位置打下基础。

32

4 基于投影曲线高通滤波的环形伪影校正

图 4.5 一阶微分后投影曲线 Fig.4.5 Projection profile of first order differential

4.1.4 均值滤波
从图 4.5 中看出微分锐化后, 部分非异常的探测器数据受工件材料和形状或探 测器系统噪声的影响,也表现为一定的突变,为了更好的选取阈值,判别出异常 和正常的投影数据,对微分锐化的投影曲线再按一定的方式进行平滑处理。空间 域中常用的平滑滤波器有平滑线性滤波器和统计排序滤波器。 平滑线性空间滤波器的输出是包含在滤波掩模邻域内信号的简单平局值。因 此,这种滤波器通常也称为均值滤波器。它的概念非常直观,就是用滤波掩模确 定的邻域内信号的平均值代替原位置的信号值,这种处理减小了信号曲线的“锐 化”变化。 统计滤波器是一种非线性的空间滤波器,它的响应基于掩模包围区域信号值 的排序,然后用统计排序结果决定的值代替中心点的值。统计滤波器中最常见的 例子是中值滤波器,正如其名,它是将邻域内的信号中值代替该信号值。该滤波 器对处理脉冲噪声非常有效。 本文对投影曲线进行平滑是为了平滑由被测物质形状和材料引起的非伪影信 号,所以采用均值滤波器能满足要求。根据 2.3 节环形伪影的特征分析可知,虽然 正弦图中所示的一条灰色伪影条带对应着一个异常的物理探测器,但它却表现为 多个像素的宽度,且一个伪影的像素宽度与高频率采样插值的次数相等。所以, 需按照系统扫描采样的插值次数来对投影曲线进行均值滤波。采用 1? W 的掩模,

W ? I num ,将高频的投影曲线压缩还原为原始投影曲线,其计算公式为:

33

重庆大学硕士学位论文

m(t ) ?

1 I num

? diff (t ? I
j ?1

I num

num

? j)

(4.13)

式中,t ? 1?n0 , n0 为实际物理探测器数目,且 n0 ? n / I num , n 为投影数据的 列数,即重建图像的大小, I num 为高频率采样插值次数。 按式 4.13 对一阶微分后的投影数据进行均值滤波, 滤波后的投影曲线如图 4.6 所示,从图中可以看出,伪影的信息得到了进一步增强,而且工件边缘和细节处 的信息得到抑制,投影曲线更加平滑。与图 4.5 一阶微分后的投影曲线进行对比, 可以发现曲线中高频更加突出,这样更容易确定阈值,找到异常的物理探测器位 置,也就是投影数据中异常数据的列位置。

图 4.6 均值滤波后投影值曲线 Fig.4.6 Projection profile after mean filtered

4.1.5 阈值选择
CT 图像噪声主要以光子波动引起的随机误差为主, 即噪声符合正态分布规律, 而环形伪影或条状伪影均属于粗大误差,所以可以根据粗大误差判别准则来确定 投影正弦图中条状伪影的位置,超过一定的置信区间的可看作是粗大误差。 正态分布是概率论中最重要的一种分布,也是自然界最常见的一种分布。该 分布由两个参数——平均值和方差决定,概率密度函数曲线以均值为对称中线, 方差越小,分布越集中在均值附近。若已知的密度函数(频率曲线)为正态函数 (曲线)则称已知曲线服从正态分布。其中 μ 、 σ 是两个不确定常数,是正态分布 的参数,不同的 μ 、 σ 对应不同的正态分布。正态函数的一般表达式为:

34

4 基于投影曲线高通滤波的环形伪影校正

f ( x) ?

? 1 e 2π σ

( x? μ )2 2σ 2

(4.14)

μ 是正态分布的位置参数,描述正态分布的集中趋势位置。正态分布以 X ? μ
为对称轴,左右完全对称。正态分布的均数、中位数、众数相同,均等于 μ 。 σ 描 述正态分布资料数据分布的离散程度, σ 越大,数据分布越分散, σ 越小,数据分 布越集中。

Confident interval

(a)投影值直方图 (a) Histogram of projection value

(b)拟合的正态曲线 (b) Fitted normal curve

图 4.7 均值滤波后数据的直方图与正态曲线 Fig.4.7 Histogram and normal curve of averaged projection value after mean filtered

图 4.7(a)为均值滤波后的投影数据的直方图统计,从它的直方图分布规律,可 以发现该分布可视为正态分布,根据正态分布的特性分析,可知 均值 μ 为
μ? m(1) ? m(2) ? ? ? m(n0 ) n0
1 n0 ? (m(t ) ? μ) 2 n0 ? 1 t ?1

(4.15)

均方差 σ 为
σ?

(4.16)

所以,根据粗大误差的判断准则,伪影阈值的选取公式为:

threshold ? μ ? a ? σ

(4.17)

其中, a 为置信系数,图像中伪影越多, a 越小,置信区间越小。它的概率密 度曲线如图 4.7(b),阴影区域为置信区间。实际应用下,当 a ? 2 时,伪影处的投 影值一般都能排除在置信区间之外。

35

重庆大学硕士学位论文

4.2 单环伪影投影数据校正
4.2.1 传统的线性插值
给定一组有序的数据点 pi (i ? 0,1,2, ?, n) ,要求构造一条曲线,顺序通过这些 数据点,称为对这些数据点进行插值。目前,常用的数据插值方法有邻域插值、 线性插值和多项式拟合插值。因为邻域插值对单个数据插值效果较好,但对多个 无值数据点的周围邻近点仍然无数据,无法进行插值。而线性插值和多项式拟合 插值都有效,但是线性插值比多项式拟合算法更简便,费时更少,所以本文采用 线性插值算法。下面对线性插值算法做个简单介绍。 假设已知直线上两点的坐标 x1 、y1 和 x2 、y 2 , 要得到 ?x1 , x2 ? 区间内某一位置 x 在直线上的值,采用线性插值可以得到: y ? y1 y2 ? y1 ? x ? x1 x2 ? x1 它的变换公式为: (4.18)

y ? y1 ?

y2 ? y1 ( x ? x1 ) x2 ? x1

(4.19)

假设投影正弦图矩阵大小为 M×N,其中第 N1 列到第 N2 列数据为伪影异常 数据, 利用第 N1-1 列和第 N2+1 列上的数据, 垂直于这两条直线的 M 条线性函数, 根据线性插值公式对第 N1 列到第 N2 列数据的点进行插值,即可得到伪影处数据 的近似真实值,从而达到抑制投影正弦图的直线伪影的目的。

4.2.2 存在的问题
根据 4.2 节所述的环形伪影所述的算法,定位出投影正弦图中每个竖直方向伪 影最后一个像素位置,从而找到对应的导致环形伪影的异常探测器位置。假设定 位出其中一条伪影的位置为 j,根据一个异常探测器所致的环形伪影的宽度与高频 采样插值次数相等的特征, 可知相对应的 j ? ( I num ? 1), j ? ( I num ? 2),?, j 处的投影数 据都可能是异常的。所以,该处伪影异常的投影数据修正公式为 P(k , j ? τ ) ? d ? P(k , j ? I num ) ? (1 ? d ) ? P(k , j ? 1) (4.20)

式中, P(k , j ? τ ) 为投影数据中第 k 行,第 j ? τ 列数据。其中, k ? 0,1, ?, m 为 探测器扫描视角, m 为投影数据矩阵的行数, τ ? 0,1,?, ( I num ? 1) , I num 为探测器系 统高频采样插值次数,而权值 d 为:

d?

τ ?1 I num ? 1

(4.21)

图 4.8 为某工件采用传统的线性插值算法校正前后的投影正弦图, 通过对比可 以发现,虽然正弦图中条状伪影得到了很大地抑制,但是校正不够彻底,投影正 弦图中依然存在条带伪影伪影。

36

4 基于投影曲线高通滤波的环形伪影校正

(a)校正前 (a) Before correction

(b)校正后 (b) After correction

图 4.8 线性插值校正前后正弦图 Fig.4.8 Sinogram before and after line interpolation correction

4.2.3 线性外推插值校正
图 4.9(a)为某一个视角下原始投影数据曲线的局部放大图,假设扫描该投影数 据时,采用的插值次数 I num 为 2,根据 4.2 节所述的算法定位出图中 j 、 j ? 1 为其 中的一个伪影位置,从曲线可以看出,该处投影数据值与周围的数据值相比,数 据偏低,从而导致了投影正弦图的暗条和重建图像中的暗环。 如果采用传统的线性插值对该处投影数据进行校正,可以发现虽然投影数据 值变大,得到了一定的修正,但是与周围的曲线趋势相比,校正后的数据值仍然 偏低,校正不彻底。这是因为探测器在扫描时通常会发生光信号或电信号串扰, 虽然 j 、 j ? 1 处是异常探测器导致的异常数据位置, 但是它左右相邻的投影数据值

j ? I num 、 j ? 1 也会因为该异常探测器串扰的影响,使得它们的数据值较真实值偏
低。所以,如果采用左右相邻的且已被串扰了的投影数据值来对伪影值进行线性 插值校正,必然会导致校正后投影曲线仍然存在拐点,校正不彻底。 综上所述,在对 j 、 j ? 1 处伪影数据进行校正前,应先对左右相邻的 j ? I num 、

j ? 1 处投影数据值进行数据外推校正,以消除它们受异常探测器采集数据的影响,
近似恢复它的真实值。根据投影曲线平滑且近似线性的特性,可以得到 P( j ? 1) ? 2 ? P( j ? 2) ? P( j ? 3) P( j ? I num ) ? 2 ? P( j ? I num ? 1) ? P( j ? I num ? 2)

(4.22)

从图 4.9(a)校正后的投影曲线可以看出,经过线性外推插值校正后,异常数据 得到进一步的修正,与传统的线性插值相比,校正后投影数据曲线更加平滑。对 比传统线性插值和线性外推插值校正后的投影正弦图 4.9(b)和 4.9(c), 也可以发现, 采用线性外推插值伪影校正更彻底。
37

重庆大学硕士学位论文

j-Inum-2

j-1 j

j+3

(a)某一视角下的投影数据曲线 (a) Projection profile on a certain view

(b)传统的线性插值校正 (b) Traditional linear interpolation correction

(c)线性外推插值校正 (c) Linear extrapolating interpolation correction

图 4.9 线性插值和线性外推插值对比图 Fig.4.9 Contrast image with traditional linear interpolation and linear extrapolating interpolation

4.3 实验结果与分析
4.3.1 伪影识别结果
本文对一幅含有环形伪影的 CT 图像进行伪影校正试验, 如图 4.10 所示, 图(a) 为 CT 图像的原始投影数据正弦图(为了清晰地现出正弦图的伪影信息,此处采用 的是 S-L 滤波后的投影正弦图),采用本节所提的基于投影曲线高通滤波的投影正 弦图中条状伪影识别的算法,检测出投影曲线中异常探测器所在的位置如图 (b)所 示。虚线为按照粗大误差准则确定的阈值分割线,以上部分就是基于投影曲线的 高通滤波伪影识别的结果。

38

4 基于投影曲线高通滤波的环形伪影校正

(a)原投影数据正弦图 (a) Original projection sinogram

(b)本节方法伪影识别图 (b) Artifacts recognized result by suggested method 图 4.10 不同方法伪影识别结果图 Fig.4.10 Artifacts recognized result by different methods

从图中可以看出,本节所提的基于投影曲线高通的环形伪影识别算法共识别 出了 7 条伪影,基本能检测出投影正弦图中的伪影,而且没有误检情况,但是也 存在一定的不足。本文所提的检测算法是:如果某一个探测器扫描时在某一时间 下出现了不一致,则对该探测器扫描一周的所有投影数据都进行校正;如算法流 程第③步采用了一阶微分的差分处理,必然对两个相连的异常探测器所致伪影识 别具有一定的局限性,这也局限本节所提的伪影识别算法只用于分散单环的环形 伪影识别有效。

4.3.2 伪影校正结果
① 伪影校正效果评价标准 实际应用中,对伪影的去除效果进行评价主要有两种途径:主观评价和客观 评价。主观评价就是通过目视对比伪影校正前后图像,属于定性评价;客观评价

39

重庆大学硕士学位论文

就是通过计算一定的参数指标,定量评价校正方法在噪声消除及细节保持等方面 的能力。 通常采用图像的标准差和信噪比作为图像的质量评价指标,若图像的大小为

M ? N ,则其标准差和信噪比分别为:
?? σ 1 M ?N

?? ( X
i ?1 j ?1

M

N

2 ij

? μ)

(4.23) (4.24)

μ SNR ? 20 lg( ) ? σ

其中, X ij 为图像中像素点 (i, j ) 的灰度值, μ 为图像灰度均值: 1 M N μ? ?? X ij M ? N i ?1 j ?1

(4.25)

图像标准差代表了图像区域中所有点偏离均值的程度, 它反映了图像的不均匀 性。图像伪影越多,图像的灰度值变化越大,图像方差也就越大。同时,如果图 像的信噪比越高,则代表图像质量越好;若校正后图像的信噪比提高越多,则代 表此方法的伪影校正能力越好。 ② 实验结果 本文对含有环形伪影的原始工业 CT 图像进行校正试验, 图像大小为 797×797, 投影正弦图如 4.10(a)所示, 截取出重建图像中部分环形伪影所在的区域 300×300, 如图 4.11(a)所示。而图(b)、(c)、(d)分别为采用本节所提的方法、极坐标变换滤波 法、小波—傅里叶变换滤波法对(a)进行环形伪影校正的结果。 比较应用各种算法对图 4.11(a)伪影校正的结果图像, 从视觉上可以看出极坐标 变换滤波校正后图像视觉模糊,细节失真最严重且伪影校正程度小;小波—傅里 叶变换滤波法校正后,虽然图像中的环形伪影得到一定的抑制,且图像边缘细节 保持效果好,但是同时给重建图像带来了新的环形伪影,降低了图像的质量;而 本节所提的校正算法则几乎完全去除 CT 图像中的环形伪影, 并且较好地保持了图 像细节纹理信息,校正效果比较理想。

40

4 基于投影曲线高通滤波的环形伪影校正







(a) 原图局部放大图 (a) Local image

(b) 本节方法 (b) Suggested method

(c) 极坐标变换滤波法 (c) Polar transformation filtering method

(d) 小波—傅里叶滤波法 (d) Wavelet—FFT filtering method

图 4.11 实验图像 Fig.4.11 Experiment images

下面通过一定的参数指标对各方法的校正结果进行评价,表 4.1 为计算图 4.11(a)、(b)、(c)、(d)中如图(a)白色虚框①所示的工件均匀区域的标准差和信噪比, 以及②所示的细节区域标准差。同时采用图 4.11(a)垂直工件边缘处 3 条白色虚线 沿梯度方向所对应位置的灰度均值,来计算工件边缘梯度变化情况,各校正方法 校正结果图像的工件边缘梯度曲线如图 4.12 所示。 通常对于细节区域来讲,标准差既反映了图像被环形伪影污染的程度,同时 也反映了图像细节信息的丰富度,所以校正前后图像在此区域的标准差应相差不 大。而对于均匀区域来讲,由于基本没有任何细节纹理的变化,标准差主要反映 了图像被环形伪影污染的程度,其值越小越好;而信噪比则代表图像中有用信号

41

重庆大学硕士学位论文

与伪影噪声的比值,其值越大越好。
表 4.1 实验数据评价指标 Table4.1 Experimental data of evaluation indexes 灰度图像 原局部图 本节方法 极坐标变换滤波花法 小波—傅里叶变换滤波法 细节区域标准差 均匀区域标准差 均匀区域信噪比 SNR(db) 11.745 14.265 19.042 10.945

? σ
67.902 69.190 64.999 53.387

? σ
3.943 2.714 1.687 6.156

图 4.12 各校正方法结果图像的工件边缘梯度曲线 Fig.4.12 Work piece edge gradient profiles of resulting images by different correction methods

从表 4.1 可以看出,小波—傅里叶变换滤波法在细节区域标准差变化最大,均 匀区域中标准差增大,信噪比降低;而极坐标变换滤波法和本节方法在细节区域 的标准差变化都不大,且在均匀区域标准差都有所减小,信噪比都得到一定的增 大。这是因为小波—傅里叶变换滤波法在去除图像环形伪影的同时,给图像带来 非常明显的新环形伪影。但是,对比图 4.12 可以看出,本节方法和小波—傅里叶 变换滤波法校正后图像与原图边缘梯度曲线几乎重合,而极坐标变换滤波法边缘 梯度变小,图像模糊。综合考虑,本节所提的校正算法不但很好地校正了图像中 的环形伪影,提高了图像的信噪比,且很好地保持图像边缘细节和分辨率。

4.3.3 分析与总结
本节中提出了基于投影曲线高通滤波的环形伪影特征识别算法,将二维空间

42

4 基于投影曲线高通滤波的环形伪影校正

线特征识别问题转化为一维空间点特征查找,实验证明,该算法实现简单,能有 效地识别投影正弦图中单个分散的环形伪影,同时解决了探测器性能不稳定的问 题。另外,分析了传统线性插值对单个分散的环形伪影校正不彻底原因,提出线 性外推插值算法,使环形伪影校正更彻底。最重要的是,本节所提的环形伪影校 正算法是仅对定位后伪影处的异常数据进行校正,避免了滤波法对非伪影处的处 理,去除伪影的同时,也很好地保护图像边缘边界。 虽然本节所提的环形伪影校正算法实现简单,且对环形伪影校正有较好的效 果,但是也存在一定的不足。在对投影正弦图伪影进行识别时,采用了差分处理, 显然对两个或者多个相连的异常探测器所致伪影识别具有一定的局限性;而且对 伪影数据进行校正时,对投影数据的每一行(时间方向)进行插值处理,如果对 多个相连的环带进行校正时,就会导致严重的马赛克效应。这些因素都局限本节 所提的伪影校正算法对分散单环的环形伪影的识别,而对多个连续的环形伪影的 校正效果较差,需待改进。

4.4 本章小结
本节提出一种基于投影曲线高通滤波算法,对投影正弦图中的环形伪影进行 识别,并分析传统的线性插值算法校正不彻底的原因,提出线性外推插值算法对 伪影数据进行校正;给出了本节算法的实验结果,并与其他学者所提的算法进行 了对比。

43

重庆大学硕士学位论文

44

5 基于改进 Canny 算法的环形伪影校正

5 基于改进的 Canny 算法的环形伪影校正
根据上节所提算法对多个连续的环形伪影校正的局限性,本节按照数字图像 处理中图像边缘检测的思想,来识别出投影正弦图中的竖直条状伪影位置,即找 到不一致的探测器通道的位置。通常正弦图中主要的边缘细节都是重建图像工件 边缘与细节的正弦线,而环形伪影只表现为竖直方向的弱边缘,如何实现投影正 弦图中竖直方向弱边缘的检测,是按照这种思想对投影正弦图中的条状伪影校正 的难点之一;若重建图像中存在多个相连的环带,采用低通滤波或上节所述的伪 影数据校正方法,必然会使校正后的投影正弦图中出现严重的马赛克现象,重建 图像细节损失,如何校正图像中多个相连的环带,又不损失图像细节,这是本节 的另一个难点所在。

5.1 基于改进的 Canny 算法的条状伪影特征识别
5.1.1 传统 Canny 算法的准则与实现
1986 年,John Canny 提出一个优良的边缘检测算子应满足以下三个判断标准
[39~40]



① 高信噪比准则:即输出信号的信噪比最大,使得非边缘点判为边缘点或边 缘点判为边缘点的概率最小。信噪比的数学表达式为
??

SNR ?

??

? G(? x) f ( x)dx
??

(5.1)
2

σ

??

?f

( x)dx

其中, f ( x) 为边界是 ?? ω,?ω? 的滤波器脉冲响应; G(? x) 为边缘函数; σ 为高 斯噪声的均方根。信噪比越大,误检率越低,提取的边缘质量越好。 ② 高定位精度准则:即检出的边缘点与实际边缘点距离最小,使得定位精度 最高。定位精度的数学表达式为
??

Localizati on ?

??

? G?(? x) f ?( x)dx
??

(5.2)
2

σ

??

? f ? ( x)dx

其中,G?( x) 和 f ?( x) 分别是 G(? x) 和 f ( x) 的一阶导数,Localizati on 的值越大, 定位精度越高。其目标是求一个函数 f ( x) ,使得下面的这个式子达到最大值。

45

重庆大学硕士学位论文

J( f ) ?

?

?w

?w

G (? x) f ( x)dx

?

?w

?w

G?(? x) f ?( x)dx
?w

?w n0 ? ? f 2 ( x)dx ? ? ? ? ?w ?

n0 ? f ? 2 ( x)dx
?w

(5.3)

③ 单一边缘响应准则:即在理想情况下,用滤波器对噪声响应的两个峰值间 的距离来近似用滤波器响应一个边缘点的长度。输出信号中相邻两个极大值的距 离是相邻两个零交叉点距离的两倍,每个边缘点有且只有一个响应,最大限度抑 制伪边缘出现。要保证单一边缘响应,检测算子脉冲响应导数的零交叉点评价距 离 D?( f ) 应满足
? ?? f ? 2 ( x)dx ? 2 ??? ? D?( f ) ? π ? ? ? f ( x) ?? ? ? f ?? 2 ( x)dx ? ? ?? ?
1

(5.4)

式中 f ??( x) 为 f ( x) 的二阶导数。 以上述指标和准则为基础,利用函数求导的方法导出一个由边缘定位精度和 信噪比乘积组成的表达式,这个表达式近似于高斯滤波的一阶导数[41], 2 2 exp( ? ( x ? y ) 2 ) 2σ G? 2πσ 2 设 n 为任意方向,则高斯函数在这个方向上的一阶导数 G n 为 ?G Gn ? ? n ? ?G ?n

(5.5)

(5.6)

事实上,方向 n 应取与边缘垂直的方向,可根据下式给出很好的估计。 ?(G ? I ) (5.7) n? ?(G ? I ) 符号 ? 表示卷积, 边缘点定义为 G n 作用于输入图像 I 后的局部极大值点, 在这 个点上有
G(? ?x)G ? I ? 0

(5.8)

基于前面所述的三个判断准则, Canny 提出了一种最佳的边缘检测算子, 称为 Canny 算法,它的具体实现过程是[42~43]: ① 用高斯滤波器平滑图像,从而达到滤除噪声的目的。 一维高斯平滑函数为

Gn ?

1 ? x2 exp( ) 2πσ 2 2σ 2

(5.9)

用高斯滤波器平滑图像的过程实际就是分别按行和列对原始图像 f ( x, y) 进行 卷积的过程,最后得到平滑图像 K ( x, y) 。这里高斯滤波函数的标准差 σ 控制着平 滑程度,当 σ 较小时,滤波器也较短,卷积运算量小,定位精度高,但信噪比低; 当 σ 较大时, 情况则恰好相反。 因此, 要根据实际需要适当选取高斯滤波器参数 σ 。

46

5 基于改进 Canny 算法的环形伪影校正

② 用一阶偏导的有限差分来计算梯度的幅值和方向。 采用 2 ? 2 邻域一阶偏导的有限差分计算平滑后图像 K ( x, y) 的梯度幅值

M ( x, y) 和梯度方向 θ ( x, y) 。
M ( x, y) ? Gx ( x, y) 2 ? G y ( x, y) 2
) E x ( x, y ) 0.5 ? ?? 0.5 0.5? ? 0.5 , fy ? ? fx ? ? ? ? ?? 0.5 0.5? ?? 0.5 ? 0.5? θ ( x, y ) ? arctan( E y ( x, y )

(5.10) (5.11) (5.12)

其中, E x 和 E y 分别为原图像 K ( x, y) 被滤波器 f x 和 f y 沿着行和列影响的结果。 ③ 对梯度幅值进行非极大值抑制。
方向 方向3方向3 方向 2 2
M[x-1,y-1] (3) M[x,y-1]
[i-1,j-1] M[x-1,y]

方向4 方向 4

M

(3)

(2)

M[i-1,j]

M[x-1,y+1]
M[i,j]

(2)

(4)

M[i-1,j+1] (4) M[i,j+1]

M[i,j-1]
(1)

(1)

(1)

M[x,y]
M[i+1,j-1](4)

M[x,y+1]
M[i+1,j] (2)

(1)

M[i+1,j+1] (3)

方向1 方向4 方向2

方向1

M[x+1,y-1](4)

M[x+1,y]

(2)

M[x+1,y+1] (3)

方向1
方向3

方向1 方向2 (b)梯度方向示意图 方向3

(a)非极大值抑制图 (a) Non-maximum suppression diagram

方向4

(b) Gradient direction schematic diagram

图 5.1 对梯度幅值进行非极大值抑制 Fig.5.1 Non-maximum suppression in gradient magnitude

非极大值抑制的具体做法见图 5.1(a)。图中每个单元格内 M 函数表示该像素 的梯度幅值,括号中的数字表示梯度方向代号,梯度方向代号相同的单元格表示 在同一梯度方向上。图(b)表示将像素 ( x, y) 的梯度方向角的变化范围缩减到 4 个方 向之一,4 个方向如(b)所示。 在同一方向上将像素 ( x, y) 所具有的梯度幅值 M ( x, y) 与其相邻像素的梯度幅值进行比较,若 M ( x, y) 非局部极大值,就将其设为 0,即 不是边缘。 ④用双阈值算法检测和连接边缘。 对图像进行双阈值操作和边缘连接。确定上、下限阈值 T1 和 T2 ,对于经过非 极大值抑制处理的图像,如果像素点的梯度值 M ( x, y) ? T2 ,则直接把此像素点标 记为边缘像素点;如果像素点的梯度值 M ( x, y) ? T1 ,则把此像素点标记为非边缘 像素点;如果像素点的梯度值 T1 ? M ( x, y) ? T2 ,则标记此像素点位“准像素点” 。 在双阈值标记完后,搜索图像中的“准像素点” ,并选择其邻域点的位置寻找

47

重庆大学硕士学位论文

是否有梯度值 M ( x, y) ? T2 的点存在,如存在,则标记此像素点为边缘,否则标记 为此像素点为非边缘像素点。Canny 算法的方框图如图 5.2 所示。
输入 图像 输出边缘 图像

高斯 平滑

梯度 计算

非极大 值抑制

双阈值 检测

图 5.2 传统的 Canny 算法框图 Fig.5.2 The flowchart of traditional Canny algorithm

5.1.2 存在的问题
目前,应用中 Canny 算子与 Roberts 算子、Sobel 算子、Laplace 算子等边缘检 测算法相比,存在很多优点。但是,将 Canny 算子应用于投影正弦图中条状伪影 的检测,仍然存在很多需改进的地方: ① 传统的 Canny 算子采用高斯滤波来对输入图像进行预处理,虽然有利于图 像噪声平滑,但是也会使得图像边缘更加模糊,从而使后续的检测过程变得困难, 同时降低了图像边缘的定位精度。 如果高斯滤波中 σ 取值小时,虽然边缘定位精度 高,但是图像平滑能力弱,抑制噪声作用小;如果要追求好的噪声抑制效果,必 须增大 σ 的取值,但又导致模板增大,从而导致最后定位的边缘精度低,同时计算 过程的运算量也急剧增加。通常环形伪影在投影正弦图中表现为竖直方向的弱边 缘,所以 σ 值的选择至关重要。 ② 传统Canny算法通过在 2 ? 2 邻域内求有限差分均值计算梯度幅值,对噪声 过于敏感,容易检测出孤立边缘点和伪边缘。 ③ 传统 Canny 算子对 4 个方向进行非极大值抑制,梯度方向被定义为 0 o 、

45 o 、 90 o 、 135 o 四个方向之一。而本文是对投影正弦图中条状伪影进行检测,即
对竖直方向进行检测,只需对梯度方向近似水平方向进行抑制。 ④ 采用梯度双阈值方法定位边缘,难以选择合理的下限阈值。下限阈值选择 过大,则无法定位全部边缘,边缘出现间断;而下限阈值选择过小,又可能出现 大量非伪影边缘。 所以,通过以上对Canny算子存在问题的分析,本文在采用Canny算子对投影 正弦图中条状伪影进行识别时,提出了改进的方案,既有效抑制了噪声,又尽可 能地保护了边缘细节,较好地解决了上述问题,进一步提高了检测精度。

5.1.3 改进的 Canny 算法
① 图像预处理 由于工业 CT 投影数据动态范围偏小,大部分 CT 图像偏暗,条状伪影信息在

48

5 基于改进 Canny 算法的环形伪影校正

原始投影正弦图中表现不明显或完全被覆盖。如果以原始正弦图作为输入图像来 检测与识别伪影,肯定无法达到预期的效果。所以,首先对投影正弦图进行增强 处理,以高通滤波后的投影正弦图作为检测的输入图像。增强伪影的算法与 4.1.1 节完全相同,这里不再重复介绍。 ② 梯度幅值的计算 由于图像的离散特性,其梯度值不能完全由 2 ? 2 邻域有限差分确定,这种做 法容易受到噪声影响。Canny算法先对图像做平滑处理,然后计算梯度;而小波变 换用高斯函数的一阶导数作为小波函数,直接对原图像进行小波变换计算梯度, 基于此,结合二者的优点,对原来的 2 ? 2 滤波器进行了改进,采用二维高斯函数 的一阶偏导数[41]

G?i

exp( ? (i 2 ? j 2 ) 2σ 2 ) πσ 2

(5.13)

构造滤波器计算梯度,从而得到梯度值的更好近似。 由式5.13构造的 i 方向滤波器( σ ? 1 )为 ? 0.0001 0.0010 0.0021 0 ? 0.0021 ?0.0014 0.0117 0.0261 0 ? 0.0261 ? ?0.0064 0.0523 0.1171 0 ? 0.1171 ? f i ? ?0.0106 0.0862 0.1931 0 ? 0.1931 ?0.0064 0.0523 0.1171 0 ? 0.1171 ? ?0.0014 0.0117 0.0261 0 ? 0.0261 ? 0.0001 0.0010 0.0021 0 ? 0.0021 ? ③ 非极大值抑制 图 5.3 中每个单元格 I 函数表示被检测的投影正弦图的灰度值,已知连续多个 探测器异常导致的环形伪影,在投影正弦图中表现为竖直方向多个相连的条带伪 影。如果像素 (i, j ) 对应其中一个条带的边界,像素 (i ? 1, j ) 和 (i ? 1, j ) 肯定都是处 于边界位置,他们的灰度值几乎相同。那么根据 Canny 算子像素导数的计算方法, 像素 (i, j ) 的导数更加“趋向” j 分量,而沿 i 方向的导数分量非常小。
I(i-1,j-1) I(i,j-1) I(i+1,j-1) I(i-1,j) I(i,j) I(i+1,j) I(i-1,j+1) I(i,j+1) I(i+1,j+1)

? 0.0010 ? 0.0001? ? 0.0117 ? 0.0014? ? ? 0.0523 ? 0.0064? ? ? 0.0862 ? 0.0106? ? 0.0523 ? 0.0064? ? ? 0.0117 ? 0.0014? ? 0.0010 ? 0.0001? ?

(5.14)

图 5.3 图像灰度表示图 Fig.5.3 The schematic diagram of the image gray

所以,为了实现投影正弦图中竖直方向条状伪影检测,只用考虑方向导数 j 分
49

重庆大学硕士学位论文

量比 i 分量大的情况,而对导数 i 分量比 j 分量大的不再判断。假设 M 为当前像素 梯度幅值,那么它与 M 1 ~ M 4 的位置对应着 2 种关系:当像素 (i, j ) 沿着 i 方向的方 向导数与沿着 j 方向的方向导数的符号相同时,如图 5.4(a)所示;当它们的符号相 异时,如图(b)所示。
M3 M3 M4 M4 M M M2 M2 M1 M1 M1 M1 M3 M3 M M M2 M2 M4 M4

(a)符号相同 (a) Identical sign

(b)符号相反 (b) Opposite sign

图 5.4 i, j 两个方向的方向导数符号关系 Fig.5.4 Directional derivative sign relations of i and j direction

传统的 Canny 算子中, 关系 1 和关系 2 梯度方向 θ 的绝对值范围都为 0 o ~ 45 o , 它包括了竖直方向和斜角方向的边缘检测。而对投影正弦图中条状伪影的识别, 只需对图中竖直方向的边缘进行检测。所以,在对梯度进行非极大值抑制时,还 需设定角度阈值来对梯度方向进行限制,来排除斜角方向边缘的检测。设定角度 阈值 Tang ,判断方向梯度 θ 的绝对值是否小于或等于 Tang 。角度阈值 Tang 越大,检 测斜角范围越大,检测的斜角边缘像素点越多;反之,检测的斜角边缘像素点越 少。 ④ 边缘检测与连接 CT 系统在扫描成像时,受光子波动的影响,采集的投影数据不可避免地会引 入随机噪声。采用 S-L 滤波器对投影正弦图进行增强预处理,图像在增强伪影的 同时,正弦图中其他的噪声也得到放大。从滤波后的投影正弦图中可以看出,通 过 S-L 滤波器高通滤波后的图像,光子波动引起的随机噪声表现为孤立的噪声点, 被测工件外轮廓的正弦线表现为强边缘,条状伪影边界表现为弱边缘。条状伪影 与正弦图中的随机噪声类似,只是在列方向上检测的边缘比其他噪声要密集。所 以采用 Canny 算子对投影正弦图中条状伪影进行识别,梯度方向和连接方法可表 述为: 1) 确定梯度阈值 T ,对于经过非极大值抑制处理的图像,如果像素点的梯度值

M (i, j ) ? T ,则直接把此像素点标记为边缘起始点;
2) 去除边缘起始点中工件外轮廓正弦线的强边缘像素点; 3) 设定链长度 T0 , 对每列的边缘起始点连通的局部模极大值与 T0 进行比较, 若大 于 T0 ,则保留;若小于 T0 ,则该列的边缘像素点全部去除;

50

5 基于改进 Canny 算法的环形伪影校正

4) 按照高频采样的插值次数合并、细化边缘像素点; 5) 对每列中边缘像素点之间连续间断的局部模与 T0 进行比较, 若小于 T0 , 则连通; 若大于 T0 ,则不作处理; 6) 设定链长度 T1 ,重复步骤(3)和步骤(5),获得条状伪影边缘检测图像。 在实际应用中,通常 T1 设定为 3 T0 ,所以,它与传统 Canny 算子相似,也是 采用了双阈值的处理方法。但是,这种双阈值的方法与传统 Canny 算法相比,对 于弱边缘检测更具有优势,而且它既可以检测出伪影边界的弱边缘,又可以避免 边缘长度较短的随机噪声的影响。

5.1.4 算法的实现
改进后 Canny 算法投影正弦图条状伪影检测算法如图 5.5 所示。
采用二维高斯滤 波的一阶偏导计 算梯度和幅值

读取投影 正弦图

使用S-L滤波 器增强伪影

高斯滤波

显示图像 边缘

梯度阈值T和链 长度阈值T0检 测和连接边缘

趋向i分量方向 的局部像素非极 大值抑制

限制梯度方向角 abs(θ)<± Tangle

图 5.5 改进的 Canny 算法流程图 Fig.5.5 The flowchart of improved Canny algorithm

5.2 环带伪影投影数据校正
5.2.1 问题分析
如果采用单个探测器异常的环形伪影校正方法,来对连续多个探测器异常的 环形伪影数据的校正,即分步地对每一时间分度下的数据进行校正。无论是采用 插值算法还是拟合算法,都必然会导致图像严重的马赛克现象,使得图像细节损 失严重。如图 5.6 所示,假设采用 5.2 所述的伪影识别算法已找到图(a)原正弦图中 伪影的位置,如白线框所示,它一共包括 5 个相连的伪影条带。若已知了该投影 数据高频采样的插值次数为 8,即表明每个伪影条带的宽度为 8,说明这里我们必 须完成 40 个像素宽度的环带校正。如果采用如第 4.3 节所述线性插值校正的方法 对这个环带伪影进行校正,校正后投影正弦图如图 (b)中白色虚线框内所示,可以 发现经过白色虚线框的两条弦线带部分被模糊。而且对比校正前后(a)和(b)的重建 图像也可以看出,直接影响校正后重建图像中其中一个孔位置的图像灰度对比度 和边缘细节信息。图(c)为经过校正前后重建图像的小孔灰度曲线对比图,从图中

51

重庆大学硕士学位论文

可以看出,校正后的图像灰度曲线梯度小于校正前灰度值梯度,进一步证明线性 插值校正使得图像边缘模糊,对比度降低。如果连续多探测器导致的伪影经过工 件外轮廓的位置,这种校正方法甚至还会导致重建后图像中工件边缘出现明显的 阶梯效应。

(a)原正弦图和重建图像

(d)校正后正弦图和重建图像

(a) Original sonogram and CT image (b) Original sonogram and CT image after correction

(c)校正前后小孔位置局部梯度对比图 (c) Local gradient comparison diagram of keyhole position before and after correction 图 5.6 线性插值校正 Fig.5.6 Linear interpolation correction
52

5 基于改进 Canny 算法的环形伪影校正

5.2.2 分段 B-样条拟合校正
根据投影数据中每行插值校正算法存在的问题,本节考虑多个连续环形伪影 在投影正弦图中校正特点,对条带伪影的校正方法进行研究。 参考 Gadallah[44]和刘正军[45]采用矩匹配去除光谱成像仪条带噪声的思想,通 过将各列的均值匹配到参考投影数据列的列均值上来恢复图像的真实的列均值分 布,只是对灰度分布不均匀的正弦图来说,其真实的列均值分布并不是一条直线, 矩匹配方法恢复出的列均值曲线为一条值为参考值的曲线,从而严重违背了图像 的真实灰度分布。根据上节改进的 Canny 算法已识别出投影正弦图中伪影的位置, 只对伪影的位置进行局部拟合校正,获得列均值曲线的趋势曲线,以作为真实的 列均值分布曲线。根据列均值与真实的列均值的比值,对伪影处的投影数据进行 校正。 光谱成像仪中条带噪声是均匀的,而投影正弦图中条状伪影是不均匀的,射 线通过的路径不同,造成伪影的严重程度也不同,如果采用同一参数值对同一探 测器在不同时间下采集的同一列投影数据进行校正,必然导致严重的伪影段校正 不彻底,而不严重的伪影段校正过度。所以,本文将检测的每一条伪影分成 n 段, 并将每一段看成均匀的伪影来求均值曲线进行校正。 前面所述已完成了投影正弦图中条状伪影的识别,进一步伪影校正的关键就 在于如何对列均值曲线进行伪影校正,较好地恢复出真实的列均值分布。目前, 对离散数据谱图的绘制主要有插值法、无理模型法以及曲线拟合法。插值法对多 值数据曲线的拟合存在着很多问题,例如有龙格现象、插值精度不够和不能处理 大斜率问题等;无理模型方法是新近发展起来的基于神经网络和遗传算法等无理 模型的曲线拟合法,在实际使用时,由于得不到统一的曲线表达式,无理模型方 法将使数据曲线的定量分析变得更加复杂;曲线拟合法是用规则的曲线来近似表 述数据的特性,不要求函数严格的经过每个插值点,但要求函数在插值点处与插 值点的拟合误差矢量按某种范数达到最小值[46]。常用的曲线拟合法主要有多项式 最小二乘法、样条函数拟合法等。 ① 最小二乘法 参考 n 次参数(n<m)多项式曲线

P(t ) ? ? ? i?i (t )
i ?0

n

(5.15)

其中:?i (t )(i ? 0,1, ?, n) 为 n 次多项式的一组基;? i ? [ xi , yi , zi ](i ? 0,1, ?, n) 为 待定的系数矢量;给定的测量数据点为 pk ? [ xk , yk , zk ](k ? 0,1, ?, m) ,相应参数分 割是 ?t:t0 ? t1 ? ? ? t m 。 如果当作插值问题来处理,将有如下线性方程组

53

重庆大学硕士学位论文

? ? ? (t ) ? P
i ?0 i i

n

k

k ? 0,1, ?, m

(5.16) (5.17)

写成矩阵形式 其中系数矩阵
?? 0 (t 0 ) ?1 (t 0 ) ? ? (t ) ? (t ) 1 1 ??? 0 1 ? ? ? ? ?? 0 (t m ) ?1 (t m )

?A ? P

?? 0 ? ?? ? A ? ? 1? , ?? ? ? ? ?? n ?

? ? n (t 0 ) ? ? ? n (t1 ) ? ? 为 (m ? 1) ? (n ? 1) 阶矩阵, ? ? ? ? ? ? n (t m )? ? P0 ? ?P ? (5.18) P?? 1? ??? ? ? ? PM ?

这里矢量方程个数 m+1,超过了未知矢量个数 n+1,这样的方程组,在一般情 况下是无解的,即一般地不存在严格依次通过这些数据点的插值曲线 P(t ) ,这时, 只能寻求在某种意义下最接近这些数据点的参数多项式曲线 P(t ) 作为逼近曲线
[47]

。 实际中,通常采用逼近曲线 P(t ) 上参数值 t k 的点 P(t k ) 与数据点 Pk 间的距离平

方和来度量逼近的程度。
J ? ? P(t k ) ? Pk
k ?0
m 2

m

2

? Jx ? Jy ? Jz
2

(5.19)

m ?n ? ?n ? 其中: J x ? ? ?? xi (t k ) ? xk ? , J y ? ? ?? yi (t k ) ? y k ? 达到最小,称为最小 k ?0 ? i ?0 k ?0 ? i ?0 ? ? 二乘逼近。 欲使 J 为最小,就应使 J x , J y , J z 都最小,即必须使下列偏导数为零

?J y ?J x ?J ? 0, ? 0, z ? 0 ?x j ?y j ?z j 由上式可推出高斯(Gaussian)正交方程组 ? T ?A ? ? T P

j ? 0,1, ?, n

(5.20)

(5.21)

最后, 采用正交化方法即可求出未知的系数阵 A, 相当于定义了逼近曲线 P(t ) 。 但是,当 m ? n 时,这里的逼近曲线也就是插值曲线。 ② 样条拟合法 区间 [a,b] 上具有分化 n 各小区间 ?:a ? x0 ? x1 ? ? ? xn ? b 的分段 k 次式

S k ( x) 为 k 次样条,在每个节点 xi (1 ? i ? n ? 1) 上具有直到 k ? 1 阶连续导数;当分段 S k ( x) 的函数为三次多项式, 在[a,b]上具有连续的二阶导数, 称为三次样条曲线[48]。
n ( xi ) ? M i (i ? 0,1,?, n) 来作为参考数据,构造函 假设以节点上的二阶导数值 S3

数 S3 ( x) 满足系列条件:

54

5 基于改进 Canny 算法的环形伪影校正

1) S3 ( x) 在[a,b]上具有连续的一阶、二阶导数;

2)在每个小区间 ?xi , xi?1 ?上, S3 ( x) 都是 n 次多项式; 3) S3 ( xi ) ? yi (i ? 0,1,?, n) 表示这个函数代表曲线上所有各点。 则在分段 ?xi?1 , xi ? 上有:

S3 ( xi?1 ) ? yi?1 , S3 ( xi ) ? yi
S ( xi?1 ) ? M i?1 , S ( xi ) ? M i
n 3 n 3

(5.22) (5.23)

由 Lagrange 插值有

xi ? x x ? xi?1 M i?1 ? Mi hi hi 式中 hi ? xi ? xi?1 ,对上式积分两次得: S3n ( xi ) ?
S3 ( x) ?

(5.24)

( xi ? x)3 ( x ? xi?1 )3 h2 x ?x h2 x ? xi?1 M i?1 ? M i ? ( yi?1 ? i M i?1 ) i ? ( yi ? i M i ) 6hi 6hi 6 hi 6 hi (5.25)

再对上式微分一次得: ( x ? x) 2 ( x ? xi?1 ) 2 h y ? yi?1 ' S3 ( x) ? ? i M i?1 ? M i ? i ( M i ? M i?1 ) ? i 2hi 2hi 6 hi

(5.26)

使得节点处的连接条件, 即:i 点在左区间 ?xi?1 , xi ? 和右区间 ?xi , xi?1 ?上的插值在

i 点处的一阶导数相等,得:
' ' S3 ( xi ? 0) ? S3 ( xi ? 0)

(5.27)

代入式 5.27,整理得: hi?1 ? h h h y ? y y ? yi?1 M i ? i?1 M i?1 ? i M i?1 ? i?1 i ? i 3 6 6 hi?1 hi 上式可表示为: (1 ? ai )M i?1 ? 2M i ? ai M i?1 ? gi 其中
(i ? 1,?, n ? 1)

(5.28)

(5.29)

ai ?
6( gi ?

hi?1 hi ? hi?1

(i ? 1,?, n ? 1)

(5.30)

yi?1 ? yi yi ? yi?1 ? ) hi?1 hi hi ? hi?1

(i ? 1,?, n ? 1)

(5.31)

边界条件的确定:对于起点的左边界 ?x0 , x1 ? ,当 x ? x0 时,代入式 5.28,可得: ? 6 ?y ? y (5.32) 2M 0 ? M1 ? ? 1 0 ? S ' (0)? h1 ? h1 ?

S3n (0) 若为已知时,由式 5.24,并取 x ? x0 时,得:
55

重庆大学硕士学位论文

n 2 M 0 ? 0 ? M 1 ? 2 ? S3 ( x0 )

(5.33)

于是起点的边界条件可写成: 2M 0 ? a0 M1 ? g0

(5.34) ? 6 ? y ? y0 ?? 时, ?; ? 时, 其中: 给定 y 0 给定 y 0 a0 ? 0 , g 0 ? 2 y0 a0 ? 1 , g0 ? ? 1 ? S ?(0)? h1 ? h1 ? 同理终点的边界条件可以写成:an M n?1 ? 2M n ? g 0 (5.35) ? 6 ? y ? yn ?; ? 时, ?? 时, ?? an ? 0 , g n ? 2 yn an ? 1 , 其中: 给定 y n 给定 y n gn ? ? ? n ? yn hn ? hn ? 构成定解问题:
a0 ? 2 ?1 ? a 2 1 ? ? ? ? 1 ? a n ?1 ? ? ? 0 0 ? ?M 0 ? ? g 0 ? ?? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ? a n ?1 ? ? ? ? ? ? ? 2 ? ?? ?M n ? ? ? ?gn ? ?

? 2 an

(5.36)

利用三次多项式拟合和三次 B 样条拟合对一投影数的其中一组列均值数据进 行拟合, 如图 5.7 所示。 通过对比发现, 虽然两者都拟合出了均值数据的趋势曲线, 拟合曲线都有较好的连续性和光滑性, 但是对比非伪影处数据的拟合效果, B 样条 拟合比多项式拟合更逼近列均值曲线,具有更好的保形功能。而且多项式拟合时, 拟合曲线容易出现振荡。所以本节考虑采用分段 B 样条拟合来获取列均值数据曲 线的趋势曲线,以作为真实的列均值数据。

图 5.7 不同拟合方法比较 Fig.5.7 Compare different fitting methods

56

5 基于改进 Canny 算法的环形伪影校正

5.2.3 算法的实现
根据前面的分析,本节采用了分段的 B 样条拟合的方法来实现 CT 图像中多 个连续的环带伪影的校正,假设投影数据的大小为 M ? N ,它的具体实现步骤如 下: ① 按照 5.2 节伪影识别的二值图,将相连的伪影投影数据按每段 m 行(时间 分度)分成均匀的 K 段来实现校正。通过分段将每一段看作均匀伪影,视其伪影 的严重程度相同。假设,需对第 k 段第 j 列的伪影数据进行校正,那么该段的列均 值为
P( k , j ) ? 1 m?k ? P(i, j) m i ?m?( k ?1)?1

(5.37)

式中,k ? 1,2,?, K , j ? 1,2,? N ,P(i, j ) 为第 i 行第 j 列处待校正的投影数据,
P(k , j ) 为第 j 列第 k 段的列均值数据。如果图像中同一探测通过伪影越均匀, 则m

越大,且越接近 M 值。 ② 排除伪影位置的列均值数据, 对剩余列均值数据 P(k , j ) 曲线进行三次 B 样 条拟合,以拟合获得的趋势曲线数据 Ptrend (k , j ) 作为真实的均值曲线。 ③ 计算第 k 段伪影的增益因子为

Ptrend (k , j ) P(k , j ) ④ 计算第 k 段伪影的真实投影数据为 c(k , j ) ?

(5.37)

P?(i, j ) ? P(i, j ) ? c(k , j )
其中, i ? n ? (k ? 1) ? 1, n ? (k ? 1) ? 2,?, k ? n 。

(5.38)

图 5.8 校正前后的投影数据值曲线 Fig.5.8 Projection profile before and after corrected
57

重庆大学硕士学位论文

根据图 5.7 的列均值曲线,B 样条拟合得到它的真实的列均值,以它们的比值 作为校正系数对这段伪影的某一行投影数据进行校正,它其中一组校正前后数据 如图 5.8 所示。对比校正前后的投影数据,校正前伪影处数据偏高,校正后数据幅 值下降,且保持了峰值细节信息。

5.3 实验结果与分析
5.3.1 伪影识别结果
本文对含有环形伪影的 CT 图像进行伪影校正实验,图 5.9(a)为它 S-L 滤波增 强后的投影正弦图,图(b)为采用本节所述的基于 Canny 算子条状伪影特征识别结 果。从图中可以看出,图(a)有较多的环形伪影,其中有单个分散的,也有多个相 连的,这些伪影在投影正弦图中表现为竖直方向的弱边缘。采用本节所述的基于 Canny 算子条状伪影特征识别结果如图(b)所示,图中白色竖线为所识别出伪影的 位置,可以发现正弦图中伪影信息基本都被识别,且定位在每个伪影右边界的像 素位置,定位精确,而且识别的边缘都是单级像素。

(a)含有伪影的原投影正弦图 (a)Original sinogram contained by artifacts

(b)伪影识别结果的投影正弦图 (b) Sinogram signed by edge with recognized artifacts

58

5 基于改进 Canny 算法的环形伪影校正

(c)分段 B 样条拟合校正后投影正弦图 (c) Sinogram after corrected by piecewise B-spline fitting 图 5.9 投影正弦图中伪影的识别与校正 Fig.5.9 Artifacts recognition and correction in sonogram

5.3.2 伪影校正结果
本节试验的原始工业 CT 图像大小为 753×753,如图 5.10(a)所示。图(b)是为 使用本节所提的方法进行环形伪影校正的结果图像,图(c)和(d)分别为图(a)白线框 区域①和②的校正前后局部放大图。其中,(a)和(b)重建图像对应的投影正弦图如 图 5.9(a)和(c)所示。



① ② ①

(a)校正前 CT 图像 (a) CT image before correction

(b)校正后 CT 图像 (b) CT image after correction

59

重庆大学硕士学位论文

(c)局部放大图① (c) Partial enlarged image of ①

(d)局部放大图② (d) Partial enlarged image of ②

(e)校正前后小孔位置局部梯度对比图 (e) Local gradient comparison diagram of keyhole position before and after correction 图 5.10 本节方法校正结果图 Fig.5.10 The result image of suggested method

对比图 5.9 校正前后的投影正弦图(a)和(c)可以看出,环形伪影在投影正弦图 中所表现的条状伪影基本都校正干净,且图中弦线细节保留较好;对图 5.9(c)与图 5.6(b)进行比较,可以发现采用分段 B 样条拟合校正法,条带伪影处的正弦线细节 得到保留,避免了各行线性插值校正时的马赛克现象;从图 5.10(c)和(d)的局部放 大图中也可以看出, 校正后工件的细节信息完全保留, 图像分辨率不会降低; 图(e) 为不同方法校正前后图(c)的小孔灰度曲线对比图,可以看出本节所提的方法与线 性插值算法相比,在图像细节和对比度保持方面都要好。这些都证明了,本节所 提的算法既能够去除多个连续相连的环形伪影,且对单个分散的环形伪影校正有 效。 根据公式 4.23 和 4.24,计算图 5.10(b)中白色虚框区域所示的工件细节区域① 的标准差,以及均匀区域②的标准差和信噪比,如表 5.1 所示。从表中可以看出, 采用本节所提的校正算法,校正前后重建图像细节区域的标准差基本没有变化,
60

5 基于改进 Canny 算法的环形伪影校正

说明采样分段 B 样条拟合校正算法不但可以校正多个相连的环形伪影,而且校正 后图像细节和边缘保持良好;对于均匀区域来讲,标准差有微量减小,且信噪比 得到提高。
表 5.1 实验数据评价指标 Table5.1 Experimental data of evaluation indexes 灰度图像 原局部图 本节方法 细节区域标准差 均匀区域标准差 均匀区域信噪比 SNR(db) 19.691 21.873

? σ
16.417 16.491

? σ
17.663 16.181

5.3.3 分析与总结
本节所提的基于 Canny 算法的环形伪影校正算法有以下优点:提出改进的 Canny 算子对投影正弦图中竖直方向线状或带状伪影弱边缘检测,实验结果证明, 提取的弱边缘线型连接程度好,边缘完整,位置准确,同时对于噪声有很好的抑 制作用, 具有较好的边缘定位精度; B 样条函数拟合逼近不但能够保证拟合的曲线 充分逼近原投影数据,且能保证拟合曲线的光滑性,对列均值曲线进行拟合,获 得各个异常探测器投影数据的增益因子,实现多个连续的环形伪影校正;采用分 段拟合校正思想,解决探测器不稳定性致同一个环形伪影亮暗不一致的问题,投 影数据中同一列增益因子不一致的问题;其中,最重要的一点是本节所提的校正 方法不但较好地校正了图像中多个相连的环形伪影,提高信噪比,且校正后图像 细节和边缘保持很好,同时它对于单个分散的环形伪影校正同样有效。 当然,本节所提的环形伪影校正算法也存在不足之处。由于它是从投影正弦 图中识别出不一致探测器所致的环形伪影位置,它的检测效果依赖于伪影的强度, 通常检测和校正了较强的环形伪影后,较弱的环形伪影又表现出来;而且伪影识 别时采用了链长度阈值来对弱边缘进行检测和链接,所以本算法对于最大连接长 度小于链长度阈值的伪影检测无效。 上节所提的基于投影曲线高通滤波的环形伪影校正算法对于单个分散的环形 伪影校正有效,但是无法对多个连续的环形伪影进行校正;而基于改进的 Canny 算法的校正方法实现了多个连续的环形伪影的校正,但是对于最大连接长度小于 链长度阈值 T0 的伪影检测受限。它们各有特点,互相互补,但也不能完全取代对 方。

61

重庆大学硕士学位论文

5.4 本章小结
本节提出一种基于 Canny 算法的环形伪影识别算法,实现投影正弦图中条状 伪影的检测;分析了各行线性插值算法校正不彻底的原因,提出分段 B 样条拟合 校正算法,实现了多个连续的环形伪影校正。

62

6 结论和展望

6 结论和展望
6.1 论文工作总结
根据重建图像环形伪影在投影正弦图中表现为竖直方向直线的特征,本文从 投影正弦图的角度提出两种环形伪影的校正方法。主要工作总结如下: ① 提出了基于投影曲线高通滤波的环形伪影特征识别算法,将二维空间线特 征识别问题转化为一维空间点特征查找。实验证明,该算法实现简单,且对单个 分散的环形伪影在投影正弦图中的识别有效。 ② 分析了传统线性插值对单个分散的环形伪影校正不彻底原因,采用线性外 推插值算法,使环形伪影校正更彻底,同时很好地保持图像边缘细节和分辨率。 ③ 结合传统 Canny 算法对图像边缘检测的优点, 改进了滤波、 梯度幅值计算、 非极大值抑制和边缘检测和链接过程,实现了投影正弦图中条带伪影边缘的检测, 它既适用于多个连续的环形伪影的识别,也对部分单个分散的环形伪影识别有效。 ④ 分析了每行线性插值对多个连续的环形伪影校正的问题, 采用分段 B 样条 拟合算法获得校正增益因子,对伪影处投影数据进行校正。实验证明,该算法对 多个连续的环形伪影校正有效,且很好地保持了图像边缘和细节。

6.2 后续工作展望
本文主要从投影正弦图的角度来对 CT 图像的环形伪影进行校正, 通过实验验 证,这两种算法都能有效地校正环形伪影的校正,且很好地保持图像细节和分辨 率。但是,在研究过程中也发现了该算法存在一些不足之处: ① 基于投影曲线高通滤波的环形伪影特征识别算法,将二维空间线特征识别 问题转化为一维空间点特征查找,使问题得到简化。但是,后续的投影曲线滤波 过程中,采用一阶微分对数据进行锐化处理,致使对多个连续的环形伪影识别受 到限制,有待改进。 ② 采用分段 B 样条拟合算法获得列均值增益因子, 对伪影处投影数据进行处 理时,如果是对多个相连、且刚经过弦线轮廓边界的伪影进行校正时,可能导致 重建图像的工件边界出现“台阶”效应。所以,对于边界处伪影的校正还有待研 究。 ③ 由于本文所提的两种环形校正方法,都是从投影正弦图的角度进行识别和 校正。它的识别效果依赖于伪影的强度,通常检测和校正了较强的环形伪影后, 较弱的环形伪影可能又表现出来。在以后的研究中,可以重点从 CT 系统的设计及 物理原理上来对探测器系统所造成的环形伪影进行校正。
63

重庆大学硕士学位论文

64









在这篇毕业论文完成之时,也意味着我的研究生阶段即将结束。在这三年之 中,我的每一份收获都离不开老师、同学、朋友和亲人的悉心指导和热心帮助, 对于他们的恩情,无以为报,借此处的只言片语来表达我的感激之情。 首先,感谢我的导师王珏教授,感谢他自本科毕业设计以来,对我的精心指 导和悉心关怀,我的课题和论文研究工作无不倾注着王老师的辛勤汗水。这三年 学习和生活中,王老师不但让我学到了更多的专业知识、严谨的治学态度,更学 到了很多做人的道理,在此我向王老师表示深深的敬意和感谢! 感谢我的指导教师蔡玉芳老师和邹永宁老师,感谢他们在我课题的研究过程 中给予的悉心指导和热情帮助,在此对蔡老师和邹老师表示由衷的感谢! 感谢自动化学院的领导和老师,感谢他们在这三年提供给我的学习环境,尤 其感谢年级辅导员温勤老师对我生活上的关心,在此表示诚挚的谢意! 感谢 ICT 研究中心提供的实验平台和学习环境,感谢中心领导、老师、员工 给我的关怀与帮助! 感谢我所有的同门师兄姐弟妹,特别是 ICT 研究中心 401 学习室的所有同学 在三年学习旅途中的陪伴、关心和帮助!同时,感谢我的朋友安传利在学习和生 活上的支持! 感谢我的父母及兄弟姐妹,感谢他们在求学之路上给予的支持和鼓励,因为 他们无私奉献,我才能顺利地完成二十年的学业,在此对我父母的养育和亲戚的 关心表示衷心的感谢! 最后,衷心地感谢在百忙之中评阅论文和参加答辩的各位专家、教授!

黄苏红
二 O 一一年四月 于重庆

65

重庆大学硕士学位论文

66

参考文献

参考文献
[1] [2] [3] [4] 张朝宗, 郭志平. 工业 CT 技术和原理[M]. 北京: 科学出版社, 2009. 谢强. 计算机断层成像技术[M]. 北京: 科学出版社, 2005. 张俊哲. 无损检测技术及其应用[M]. 北京: 科学出版社, 2010. 王珏, 黄苏红, 蔡玉芳. 工业 CT 材料密度测量方法研究[J]. 计算机工程与应用, 2010, 46(2): 203-205. [5] Jiang Hesieh. Reconstruction bias resulting from weighted projection and iso-center misalignment[J]. SPIE, 1999, 3661: 442-449. [6] G. R. Davis, J. C. Elliott. X-ray microtomography scanner using time-delay integration for elimination of ring artifacts in the reconstructed image[J]. Nuclear Instruments And Methods in Physics Research A, 1997, 394:157-162. [7] 李俊江, 路宏年, 李保磊. X 射线图像增强器像元响应不一致性的分析及校正[J]. 光学技 术, 2006, 32(5): 779-784. [8] 王苦愚, 张定华, 黄魁东, 等. 一种锥束 CT 中平板探测器输出图像校正方法[J]. 计算机 辅助设计与图形学学报. 2009, 21(7): 954-961. [9] 梁丽红, 路宏年. 射线面阵探测器成像系统校正研究[J]. 光子学报, 2004, 33(10): 1277– 1280. [10] [11] 夏雄军. 辐射图像噪声及 CT 环形伪影消除[D]. 北京: 清华大学, 2007. 李俊江, 胡少兴 , 李保磊, 等. CT 图像环状伪影校正方法 [J]. 北京航空航天大学学报 , 2007, 33(11): 1382-1398. [12] Jan Sijbers, Andrei Postnov. Reduction of ring artifacts in high resolution micro-CT reconstructions[J]. Physics in Medicine and Biology, 2004, 49(14): 1-8. [13] Yang Jun, Zhen Xin, Zhou Linghong, et al.. Geometric correction for cone-beam CT reconstruction and artifacts reduction[C]. 2nd International Conference on Bioinformatics and Biomedical Engineering, shanghai, P.R. China: IEEE, 2008: 2386-2389. [14] Chen Yenwei, Duan Guifen. Independent component analysis based ring artifact reduction in cone-beam CT images[C]. IEEE International Conference on Image Processing, November 7-12. 2009: 4189-4192. [15] 王珏 , 黄苏红 , 蔡玉芳 . 工业 CT 图像环形伪影校正 [J]. 光学精密工程 , 2010, 18(5): 1226-1233. [16] Kowalski G. Suppression of ring artifacts in CT fan-beam scanners[J]. IEEE Trans On Nucl Sci, 1978, NS 25(5):1111-1116.
67

重庆大学硕士学位论文

[17]

Sun Haining, Qiu Shaokun, Lou SHanshan, et al.. A correction method for nonlinear artifacts in CT imaging[C]. Proceedings of the 26th Annual International Conference of the IEEE EMBS, San Francisco, P.R.USA: IEEE, 2004, 26 II:1290–1293.

[18]

A. Matani, K. Terakawa. Artifact reduction filtering method for CT images[J]. Proceedings of the Second Joint EMBSBMES Conference. HOUSIM, TX. USA October 23-26. 2002: 1035-1036.

[19]

Carsten Raven. Numerical removal of ring artifacts in microtomography [J]. Review of Scientific Instruments, 1998, 69(8):2978-2980.

[20]

Beat Mü nch, Pavel Trtik, Federica Marone, et al.. Stripe and ring artifact removal with combined wavelet-Fourier filtering [J]. Optics Express, 2009, 17(10):8567-8591.

[21]

Rainer Raupach. Method for removing rings and partial rings in computed tomography images [P]. United States Patent: US 6819734 B2, Nov 16, 2004.

[22]

余晓锷 , 罗君芳 , 陈武凡 . 基于弦图的 CT 图像环形伪影校正 [J]. 第四军医大学学报 , 2009, 30(3):207-209.

[23] [24]

罗君方. 去除 CT 图像环形伪影算法研究[D]. 广州: 南方医科大学, 2009. Sun Shaohua, Zhang Li, Gao Wenhuan. A correction algorithm for detector nonconformity in computed tomography [C]. Nuclear Instruments & Methods in Physics Research A, 2003, 505:552-555.

[25]

李 保 磊 , 杨 民 , 傅 健 , 等 . 两 种 CT 成 像 环 状 伪 影 校 正 方 法 [J]. 光 学 学 报 , 2009, 29(7):1849-1853.

[26] [27] [28] [29] [30] [31] [32] [33]

庄天戈. CT 原理和算法[M]. 第一版.上海: 上海交通大学出版社, 1992:1-99. 曹思远. 锥束 CT 三维重建算法加速技术研究[D]. 重庆: 重庆大学, 2009. 谭采云. 第三代 X 线 CT 环形伪影分析[J]. CT 理论与应用研究, 1992, 1(3): 37-40. 朱险峰, 何希儒. CT 环形伪影原因分析与处理[J]. 牡丹江医学院学报, 2001, 22(2): 56-57. 庞彦伟, 王召巴, 林敏. CT 图像环形伪影分析[J]. 华北工学院学报, 2001, 22(1): 16-19. 谭采云. 第三代 CT 环形伪影的分析[J]. 中国医学装备, 2005, 2(5): 49-50. 刘宾. 高分辨率 ICT 重建技术研究[D]. 太原: 中北大学, 2006. 刘宾, 潘晋孝, 王明泉. 基于错位检测的高空间分辨率 ICT 重建技术研究[J]. 科技情报 开发与经济, 2006, 26(4): 173-174.

[34]

黄战华, 蔡怀宇, 张以谟, 等. 基于半像素错位的多幅图像重建高分辨率图像技术研究 [J]. 光学技术, 2002,28(5): 389-394.

[35]

Klaus J. Engel, Lothar Spies, Gercon Vogtmeier. Impact of CT detector pixel-to-pixel crosstalk on image quality[J]. Progress in Biomedical Optics and Imaging - Proceedings of SPIE, v 6142 II, 2006.
68

参考文献

[36]

傅健 , 路宏年 . 扇束 X 射线 ICT 中环状伪影的一种校正方法 [J]. 光学精密工程 , 2002,10(6): 542-546.

[37]

Sanjit K.M. Digital signal processing: A computer-based approach[C]. 2nd ed NewYork : MeGraw-Hill.2001.

[38] [39]

Rafael C. Gonzalez, Richard E. Woods. 数字图像处理[M]. 北京: 电子工业出版社, 2008. Canny J.A Computational Approach to Edge Detection[J].IEEE Tram PAMI, 1986, 8(6): 679-698.

[40] [41]

戴艳. 图像边缘检测与应用[D]. 西安: 西安科技大学, 2009. 张震, 马驷良, 张忠波, 等. 一种改进的基于 Canny 算子的图像边缘提取算法[J]. 吉林大 学学报, 2007, 45(2): 244-248.

[42]

Jun Li, Sheng Ding. A Research on Improved Canny Edge Detection Algorithm[C]. 2010 3rd International Conference on Computational Intelligence and Industrial Application (PACIIA):68-70.

[43] [44]

吴云. 光纤面板传像暗影的自动检测[D]. 大庆: 大庆石油学院, 2008. Gadallah F L, Csillag F, Smith E J M. Destriping multisensor imagery with moment matching[J]. Int. J. Remote Sentsing, 2000, 21(12): 2505-2511.

[45]

刘正军, 王长耀, 王成. 成像光谱仪条带噪声去除的改进矩匹配方法[J]. 遥感学报, 2002, (4): 280-284.

[46]

王先培,张爱菊,李少雄. 基于非均匀 B 样条曲线的红外数据的精确拟合及校正[J]. 光 谱学与光谱分析. 2006, 26(10): 2850-2853.

[47] [48]

张玲. 基于三次样条曲线拟合公路平面线形方法研究[D]. 武汉: 武汉理工大学, 2007. 梁芳. 遗传算法的改进及其应用[D]. 武汉: 武汉理工大学, 2008.

69

重庆大学硕士学位论文

70









A 作者在攻读硕士学位期间发表录用及投稿的论文目录
[1] 王珏, 黄苏红, 蔡玉芳. 工业 CT 材料密度测量方法研究. 计算机工程与应用, 2010, Vol.46, No.2, 2010. (CSCD 检索,已发表) [2] 王珏, 黄苏红, 蔡玉芳. 工业 CT 图像环形伪影校正. 光学精密工程, Vol.18, No.5, 2010. (EI 检索,已发表) [3] 王珏, 黄苏红, 蔡玉芳. 基于改进的 Canny 算法的 CT 图像环形伪影校正. 光学精密工程. (审稿中)

B 作者在攻读硕士学位期间参加的课题与基金项目
[1] 黄苏红,王慧倩,陶李,伍立芬,徐维. 基于 CT 图像几何尺寸定量测量的方法研究. 重 庆大学研究生科技创新基金(CDJXS11171159).

71

重庆大学硕士学位论文

72

CT图像环形伪影校正方法研究
作者: 学位授予单位: 黄苏红 重庆大学

引用本文格式:黄苏红 CT图像环形伪影校正方法研究[学位论文]硕士 2011


相关文章:
CT图像伪影的分析和处理
CT图像伪影的分析和处理_医药卫生_专业资料。CT 图像伪影的分析和处理 无论哪一代 CT 都是由硬件和软件两部分组成。硬件是躯壳,软件是灵魂, 硬件的工作必须依靠...
锥束CT图像伪影和校正方法综述
一种快速校正CT环形伪影... 4页 4下载券 工业CT图像环形伪影校正 8页 免费 ...[8] 魏英,李春云.CT 图像“条状”伪影矫正方法研究[J].2007,33(1) [9] ...
CT图像的运动伪影校正
CT图像的运动伪影校正_其它考试_资格考试/认证_教育专区。Sample 20mm 30mm 光固化快速成型过程零件变形的数值模拟武殿梁,丁玉成,卢秉恒 西安交通大学先进制造技术研究...
三代CT环形伪影分析
一种快速校正CT环形伪影... 4页 4下载券 CT 图像伪影分析 59页 1下载券 多层...参考文献: 1)《CT 理论与应用研究》第一卷第三期 作者简介: 谭彩云: 北京...
CT图像伪影的分析与处理
【关键词】 ct 图像;伪影;处理 pseudomorph analysis and disposal of ct ...1.2 分析与处理 利用计算机将存储于磁带上的标准水模的原始数据和空气校准数据...
CT图像伪影的分析与处理
CT图像伪影的分析与处理_基础医学_医药卫生_专业资料...水模的原始数据和空气校准数据重建(相当于 正常采集...利用排除法是判定故障的有效方法。 综上所述,由于...
影工影像设备学06-07-2第二套
校正探测器通道的非线性和 X 线硬化效应引起的非效应,消除环形伪 影,提高图像...二、填空题(每空 1 分,共 10 分) 1.CT 环行伪影常见与第___代 CT,主要...
CT图像去伪影和简单重建
(1) 掌握 CT 图像伪影形成原因和特点; (2) 掌握去伪影和重建理论方法研究与...基于插值方法去除CT图像... 4页 免费 工业CT图像环形伪影校正 8页 免费 CT ...
CT设备原理
1917 年,奥地利数学家 Radon 提出了图像重建理论的数学方法 。他指出对二维或...并从软件上解决了环形伪影的产生条件和配置了软件校正措施,第四代 CT 则因探测...
影工影像设备学07-08-2第二套
2.GE8800 CT 主、副可控硅触发电路分析? 答:A1、A2 117V 交流电→59T3→...(6)环形校正 消除环形伪影 5.画图和文字说明永磁射频平面板线圈产生射频磁场的...
更多相关标签:
ct图像伪影 | ct图像伪影分类 | ct伪影 | ct金属伪影 | 图像伪影 | 图像伪影的滤除 | 图像环形内存 | 图像畸变校正 |