体层合成的各向异性扩散滤波重建算法

俞龙江，戴道清，邹鲁民

（1.中山大学 数学与计算科学学院，广东 广州510275；2.珠海友通科技有限公司 医疗信息系统部，广东 珠海519085）

0 引 言

体层合成 （tomosynthesis）［1］是一种有限角度图像重建技术，其成像角度通常在60度以内，与传统CT 相比剂量明显减小，能够极大减少病人所受的辐射剂量和检查所需的时间，降低医疗设备成本，提高了医生的诊断效率，具有广泛的应用前景，能产生巨大的经济效益并具有良好的社会效益。目前体层合成重建技术已经开始用于国外的高端X 线设备上，该技术仅被少数公司所掌握，尚未普遍应用，但是这项技术可以应用的范围非常广泛，在牙科检查［2］、关节 检 查［3］、乳 腺 癌 早 期 筛 查［4］、胸 部 检 查［5］、放疗定位［6］等方面都可发挥CT 无法替代的作用。然而，这种有限角度的成像结构会产生切片间的伪影，这种伪影会影响阅片时对病灶的判断，从而影响诊断的正确率。传统的滤波后投影重建算法无法消除切片间伪影，而目前迭代类型的重建算法难以保证实时性，尚未得到普遍应用。

1 各向异性扩散滤波伪影消除算法

1.1 滤波后投影重建算法

自从1972年体层合成由Grant提出以来，针对这种成像结构提出的各种重建算法，主要分为滤波反投影算法和迭代算法两大类［7］。滤波反投影算法主要是在平移－叠加的后投影操作基础上，再利用频率域斜坡滤波来提高重建图像的对比度［8］。然而，频率域斜坡滤波在放大高频信号以达到提高图像对比度的同时，对有限角度投影产生的伪影也进行了放大，特别是在体层合成应用中，由于投影角度的数目与360度相比非常小 （通常在60度以内），经滤波后投影重建后产生的伪影尤为明显，严重影响了重建图像的质量。而迭代重建算法收敛速度较慢，需要上百次迭代才能达到较好的重建效果，不能满足体层合成重建应用的实时性要求，故目前滤波后投影算法依然是提供体层合成功能的医疗设备主流厂商的选择，如Siemens［9］、GE［10］、Philips［11］、岛津［12］等。如何抑制体层合成产生的伪影，是各大厂家秘而不宣的申请专利保护的技术，这固然保护了厂家的利益，但也阻碍了体层合成技术的应用与推广。

从实际应用角度考虑，本文选择滤波后投影类型的算法进行图像重建，通过将频率域斜坡滤波器替换为可抑制体层合成产生的伪影的滤波器，从而满足体层合成重建的图像质量要求。接下来对本文提出的重建算法进行具体阐述。

重建算法的具体构造取决于成像系统几何结构，下面分别给出锥束CT 与体层合成的成像几何结构图进行对比，通过对比来演示体层合成的成像特点及与锥束CT 的异同点。图1是锥束CT 的成像几何结构。

图1 锥束CT 成像几何结构

图2是体层合成的成像几何结构。

图2 （a）是体层合成系统图像采集方式的原理演示图，图2 （b）是对应的运动坐标系演示图。

由图1与图2对比可见，图1所示的锥束CT 成像结构，成像物体绕轴旋转，而从成像物体坐标系来看，可把X 线光源与探测器视为绕轴旋转，旋转轴与运动轨道平面垂直。图2 （a）所示的是一种常见的体层合成的成像结构，成像物体静止，球管与探测器绕物体中心轴旋转，旋转范围从θmax到－θmax，从图2 （b）所示的成像物体坐标系来看，旋转轴与运动轨道平面也是垂直的。因此，从成像几何结构角度来看，可以将体层合成视为有限投影角度范围（2倍θmax）锥束CT，且旋转轴方向从图1中沿纸面向上变换到图2 （b）中从纸内指向读者。

图2 一种常见的体层合成的成像几何结构

根据上述分析得出的结论，即体层合成可看作有限投影角度范围的锥束CT，故锥束CT 上的滤波后投影算法同样适用于体层合成图像重建。根据体层合成的成像结构稍作调整，得到的滤波后投影算法具体如下：

（1）旋转轴变换。旋转轴方向从图1中沿纸面向上变换到图2中从纸内指向读者。

（2）后投影。采用后投影算法所遵循的锥束CT 几何结构进行后投影操作。

（3）滤波操作。由于滤波后投影算法中滤波与后投影操作的顺序是可以调换的，故可采用先进行后投影再滤波的重建算法，这里将采用2个滤波器，分别是各向异性扩散滤波器和切片厚度滤波器，前者用于抑制体层合成由于投影角度不足而产生的伪影；后者用于调整切片厚度，使切片切换时呈现自然的视觉过渡。

其中，切片厚度滤波器操作在垂直于切片平面的坐标轴方向，目的是对切片厚度进行调整，使感兴趣区域清晰呈现于特定切片内，且切片间的过渡要足够平滑，不至于使人眼在阅片时感觉变化太突然。切片厚度滤波器可以选择窗函数进行，例如Siemens方法的汉宁 （Hanning）窗函数［9］，以及岛津方法的高斯窗函数［12］。其中，高斯窗函数参数设置简单，且能提供切片间的平滑过渡，故本文选择高斯窗函数进行切片厚度调整。

各向异性扩散滤波器是本文提出的算法核心，下面进行具体阐述。

1.2 各向异性扩散滤波算法

在体层合成重建的三维图像中，感兴趣的对象应能在其对应的某一个或某几个切片图像中清晰可见，而在其它切片图像中模糊不清甚至是完全不可见。然而，体层合成具有的有限角度投影的特点，由其本身的图像采集几何结构决定了切片间的伪影不可能完全消除，即在某一切片图像中应该不可见的对象，在该切片图像中产生伪影。这种伪影会遮挡该切片图像中本应清晰呈现的对象，造成模糊，还会使本应在其它切片图像中的对象的轮廓呈现在该切片图像中，影响对该对象的准确位置的判断。切片间伪影来源于体层合成的图像采集结构，无论采用滤波反投影重建算法还是迭代重建算法，都会在重建结果中出现这种伪影。

为消除这种切片间伪影，本文提出采用三维各向异性扩散滤波器对重建体数据进行伪影消除。使用该滤波器的目的是消除体层合成的切片间伪影，并在滤波操作后保持对象的对比度不会降低。各向异性扩散滤波器的扩散方程通 常 为［13］

式中：u——三维体素空间，t——扩散时间，D——扩散张量。对式 （1）的扩散方程进行迭代前向差分近似得到

式中：k——迭代次数。扩散张量D 是对图像进行结构张量的特征分解得到的，结构张量J 由下式得到

式中：Kρ——一个高斯加权函数，ρ——该高斯函数的参数σ 值，在使用中一般将其设置为1。对结构张量J 的特征分解得到

式中：λ——特征分解得到的特征值，v——特征值λ 对应的特征向量。从而根据式 （4）得到扩散张量D 为

式中：D——对称矩阵，满足

将式 （6）代入式 （1），可得到等式右边的散度运算为

式中：ix、iy、iz——x、y、z方向上的通量

体层合成图像重建的要求是消除体层合成设备在图像采集过程中产生的切片间伪影，这意味着滤波器要有一定平滑能力，而同时要对滤波操作造成的平滑进行增强，这又意味着滤波器要有一定的增强能力。综上，选择边缘增强型各向异性扩散滤波器，其特征值满足［14］

由于切片间伪影在切片图像中表现为一种模糊的伪影，并随着切片图像数目增加而缓慢衰减，因此各向异性扩散滤波器根据其滤波器原理可通过扩散方式，将这种模糊伪影通过切片图像内和切片图像间的扩散，将其逐渐减弱，从而达到消除切片间伪影的目的。同时由于选择的各向异性扩散滤波器是边缘增强型的，故不会因为滤波操作而造成对象边缘的平滑，保持了滤波后图像的对比度。

2 仿真实验

在X 线体层合成图像重建系统设计和开发过程中，通常要进行一系列试验来测量系统各项性能指标，衡量各种影响图像重建质量的因素并验证设计的有效性。由于影响X 线体层合成图像重建的因素较多，整个试验过程非常耗时且成本极高，因此图像重建算法的开发阶段采用计算机仿真技术，可避免各种扰动因素的叠加，从而降低处理难度，加速研究工作的进展。这种仿真技术是根据体层合成成像系统的物理模型，并对其加入各种影响成像质量的因素模型进行数字仿真，例如扫描机构的机械运动不稳定、采集系统几何结构的偏差、散射噪声、射束硬化伪影、金属伪影等。本文在实验中采用珠海友通公司开发的体层合成数字仿真平台，生成符合体层合成图像采集结构的投影数据，进行重建算法的设计。仿真实验中的体层合成系统的图像采集及重建参数见表1。

仿真实验采用一个数字生成的体模，模拟人胸腔及其内部的肺部结节，如图3所示。

表1 体层合成仿真实验参数设置

图3 人胸腔及肺部结节的数字仿真体模

对人胸腔及肺部结节的数字仿真体模按表1所示的采集及重建参数进行投影，对生成的投影数据进行图像重建，重建算法选择常用的滤波后投影算法，得到如图4所示的重建结果。

图4 滤波后投影重建得到的切片图像

图4是重建结果中的一张切片图像。从图4 中可见，该切片图像中呈现出3个圆拱形的伪影，是体模中的肋骨在该切片图像上形成的切片间伪影。在没有先验知识的情况下，无法知道这是伪影，还是组织本来就应呈现的影像。

在使用本文提出的算法进行实验时，需要对算法参数对性能的影响进行分析，以便正确设置算法参数。

各向异性扩散滤波算法的参数主要包括式 （1）中的扩散时间t和式 （9）中的控制增强对比度的参数λe，参数设置及原则如下：

扩散时间t是控制各向异性扩散滤波器的平滑程度的参数，随着扩散时间的增加，图像将变得越来越平滑，对象与背景间的对比度会越来越低。根据图像总体对比度选择一个合适的扩散时间，使扩散操作后保证组织对比度满足正常的阅片要求。

参数λe是用于平滑切片间伪影的门限值，梯度小于这个门限值的图像区域，各向异性扩散滤波器的平滑程度随梯度值减小而增大，即这些区域的梯度值越小，平滑程度就越大，而梯度大于这个门限值的图像区域，各向异性扩散滤波器的平滑程度随梯度值增大而减小，即这些区域的梯度值越大，平滑程度就越小，通常存在切片间伪影区域的梯度小于组织区域的梯度，为消除切片间伪影，设置平滑梯度门限值参数λe在切片间伪影区域的梯度和组织区域的梯度数值之间。如果设置该门限值越接近切片间伪影区域的梯度，组织区域受滤波器平滑操作的影响就越小，即正常区域的对比度保持的就越好；如果设置该门限值越接近组织区域的梯度，切片间伪影受滤波器平滑操作的影响就越大，即切片间伪影的消除程度就越大，图像就越平滑。一般选择以保持图像内组织的对比度为前提的情况下进行切片间伪影的平滑，例如在切片图像总体梯度值范围归一到0到1之间后，切片间伪影区域的梯度数值处于0.02以下的范围，则可把该门限值设置为0.05，略大于切片间伪影区域的梯度数值即可。

图5是采用本文提出的算法重建得到的结果。图5是与图4所示的相同位置的切片图像。在实验中，取扩散时间t＝50，梯度门限值λe＝0.02。从图5中可见，本文提出的算法处理后得到的切片图像，与图4所示的传统的滤波后投影重建结果相比，切片间伪影得到显著的消除，图像中的对象与周围背景过渡自然，从而提高了重建图像质量。

图5 各向异性滤波重建结果

为定量地对比切片间伪影的消除性能，这里定义目标可见度这一指标来衡量切片间伪影的消除结果。设切片中某一目标区域内的平均灰度值为O，其周围背景的平均灰度值为B，且满足0＜O＜B＜1，则目标可见度D 表示为

根据体层合成重建原理，在理想情况下，某一目标应在某一特定切片最为清晰，可设在这一切片内的目标可见度为1，随着切片数目增加，该目标在其它切片内应迅速变得模糊，即O 显著变大，B 与O 在数值上迅速接近，使该目标可见度D 迅速衰减为0。而实际重建中，由于存在切片间伪影，该目标可见度D 衰减缓慢，在其它切片中仍可见到该目标的影子，D 的数值越低，衰减越慢，对应的切片间伪影就越强。

下面通过实例来演示目标可见度是如何计算得到的。图6是滤波后投影重建算法和本文提出的各向异性滤波重建算法得到的第60张切片：图6 （a）滤波后投影重建得到的切片，图6 （b）是各向异性滤波重建得到的切片。图6（a）和图6 （b）中位于图像下方用矩形框标出的圆形黑点，是体模中对人肺部结节的模拟，选择该目标进行目标可见度计算，根据体模的预先设定，该目标应只存在于第60张切片内，故设在第60张切片内的该目标可见度为1。

图6 2个算法重建得到的第60张切片

图7是滤波后投影重建算法和本文提出的各向异性滤波重建算法得到的第100张切片：图7 （a）滤波后投影重建得到的切片，图7 （b）是各向异性滤波重建得到的切片。图7 （a）和图7 （b）中位于图像下方用矩形框标出的区域，是选择的肺结节目标对应的伪影区域，在这2个矩形区域中分别按式 （10）计算目标可见度D，同理，分别对2个算法重建得到的其它切片计算目标可见度，得到该目标从第60张切片到第100张切片对应的目标可见度结果，如图8所示。

图7 2个算法重建得到的第100张切片

图8是按式 （10）分别计算滤波后投影算法与本文提出的各向异性扩散滤波算法的目标可见度，以对比2个算法对切片间伪影的消除结果。

图8 伪影可见性对比结果

从图8中可见，横坐标是从第60张切片到第100张切片的切片序数，纵坐标是目标的可见度，该目标在重建结果中应只存在于第60张切片中，因此在第60张切片的目标可见度为1，随着切片数的增加，该目标的可见度应迅速衰减为0。图6中实线为滤波后投影重建算法得到的结果，从图6中可见，随着切片数目的增加，目标可见度的衰减缓慢，直到第100张切片时，该目标的可见度仍有0.3，从而造成该目标在其它切片内仍依稀可见，存在明显的切片间伪影。相比而言，虚线所示的本文提出的算法结果，可随切片数目增加显著地降低该目标的可见度，到100张切片时，目标可见度已衰减到接近0，从而显著减少了切片间伪影。

3 结束语

本文提出用于体层合成的各向异性扩散滤波重建算法，经过实验验证，可显著消除切片间伪影，与传统的滤波后投影重建算法相比，图像质量得到明显改善，从而确保了更好的阅片效果。

目前实验在数字仿真的体模上进行，有利于加入各种影响伪影消除的因素比如噪声、运动模糊、高衰减对象等，来验证算法性能，数字仿真实验同时可对重建算法进行不断完善，使其应用于实际采集图像成为可能。

［1］Karellas A，Lo JY，Orton CG.Point／counterpoint：Cone beam x－ray CT will be superior to digital X－ray tomosynthesis in imaging the breast and delineating cancer ［J］.Medical Physics，2008，35 （2）：409－411.

［2］Cho MK，Kim HK，Kim SS，et al.Development of dental tomosynthesis system ［C］／／IEEE Nuclear Science Symposium，2007：3788－3791.

［3］Kalinosky B，Sabol JM，Piacsek K，et al.Quantifying the tibiofemoral joint space using X－ray tomosynthesis［J］.Medical Physics，2011，38 （12）：6672－6682.

［4］Sahiner B，Chan HP，Hadjiiski LM，et al.Computer－aided detection of clustered microcalcifications in digital breast tomosynthesis：A 3D approach ［J］.Medical Physics，2012，39（1）：28－39.

［5］Wang J，Dobbins JT III，Li Q.Automated lung segmentation in digital chest tomosynthesis［J］.Medical Physics，2012，39（2）：732－741.

［6］Brunet－Benkhoucha M，Verhaegen F，Lassalle S，et al.Clinical implementation of a digital tomosynthesis－based seed reconstruction algorithm for intraoperative postimplant dose evaluation in low dose rate prostate brachytherapy ［J］.Medical Physics，2009，36 （11）：5235－5244.

［7］Baker JA，Lo J.Y.Breast tomosynthesis：State－of－the－art and review of the literature ［J］.Academic Radiology，2011，18 （10）：1298－1310.

［8］Gengsheng Lawrence Zeng.Medical Image Reconstruction：A conceptual tutorial［M］.NY：Springer，2010.

［9］Zhao B，Zhao W.Three－dimensional linear system analysis for breast tomosynthesis［J］.Medical Physics，2008，35 （12）：5219－5232.

［10］Deller T，Jabri KN，Sabol JM，et al.Effect of acquisition parameters on image quality in digital tomosynthesis ［C］／／Proc SPIE，2007：1－11.

［11］Erhard K，Grass M，Hitziger S，et al.Generalized filtered back－projection for digital breast tomosynthesis reconstruction［C］／／Proc SPIE，2012：1－7.

［12］Gomi T，Hirano H，Umeda T.Evaluation of the X－ray digital linear tomosynthesis reconstruction processing method for metal artifact reduction ［J］.Computerized Medical Imaging and Graphics，2009，33 （4）：267－274.

［13］WANG Dakai，HOU Yuqing.Image processing based on partial differential equations［M］.Science Press，2008 （in Chinese）. ［王大凯，侯榆青.图像处理的偏微分方程方法［M］.科学出版社，2008.］

［14］Mendrik A，Vonken E，Rutten A，et al.Noise reduction in computed tomography scans using 3－D anisotropic hybrid diffusion with continuous switch ［J］.IEEE Transactions on Medical Imaging，2009，28 （10）：1585－1594.