基于顶帽运算和各向异性扩散的CT图像金属伪影去除

牛晓伟, 孔慧华, 邸云霞

(中北大学 数学学院 山西省信息探测与处理重点实验室, 山西 太原 030051)

0 引言

自1895年伦琴发现X射线以来,计算机断层扫描技术(computed tomography,CT)已广泛用于临床诊断.现代CT扫描仪可以提供具有可靠解剖信息的高分辨率图像,然而,如果扫描场中存在金属植入物(如牙科填充物或修复体),金属伪影往往会影响CT图像质量并降低其诊断价值[1].导致金属伪影产生的原因可归纳为射束硬化、部分容积效应、噪声、散射和运动等[2],这些问题会导致检测器接收到的投影数据损坏或丢失.在图像重建中,损坏的数据表现为或亮或暗的伪影,除此之外,金属植入物的大小、形状、位置也会影响金属伪影的情况,这会降低CT图像的质量,影响金属植入物周围组织结构的诊断信息.

近年来,越来越多的金属伪影去除(metal artifact reduction,MAR)方法用来校正金属伪影,它们大致可分为迭代重建、正弦图校正和混合方法[3].迭代重建和混合方法在抑制金属伪影方面表现出良好的性能,其中最大似然和代数重建技术是两种最常见的方法.如Kubo等[4]比较了基于模型的迭代重建和混合迭代重建与金属伪影去除算法相结合时改善牙齿金属伪影CT图像的能力,发现混合迭代重建算法显著减少了颈部CT图像的金属伪影.但由于迭代算法在临床应用中受到限制,计算成本高,耗时严重等问题,因此,在临床实践中应用有限.

在正弦图校正方法中,大致可分为两组:基于插值的方法和基于先验的方法.基于插值的方法使用插值技术来估计金属投影区域的替代数据,例如多项式插值、小波插值或线性插值等[5-7].线性插值(linear interpolation,LI)是最常用的方法,Kalender等[6]在1987年首次提出了关于线性插值的金属伪影去除算法,使用金属轨迹两侧未损坏的数据对金属投影区域进行线性插值,这样往往会造成组织信息的丢失,而且可能会产生严重的二次伪影.为了克服基于插值方法的缺点,基于先验的方法使用来自先验图像的前向投影来估计替代数据,先验图像的结构信息可以添加到金属投影区域,这样可以获得更好的估计[8-15].Wang等[8]提出了一种基于融合先验的MAR(fusion prior-based MAR,FP-MAR),将滤波后的线性插值图像和去除金属部分的原始重建图像融合生成一个先验图像,保留了更多的组织信息,但会残留白色伪影.Meyer等[9]基于初始图像的多阈值分割获得的三值图像的正向投影来执行归一化,提出了归一化金属伪影去除(normalized metal artifact reduction,NMAR)算法,以便更好地保持空气和水等物体之间的对比度,对于带有金属植入物的骨骼结构有相当大的改进.Yakdiel等[12]利用Mumford-Shah图像模型和L0梯度最小化方法,通过高斯滤波、Parisotto和Schoenlieb修复方法三种变换获得先验图像,可以有效减少条纹金属伪影,避免金属植入物周围出现新伪影.

由于一般基于先验的方法通常会在最终校正图像中造成伪影残留和金属附近骨骼结构丢失,如基于线性插值的NMAR算法等,无法取得理想的校正效果.因此,本文提出一种基于顶帽运算和各向异性扩散[16]的归一化金属伪影去除算法(top hat and anisotropic diffusion-NMAR,THAD-NMAR).该算法对CT图像执行顶帽运算并与线性插值图像进行融合,可以保留大部分骨骼信息和细节特征,且各向异性扩散滤波[16]具有保留图像边缘细节同时减少噪声的作用,由此可以获得质量更好的先验图像.临床图像测试结果表明,与现有的常规图像域金属伪影去除算法相比,本文所提出的THAD-NMAR方法在抑制伪影和组织特征保留方面表现更好.

1 理论与方法

含金属物体的重建图像常常表现为强烈的黑色或白色金属伪影,或像骨骼一样明亮,或像空气一样暗,经常会出现分割错误的情况,严重的伪影会覆盖感兴趣区域,干扰正确诊断.而先验图像可以在一定程度上弥补常规算法的不足,因此先验图像的质量在基于先验的MAR中至关重要,为了克服由线性插值引起的伪影残留和组织结构丢失的问题,本文使用了一种数学形态学算法.主要思想是通过对CT图像执行顶帽运算来恢复因线性插值方法丢失掉的骨骼信息,然后使用各向异性扩散滤波去除残留伪影,从而得到质量更好的先验图像.本文算法主要分为以下四个步骤.

1.1 金属分割

当被扫描物体中包含多个金属材料时,会导致重建的CT图像中出现严重的金属伪影,如图1(a)所示,黄色框为黑色伪影,红色框为白色伪影.首先需要预先从图像中分割出这些金属来得到金属轨迹,而阈值分割法具有速度快的优点,能准确识别出金属部分,可以满足大多数CT图像的金属分割.由于CT图像中不同的组织具有不同的CT值,例如空气为-1 000 HU左右、脂肪为-90～-70 HU、水为0 HU左右、软组织为20～40 HU,骨骼为300～1 000 HU,而金属的CT值远远大于2 000 HU,因此将分割的阈值选择为2 500 HU,可以得到CT图像的金属部分,对金属部分前向投影得到金属轨迹.

图1 黑色伪影校正图像

1.2 分类线性插值

在投影域中对原始图像投影所有角度下的金属轨迹部分使用线性插值,公式为:

(1)

式(1)中:gβ为角度β下校正后的投影值,[mβ,nβ]为该角度下的金属投影所在的区间,即插值区间,坐标α为对应的插值点,且α∈[mβ,nβ].然后对线性插值的投影重建得到线性插值图像Ili,线性插值图像虽然可以有效去除白色和黑色金属伪影,但会引入严重的二次伪影,造成金属植入物骨骼丢失的情况,如图1(b)所示.通过下面公式融合可以得到无黑色伪影图像:

Ino_black=Iorign+Iblack

(2)

其中,

(3)

式(3)中:Iblack是黑色伪影的补偿信号,Iorign是原始图像,通过式(2)可以得到无黑色伪影图像Ino-black.无黑色伪影图像通过添加补偿信号只是去除了原始图像的黑色伪影,并没有校正白色伪影,可以保留原始图像中的所有细节特征,如图1(d)所示.

1.3 先验图像生成

1.3.1 基于顶帽运算的骨骼修复

MAR中的插值算法常常会导致最终校正图像中金属附近的组织结构丢失.数学形态学是图像处理过程中一种非常重要的方法,可以提取图像中的分量信息,包括膨胀、腐蚀、开运算和闭运算四种基本运算,开运算可以去除相对于结构元素较小的明亮细节,保持整体的灰度和较大的明亮区域不变,顶帽运算正是用来提取这些被移除的部分,选择不同的结构元素会导致不同的分割结果,即提取出不同的特征.对无黑色伪影图像Ino-black执行顶帽运算,可以恢复丢失掉的骨结构.顶帽运算的定义为:

假设原始灰度图像为I,使用半径大小为d的圆盘结构元素z进行顶帽运算,公式如下:

Itophat=I-(I∘z)

(4)

其中,

I∘z=(I⊖z)⊕z

(5)

式(4)中:Itophat是顶帽运算提取到的信息,⊖表示腐蚀,⊕表示膨胀,结构元素z的大小d依据骨骼信息的情况而定.式(5)为开运算,用∘表示,开运算是对图像先腐蚀,再膨胀.腐蚀的公式为:

I⊖z=min[I(x+i,y+j)-z(i,j)]

(6)

式(6)中:(x+i,y+j)是定义在目标I上的坐标,(i,j)是定义在结构元素z上的坐标.膨胀的公式为:

I⊕z=max[I(x+i,y+j)+z(i,j)]

(7)

先验图像的质量在基于先验的MAR中起着至关重要的作用,而传统的插值方法往往会由于分割错误,在插值过程中丢失骨骼等组织结构,从而影响到图像质量.通过顶帽运算可以恢复丢失的骨骼信息,从而获得更好的先验图像.

将无黑色伪影图像减去开运算的结果,可以提取到图像的骨骼部分,这是由于伪影的CT值一般小于300 HU,而骨骼的CT值为300～1 000 HU,因此选择阈值200 HU可以去除其他细节(包含噪声).然后对线性插值图像和无黑色伪影图像进行融合,融合过程为:

(8)

式(8)中:Itophat表示CT图像顶帽运算的结果,如图2(a)所示,Ifusion表示融合后的图像,它会将线性插值丢失的骨骼信息恢复,同时去除图像中的大部分伪影,如图2(b)所示.

图2 顶帽运算骨骼恢复图像

1.3.2 各向异性扩散滤波

融合后的CT图像往往会残留一部分伪影,因此需要进一步处理.传统的线性滤波算法在去噪时通常会牺牲掉图像的边缘信息,而基于非线性偏微分方程的各向异性扩散方法会弥补这个缺陷,具有保留图像边缘细节同时减少噪声的作用,可以增强图像边缘的对比度.Perona和Malik给出的各向异性扩散方程(P-M方程)的离散表达式为[16]:

(9)

(10)

对融合后的图像Ifusion执行各向异性扩散滤波,可以得到先验图像Iprior,如图3所示.

图3 各向异性扩散滤波图像

1.4 归一化金属伪影去除

基于原始CT图像、先验图像和金属部分的投影,执行Meyer等[9]提出的归一化和反归一化操作,对受金属影响正弦图进行校正,过程如下:

(11)

式(11)中:M(·)表示投影的线性插值,R表示拉东变换,Pcorrect表示最终校正的正弦图.最后通过滤波反投影算法重建得到最终的校正图像Icorrect.

2 实验结果与分析

2.1 实验设置

为了验证本文算法在去除金属伪影方面的有效性,使用临床CT图像进行测试.本文选择了6个临床病例来测试所提出的THAD-NMAR算法在实际临床中应对不同种类金属伪影的能力.所有原始未校正金属伪影图像均来自REVISION RADIOLOGY(http://www.revisionrads.com/cgi-bin/MDT_cgi).图4中包含了三组金属结构简单的图像,图5包含了三组金属结构较为复杂、伪影严重的图像.本文选取线性插值法、NMAR、LI-NMAR、金属删除技术(metal deletion techniq-ue,MDT)和FP-MAR算法作为比较算法,其中FP-MAR算法的参数均设置为原文中的值,本文算法各向异性扩散滤波的参数都是根据经验选择,可以达到最优效果,其中迭代次数t设置为100,热传导系数K设置为15,控制扩散总体强度的λ设置为0.25,可满足所有的实验需求.图像大小均为512×512,CT图像窗位和窗宽均设为C=500/W=3 000 HU.

2.2 实验分析

图4显示的是临床数据金属植入物对比方法的测试结果,案例1是脑动脉瘤植入物,案例2和案例3是深部脑刺激器,且参数d的值均为0.图4是几个金属伪影较轻的情况,所有的方法都可以消除大部分伪影,其中线性插值法在校正后会引入二次伪影,NMAR算法会因为伪影较严重而分类错误最终残留下一部分黑色伪影,LI-NMAR算法较NMAR算法校正效果会好很多.而FP-MAR算法由于是通过融合滤波后的线性插值图像和去金属的原始伪影图像来恢复丢失掉的骨信息,因此会残留下白色伪影,如图4(f)中的案例3所示.本文的THAD-NMAR算法通过有效的伪影抑制和骨骼结构保留获得了最佳的图像质量.

图4 临床数据金属伪影去除方法结果对比

图5是更为复杂的情况,包含了多个不规则或大型金属,其中案例4是胸骨固定螺钉,案例5是单侧髋关节假体,案例6是椎弓根螺钉,d值分别为10、12、35.线性插值法、NMAR算法和LI-NMAR算法虽然去除了大部分的金属伪影,也损失了金属周围的骨骼信息,如图5(b)、(c)、(d)所示,这样的结果不利于诊断.MDT方法有很好的伪影抑制效果,但也会缺少一部分金属附近的组织结果并残留白色伪影,如图5(e)中的案例4.FP-MAR算法在很大程度上抑制了金属伪影和保留大部分的骨骼信息,却会残留下一部分明显的白色伪影,如图5(f)中的案例4和案例6所示,因此具有一定的局限性.显然,本文提出的THAD-NMAR算法在抑制伪影方面优于其它MAR算法,在恢复那些靠近金属物体的组织结构方面也比MDT和FP-MAR算法表现更好.

图5 临床数据金属伪影去除方法结果对比

金属伪影或像骨骼一样明亮,或像空气一样暗,就经常会出现分割错误的情况,最终会残留下一部分黑白伪影,因此拥有更好的先验图像显得尤为重要.LI-NMAR将线性插值图像作为先验图像,由于在金属较多时线性插值图像会造成明显的残留伪影和骨骼缺失,最终的校正图像也会出现骨骼缺失和伪影残留的情况.FP-MAR算法也有一定的局限性,例如它可能会将金属物体周围的一些明亮伪影错误地合并到先验图像中作为骨骼结构,因此可能导致最终校正后残留下白色伪影.而本文的THAD-NMAR算法恰好弥补了这一缺陷,它可以准确地提取出原始图像中的骨骼信息.

由于图4中的MAR算法不会出现骨骼缺失的情况,因此参数d的值均设为零,而图5中的CT图像有较为严重的金属伪影,所以在常规MAR算法校正后会出现伪影残留和骨骼信息丢失的情况,依据不同的情况需要设置不同的d值.图6是从无黑色伪影图像中提取到的骨骼信息,可以看到,不同的d值会提取到不同的信息,当d值为10和15时可以提取到程度相同的骨骼和金属部分.d值过小,提取到的骨骼信息可能会不完整,造成组织信息丢失;d值过大,会提取到原始CT图像中的白色伪影和其他信息,进而造成伪影残留.因此需要手动设置选取合适的d值,当骨骼少或较细时d值应小一点,当骨骼较多或复杂时d值要大一点.

图6 案例4顶帽运算不同d值提取到的骨骼信息

图7和图8显示了沿案例4原始CT图像中第212行(红色线条)和296行(黄色线条)得到的各种算法校正图像的强度分布曲线,用来对比黑色伪影和白色伪影的校正情况.从图7可以看到,当前黑色伪影较为严重,每种算法都可以较好的去除黑色伪影,保持正常的CT值.

从图8可以看到,当前白色伪影较为严重,CT值远远大于正常值,除FP-MAR算法外都可以很好地校正骨骼附近的白色伪影,但LI算法、NMAR算法、LI-NMAR算法会造成骨骼丢失的情况,如图8(a)中的红色虚线框所示,而THAD-NMAR算法不仅将CT值恢复到了正常的CT值附近,而且可以保留骨骼等重要信息,因此本文算法在去除白色伪影和保留骨骼信息方面要优于其他算法.

图7 案例4第286行校正CT图像的强度分布

图8 案例4第212行校正CT图像的强度分布

3 结论

CT图像中的金属伪影一直是影响CT图像质量的关键因素之一,这些金属伪影不仅模糊了CT图像,而且影响了疾病的诊断效果,甚至影响了CT成像系统的进一步发展.而一般基于先验的金属伪影去除方法通常会在最终校正图像中造成伪影残留和金属附近骨骼结构丢失,为了解决这个问题,本文提出了一种基于顶帽运算和各向异性扩散滤波的归一化金属伪影去除算法.

通过顶帽运算来恢复因线性插值引起的组织信息丢失问题,然后通过各向异性扩散滤波可以得到先验图像,最后基于原始CT图像投影、先验图像投影和金属轨迹执行归一化和反归一化操作得到最终的校正图像.顶帽运算通过调整参数d可以恢复丢失掉的骨结构,而各向异性扩散滤波具有保留图像边缘细节同时减少噪声的作用.

与其他经典的MAR算法相比,在存在多个或较为复杂的金属植入物时,本文算法可以去除CT图像中强烈的暗和亮的金属伪影,不会引入新的伪影,并且可以完整地保留金属周围的组织结构,适合较大金属植入物周围存在骨骼等结构的情况.