0 引 言

1 方法

1.1 胸片图像数据库

1.2 基于解剖结构的人工神经网络

为了能够抑制骨骼，MTANN^[17]采用了线性输出多尺度ANN的回归模型，它可以直接对像素数据进行操作，输出处理后的像素值^[18]。这个模型采用线性函数代替S型函数作为输出层的激励函数。使用线性函数是因为当人工神经网络将线性函数应用于连续映射的值的图像处理中，它的性能会得到显著的提高^[19-20]。其他的机器学习回归模型也能被用于MTANN框架(也叫做基于像素的机器学习)中，如支持向量回归机和非线性高斯处理回归模型^[20]，输出为连续值。

为了分割出骨骼，将标准的胸片图像按像素分成大量重叠的子区域。从训练图像中提取的单个像素相当于输入子区域，作为训练值。通过大量子区域和对应骨质图像的像素值训练线性输出的ANN。一旦ANN训练好后，模型将不再需要骨质图像。用训练好后的ANN处理其他标准胸片，输出虚拟的骨质图像。MTANN的输出是连续的灰度值，与子区域中心点的像素有关，其表达式为

$o\left( {x,y} \right) = ML\left\{ {I\left( {x - i,y - j} \right)|\left( {i,j} \right) \in {\bf{R}}} \right\}$

(1)

式中，ML(·)是回归模型的输出，其输出值作为子区域中心点新的像素值。I(x,y)是输入图像的像素值。训练ANN时，最小化的误差可表示为

$E = \frac{1}{P}\sum\limits_c {\sum\limits_{\left( {x,y} \right) \in {R_T}} {{{\left\{ {{T_c}\left( {x,y} \right) - {O_c}\left( {x,y} \right)} \right\}}^2}} } $

(2)

式中，c是训练的样本数，O_c是第c个样本的输出值,T_c是第c个样本的学习值，P是训练区域中像素的总数目。

胸片中，不同位置的骨头方向、宽度、对比度等也不同，由于单个MTANN性能有限，为了提高MTANN抑制不同位置的骨骼的性能,对单个MTANN进行了扩展，用8组MTANN对特定解剖结构的骨骼进行抑制，如图 1(a)所示。对于每个特定解剖结构，采用正常的图像样本与在对应解剖结构区域含有节结的图像样本对MTANN进行训练。本文将肺区分为8个部分：左侧肋骨和锁骨的左侧上部区域、肺门的左肺门区域、肺区中间骨头的左侧中间区域、肺下叶区域肋骨的左侧下部区域、右侧肋骨和锁骨的右侧上部区域、肺门区域骨头的右肺门区域、肺区中间骨头的右侧中间区域以及肺下叶区域肋骨的右侧下部区域。对于每个特定解剖结构的MTANN，用其对应的分割图像提取训练样本。

图 1 特定解剖结构的多个大规模训练人工神经网络模型的结构

Fig. 1 Architecture and training of new anatomically specific MTANNs ((a)training phase；(b)execution phase)

训练完成后，用训练好后的每个特定解剖结构的MTANN处理其他胸片图像，输出为各个解剖结构对应区域的骨质图像。如图 2所示，MTANN 1输出左侧上部有锁骨的肺部区域，MTANN 2输出左肺门区域等。

图 2 特定解剖结构的MTANN训练后的输出图像

Fig. 2 Output bone images of the trained anatomically specific anatomical MTANNs ((a) the hilar region；(b) the lower region； (c) the middle region； (d) the upper region)

然后将8个区域的骨质图像按照8个解剖分区的模板合成一幅完整的虚拟双能量骨头图像。为了使得混合的子图像边缘更加平滑，采用高斯滤波器对分割模板进行平滑，即

${f_b}\left( {x,y} \right) = \sum\limits_{i = 1}^8 {{O_i}\left( {x,y} \right)} * {f_G}\left[ {{M_i}\left( {x,y} \right)} \right]$

(3)

式中，f_b(x,y)是合成的骨头图像，O_i是第i个ANN，f_G是高斯滤波器，M_i是第i个解剖结构模板。

1.3 训练解剖结构分区的人工神经网络

1.4 解剖结构分区的自动分割

为了将训练好后的ANNs自动地用于图像的骨质结构抑制，需要自动分割胸片肺区中各解剖结构。因为每个解剖结构的ANN所对应的分区在所有的分区中信噪比是最高的，如图 2所示。自动合并所有的分区将产生一幅完整的虚拟骨质图像，它的信噪比要高于所有的单个分区。

为了确定8个解剖结构的分区，提出了一种基于活动形状模型^[21](ASMs)的自动解剖结构的分割方法。采用一种适用于每个分区的多分区ASMs^[22-23](M-ASM)将肺区自动分割，如图 3所示。由于传统的ASM方法中节点沿着整个肺的边缘均匀分布，这样有些曲率很大的部分就无法适应。本文方法通过修正在特殊结构边界的节点(即移动地标)来对模型进行了改进。移动地标可以识别从一种边界类型(如肺野与心脏的边界)到另一种边界类型(如心脏与膈的边界)的改变。这样就可以对肺区边界进行了多处分段并且每个分区对应了一种特有的边界类型(如心、主动脉、胸腔、隔膜等)。特定节点的ASM对于每个分区边界采用固定数量的节点并且等距分布。对整个肺边界采用了50个节点的非等距分布的M-ASM，每个分区边界都分配固定数量的节点，并沿着边界等距分布。图 3中，无论实际肺区边界如何分布，每幅图像中左肺野与心脏之间的边界包含了11个点。这使得节点可以适应特定分区的边界而不是整个肺区，从而显著提高分割的准确程度。在训练图像中，每段分区中节点在空间上的相对关系将会被学习，以形成了形状模型。这些节点可形成向量x并投射到主分量形状空间，即

图 3 M-ASM自动分割的结果

Fig. 3 Results of M-ASM automatic segmentation

$\boldsymbol{b} = {\boldsymbol{V}^{\rm{T}}}\left( {\boldsymbol{x} - \boldsymbol{\bar x}} \right)$

(4)

式中，V=[V₁，V₂，V₃，…,V_m]是前m个特征向量的协方差矩阵，b=[b₁，b₂，b₃，…,b_m]^T是主要数轴形状系数的矢量。形状系数取值范围为$ \pm m\sqrt {{\lambda _i}} $，且只产生唯一的合理的形状并投射回节点坐标，表示为

$\boldsymbol{x} = \boldsymbol{\bar x} + \boldsymbol{Vb}$

(5)

m通常取值为2～3^[24]，实验中将其设置为2.5。

肺部分割后，通过边界类型和标记将它们自动分割为8个解剖学分区。文中使用标记点，可获得上部区域，下部区域和肺门区域，如图 3所示。由多个ANN输出的8幅分区的骨质图像组成1幅完整的VDE骨头图像，即

${f_{\rm{b}}}\left( {x,y} \right) = \sum\limits_{i = 1}^8 {f_{\rm{b}}^i\left( {x,y} \right) * m_{\rm{b}}^i\left( {x,y} \right)} $

(6)

式中，f_bⁱ(x,y)是第i个ANN的输出图像，m_bⁱ(x,y)是第i个ANN的分割模板。分割模板经过高斯滤波器平滑后不同分区间的不连续性在合并后就消失了。

1.5 虚拟软组织图像的生成

在获取了VDE骨图像之后(如图 4)，骨骼图像含有大量噪声信息，在将其与原图相减获得软组织图像之前，需要对其进行平滑去噪处理。引入能够保留骨骼图像边缘信息的总变分最小化方法平滑骨骼图像，最终生成虚拟软组织图像。

图 4 从骨骼图像中获取软组织图像的方法

Fig. 4 Method for obtaining a soft-tissue image from a bone image

通过与原始图像进行去除法处理获得软组织图像。在分割骨骼后，采用高斯滤波用于平滑分割区域的边缘，产生的图像为掩膜图像m(x,y)，掩膜图像的标准值域为[0, 1]。

为了抑制标准胸片中的肋骨等结构，需要在原图的肺区m(x,y)中减去用多尺度MTANN 产生的虚拟骨质图像，即

$\begin{array}{l} {f_s}\left( {x,y} \right) = g\left( {x,y} \right) - {w_c} \times \\ \;\;\;\;\;\;\;{f_b}\left( {x,y} \right) \times m\left( {x,y} \right) \end{array}$

(7)

式中，w_c是决定肋骨对比度的权重变量。通过改变w_c，可以获得不同肋骨对比度的虚拟软组织胸片。

VDE骨图像中存在噪声，使用高斯滤波的方法虽然可以滤掉VDE骨头图像中的噪声，但也可能模糊骨头的边缘。提出一种总变分(TV)最小化噪声去除方法，TV最小化问题第一次是被Rudin用来进行图像的滤波^[25]。该方法可以在滤除VDE骨头图像的噪声的同时保留骨头的边缘信息(如图 4)，这是因为TV准则的分段滤波标准属性。假设VDE骨头图像的噪声为高斯白噪声，即

$z\left( {x,y} \right) = u\left( {x,y} \right) + \eta \left( {x,y} \right)$

(8)

式中，u(x,y)是未知的2维分段连续函数，它表示原始无噪声图像，z(x,y)是u(x,y)的噪声图像，η(x,y)是高斯白噪声。传统的噪声抑制使用的是Wiener滤波，它是将能量函数最小化，即

$T\left( u \right) = \frac{1}{2}{\left\| {u - z} \right\|^2} + aJ\left( u \right)$

(9)

J的常见值为

$J\left( u \right) = \int {{u^2}{\rm{d}}x} $

(10)

但当u不连续时，式(9)常会使图像产生模糊以及假震荡。

因此，考虑非线性的TV能量函数

${J_{{\rm{TV}}}}\left( u \right) = \int_\Omega {\left| {\nabla u} \right|} {\rm{d}}x$

(11)

式中，▽u表示u的梯度，$\nabla u = \left( {\frac{{\partial u}}{{\partial x}},\frac{{\partial u}}{{\partial y}}} \right)$，并且不要求u连续。

但是，Euclidean 准则在0处是不可导的。为了避免不可导引起的麻烦，将式(11)修改为

${J_\beta }\left( u \right) = \int_\Omega {\sqrt {\left| {\nabla u} \right| + {\beta ^2}} {\rm{d}}x} $

最后简化的能量函数为

$\begin{array}{l} T\left( u \right) = \frac{1}{2}{\left\| {u - z} \right\|^2} + \\ \;\;\;a\int_\Omega {\sqrt {\left| {\nabla u} \right| + {\beta ^2}} {\rm{d}}x} \end{array}$

(12)

将欧拉—拉格朗日等式与式(12)合并可得到一个椭圆型偏微分方程等式，即

$\begin{array}{l} u + \alpha L\left( u \right)u = z,x \in \Omega \\ \;\;\;\;\;\frac{{\partial u}}{{\partial n}} = 0,x \in \partial \Omega \end{array}$

(13)

式中，L(u)是微分算子，u的计算公式为

$L\left( u \right)u = - \nabla \left( {\frac{1}{{\sqrt {\left| {\nabla u} \right| + {\beta ^2}} }}\nabla u} \right)$

(14)

通过式(14)可以发现，当梯度强度增加时，滤波能力就会降低，到达边缘时停止滤波。

为了解决TV最小化的问题，有很多标准的数字优化技术，如共轭梯度法，但这些标准方法不能很好的解决骨头边缘被模糊的问题。采用非线性多网格方法来处理这个问题。与传统方法不同，多网格算法可以解决含可变系数的非线性椭圆型偏微分方程问题，并且只需求解粗网格上的方程，而不会降低处理的效率。

假设将式(13)离散到一个均匀的网格大小为h的网格上，即

${T_h}\left( {{u_h}} \right) = {z_h} = {u_h} + \alpha {L_h}\left( {{u_h}} \right){u_h}$

(15)

以${\tilde u_h}$表示某些近似解和式(15)的确定解。更正为${v_h} = {u_h} - {{\tilde u}_h}$，则

$\begin{array}{l} {T_h}\left( {{{\tilde u}_h} + {v_h}} \right) - {T_h}\left( {{{\tilde u}_h}} \right) = \\ \;\;\;{f_h} - {T_h}\left( {{{\tilde u}_h}} \right) = - {d_h} \end{array}$

(16)

现在构造一个网格大小为H(通常取H=2h)的粗网格T_h近似值T_H。剩余等式(16)近似为

${T_H}\left( {{u_H}} \right) - {T_H}\left( {{{\tilde u}_H}} \right) = - {d_H}$

(17)

因为T_H有更小的维数，则该等式将更容易求解。为了定义粗网格上的${{{\tilde u}_H}}$和d_H，需要给粗网格一个约束参数H来约束${{{\tilde u}_H}}$和d_h，也就是在粗网格上求解

${T_H}\left( {{u_H}} \right) = {T_H}\left( {R{{\tilde u}_h}} \right) - R{d_h}$

(18)

然后将粗网格更改为${\tilde v_H} = {u_H} - R{\tilde u_h}$，有了粗网格${\tilde v_H}$的解之后，需要一个延伸参数P将更正值插入到细网格中，得到${\tilde v_h} = P{\tilde v_H}$。因此得

$\tilde u_h^{{\rm{new}}} = {\tilde u_h} + P{\tilde v_H}$

(19)

式中，$\tilde u_h^{{\rm{new}}}$是当前多网格尺度上，修正后的解。

2 实验结果

2.1 肺区分割

肺区的分割对于骨骼的抑制特别重要，不精确的分割会使得胸片图像中解剖结构的区域和特定解剖结构的ANN所训练的解剖结构子区域不对应，这会导致VDE软组织图像中的骨骼结构不能受到很好地抑制。图 5就是一个肺野分割不精确导致骨骼分割失败的案例。从图 5中可以看到，虽然有些骨头被抑制了但是锁骨还是存在。当把肺野手动地分解成8个解剖段时，对锁骨的抑制效果要比自动分割好得多。

图 5 分割不完整导致骨骼抑制失败

Fig. 5 Illustration of incomplete suppression caused by a lung segmentation failure ((a)an original image; (b)lung field segmentation; (c)bone suppression within the segmented lung fields,the right clavicle in (c)is not suppressed)

选取了93幅正常的胸片图像用来训练M-ASM。使用重叠度Ω来计算分割的精度，即

$\Omega = \frac{{T{P_{{\rm{seg}}}}}}{{T{P_{{\rm{seg}}}} + F{P_{{\rm{seg}}}} + F{N_{{\rm{seg}}}}}}$

(20)

式中，TP_seg是正确划分的肺野区域，FP_seg是错误划分的肺野区域，FN_seg是错误划分的背景。实验中采用的118个正常的胸片图像，由于并没有真实的肺野与M-ASM分割的结果作比较，所以只给出了可视化的分割结果的评估。其中有10例的分割结果不是很理想，这是因为在训练M-ASM时，选取的正常的胸片图像都是从胶片中获取的数字化影像，但是这些数字X线光片是从CR设备中获取的。肺的分割性能在以后的工作中还有待改善。

2.2 VDE骨骼图像的平滑

采用不同的方法来平滑VDE骨骼图像，通过平滑结果的对比来证明TV最小化平滑方法的有效性。图 6(d)中肋骨的边缘和骨骼的结构都消失了，而图 6(c)中是经过高斯平滑处理的结果，其肋骨边缘与骨骼结构可看得非常清楚。

图 6 使用特定解剖结构的ANN处理后的图像

Fig. 6 Illustration of images which all by use of the anatomically specific multiple MTANNs((a) a VDE bone image with Gaussian smoothing；(b) a VDE bone image with TV-minimization-based smoothing；(c) a VDE soft-tissue image corresponding to (a)；(d) a VDE soft-tissue image corresponding to (b))

实验中改善处理的VDE骨骼图像所选用的平滑参数比原始的VDE骨骼图像的平滑参数要大。这是考虑了经过改善的骨骼抑制方法，每一组特定解剖结构的ANN只需要用一种简单的模式处理一个解剖段就行。这样能得到高于整个肺野的信噪比。从原始的没有经过平滑处理的胸片图像中分别减去原始的VDE骨骼图像和处理后的VDE骨骼图像，得到的软组织图像中，处理过的VDE软组织图像效果要优于原始的VDE软组织图像。

相比于高斯平滑的方法，TV最小化的处理时间只有1 s，因为它使用了多梯度算法。

2.3 评估

使用110个含有结节的胸片图像对提出的基于解剖结构的回归模型进行了测试。采用绝对误差对胸片图像中骨骼的抑制性能进行了定量评估,即

${E_N} = \frac{{\sum\limits_{x,y \in {R_L}} {\left| {b\left( {x,y} \right) - {f_b}\left( {x,y} \right)} \right|} }}{{{N_L}\left( {{b_{\max }} - {b_{\min }}} \right)}}$

(21)

式中，f_b(x,y)是VDE骨骼图像，b(x,y)是对应的双能量图像，R_L是肺区域，N_L是区域中像素的数目，b_max和b_min分别为双能量骨骼图像肺区R_L中的最大值和最小值。这110幅胸片的平均绝对误差为0.072，标准差为0.012。

图 7第1行展示的是正常胸片中骨骼抑制的效果，和传统MTANN获得的VDE软组织相比，本文方法中肋骨边缘、锁骨以及近肺壁侧肋骨可以得到很好的抑制，同时也能保证软组织如血管的可见性。其E_N值为0.031。

图 7 实验结果

Fig. 7 Experimental result ((a) original chest images；(b) a VDE soft-tissue image obtain by use of our new MTANN technique；(c) a VDE soft-tissue image obtained by use of our new MTANN technique; (d) the corresponding gold-standard dual-energy soft-tissue image)

图 7第2行所示的胸片图像中，结节不仅被肋骨所遮挡还靠近肺壁。在传统的方法中，靠近肺壁的肋骨分割的都不是很成功，并且此区域中结节的对比度和原始的胸片非常相近。采用新提出的方法，在抑制结节周围肋骨的同时结节也可以得到保留，而且保留的结节比原始胸片中的结节要清楚的多，其E_N值为0.050。

图 7第3行的胸片图像中，结节位于肺门区，其E_N值为0.029。相比于方法中出现的肺结节的消失和边界的模糊，本文方法能很好地保留结节的对比度和形状。

3 讨 论

4 结 论

参考文献

[1] Murray C J L, Lopez A D. Mortality by cause for eight regions of the world:global burden of disease study[J]. The Lancet , 1997, 349 (9061) : 1269–1276. DOI:10.1016/S0140-6736(96)07493-4
[2] Holleb A I, Fink D J, Murphy G P. American Cancer Society Textbook of Clinical Oncology[M]. New York: American Cancer Society, 1991 .
[3] Austin J H, Romney B M, Goldsmith L S. Missed bronchogenic carcinoma:radiographic findings in 27 patients with a potentially resectable lesion evident in retrospect[J]. Radiology , 1992, 182 (1) : 115–122. DOI:10.1148/radiology.182.1.1727272
[4] Shah P K, Austin J H M, White C S, et al. Missed non-small cell lung cancer:radiographic findings of potentially resectable lesions evident only in retrospect[J]. Radiology , 2003, 226 (1) : 235–241. DOI:10.1148/radiol.2261011924
[5] Li F, Engelmann R, Doi K, et al. Improved detection of small lung cancers with dual-energy subtraction chest radiography[J]. American Journal of Roentgenology , 2008, 190 (4) : 886–891. DOI:10.2214/AJR.07.2875
[6] Kido S, Kuriyama K, Kuroda C, et al. Detection of simulated pulmonary nodules by single-exposure dual-energy computed radiography of the chest:effect of a computer-aided diagnosis system (Part 2)[J]. European Journal of Radiology , 2002, 44 (3) : 205–209. DOI:10.1016/S0720-048X(02)00269-3
[7] Kido S, Nakamura H, Ito W, et al. Computerized detection of pulmonary nodules by single-exposure dual-energy computed radiography of the chest (part 1)[J]. European Journal of Radiology , 2002, 44 (3) : 198–204. DOI:10.1016/S0720-048X(02)00268-1
[8] Luo Q, Zhang G H. Clinical evaluation of pulmonary nodules and Calification in nodules with dual-energy subtraction digital radiography[J]. Journal of Practical Radiology , 2006, 22 (1) : 30–32. [ 罗清, 张国华. 数字化双能量减影技术在胸部肺结节及结节内钙化诊断中的应用价值[J]. 实用放射学杂志 , 2006, 22 (1) : 30–32. DOI:10.3321/j.issn:0253-2239.2008.03.009 ]
[9] Glocker R, Frohnmayer W. Vberödie r ntgenspektroskopische bestimmung des gewichtsanteiles eines elementes in gemengen und verbindungen[J]. Annalen der Physik , 1925, 381 (4) : 369–395. DOI:10.1002/andp.19253810403
[10] Jacobson B, Mackay R S. Radiological contrast enhancing methods[J]. Advances in Biological and Medical Physics , 1959, 6 : 201–261. DOI:10.1016/B978-1-4832-3112-9.50010-4
[11] Suzuki K, Abe H, MacMahon H, et al. Image-processing technique for suppressing ribs in chest radiographs by means of massive training artificial neural network[J]. IEEE Transactions on Medical Imaging , 2006, 25 (4) : 406–416. DOI:10.1109/TMI.2006.871549
[12] Suzuki K,Abe H,Li F,et al.Suppression of the contrast of ribs in chest radiographs by means of massive training artificial neural network[C]//Proceedings of SPIE 5370,Medical Imaging 2004:Image Processing.San Diego,CA:SPIE,2004.[DOI:10.1117/12.536436]
[13] Suzuki K. Pixel-based machine learning in medical imaging[J]. International Journal of Biomedical Imaging , 2012, 2012 : #792079. DOI:10.1155/2012/792079
[14] Oda S, Awai K, Suzuki K, et al. Performance of radiologists in detection of small pulmonary nodules on chest radiographs:effect of rib suppression with a massive-training artificial neural network[J]. American Journal of Roentgenology , 2009, 193 (5) : W397–W402. DOI:10.2214/AJR.09.2431
[15] Ahmed B,Rasheed T,Khan M A U,et al.Rib suppression for enhancing frontal chest radiographs using independent component analysis[C]//Proceedings of the 8th International Conference.Berlin Heidelberg:Springer,2007:300-308.[DOI:10.1007/978-3-540-71629-7_34]
[16] Loog M, van Ginneken B, Schilham A M R. Filter learning:application to suppression of bony structures from chest radiographs[J]. Medical Image Analysis , 2006, 10 (6) : 826–840. DOI:10.1016/j.media.2006.06.002
[17] Suzuki K, Armato Ⅲ S G, Li F, et al. Massive training artificial neural network (MTANN) for reduction of false positives in computerized detection of lung nodules in low-dose computed tomography[J]. Medical Physics , 2003, 30 (7) : #1602. DOI:10.1118/1.1580485
[18] Suzuki K, Horiba I, Ikegaya K, et al. Recognition of coronary arterial stenosis using neural network on DSA system[J]. Systems and Computers in Japan , 1995, 26 (8) : 66–74. DOI:10.1002/scj.4690260807
[19] Suzuki K, Horiba I, Sugie N. Neural edge enhancer for supervised edge enhancement from noisy images[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence , 2003, 25 (12) : 1582–1596. DOI:10.1109/TPAMI.2003.1251151
[20] Xu J W, Suzuki K. Massive-training support vector regression and Gaussian process for false-positive reduction in computer-aided detection of polyps in CT colonography[J]. Medical Physics , 2011, 38 (4) : #1888. DOI:10.1118/1.3562898
[21] Van Ginneken B, Stegmann M B, Loog M. Segmentation of anatomical structures in chest radiographs using supervised methods:a comparative study on a public database[J]. Medical Image Analysis , 2006, 10 (1) : 19–40. DOI:10.1016/j.media.2005.02.002
[22] Chen S, Suzuki K, MacMahon H. Development and evaluation of a computer-aided diagnostic scheme for lung nodule detection in chest radiographs by means of two-stage nodule enhancement with support vector classification[J]. Medical Physics , 2011, 38 (4) : 1844–1858. DOI:10.1118/1.3561504
[23] Chen S, Suzuki K. Computerized detection of lung nodules by means of "virtual dual-energy" radiography[J]. IEEE Transactions on Biomedical Engineering , 2013, 60 (2) : 369–378. DOI:10.1109/TBME.2012.2226583
[24] van Ginneken B, Frangi A F, Staal J J, et al. Active shape model segmentation with optimal features[J]. IEEE Transactions on Medical Imaging , 2002, 21 (8) : 924–933. DOI:10.1109/TMI.2002.803121
[25] Rudin L I, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms[J]. Physica D:Nonlinear Phenomena , 1992, 60 (1-4) : 259–268. DOI:10.1016/0167-2789(92)90242-F