乳腺癌是女性最常见的恶性肿瘤。据世界卫生组织国际癌症中心统计, 2012年全球女性乳腺癌新发病例达167万, 同期约52万女性因乳腺癌死亡, 占所有女性恶性肿瘤死亡的14.7%, 乳腺癌的发病率及死亡率均位居女性恶性肿瘤中的首位^[1]。据统计, 我国每年乳腺癌新发病例为16.9万, 是女性第2位常见恶性肿瘤, 死亡人数约4.5万人/年, 在恶性肿瘤致死率中排第6位^[2]。研究表明, 乳腺癌的早期发现和治疗, 能有效预防控制肿瘤恶化, 提高患者5年生存期。

全数字化X射线乳腺摄影(FFDM)作为当前临床最常用的乳腺癌筛查技术^[3], 具有辐射剂量低, 筛查成本低, 病人吞吐量大等优点。但是, FFDM为平片成像, 具有乳腺组织的堆叠效应, 致使其病灶检出率降低, 易造成误诊、漏诊等临床事故。数字乳腺层析成像(DBT)技术以传统的数字乳腺摄影为基础, 在有限角度内, 对乳腺组织进行数次曝光, 从而对乳腺进行断层成像。相比于FFDM, DBT成像消除了乳腺组织重叠的干扰, 进而提高了乳腺病灶的检出率^[4-5]。

针对二维平面成像的缺陷, 国外对三维DBT系统进行大量研究的基础上, 推出可实现乳腺3D成像的DBT系统^[6]。该系统不仅可以产生二维数字化乳腺图像, 还可以重建乳腺的断层图像, 能够更加清晰的看到乳腺组织和病灶结构。但目前国内对于DBT成像系统的研究多局限于某一项技术, 比如图像重建、剂量控制、图像质量评估, 但缺乏对整个DBT成像系统流程的研究, 没有考虑到各个环节的处理对最终图像质量的影响; 同时也缺乏对几何校正、滤波调制、图像增强等关键技术的研究。

针对上述问题, 本研究团队搭建了DBT硬件成像平台, 并深入研究DBT成像系统流程的设计以及成像关键技术研究, 并结合三维断层成像, 实现乳腺层析成像。该DBT系统主要包括X射线曝光控制模块、机械控制模块、数据采集模块、数据校正与图像重建模块。其中, 投影增强、几何校正、滤波调制及图像重建是实现优质成像的关键技术。实验中, 使用该原型机系统扫描多种体模(ACR体模和高分辨率体模)并重建, 评估系统的稳定性及成像性能。

1 系统与技术 1.1 DBT整机系统简介

DBT是以传统体层摄影几何原理为基础, 结合数字影像处理技术开发的新型体层成像技术, 运用与传统数字乳腺成像设备类似的压迫板, 将乳腺组织挤压到探测器上方的碳纤维板上, 减少组织结构重叠^[7]。在数据采集过程中, 探测器保持静止或者相对运动, X线源围绕旋转中心以某一个限定角度做弧形运动, 采集有限角度内的投影数据, 重建与探测器平面平行的任意断层层面的图像。考虑系统机械稳定性, 本文选择探测器静止策略, 图 1 (A、B)分别显示了DBT成像几何示意图及原型机系统。该硬件系统主要由X射线曝光控制模块、机械控制模块、数据采集模块等组成。

1.1.1 X射线曝光控制模块

该模块主要包括X射线源和高压发生器。X射源采用旋转的钨阳极射线球管, 其热容量为150 KJ, 阳极倾角为15°。射源焦点为双焦点, 分别为0.1/0.3 mm。该系统采用乳腺成像专用的Spellman VMX高压发生器, 其功率为5 KW, 标称电压为20~40 KV, 最大电流为200 mA。DBT系统的曝光控制模块控制高压发生器在短时间内进行连续高压切换, 控制球管曝光, 产生X射线。数据采集时, 球管的旋转角度为±20°, 采集21帧投影, 采集帧率为2 fps, 射源锥角范围为25°。

1.1.2 机械控制模块

机械控制模块包括机架升降控制、C臂旋转控制、压迫器运动控制等。该系统机架呈站立式, 其主要作用为保持扫描系统的稳定性, 以获得特定扫描方式下的乳腺投影数据。在采集数据的过程中, 探测器静止不动, C臂和高压球管绕着探测器进行有限角度的曝光扫描。系统采用自动/手动的方式控制压迫器按压乳腺, 最大压力为200 N, 压迫板运动范围5~245 mm, 具有快速解压装置。压迫器按压乳腺能减小乳腺厚度, 有利于X射线穿透, 降低几何模糊, 提高空间分辨率。

1.1.3 数据采集模块

系统数据采集采用第二代非晶硒数字乳腺平板探测器(型号:Analogic Canda AXS- 2430), 将接收的X射线光子直接转换成电荷, 再由晶体管阵列将电信号读出并数字化成投影数据, 具有转换效率高, 动态范围广等优点。探测器实际尺寸为239.36 cm×304.65 cm, 维度为2816×3584, 像素尺寸为0.085 mm×0.085 mm, 空间分辨率为5.8l p/mm。控制台将数据传输到工作站, 进行数据处理及图像重建。

1.2 数据校正与图像重建模块

数据校正与图像重建模块具体包括投影数据校正、投影增强、系统几何校正、滤波调制以及图像重建, 该模块为DBT成像系统的关键技术部分, 其技术处理流程如图 2所示。

1.2.1 投影数据校正

因探测器生产工艺、制造批次等问题, 造成探测器单元响应不一致以及无曝光情况下出现暗电流信号^[8]; 因锥角范围内X射线强度的不均性造成的光场效应等因素^[9-10], 均会导致采集的投影数据质量下降, 因此需要对投影数据进行暗场校正及增益校正^{[5, 11]}。

对于暗场校正最有效的消除方法就是用投影数据与暗场数据均值进行相减^[8], 公式如下:

$ P'\left( {x, y} \right) = {P_0}\left( {x, y} \right) - {P_{dark}} $

(1)

其中, (x, y)表示投影数据的索引, P₀、P'分别为暗场校正前后投影数据, ${P_{dark}} = \frac{1}{n}\sum\limits_{i = 1}^n {{D_i}\left( {x, y} \right)} $为暗场数据的均值, D_i为第i个暗场数据, n为暗场数据的个数。

增益校正主要是使每个探测单元对相同X射线强度的响应一致, 减少图像中的固有噪声^[7]。其做法是在相同的曝光条件下, 采集多组空扫描光场数据, 并将相同角度下的光场数据取平均, 获得对应的增益场数据。利用上述的暗场数据, 求取增益系数, 进行投影的增益校正, 其公式如下:

$ P\left( {x, y} \right) = \frac{{P'\left( {x, y} \right)}}{{G\left( {x, y} \right) - {P_{dark}}}} \times {\left( {G\left( {x, y} \right) - {P_{dark}}} \right)_{median}} $

(2)

其中, $G\left( {x, y} \right){\rm{ = }}\frac{1}{k}\sum\limits_{j = 1}^k {{R_j}\left( {x, y} \right)} $为增益场数据, R_j为空扫数据, k为空扫次数, P'为暗场校正后的投影数据, P为增益校正后的数据, (.)median用于求取括号内数据中值。上述数据校正为传统的探测器校正技术, 已集成到DBT数据采集的过程中, 校正后的数据用于图像重建。

1.2.2 投影增强

由于乳腺组织之间衰减系数接近, 导致DBT采集的投影数据对比度较低, 需要在重建前对原始的投影数据进行增强。本文基于限制对比度直方图均衡算法(CLAHE)进行投影增强, 通过限制局部直方图的高度来控制局部的对比度, 从而防止噪声的放大以及局部的对比度过度增强。但传统的CLAHE算法会导致增强后的图像有块状伪影^[12-16], 这是由于该算法对局部的图像块进行不同程度的增强造成的。

本文改进了CLAHE算法的取块方式, 让图像块之间相互重叠, 而不是之前的无重叠, 这样增强后的图像块之间的边界能够很好的融合, 从而消除了块状伪影。为了抑制噪声获得更好的增强效果, 本算法还结合了图像灰度拉伸以及引导滤波去噪。具体步骤如下:(1)利用最大类间方差, 分割出投影中的空气和乳腺区域; (2)统计乳腺区域的灰度直方图; (3)利用致密腺体组织在直方图中的分布, 分割致密腺体组织; (4)对整个乳腺区域进行线性灰度拉伸, 再对分割出的致密腺体组织进行拉伸, 生成归一化的系数; (5)采用改进的CLAHE算法增强拉伸后的投影数据; (6)将增强后的投影数据与归一化的系数相乘; (7)采用引导滤波对增强后的投影进行保边去噪, 得到滤波增强后的投影数据。由于DBT的小角度成像, 无法重建准确的CT值, 所以才能在步骤(4)中对投影数据进行灰度拉伸, 算法处理流程图(图 3)。

1.2.3 系统几何校正

在硬件加工和搭建过程中, 因机械抖动、重力等因素导致实际成像的几何与设计的几何之间存在偏差。本文采用基于钢珠点体模的方法进行几何校正, 利用钢珠点在世界坐标系中球心坐标和投影坐标系中圆心坐标的映射关系, 求解系统几何参数^[17-21]。射源S, 体模中钢珠点的空间坐标(x, y, z)及投影中钢珠的圆心坐标(u, v)的关系如图 4 (A)所示, 其中(u_d, v_d)表示探测器的中心坐标。几何校正体模实物图如图 4 (B)所示, 该体模包含有84个钢珠点, 分布在上下两层塑料板中, 每层42个钢珠点。体模上钢珠点球心坐标和探测器上圆心坐标之间的映射关系可用一个矩阵形式表示:

$ [uw, vw, w] = P{[x, y, z, 1]^{\rm{T}}} $

(3)

P是一个3 × 4的投影矩阵; w是距离权重因子。该矩阵包含系统所有的几何信息, 可以被分解为:

$ P = K[R|t] $

(4)

R是一个3 × 3的旋转矩阵^[20], t是一个3 × 1的变换向量。K是一个3 × 3系统内部矩阵, 包含了成像系统的几何信息。

$ K = \left( {\begin{array}{*{20}{c}} {\frac{f}{{{p_u}}}}&{\frac{f}{{{p_u}\tan \alpha }}}&{{u_0}}\\ 0&{\frac{f}{{{p_v}\sin \alpha }}}&{{v_0}}\\ 0&0&1 \end{array}} \right) $

(5)

f是射源到探测器的距离, p_u和p_v表示探测器单元的宽和高, u₀和v₀是射源中心射束与探测器交点的坐标。对于平板探测器, α等于90 °。公式(4)中的R和t包含系统的外部参数, 表示物体上的世界坐标系到平板探测器的坐标变换。R可以进一步用三个欧拉角表示:

$ \begin{array}{l} R = \left( {\begin{array}{*{20}{c}} {\cos {\theta _z}\cos {\theta _y}}&{ - \sin {\theta _z}\cos {\theta _x} + \cos {\theta _z}\sin {\theta _y}\sin {\theta _x}}\\ {\sin {\theta _z}\cos {\theta _y}}&{\cos {\theta _z}\cos {\theta _x} + \sin {\theta _z}\sin {\theta _y}\sin {\theta _x}}\\ { - \sin {\theta _y}}&{\cos {\theta _y}\sin {\theta _x}} \end{array}} \right)\\ \;\;\;\;\;\;\left( \begin{array}{l} \sin {\theta _z}\sin {\theta _x} + \cos {\theta _z}\sin {\theta _y}\cos {\theta _x}\\ - \cos {\theta _z}\sin {\theta _x} + \sin {\theta _z}\sin {\theta _y}\cos {\theta _x}\\ \cos {\theta _y}\cos {\theta _x} \end{array} \right) \end{array} $

(6)

其中, θ_x, θ_y, θ_z表示探测器绕着x轴, y轴, z轴的偏转角度。通过体模钢珠点坐标(x, y, z)与投影坐标(u, v)之间的映射关系, 可解析出投影矩阵P, 进而可求解出成像系统的几何参数, 用于图像的精确重建^[17]。

1.2.4 FDK重建及滤波核函数调制

本文采用三维的滤波反投影(FDK)算法重建乳腺图像^[22-27]。由于乳腺组织之间的衰减系数接近, 使用传统的"Ramp"滤波核所重建的图像对比度低, 噪声强, 不适合乳腺组织成像。因此, 对滤波核进行调制^[28], 在傅里叶域中同时对两个方向进行滤波, 即分别沿着探测器行方向(球管运动方向)和图像层厚方向, 以增强图像的软组织对比度及抑制图像噪声。滤波核调制前后的波形如图 5所示。

调制的滤波核函数H_filter为:

$ \begin{array}{l} {H_{filter}}\left( {{\omega _y}, {\omega _z}} \right) = {H_{sprec\;{\rm{tr}}\;um}}\left( {\omega y} \right) \cdot {H_{profile}}\left( {{\omega _y}, {\omega _z}} \right)\\ \cdot {H_{in\;{\rm{vers}}\;e}}\left( {{\omega _y}, {\omega _z}} \right) \end{array} $

(7)

其中, H_{in vers e}是调制传输函数的反转, 并且与Ramp滤波核成正比。

$ {H_{in\;{\rm{vers}}\;e}}\left( {{\omega _y}, {\omega _z}} \right) = 2\alpha \sqrt {\omega _y^2 + \omega _z^2} $

(8)

α表示比例系数, ω_y、ω_z分别表示y, z方向的频率分量。H_{spec tr um}滤波核和H_profile滤波核可以表示为:

$ {H_{spec\;{\rm{tr}}\;um}}\left( {{\omega _y}} \right) = \left\{ \begin{array}{l} 0.5 \cdot \left( {1 + \cos \frac{{\pi \omega }}{{{D_y}}}} \right)\\ \;\;\;for\left| {{\omega _y}} \right| < {D_y}\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;elsewhere \end{array} \right. $

(9)

$ {H_{profile}}\left( {{\omega _z}} \right) = \left\{ \begin{array}{l} 0.5 \cdot \left( {1 + \cos \frac{{\pi \omega }}{{{D_z}}}} \right)\;for\left| {{\omega _z}} \right|\\ < {D_z}\;\;and\;\left| {{\omega _z}} \right| < \tan \;\alpha \left| {{\omega _y}} \right|\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;elsewhere \end{array} \right. $

(10)

D_y和D_z分别表示在ω_y、ω_z方向的截止频率, ω是三维傅里叶空间的频率。

1.3 系统评估 1.3.1 实验数据

实验采用体模数据和病人数据, 其中体模数据包括高分辨率体模和ACR标准体模。高分辨率体模用于检验DBT成像系统的几何精度以及图像的分辨率, 如图 6 (A、B)所示, 其中不同的数字对应不同的分辨率, 数字越大, 分辨率越高。测试系统分辨率时, 体模倾斜45°摆放可以避免体模中的线对与X射线形成缝隙干涉, 便于读数, 且能综合性评估系统的分辨率。ACR体模是DBT的标准测试体模^[28], 其内部包含与乳腺衰减系数相似的物质, 主要用于仿真乳腺的囊肿、腺体、钙化点、纤维组织等。其中1~6为尼龙纤维, 用于仿真乳腺中的纤维组织; 7~11为点簇状的三氧化二铝, 用于仿真乳腺中的钙化点; 12~16为空囊性结构, 用于仿真乳腺中的囊肿组织, 体模的材质为有机玻璃, 腺体结构仿真区域为蜡制材料, 体模的实物图及结构信息图 6 (C、D)所示。此外还使用商用DBT设备对8位患者进行扫描, 获取患者投影数据, 并采用文中的数据处理方法进行图像重建, 将之与商用设备重建图像做比较, 以说明本文方法的有效性。

1.3.2 图像质量评估

本文采用对比度噪声比[5(] CNR)对图像进行定量分析。其中, CNR表示图像两个不同区域信号强度的差值与背景噪声强度之间的比值, 该指标综合考虑了对比度和信噪比。

此外, 采用专家评分的方式对病人数据进行定性评估。由三位资深的影像科医师进行回顾性双盲法阅片, 参照美国放射学会的乳腺影像报告和数据系统, 以每侧乳腺为单位, 对所有乳腺数据DBT图像进行分析, 分别记录阅片结果。观察内容包括:病变大小、形态、边缘、毛刺、结构扭曲、血管穿入以及其他特征。综合评分范围为0~5分。

2 实验分析 2.1 系统几何校正对重建结果的影响

图 7中(A、B)分别为使用几何校正前及几何校正后的参数所重建高分辨率体模图像。从相互垂直的两组线对的局部放大图可以看到, 使用未校正的几何参数重建的图像中包含严重的伪影, 线对和数字模糊; 使用校正后的几何参数进行重建, 体模中两组线对及数字可清晰显示, 最高线对分辨率可达到7l p/mm。

为评估DBT系统的机械稳定性, 在相同扫描剂量下(28 kV, 80 mA), 连续采集7组几何校正体模的数据。利用7组几何校正数据, 求出对应的几何参数, 并求取几何参数的均值和方差。其中, DSD表示射源到探测器的距离, DSO表示射源到旋转中心的距离, airgap表示碳纤维板到探测器表面的距离, angle表示射源的偏转角度如图 1 (A), ΔU₀, ΔV₀分别表示射源中心射束和探测器的交点与探测器中心的偏移量(ΔU₀=u_d-u₀; ΔV₀=v_d-v₀); 如图 4(A), 除angle之外, 其他参数的单位均为mm, 结果如表 1所示, 系统的几何参数波动小, 扫描系统运行稳定。由于DBT系统采集数据时, 探测器静止不动, 导致每个偏转角度下的DSD都是随角度的变化而变化的, 所以表中没有给出DSD的均值和方差。

2.2 滤波核函数调制对重建结果的影响

图 8所示为滤波调制前后ACR体模的重建结果对比。其中, (A)为传统的Ramp滤波核重建的结果, (B)为本文方法重建结果。从图中可以看出, 滤波调制后重建图像的对比度更高, 囊肿和纤维束的显示更清晰, 噪声水平更低(如图中红色箭头指示区域)。从两幅重建的ACR体模图像中, 选择3组ROI区域(红色区域为囊肿区域; 黄色区域为背景区域), 从上往下依次为ROI1- ROI3, 分别计算所选区域的CNR值, 如表 2所示。可以看出, 经滤波调制后, 图像的CNR得到显著提升。

2.3 投影增强对重建结果的影响

图 9为乳腺投影数据增强前后的对比图。从图中可以看出, 增强后的投影中乳腺血管、纤维组织之间的对比度更高。为了定量的评估图像的对比度, 我们在乳腺投影中选取5个ROI, 从上到下依次为ROI1-ROI5, 计算ROI的对比度。对比度的计算, 需要测量图像中的最亮点和最暗点之间的灰度级数。根据对比度的计算公式为: $Con = \sum\limits_\delta {\delta {{(i, j)}^2}{P_\delta }(i, j)} $, 其中, i, j表示像素索引, δ(i, j)=|i-j|表示相邻像素间的灰度差, P_δ(i, j)表示相邻像素间的灰度差为δ的像素分布概率^[30]。计算结果如表 3所示, 增强后投影的ROI对比度更高。

2.4 病人数据成像

为了评估本文DBT系统成像影像链的处理效果, 采用临床病人的数据进行重建, 并与商用DBT设备的结果进行对比。影像专家将根据观察乳腺的病变大小、形态、边缘、毛刺、结构扭曲、血管穿入以及其他特征综合评分, 满分为5分。实验重建了8组病人数据, 经统计, 本文DBT成像系统重建结果平均分为4.52分, 最高分为4.7分; 商用DBT设备重建图像平均分为4.57分, 最高分为4.9分。从评分结果可以看出, 本文DBT系统成像影像链重建的图像质量接近目前商用DBT设备重建效果。

本文选取1组病人重建图像进行显示, 如图 10所示, 其中, 左图为商用DBT设备重建结果, 右图为本文DBT成像系统重建结果。为了比较乳腺中的病灶区域的显示, 我们选取图像中的肿瘤和微钙化点簇进行放大比较, 如局部放大图所示, 肿瘤边界显示清楚, 与周围组织对比明显; 钙化点显示清晰, 数目准确。重建图像病灶显示清楚, 乳腺细节明显, 成像效果与商用DBT系统相当。

3 讨论

对于乳腺癌的筛查, DBT成像相比于传统的FFDM能够消除组织重叠, 提高病灶的检出率, 获取乳腺的断层图像。但因DBT系统采用低剂量的限角扫描, 投影数据少, 且其对于扫描系统的精度和稳定性要求更高; 若要获取高质量的乳腺断层图像, 需要根据DBT成像系统的特性以及乳腺数据特点, 针对性的进行投影数据校正、几何校正、投影增强、滤波调制及图像重建。鉴于此, 作者所在实验室研发了DBT整机系统, 实现其成像影像链的设计与优化。

本文简单介绍了DBT系统的成像模式以及硬件控制系统, 着重描述了DBT成像的相关关键技术, 主要包括X射线曝光控制模块、机械控制模块、数据采集模块、数据校正模块、几何校正模块、投影增强模块、滤波调制及图像重建模块等7大模块。在实验中, 通过对探测器的暗场效应以及响应不一致进行数据校正, 并对乳腺组织间的对比度不足进行投影增强, 增加其对比度; 此外, 利用钢珠点校正体模对系统进行几何校正, 测试系统的稳定性, 并结合FDK重建算法, 调制适合DBT采集模式的滤波器, 抑制图像中的噪声和几何伪影, 保持图像的细节信息, 获取高品质的乳腺断层图像。

虽本文部分工作引用现有算法, 但并不是简单的算法组合, 而是结合整机扫描系统及数据特点, 多次试验验证的结果。如在设计投影增强时, 把一般用于图像后处理的增强技术, 用于乳腺投影的数据增强, 这样改进后, 重建的图像对比度更好。此外, 本文采用的是改进后的CLAHE增强算法, 消除了原增强算法引入的块状伪影, 增强效果更好。所以本文工作是具有一定的创新性和新颖性的, 不仅有技术点的创新, 也有系统流程设计的创新。

除却上述优点, 本文所介绍的DBT成像系统还存在一些不足。比如, 乳腺增强中, 小波变换^[31]可通过构造滤波器刻意提升或者压低变换域的某些系数, 再通过反变换回到图像域从而达到乳腺全局增强的目的, 但对局部细节增强效果不好。而相比于小波变换, 本文增强方法是利用乳腺空间域的局部信息进行增强, 未充分使用乳腺的全局信息, 可能会出现局部增强过度的问题。此外, 本文还未对DBT成像中的层间伪影^[28]和金属伪影^[32]进行校正, 影像链的设计还不完善。在未来工作中, 将利用乳腺的局部信息与全局信息相结合的方式, 进行乳腺增强; 还将针对DBT成像中的层间伪影和射束硬化伪影进行校正, 进一步提升图像质量。

综上所述, 本文所介绍的DBT系统能够对乳腺进行优质的断层成像, 增加组织间的对比度, 提高病灶的检出率, 为乳腺疾病的早期诊断和筛查提供了重要的影像信息。该系统的成功研制为国内DBT系统的研究提供了技术支持, 推动了国产DBT设备的研发进程。