一种乳腺肿块超声图像超声区域分离方法与流程

本发明属于图像处理技术领域，特别涉及一种医学超声图像自动分割方法，是一种乳腺肿块超声图像超声区域分离方法。

背景技术：

乳腺肿块是妇女最常见的疾病之一，目前最有效的检测与诊断方式是乳腺x线摄影。但由于乳腺x线摄影的特异度较低，导致了大量不必要的活检，不仅给患者带来痛苦，而且增加了成本。另外，乳腺x线摄影产生的电离辐射也会给患者和医生带来健康风险。超声影像具有低成本、无侵性、实时性等特点，已成为乳腺肿块检测的重要手段之一。但乳腺超声图像具有低对比度、斑点噪声以及与组织相关的纹理结构，会给医生的诊断带来一定的困难；且不同的医生之间存在对乳腺超声的理解和诊断方面的差异。因此，乳腺超声的计算机辅助诊断系统具有很大的意义。有研究表明，乳腺超声计算机辅助诊断系统的输出结果可以为医生在乳腺肿块检测、定征和诊断决策时提供一种参考。

图像分割是乳腺超声计算机辅助诊断系统的重要环节。近几十年来，国内外在乳腺超声图像分割方面已经提出许多方法，包括阈值分割、区域生长、基于模型的方法(包括活动轮廓模型，水平集，马尔科夫随机场等)、基于图论的方法(包括图割，归一化割等)、基于机器学习的方法、基于聚类的方法(包括k均值，模糊c均值等)。这些方法存在一些主要缺点：(1)多数方法需要手工交互，诸如手工选择种子点或初始轮廓；(2)多数方法对斑点噪声敏感，且由于超声图像对比度低，具有与组织相关的纹理，故准确分割较为困难；(3)多数方法算法复杂度高，处理时间长，难以满足临床要求。

准确、快速、自动化分割是图像分割的目标。但分割的准确性、效率以及自动化水平也是相互矛盾的。近些年来，也有一些乳腺超声图像自动分割方法相继被提出，这些方法一般思路为：(1)通过某种手段自动找到肿块区域的种子点或者找到能够将肿块涵括在内的感兴趣区域；(2)利用种子点或者感兴趣区域内的某一点作为区域生长算法的种子点，进行区域生长粗分割；(3)将粗分割得到的边界作为活动轮廓模型的初始轮廓，进行活动轮廓模型精分割。但这些方法存在以下问题：(1)自动查找种子点或感兴趣区域的算法往往鲁棒性不高，仅对部分图像有效，影响了自动分割的准确性；(2)算法复杂度高，处理费时，难以达到临床上对自动分割的速度要求。分水岭算法是图像分割中应用较广的一种快速分割方法，但由于超声图像的固有特性，在淹没过程会产生过多的区域，会导致分割较慢并且准确性低的问题。

技术实现要素：

本发明旨在提出一种乳腺肿块超声图像的准确、快速、自动化分离方法。本发明将分水岭算法与均值漂移算法结合在一起，并进行了改进。首先采用金字塔均值漂移算法对乳腺超声图像进行滤波，大幅提高了图像的同质性，使其更加平滑，再使用分水岭算法对滤波后的图像进行分割处理，从而能够大量减少生成的区域个数，提高了分割的速度和精度。分水岭分割得到的结果图像中包含了若干区域，每个区域内部的灰度相同，但各区域之间灰度不同。此外，本发明还结合了临床医生的经验知识，即乳腺肿块一般位于图像的中部或上部且乳腺肿块内部平均回声强度较低，提取分水岭分割结果图像中特定感兴趣区域内的最小灰度，作为目标肿块区域的灰度，再遍历整个分水岭分割结果图像，将灰度等于最小灰度的像素赋为前景，否则赋为背景，即得到最终的肿块分割结果图像。

具体技术内容如下：

本发明的一种乳腺肿块超声图像超声区域分离方法包括以下步骤：

1.1.对于一幅乳腺肿块超声图像i，利用金字塔均值漂移算法对其进行滤波，得到滤波后的图像if；所述的金字塔均值漂移算法包括如下步骤：

(a).对乳腺肿块超声图像i进行最高层数l的高斯金字塔分解，l≥2，得到l层图像i1,…,il，图像il为金字塔的底层；

(b).对第l层图像il进行均值漂移滤波，得到结果图像(il)ms；

(c).使用拉普拉斯金字塔将第l层的均值漂移滤波结果图像(il)ms向上采样得到第l-1层图像i'l-1，再对i'l-1进行均值漂移滤波，得到(il-1)ms；

(d).重复步骤(c)，直到得到第1层均值漂移滤波结果图像(i1)ns，即金字塔均值漂移滤波后的图像if；

1.2.利用分水岭算法对图像if进行分割，得到分割后的图像is,该步骤包括：

1.2.1输入图像if，定义起始变量i＝1，i表示圆盘型结构元素bi的半径，i＝2表示将梯度图像初始划分为2个等级；

1.2.2利用矢量梯度方法计算图像if的梯度图像g；

1.2.3对梯度图像g进行开重建运算，得到第一次修正后的梯度图像grec(1)；

1.2.4对梯度图像grec(1)进行标准分水岭变换，统计分割后的区域数目为num(1)；

1.2.5自适应梯度重建，得到修正后的梯度图像gm，包括以下步骤：

(a)i加1，圆盘型结构元素bi为尺寸(2i+1)×(2i+1)的圆形结构，对grec(i-1)进行开重建，得到修正后的梯度图像grec(i)；

(b)修正梯度图像gm＝σi＝1igrec(i)×bw(g,t,i),

其中，二值化运算

t＝max(g)-min(g)l表示梯度等级间隔；

(c)若i≤l，返回步骤(a)，否则进入下一步；

(6)对gm进行标准分水岭变换，得到分割后的结果为l²，统计分割后的区域数目为num(l)；

(7)判断分割后的区域数目变化，如果num(i)＝num(i-1)，则输出最终的分割结果l²；否则，加1，返回步骤1.2.5。

1.3.计算is中特定感兴趣区域内的最小灰度mingray；

1.4.遍历is，对于灰度等于mingray的像素，将其灰度赋为前景，即为目标肿块区域，否则赋为背景，得到最终的肿块分割结果二值图像。

进一步的步骤1.1中，还可以采用自适应多尺度形态学滤波，其包括以下步骤：

(a)选取结构元素与尺度

多尺度结构元素定义为:nb＝b⊕b⊕b..⊕b；其中:n为尺度参数；⊕为形态学膨胀运算，即大尺度结构元素，是由基本元素b膨胀而来的；

(b)单尺度下滤波

设f(x，y)为输入图像，fi(x，y)为输出图像，bi(x，y)为结构元素(1≤i≤3)，则形态学开运算与闭运算分别为形态学开运算:

f(x，y)·bi(x，y)＝(f(x，y)θbi(x，y))⊕bi(x，y)；

形态学闭运算:

f(x，y)·bi(x，y)＝(f(x，y)⊕bi(x，y))θbi(x，y)；

其中，⊕和θ分别表示形态学膨胀和腐蚀运算；形态学开闭运算滤波结果为：

fi(x，y)＝(f(x，y)·bi(x，y))·bi(x，y)；

(c)多尺度自适应滤波结果

设f(x，y)为各尺度下形态滤波得到的图像加权合成，fi(x，y)为各尺度下滤波后的图像，βi为不同尺度结构元素的权值，则多尺度形态学滤波结果为

自适应权值βi的选取方法为:设b1、b2、b3可填入图像f(x，y)的次数分别为α1、α2、α3，则它们对应的自适应权重值为：

可填入次数α1、α2、α3的求值；设f1为第一次腐蚀滤波后得到的图像:

f1＝fθb1

下一次进行腐蚀滤波后的图像:

f2＝f1θb1；

...

直到{(x，y)/fn－fn－1＝0}，则α1＝n，可同理求出α2、α3。

进一步，所述自适应多尺度形态学滤波，所采用的多尺度形态学梯度算法为：

其中:bi(1≤i≤n)为结构元素，n为尺度参数。

进一步，在所述多尺度形态学滤波完成后，应用基于形态学的扩展最小变换技术进行梯度图像的标记提取，根据自适应获得的阈值参数h对梯度图像进行二值化标记，设经过h-minima变化后获得的标记图像为δimark，则:δimark＝hmin(δi，h)，

其中:hmin(·)表示h-minima变换；δimark为二值标记图像，使用“1”元素标记局部极小值区域；

其中:m0表示梯度图像的均值；m1表示梯度图像局部极小值的均值；m2表示梯度图像局部极大值的均值；m2－m1为汇水盆地的平均深度。

进一步，所述的步骤1.3中，特定感兴趣区域为矩形，矩形的左上角坐标为(λx·width，λy·height)，坐标原点为乳腺肿块超声图像i的左上角，矩形的宽度和高度分别为λw·width和λh·height，width和height分别表示i的宽度和高度，0.2≤λx≤0.3，0.05≤λy≤0.1，0.2≤λw≤0.3，0.3≤λh≤0.4。

有益效果：

1、本发明针对分水岭的过分割问题，利用自适应结构元素对梯度图像进行重建，从而克服了传统的单一结构元素导致分割精度较低的问题。

2、本发明针对现有的彩色图像分割方法性能依赖于参数设定的问题，利用自适应形态学重建理论结合分水岭分割区域的封闭性，在不需要设定任何参数的情况下，得到稳定的分割结果。

3.采用在分水岭变换前，对原始图像进行多尺度自适应形态学滤波与标记提取处理；在分水岭变换后，通过区域相似性与边界相似性的合并准则，进行区域合并，完成分水岭变换的改进。实验结果表明，本发明的算法不仅有效地抑制了过分割现象，而且图像的边缘细节也很到很好的保留，能够较好地满足医学图像的使用要求，具有较强的实用价值。

4.本发明方法在分割乳腺肿块超声图像时自动化水平高，是完全、充分的自动化分割，避免了手工交互。

5.本发明算法复杂度低，容易实现，处理速度比现有的乳腺肿块自动分割方法快一个数量级。

6.本发明算法经过验证，除少数图像边界十分模糊外，均能快速、准确、自动提取乳腺肿块边界，达到预期目的。

附图说明

本发明可以通过附图给出的非限定性实施例进一步说明；

图1：本发明方法的流程图。

具体实施方式

为了使本领域的技术人员可以更好地理解本发明，下面结合附图和实施例对本发明技术方案进一步说明。

结合附图和实际例子对提取过程进行具体的描述。下面分步进行介绍：

1.使用金字塔均值漂移算法对一幅待分割的原始乳腺肿块超声图像i进行滤波，得到滤波后的图像if。具体实施步骤如下：

(1).对乳腺肿块超声图像i进行最高层数l的高斯金字塔分解，l≥2，得到l层图像i1,…,il，图像il为金字塔的底层；

(2).对第l层图像il进行均值漂移滤波，得到结果图像(il)ms，均值漂移滤波是常用的滤波方法，其实施步骤分为两步：不连续性保持滤波和均值漂移聚类，以下进行简要介绍。

(2.1)不连续性保持滤波的实施步骤为：

设xi表示原始图像中第i个像素的灰度值，i≤n，n表示原始图像中像素总数，zi为滤波后的图像中第i个像素的灰度值，这些像素在空域-值域联合域内表达，

(a)对于第i个像素，初始化步数j＝1，yi,1＝xi；

(b)计算对应于第i个像素的第j+1步的灰度值yi,j+1，直到(yi,j+1-yi,c)≤ε，停止计算，yi,c为对应于第i个像素的收敛点的灰度值，其中，yi,j+1的计算方式为：h为核大小，g(x)＝-k'(x)，k(x)为核轮廓函数；

(c)第i个像素滤波后的灰度值zi赋为即：在处的滤波后的像素的灰度值被赋值为收敛点的灰度值其中s，r分别表示空域分量和值域分量。

(2.2)均值漂移聚类的实施步骤为：

(a)使用不连续性保持滤波，保存关于每个收敛点的灰度值yi,c的所有信息；

(b)将所有zi按照在空域用核hs和在值域用核hr聚类得到m个类别cp,p＝1,...,m；

(c)对均值漂移滤波后的每个像素的灰度值ti赋值ti＝{p|zi∈cp}，得到均值漂移滤波后的结果图像(il)ms。

(3).使用拉普拉斯金字塔将第l层的均值漂移滤波结果图像(il)ms向上采样得到第l-1层图像i'l-1，再对i'l-1进行均值漂移滤波，得到(il-1)ms。

(4).重复步骤(3)，直到得到第1层均值漂移滤波结果图像(i1)ms，即金字塔均值漂移滤波后的图像if。

本实施例中设定高斯金字塔分解最高层数l为3；均值漂移滤波的空域核大小为20，值域核大小为40。

2.使用分水岭算法对滤波后的图像if进行分割，得到分割后的图像is。分水岭算法的具体实施步骤如下：

(1)初始化：输入图像，图像尺寸为m×n，m和n分别表示图像的高度和宽度，定义起始变量i＝1，i表示圆盘型结构元素bi的半径，i＝2表示将梯度图像初始划分为2个等级；

(2)利用矢量梯度方法计算图像f的梯度图像g，

(a)利用sobel算子分别计算图像f的三个分量(r、g、b)在水平和垂直方向的梯度，其

中，

(b)计算水平和垂直方向的矢量梯度

(c)根据gxx、gyy和gxy计算方向角度

(d)由方向角度θ确定矢量梯度g

g(x,y)＝max(g1(x,y),g2(x,y))

g是g¹、g²在每个点(x,y)处的最大值。

(3)对梯度图像g进行开重建运算，得到重建后的图像grec(1)；

(a)初始化；b1是尺寸为(2i+1)×(2i+1)的圆形结构元素；i＝1，k＝1；

(b)定义掩膜图像为gmask，标记图像fmarker；

gmask＝g

fmarker＝gθbi

h1＝fmarker

(c)腐蚀重构运算；

(d)判断，如果hk+1＝hk，则得到腐蚀重建结果g^ε＝hk；否则，k加1，返回步骤(c)；

(e)变换掩膜图像gmask和标记图像fmarker；i＝1，k＝1；

g'mask＝(g^ε)c

f′marker＝(g^ε)cθbi

h′1＝f′marker

其中，(g^ε)c表示g^ε的补运算。

膨胀重建运算；

(f)判断，如果h'k+1＝h'k，则得到开重建结果grec(1)＝(h'k)c；否则，k加1，返回步骤(f)；

(4)对重建后的梯度图像进行标准分水岭变换，分割后的结果统计分割后的区域数目为num(1)；

其中，ws表示标准分水岭分割。

(5)自适应梯度重建，得到修正后的梯度图像gm；

(a)结构元素bi为尺寸(2i+1)×(2i+1)的圆形结构，2≤i≤l，对grec(i-1)进行开重建，得到重建后的梯度图像为grec(i)；

(b)修正梯度图像；

bw表示二值化运算：

表示梯度等级间隔。

(6)对修正后的梯度图像gm进行标准分水岭变换，分割结果为l²，并统计分割后的区域数目为num(i)；

l2＝ws(gm)

其中，ws表示标准分水岭分割。

(7)判断分割后的区域数目变化，如果num(l)＝num(l-1)，则输出最终的分割结果l2；否则，l加1，返回步骤5；

根据肿块一般位于图像中部或上部并且平均回声强度较低的经验知识，计算分水岭分割结果图像is中特定感兴趣区域内的最小灰度mingray，特定感兴趣区域为矩形，矩形的左上角坐标为(λx·width，λy·height)，坐标原点为乳腺肿块超声图像i的左上角，矩形的宽度和高度分别为λw·width和λh·height，width和height分别表示i的宽度和高度，0.2≤λx≤0.3，0.05≤λy≤0.1，0.2≤λw≤0.3，0.3≤λh≤0.4，本实施例中设定λx＝0.25，λy＝0.08，λw＝0.25，λh＝0.31。

遍历分水岭分割结果图像is，对于灰度等于mingray的像素，将其灰度赋为前景，即为目标肿块区域，否则赋为背景，得到最终的肿块分割结果二值图像，本实施例中设定前景为255，背景为0。

本实施例还可采用自适应多尺度形态学滤波，其主要包括以下几个步骤:

a)结构元素与尺度的选取

多尺度结构元素定义为:nb＝b⊕b⊕b..⊕b；

其中:n为尺度参数；⊕为形态学膨胀运算，即大尺度结构元素

是由基本元素b膨胀而来的。

实验结果表明，多尺度形态学滤波中选择圆形3×3的结构元素，尺度参数n＝3，更加适合医学图像的滤波。结构元素b2(5×5)和结构元素b3(7×7)分别是由结构元素b1经过1次与2次膨胀得到的。

b)单尺度下滤波结果

设f(x，y)为输入图像，fi(x，y)为输出图像，bi(x，y)为结构元素(1≤i≤3)，则形态学开运算与闭运算分别为形态学开运算:

f(x，y)·bi(x，y)＝(f(x，y)θbi(x，y))⊕bi(x，y)

形态学闭运算:

f(x，y)·bi(x，y)＝(f(x，y)⊕bi(x，y))θbi(x，y)

其中，⊕和θ分别表示形态学膨胀和腐蚀运算；形态学开闭运算滤波结果为

fi(x，y)＝(f(x，y)·bi(x，y))·bi(x，y)

c)多尺度自适应滤波结果

设f(x，y)为各尺度下形态滤波得到的图像加权合成，fi(x，y)为各尺度下滤波后的图像，βi为不同尺度结构元素的权值，则多尺度形态学滤波结果为：

自适应权值βi的选取:当不同尺度的结构元素探测一幅图像时，它们可填入的次数是不同的，尺度越小，可填入的次数越多；反之，则少。设b1、b2、b3可填入图像f(x，y)的次数分别为α1、α2、α3，则它们对应的自适应权重值为：

可填入次数α1、α2、α3的求值。设f1为第一次腐蚀滤波后得到的图像:

f1＝fθb1

下一次进行腐蚀滤波后的图像:

f2＝f1θb1；

...

直到{(x，y)/fn－fn－1＝0}，则α1＝n。同理求出α2、α3。

多尺度形态学梯度计算

为了抑制过分割现象，通过把梯度图像作为分水岭的输入图像，传统形态学梯度的性能取决于结构元素的形状和大小，为了弥补这一不足，本文采用多尺度形态学梯度算法，定义如下:

其中:bi(1≤i≤n)为结构元素，n为尺度参数。结构元素和尺度选值与多尺度滤波选值相同，选取圆形3×3的结构元素，尺度参数选取n＝3。

自适应标记提取

原始图像经过多尺度形态学滤波后，图像中仍然存在一些暗噪声以及不规则细纹。为了进一步抑制分水岭过分割现象，本文在多尺度形态学梯度图像的基础上，采用目前应用最广泛的标记提取技术———基于形态学的扩展最小变换技术进行梯度图像的标记提取。基于形态学的扩展最小变换通过人为设定阈值参数h，消除深度小于h的局部极小值，从而减少了过分割区域，实现抑制过分割现象的目的。然而参数h的选择存在噪声抑制与细节保留之间的矛盾，过小时噪声不能得到很好的抑制，导致过分割现象仍然严重；过大时图像边缘细节信息被模糊掉，不能得到很好的保留，而这些边缘细节信息是分析病情与临床诊断的关键。针对这一不足，本文采用参数h的自适应取值方法，实现自适应标记提取，定义如下:

其中:m0表示梯度图像的均值；m1表示梯度图像局部极小值的均值；m2表示梯度图像局部极大值的均值；m2－m1为汇水盆地的平均深度。

根据自适应获得的阈值参数h对梯度图像进行二值化标记，设经过h-minima变化后获得的标记图像为δimark，则:δimark＝hmin(δi，h)

其中:hmin(·)表示h-minima变换；δimark为二值标记图像，使用“1”元素标记局部极小值区域。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。