近年来随着影像设备不断更新,乳腺癌的探查技术得到显著提升,全数字化乳腺摄影系统(FFDM)以其良好的对比度和高分辨率成为乳腺癌早期检查的“金标准”[1] 。然而,传统的乳腺X线摄影依赖于二维的投影图像,将三维结构信息投照在二维平面图像上,使得病理结构重叠造成类似病灶的“假阳性”问题难以避免[2] ,同时,易降低病灶的可见性,造成漏诊及误诊现象。对于致密型乳腺腺体和非年轻患者乳腺腺体的诊断,其敏感性和特异性将大大降低。
数字乳腺层析(DBT)成像是一项基于平板探测器的高级应用技术,是在传统体层摄影的几何原理基础上结合数字影像处理技术开发的新型体层成像技术[3] 。与传统乳腺X线摄影得到的二维投影图像相比,DBT重建的断层图像很大程度上消除了重叠组织的干扰,使图像的“可见度”大大提高,尤其对突出组织边缘轮廓、检出微小钙化点等方面颇具优势[4] 。这使得致密型乳腺腺体中病灶的检出率得到提高,患者召回率明显降低[5] 。
因平板探测器具有高吸收率、高图像质量等优势,被广泛运用于DBT成像设备中进行投影数据采集。尽管平板探测器是一种高度精密的数字化仪器,但在性能方面仍存在些许问题[6-8] :因平板探测器中将电荷信号放大的TFT阵列存在,在没有X线照射情况下所获取的图像仍存在一定的输出暗电流现象;因探测器生产过程中制造工艺、制造批次等原因,无法满足数百个探测器一致而出现的响应曲线不一致现象;锥角范围内X线强度分布不均匀,到达平板探测器各像素所穿过空气厚度不一致造成的光场效应等。以上各种平板探测器性能问题均会导致投影图像不理想,最终在重建图像上易形成多种伪影[9] ,如环形伪影、条状伪影等。其结构特殊性也使得在输出图像的过程中易出现结构性失真[10] ,直接影响投影图像的质量。因此,为了提高重建图像的质量,获取真实可靠的投影数据,平板探测器及其所获取的投影数据的校正是图像重建与处理过程中必不可少的步骤。
本文中,我们结合DBT成像系统数据采集方式及成像几何特性,将实验原型机所使用的Analogic 公司AXS-2430探测器采集的投影数据进行暗场校正、光场及其增益校正等投影数据校正过程。而后使用基于全广义变分正则化的可获得惩罚加权最小二乘(PWLS)迭代重建算法,对校正好的投影数据进行图像重建,比较重建结果的差异,验证投影数据校正过程对DBT重建图像显像质量的影响。具体投影数据校正及重建过程流程如图 1所示。
由于平板探测器内部结构和周围环境温度的关系,探测器中会出现暗电流,暗电流的存在使得探测器在没有X射线入射的情况下仍有一定大小的信号输出,此时采集到的图像称为暗场图像,它主要是由薄膜晶体管的漏极电流、光电二极管的反向电流和集成放大电路的零点漂移引起的[11] 。由于每个像素单元情况不一,并且暗电流产生具有随机性,与X射线强度不具有相关性,致使在曝光过程中,这些暗电流仍会被探测器采集输出并叠加到投影图像中,从而形成噪声。
在X线源开启曝光的过程中,探测器获得的数据除了与物质发生作用后的X线能量衰减值外,还有一部分由暗场图像的灰度所贡献。为了真实的展示被照体内部结构与灰度值之间的对应关系,一般地将暗场图像叠加在输出响应上的灰度值减去,即暗场校正。对于暗场图像,最有效的消除方法就是采用多幅暗场图像叠加取平均值的经典方法[12] ,即:
${P^\prime }\left( {x,y} \right) = {P_0}\left( {x,y} \right) - \frac{1}{n}\sum\limits_{i = 1}^n {{D_i}} \left( {x,y} \right)$ | (1) |
其中,P0(x,y)、P'(x,y)分别为暗场校正前、后投影图像中点(x,y)的灰度值,Di(x,y)为第i个暗场图像点(x,y)的灰度值,
平板探测器增益系数主要是由X射线转换系数、光电二极管的光电转换系数和集成放大电路的放大系数决定[11] 。理想状况下,每个像素元具有相同的增益系数,并且每个像素元在每次输出中能保持各自的增益系数不变。但实际生产中,因探测单元表面的差异,导致各探测模块间的性能参数不一致,即使在同等射线强度的照射下,光电二极管的转换效应差异及射线的多能谱特性都将导致探测器获取图像灰度的不均匀[9] 。
增益校正主要是使每个探测单元对相同X射线强度的响应一致,减少图像中的固有噪声[9] 。一般利用一定的曝光剂量得到一组空气投影图像,即空扫描投影图像,对其在每个像素元上的增益系数做降权平均,便可得到整幅图像的平均增益系数,以此来反映不同像素元的相对增益系数大小。
假定各像素对射线的响应为线性,考虑到平板探测器的输出特点和每个探测器单元对X射线接收的良好线性性能,可将平板探测器校正前输出图像的表达式表示为:
${P_0}\left( {x,y} \right) = offset\left( {x,y} \right) + {k_{\left( {x,y} \right)}} \times N\left( {x,y} \right)$ | (2) |
其中,offset(x,y)表示在点(x,y)的平均暗电流,即1.1.1节中的
通过增益校正,将各像素的k(x,y)与对应的补偿因子相乘,可以获得统一的射线强度与响应输出之间的转换关系。增益校正的具体步骤如下:
(1)初始化探测器;(2)获取暗场图像。关闭或使用铅板遮挡射线源,采集100帧或更多帧图像,通过1.1.1节计算暗场平均值的方法得到一幅平均暗场图像offset(x,y);(3)采集光场图像。打开X线源或撤除铅板遮蔽,确保探测器上无任何物体。按照DBT成像设备数据采集的曝光条件,设计曝光参数采集25个投影角度下共20组光场图像,将相同角度下的光场图像取平均获得每个角度的增益场图像;(4)投影数据增益校正。分别对探测器每个角度下采集的投影数据进行增益校正,其表达式为:
$P\left( {x,y} \right) = \frac{{{P^\prime }\left( {x,y} \right)}}{{R\left( {x,y} \right) - offset\left( {x,y} \right)}} \times {\left( {R - offset} \right)_{median}}$ | (3) |
其中R(x,y)为点(x,y) 在增益场图像中的灰度值,(R - offset)median为经过暗场校正后的增益场图像灰度中值。
1.2 DBT图像重建我们拟基于投影数据几何校正的相关方法,利用数字乳腺层析成像几何特性,将角度约束下的投影状态校正到传统三维断层成像的成像模式,而后使用统计迭代重建方法,对经由探测器校正后的投影数据进行图像重建。基于投影数据统计特性,可获得PWLS目标函数[13-14] :
$\varphi \left( \mu \right)\mathop {\arg \min }\limits_{\mu \ge 0} {\left( {p - H\mu } \right)^T}{\sum ^{ - 1}}\left( {p - H\mu } \right) + \beta R\left( \mu \right)$ | (4) |
其中,μ是待估计物体的线性衰减系数,即重建图像,p是探测器获取的对数变换后的投影数据,H为系统投影矩阵,矩阵∑为对角矩阵[15] ,对应矩阵中的每个元素为测量数据的方差估计σi2 ,符号T表示矩阵转置。目标函数中R(μ)为带有图像先验信息的正则化项,β为权衡数据保真项和正则化项的超参数。
本文以课题组前期对PWLS重建模型的研究[16-18] 为基础,使用基于TGV正则化项的PWLS重建算法[19] ,运用校正后的投影数据进行图像重建。对(4)式使用一个修正的加权最小二乘替代原来的最小二乘模式即得到PWLS-TGV重建模型:
$\mathop {\min }\limits_{\mu \ge 0} {\left( {p - H\mu } \right)^T}{G^{ - 1}}\left( {p - H\mu } \right) + {\beta _2}TGV_a^2\left( \mu \right)$ | (5) |
其中,β1和β2是两个用于平衡数据保真项和正则化项的超参数项。当β1足够大时,可将修正加权最小二乘的(5)式等同于(4)式。G是系统矩阵的复合计算,其具体形式为
$\begin{array}{l} TGV_a^2\left( \mu \right) = \\ \sup \left\{ {{\smallint _\Omega }u{\rm{di}}{{\rm{v}}^2}v{\rm{d}}x|v \in C_c^2\left( {\Omega ,{S^{d \times d}}} \right),{{\left\| v \right\|}_\infty } \le {a_0},{{\left\| {{\rm{div}}v} \right\|}_\infty } \le {a_1}} \right\} \end{array}$ | (6) |
其中,Sd×d是d×d维矩阵空间,α = (α0,α1)是正数。在(5)式的基础上,引入一个新的变量f 可得:
$\begin{array}{l} \mathop {\min }\limits_f {\left( {p - Hf} \right)^T}{\sum ^{ - 1}}\left( {p - Hf} \right) + {\beta _1}\left\| {f - \mu } \right\|_2^2 = \\ {\left( {p - H\mu } \right)^T}{G^{ - 1}}\left( {p - H\mu } \right) \end{array}$ | (7) |
则可将PWLS-TGV的目标函数写为如下表达式:
$\mathop {\min }\limits_{f - \mu } {\left( {p - Hf} \right)^T}{\sum ^{ - 1}}\left( {p - Hf} \right) + {\beta _1}\left\| {f - \mu } \right\|_2^2 + {\beta _2}TGV_a^2\left( \mu \right)$ | (8) |
对于目标函数使用TGV正则化项的优势在于,其利用了高阶导数信息,在图像不满足分段常量的条件下,仍能够在有效去噪的同时保持图像的边缘和纹理信息,并且抑制图像伪影,消除了阶梯效应。对于(8)式的求解,我们将其拆分为两个子问题[20] ,分别求解最终获取的目标函数最优解:
$f = \mathop {\arg \min }\limits_f {\left( {p - Hf} \right)^T}{\sum ^{ - 1}}\left( {p - Hf} \right) + {\beta _1}\left\| {f - \mu } \right\|_2^2$ | (9) |
$\mu = \mathop {\arg \min }\limits_\mu {\beta _1}\left\| {f - \mu } \right\|_2^2 + {\beta _2}TGV_a^2\left( \mu \right)$ | (10) |
对于(9)式来说,由于其解是依照数据保真项求出,所以通常含有噪声,本文使用可分离抛物替代[21] 的方法对其进行求解,而对于(10)式,我们使用Chambolle-Pock的原始对偶算法[22] 进行求解。同时,将表达式(9)的解作为表达式(10)的初始值进行优化。
2 结果 2.1 实验数据获取实验数据主要来源于DBT成像原型机采集的物理体模投影数据,所使用的两套物理体模分别为低对比度体模及ACR标准体模。低对比度体模使用的是一块含有不同半径大小圆柱碘的衰减板,如图 2A所示,其中最右侧圆柱半径为2 cm,从右到左圆柱半径依次减小0.15 cm,并且碘浓度从上至下依次减少。该体模主要用于仿真乳腺的腺体和囊肿结构,实物如图 2B所示。实验中,我们将含碘衰减板放置于两块纯衰减板中间,为了仿真3 cm厚的乳腺照射过程,如图 2C所示。ACR标准体模是乳腺成像所使用的标准测试体模,其内部包含与乳腺结构相似衰减系数的物质,结构信息及实物如图 2D、E所示,主要用于仿真乳腺的腺体、纤维组织、钙化点、囊肿等。其中1~6为尼龙纤维,宽度分别为1.56 、1.12、0.89、0.75、0.54、0.40 mm,用于仿真乳腺组织中的纤维组织;7~11 为点状Al2O3,直径大小分别为0.54、0.40、0.32、0.24、0.16 mm,用于仿真乳腺组织中的钙化点;12~16 为空囊性结构,直径大小分别为2.00、1.00、0.75、0.50、0.25 mm,用于仿真乳腺组织中的囊肿类物质。整个体模为有机玻璃材质,腺体结构仿真区域为蜡制材料。
数据采集系统所使用的是Analogic 公司AXS-2430平板探测器[23] ,投影数据的采集参数为:管电压25 kVp,管电流100 mA,扫描时间大约为1000 ms。X线源以连续曝光的方式在±24°角度范围内每隔2°曝光一次,探测器以相同时序等间隔采集25 个投影数据。系统进行DBT数据采集过程时,设备X线源到旋转中心的距离为640 mm,球管旋转中心到探测器的距离为10 mm,原始投影数据大小为3584×2816,像素大小为0.085×0.085 mm2。实验过程中,我们将投影数据降维为1792×1408,相应像素大小为0.17×0.17 mm2。
2.2 结果分析 2.2.1 暗场校正结果图 3表示的是文中所述方法暗场校正前后低对比度体模投影数据的差异。暗场校正的目的主要是去除探测器在没有X射线入射的情况下所形成的暗电流,降低投影数据在此情况下产生的噪声及对重建图像质量的影响。从图 3的结果可以看出,经过暗场校正处理的投影图像与原始投影图像存在较大差异。我们选取两幅投影图像中五个匀质感兴趣(ROI)区域计算均值参数,红色方框标注的区域从左到右依次为区域1~5,每个区域大小为150*150个像素值。从表 1中的结果可以看出,五个ROI区域校正前与校正后投影数据的均值均存在明显差异,证明对投影数据的暗场校正取得了一定效果,暗电流影响得以抑制。
图 4给出了某两个X线源角度下探测器获取的投影数据。其中左侧第一列为探测器获取的原始投影图像,中间一列为经过单纯对数变换后的投影图像,从图像中可以很直接的看出,由于DBT设备主要以探测器静止不动X线源小范围角度这种非全等中心运动进行数据采集,使得在非中心角度的投影图像易产生“足跟效应”,即离探测器较近的区域获得的光子数大于离探测器较远的区域,物理因素造成整个平板探测器获取的光子数不均。在这种状态下获取的投影图像对比度有较大差异,同时,用于重建后重建图像空间分辨率较低,易在乳腺腺体边缘区域产生明显的环状伪影,影响重建图像质量。最右侧一列为使用文中光场及其增益校正方法校正后的投影图像,可以明显的看出,因光子数分布不均造成的足跟效应得到了很好的抑制,整幅图像亮暗分布较为均匀,图像噪声也有所降低,低对比度结构更清晰,投影图像质量得到很大提升。
图 5和图 6分别给出了低对比度体模和ACR标准乳腺体模投影数据校正前后某一断层的重建情况,左图展现的是未经过校正的投影数据用于DBT成像的图像重建结果,右图展现的是使用本文描述的投影数据校正方法校正后的投影数据用于图像重建的结果。可以明显的看出,未经校正的投影数据直接用于图像重建,会使得重建图像在体模的边缘区域有明显亮暗不均的环状伪影分布,易掩盖边缘区域的细节信息,两图箭头所指区域均存在这一现象。同时,易使得整幅图像受光场分布不均匀的影响产生匀质区域灰度不一致的现象,图像的对比度明显下降,所模拟的低对比度囊肿结构边界特征不明显,如图 6箭头所指的囊肿区域。而经过本文校正后的投影数据用于图像重建,体模边缘亮暗不均的环状伪影得到了明显的抑制,整体边缘轮廓更加明显。体模中匀质区域灰度分布较一致,低对比度囊肿结构信息明显。图像整体噪声水平较低,对比度较高,大幅度提升了图像的整体质量,达到了预期的效果。
我们对图 5和图 6中投影数据校正前后重建的两幅图像选择需增强显示的部分计算其对比度噪声比(CNR)。CNR是指图像两个不同区域信号强度的差值与本底噪声强度之间的比值,该指标综合考虑了对比度和信噪比,被认为是较为可靠的图像质量评价指标。CNR表达式如式11所示:
$CNR = \frac{{{{\bar I}_{ROI}} - {{\bar I}_{BG}}}}{{\sqrt {\sigma _{ROI}^2 + \sigma _{BG}^2} }}$ | (11) |
其中IROI和IBG分别为感兴趣区域内和背景区域内像素值的平均值,σROI和σBG分别为感兴趣区域和背景区域的标准差。对于图 5我们选择三个红色方框从左到右依次为ROI区域1~3,绿色方框则为所选择的背景区域;图 6 选择两个区域,图左边的红色和绿色方框依次为ROI区域1及其背景区域,图右边则为ROI区域2及其背景区域。分别计算所选区域的CNR值,可得表 2校正前后重建图像比较结果,可以看出,校正前后投影数据重建图像CNR值有明显差异,校正过程有效提高了图像的对比度,提升了图像的质量。
DBT技术是基于传统数字合成体层摄影技术开发的新型体层成像技术,这种特殊的数据采集模式极易因X线光子数量的分布不均、X线与探测器的相对角度位置等因素对探测器获取的投影数据造成影响[23] 。同时,投影数据采集过程中所使用的平板探测器,因其存在性能上的缺陷,极易导致投影数据质量的不理想,最终在重建图像上易形成多种伪影。现阶段,各领域已广泛展开对平板探测器成像校正的相关研究,相比于单纯的投影数据的校正工作[12, 24] ,本文侧重于设备的系统调试与数据校正,从设备研发过程的原型机开始,参考传统探测器的数据校正方法,结合DBT成像过程特有的数据采集模式,探讨研究了投影数据校正对数字乳腺层析成像质量的影响。通过对投影数据进行暗场校正、光场及其增益校正,消除投影图像存在的暗电流现象、响应曲线不一致及光子分布不均匀现象等,大幅提升投影数据的质量。同时,DBT图像重建过程中比较关键的一步为重建算法的设计,这种有限角度成像方式[25-26] ,由于投影数据的稀疏性易使重建图像质量受到噪声和伪影的严重影响。本文中,我们使用基于TGV正则化项的PWLS迭代重建算法进行DBT图像重建,验证投影数据校正过程对重建图像质量的影响。相比于基于全变分变换稀疏特性的图像重建来说[27] ,TGV先验信息的引入,一是为了在有效去噪的同时保持图像边缘和纹理信息,二则消除去噪过程中导致的阶梯效应和块状伪影。从实验结果来看,本文提出的方法能够有效降低投影图像因探测器获取光子数不同造成的亮暗差异影响,突出投影图像的细节信息,同时能够有效的抑制重建图像中的噪声和伪影,更加清晰的刻画乳腺腺体体模内低对比度组织的边缘轮廓信息,为组织解剖位置的确定奠定了良好的基础,也为后续对于探测器散射校正做出了良好的铺垫。
[1] | Antonio S, Giovanni M, Paolo R. Dedicated breast computed tomography: Basic aspects[J]. Med Phys, 2015, 42 (6): 2787-804. |
[2] | 李慧君. 数字乳腺三维断层合成摄影重建算法研究[D] . 南方医科大学, 2014. |
[3] | Park AM, Franken EA, Garg M. Breast tomosynthesis: prescntconsiderations and future applications[J]. Riadiographies, 2007, 27 : 231-40. DOI: 10.1148/rg.27si075511. |
[4] | 秦耿耿. 数字乳腺断层摄影(DBT)的剂量优化研究[D] . 南方医科大学, 2014. |
[5] | Lu Y, Chang H, Wei J, et al. A diffusion-based truncated projection artifact reduction method for iterative digital breast tomosynthesis reconstruction[J]. Phys Med Biol, 2013, 58 (3): 569-87. DOI: 10.1088/0031-9155/58/3/569. |
[6] | TAM K. Three-dimensional computerized tomography scanning methods and system for large objects with smaller area detector[J]. US Patient, 1995, 2 (14): 95-112. |
[7] | Schaller S, Noo F, Sauer F, et al. Exact radom rebinning algorithm for the long object problem in helical cone-beam CT[J]. IEEE Trans on Med Imaging, 2000, 19 (5): 361-75. DOI: 10.1109/42.870247. |
[8] | Defrise M, Noo F, Kudo H. A solution to the long object problem in helical cone-beam tomography[J]. Phys Med Biol, 2000, 45 (3): 623-43. DOI: 10.1088/0031-9155/45/3/305. |
[9] | 曾亚斌. 平板探测器数据采集及成像性能测试方法研究[D] . 南昌航空大学, 2015. |
[10] | 董歌, 罗守华, 陈功. Micro CT 投影图像噪声的去除[J]. 医疗卫生装备, 2009, 2 (30): 7-10. |
[11] | 周正干, 滕升华, 江巍, 等. X射线平板探测器数字成像及其图像校准[J]. 北京航空航天大学学报, 2004, 30 (8): 698-701. |
[12] | 王苦愚, 张定华, 黄魁东, 等. 一种锥束CT中平板探测器输出图像校正方法[J]. 计算机辅助设计与图形学学报, 2009, 21 (7): 954-61. |
[13] | Ma JH, Liang ZY, Fan Y, et al. Variance analysis of x-ray CT sinograms in the presence of electronic noise background[J]. Med Phys, 2012, 39 (7): 4051-65. DOI: 10.1118/1.4722751. |
[14] | Zhang H, Huang J, Ma JH, et al. Iterative reconstruction for x-ray computed tomography using prior-image induced nonlocal regularization[J]. IEEE Trans Biomed Eng, 2014, 61 (9): 2367-78. DOI: 10.1109/TBME.2013.2287244. |
[15] | Ma JH, Zhang H, Gao Y, et al. Iterative image reconstruction for cerebral perfusion CT using a pre-contrast scan induced edgepreserving prior[J]. Phys Med Biol, 2012, 57 (22): 7519-42. DOI: 10.1088/0031-9155/57/22/7519. |
[16] | Zeng D, Huang J, Zhang H, et al. Spectral CT image restoration via an average image-induced nonlocal means filter[J]. IEEE Trans. Biomed. Eng, 2015, 63 (5): 1044-57. |
[17] | Zeng D, Zhang X, Bian Z, et al. Cerebral perfusion computed tomography deconvolution via structure tensor total variation regularization[J]. Med Phys, 2016, 43 (5): 2091-107. DOI: 10.1118/1.4944866. |
[18] | Zeng D, Huang J, Bian Z, et al. A simple low-dose X-ray CT simulation from high-dose scan[J]. IEEE Trans. Nucl. Sci, 2015, 62 (5): 2226-33. DOI: 10.1109/TNS.2015.2467219. |
[19] | Niu SZ, Gao Y, Bian ZY, et al. Sparse-view x-ray CT reconstruction via total generalized variation regularization[J]. Phys Med Biol, 2014, 59 (12): 2997-3017. DOI: 10.1088/0031-9155/59/12/2997. |
[20] | 牛善洲. 基于变分正则化的低剂量CT成像方法研究[D] . 南方医科大学, 2015. |
[21] | Elbakri IA, Fessler JA. Statistical image reconstruction for polyenergetic x-ray computed tomography[J]. IEEE Trans Med Imaging, 2002, 21 (2): 89-99. DOI: 10.1109/42.993128. |
[22] | Chambolle A, Pock T. A first-order primal-dual algorithm for convex problems with applications to imaging[J]. J Math Imaging Vis, 2011, 40 (1): 120-45. DOI: 10.1007/s10851-010-0251-1. |
[23] | Lu YH. Improvements in lesion detection in X-ray breast tomosynthesis and in emission tomography[D] . Stony Brook University, 2015. |
[24] | 孝大宇, 蔡鹏飞, 李宏, 等. 基于阈值的平板探测器增益校正方法[J]. 东北大学学报: 自然科学版, 2014, 35 (2): 208-11. |
[25] | Chang M, Li L, Chen Z, et al. A Fiew-view reweighted sparsity hunting(FRESH) method for CT imaging reconstruction[J]. J Xray Sci Technol, 2013, 21 (2): 161-171. |
[26] | Sen Sharma K, Jin X, Holzner C, et al. Experimental studies on fiew-view reconstruction for high-resolution micro-CT[J]. J Xray Sci Technol, 2013, 21 (2): 25-42. |
[27] | Sidky EY and Pan X, Holzner C. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[J]. Phys Med Biol, 2008, 53 (17): 4777-807. DOI: 10.1088/0031-9155/53/17/021. |