CN104143101A - 一种基于超声图像的乳腺肿瘤区域自动识别方法 - Google Patents

Abstract

本发明公开了一种基于超声图像的乳腺肿瘤区域自动识别方法，包含以下步骤：获取乳腺超声图像，并对其进行预处理；使用图像分割算法对预处理后的超声图像进行分割，得到多个分割子区域；提取超声图像的灰度直方图、纹理特征、梯度特征、形态特征，再加上二维位置信息，得到高维度的特征向量；对高维度的特征向量，通过基于双聚类的特征排序和选择方法来选取最有效的特征子集；选择出最有效的特征子集后，使用分类器对其进行学习分类，即可自动识别乳腺肿瘤区域。本发明的方法，可以从乳腺肿瘤超声图像的分割结果中自动识别出肿瘤区域，从而提高计算机辅助诊断的自动化性能，减少临床医生的手工操作，降低临床医生的主观影响。

Description

一种基于超声图像的乳腺肿瘤区域自动识别方法

技术领域

本发明涉及计算机辅助诊断领域，特别涉及一种基于超声图像的乳腺肿瘤区域自动识别方法。

背景技术

如今，乳腺癌已成为威胁女性生命的主要疾病之一，然而人类对其具体病理还处于未知状态，早发现早治疗便成了得到有效治疗的主要方式。随着社会的发展和技术的进步，超声成像以其非侵入、成本低、简单可靠的特点成为了乳腺肿瘤的主要检测方式之一。

然而，依靠超声图像的诊断是经验依赖性的，不可避免地会受到临床医生的主观因素影响，从而会大大降低诊断结果的准确性和可靠性。因此，为了避免主观因素的影响和提高诊断结果的准确性和可靠性，计算机辅助诊断的技术被提了出来并取得了很多的研究成果。

在已提出的基于超声图像的乳腺癌计算机辅助诊断***中，对超声图像进行分割后，通常需要临床医生从分割出的各子区域中人工选择出肿瘤区域，才能继续进行后续的分析操作，这严重降低了计算机辅助诊断的自动化性能，而且引进了临床医生的主观性。

因此人们需要一种新的肿瘤区域识别方法来克服现有技术中存在的缺点。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供一种基于超声图像的乳腺肿瘤区域自动识别方法。

本发明的目的通过以下的技术方案实现：

一种基于超声图像的乳腺肿瘤区域自动识别方法，包含以下顺序的步骤：

S1.获取乳腺超声图像，并对其进行预处理；

S2.使用图像分割算法对预处理后的超声图像进行分割，得到多个分割子区域；

S3.提取超声图像的灰度直方图、纹理特征、梯度特征、形态特征，再加上二维位置信息，得到高维度的特征向量；

S4.对高维度的特征向量，通过基于双聚类的特征排序和选择方法来选取最有效的特征子集，具体包括以下步骤：

(1)提取双聚类种子点：对每一单独列内的元素做聚类，每一个聚类的结果都被看做一个双聚类的种子点；对于给定的一个R行C列的数据矩阵，应用一种凝聚分等级聚类算法到目标矩阵的每一列；

(2)启发式构造双聚类：利用一种启发式的方法对第一步中提取的双聚类种子点做扩展，根据预先设定的准则构造双聚类；

(3)消除冗余的双聚类：首先，对双聚类集合按照每个双聚类包含的列数目升序排序；然后，从包含最少列数目的双聚类开始迭代，检测它是否属于排在后面的双聚类的子集，如果是，则消除该双聚类；

(4)特征排序和选择：对数据矩阵提取出双聚类之后，需要从这些双聚类中提取出能够衡量各个特征的信息，即计算出每个特征相应的双聚类分值，分值越高，表明该特征越能表达该数据集，与其它特征的联系越紧密，然后根据所需要的特征数目选择出分值最高的特征子集，即为最有效的特征子集；

S5.选择出最有效的特征子集后，使用分类器对其进行学习分类，即可达到乳腺肿瘤区域自动识别的目的。

步骤S1中，所述的预处理采用全变分模型，通过最小化全变分来达到滤波去噪的目的，全变分定义为

$\mathrm{TV}\left[u\right]=\underset{\mathrm{\Ω}}{\∫}|\▿u|\mathrm{dx}$

其中，Ω表示连续的信号域，表示图像梯度，dx则表示Ω的一个元素。

由于要在平滑斑点噪声的同时不损害图像的边缘特性，选用了对消除斑点噪声和边缘增强都十分有效的全变分模型(Total Variation，TV)。

步骤S2中，所述的图像分割算法阈值法、聚类法、最大后验法、马尔可夫随机场法、动态规划法、贝叶斯水平集法、基于图论的分割算法中的一种。

步骤S4的分步骤(1)中，所述的一种凝聚分等级聚类算法具体如下：

a、计算所有原始数据的距离的最大值和最小值；

b、将所有原始数据每个元素初始化为单独一类；

c、建立距离数组，记录各类之间的距离，初始化为各个原始数据的距离；

d、对所有数据做以下循环：找到距离数组的最小值，如果该值小于阈值τ，则将这两个类归为一类，重新计算距离数组，直到距离数组的最小值大于或等于τ，终止循环；

凝聚分等级聚类后，将得到目标矩阵的一系列双聚类种子点：

[C_s(i,j),N_cl(j)]＝HC(j,τ),j＝1,...,C

其中，HC代表凝聚分等级聚类算法，τ是聚类算法预设的阈值，C_s(i,j)表示第j列的第i个双聚类的种子点，N_cl(j)表示第j列双聚类种子点的总数目。

步骤S4的分步骤(2)中，所述的对双聚类种子点做扩展，其方法包含以下步骤：

a、按照每个种子点包含的行数做升序排序；

b、从包含最小行数R_j的种子点开始，向其他列扩展相应的行，形成一个R_j行C列的矩阵M；

c、对于矩阵M，计算其均方残基得分：

$\mathrm{MSRS}=\frac{1}{\left|R\right|\left|C\right|}\underset{i\∈R,j\∈C}{\mathrm{\Σ}}{(e_{\mathrm{ij}}-e_{\mathrm{iC}}-e_{\mathrm{Rj}}+e_{\mathrm{RC}})}^2$

$e_{\mathrm{iC}}=\frac{1}{\left|C\right|}\underset{j\∈C}{\mathrm{\Σ}}{e}_{\mathrm{ij}},e_{\mathrm{Rj}}=\frac{1}{\left|R\right|}\underset{i\∈R}{\mathrm{\Σ}}{e}_{\mathrm{ij}},e_{\mathrm{RC}}=\frac{1}{\left|R\right|\left|C\right|}\underset{i\∈R,j\∈C}{\mathrm{\Σ}}{e}_{\mathrm{ij}}$

其中，R是矩阵的行数，C是矩阵的列数，e_ij是该矩阵第i行第j列的元素值；如果计算得出的MSRS大于设定的阈值，则对于矩阵中的每一个点，计算除去该点所在的行或者列，计算新矩阵的MSRS,以新旧矩阵MSRS得分之差来衡量该点所在行或者列对MSRS的贡献值，删除MSRS贡献值最大的点所在的行或者列；不断迭代这一过程，直到矩阵的MSRS值少于预设定的阈值，则该矩阵为最终求得的一个双聚类。

所述的步骤S4的分步骤(4)，具体如下：

A、对数据矩阵提取出双聚类之后，需要从这些双聚类中提取出能够衡量各个特征的信息：基于双聚类的特征排序方案考虑了特征的相关性和样本的分离性两个因素，其中特征的相关性用相关性得分表示，用于衡量特征子集中的特征之间联系的紧密程度，样本的分离性用分离性得分表示，用于衡量一个特征的独立程度；对于第k个特征，如果它出现在双聚类子集Z_k中的任一双聚类中，则该特征的相关性得分和分类性得分分别表示如下：

$\mathrm{correlation}\_\mathrm{score}=\underset{i=1}{\overset{n_{b,k}}{\mathrm{\Σ}}}\frac{n_{f,k}\left(i\right)}{n_c}$

$\mathrm{separability}\_\mathrm{score}=\frac{n_{s,k}}{n_r}\sqrt{\underset{i=1}{\overset{n_{b,k}}{\mathrm{\Σ}}}{({\mathrm{\μ}}_{i,k}-{\mathrm{\μ}}_{a,k})}^2/n_{b,k}}$

其中，n_b,k表示Z_k中双聚类的数目，n_f,k(i)表示Z_k中第i个双聚类包含的特征的数目，n_s,k表示Z_k中不重复的行的数目，μ_i,k表示Z_k中第i个双聚类的第k个特征的均值，μ_a,k表示μ_i,k的均值；

第k个特征的最终得分用双聚类分值表示：

$\mathrm{bicluster}\_\mathrm{score}=a\·\stackrel{\‾}{\mathrm{correlation}\_\mathrm{score}}+\stackrel{\‾}{\mathrm{separability}\_\mathrm{score}}$

其中，和分别表示correlation_score和separability_score归一化结果；a是一个常数，用于平衡correlation_score和separability_score对最终得分的贡献；

特征的bicluster_score分值越高，表明该特征越能表达该数据集，与其它特征的联系越紧密；

B、根据所需要的特征数目选择出分值最高的特征子集。

步骤S5中，所述的分类器K最近邻分类器、神经网络分类器、支持向量机分类器中的一种。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明主要工作原理为：先用机器学习方法对肿瘤区域和非肿瘤区域进行学习，在图像分割后，使用训练好的分类器对肿瘤区域和非肿瘤区域进行分类识别，自动识别出肿瘤区域；利用本发明所提供的方法，可以从乳腺肿瘤超声图像的分割结果中自动识别出肿瘤区域，从而提高计算机辅助诊断的自动化性能，减少临床医生的手工操作，降低临床医生的主观影响。

2、由于成像环境和成像原理等复杂因素，超声图像具有先天性的斑点噪声、多阴影和低对比度等特征，使得超声图像的总体质量较差。本发明对超声图像进行分割之前，对超声图像进行去噪等预处理，抑制斑点噪声所带来的不良影响。

3、本发明的特征提取选取灰度直方图、纹理特征、梯度特征、形态特征，是因为痛过研究表明这些特征对于超声图像中的乳腺肿瘤具有较强的表达能力，能够清楚描述肿瘤区域和非肿瘤区域的特性及两者之间的异同。

4、本发明进行特征选择的原因：研究表明，以适当的方法降低特征的维度，去除冗余低效的特征，选择最具有表达能力的特征组合，不仅不会降低分类识别的准确率，而且还可以在一定程度上改善分类性能。

附图说明

图1为本发明所述的一种基于超声图像的乳腺肿瘤区域自动识别方法的流程图；

图2为图1所述方法的GLCM计算例子(a)数据矩阵(b)GLCM示意图；

图3为图1所述方法的双聚类种子点扩展示意图；

图4为图1所述方法的自动识别准确率随特征集大小的变化图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

采用C语言等各类编程语言均能很好地实现本发明。

如图1所示，一种基于超声图像的乳腺肿瘤区域自动识别方法，包括以下五个步骤：

步骤(1)预处理：由于成像环境和成像原理等复杂因素，超声图像具有先天性的斑点噪声、多阴影和低对比度等特征，使得超声图像的总体质量较差。因此，对超声图像进行分割之前，需要进行去噪等预处理，抑制斑点噪声所带来的不良影响。由于要在平滑斑点噪声的同时不损害图像的边缘特性，本发明选用了对消除斑点噪声和边缘增强都十分有效的全变分模型(Total Variation，TV)，通过最小化全变分来达到滤波去噪的目的，全变分定义为

$\mathrm{TV}\left[u\right]=\underset{\mathrm{\Ω}}{\∫}|\▿u|\mathrm{dx}$

其中，Ω表示连续的信号域，表示图像梯度，dx则表示Ω的一个元素。

步骤(2)图像分割：使用图像分割算法对预处理后的超声图像进行分割，得到众多的分割子区域。可选用的图像分割算法有很多，如：阈值法、聚类法、最大后验法、马尔可夫随机场法、动态规划法、贝叶斯水平集法和基于图论的分割算法等。

步骤(3)特征提取：研究表明，灰度直方图、纹理特征、梯度特征、形态特征等对超声图像中的乳腺肿瘤具有较强的表达能力，为了描述清楚肿瘤区域和非肿瘤区域的特性及两者之间的异同，本发明选择对上述几类特征都进行提取，再加上两维的位置信息，可得到73维的特征向量。具体如下：

灰度直方图：如果将乳腺肿瘤超声图像的灰度等级划分为L级，则对于任一分割子区域的灰度直方图可以通过以下概率分布式来表述：

p_i＝n_i/N,p_i≥0,i＝1,2,...,L

其中，N为该子区域的像素点数，n_i为第i灰度级的像素点数。然而，直接从灰度直方图得到图像子区域的特征不够直观，因此计算以下几种常用的描述图像灰度直方图的一阶统计量作为子区域的量化特征：

(1)灰度均值：

$\mathrm{\μ}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}i\·p\left(i\right)$

(2)灰度标准差：

${\mathrm{\σ}}^2=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{(i-\mathrm{\μ})}^{2}p\left(i\right)$

(3)偏斜度：

$\mathrm{Skew}=\frac{1}{{\mathrm{\σ}}^2}\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{(i-\mathrm{\μ})}^{3}p\left(i\right)$

(4)峰度：

$\mathrm{Kurtosis}=\frac{1}{{\mathrm{\σ}}^4}\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{(i-\mathrm{\μ})}^{4}p\left(i\right)-3$

(5)能量：

$\mathrm{Energy}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}p{\left(i\right)}^2$

(6)熵：

$\mathrm{Entropy}=-\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}p\left(i\right){\log }_{2}p\left(i\right)$

纹理特征：灰度共生矩阵是一种较为常用的描述图像纹理特征的方法，它实际上描述的是图像中像素灰度之间的空间依赖特性。灰度共生矩阵中的每一个点记录的是互成角度θ，相互距离为d的点对在整个子区域中出现的概率，可以用以下公式来表述：

p(i,j,d,θ)＝#{(x₁,y₁)(x₂,y₂)|f(x₁,y₁)＝i,f(x₂,y₂)＝j,

|(x₁,y₁)-(x₂,y₂)|＝d,∠((x₁,y₁)(x₂,y₂))＝θ}

其中，(x₁,y₁)，(x₂,y₂)分别是点对中的两个点在子区域中的位置，f(x_n,y_n)是相应位置像素点的灰度等级，∠表示两个点之间的角度，#表示满足所述条件的点对在整个子区域中出现的总次数。图2是计算灰度共生矩阵的一个例子，其中θ＝0，d＝1。

然而，在实际应用中，很少直接将灰度共生矩阵作为图像的纹理特征，而是像求灰度直方图特征一样，计算一些统计量来表示特征。令p_ij表示灰度共生矩阵中第i行第j列元素值，L表示划分的灰度等级，则常用的统计量表达式如下：

(1)最大概率：

Max Pro＝maxp_ij

(2)差异度：

$\mathrm{Dissimilarity}=\underset{i,j=0}{\overset{L-1}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}|i-j|$

(3)对比度：

$\mathrm{Contrast}=\underset{i,j=0}{\overset{L-1}{\mathrm{\Σ}}}{p}_{i,j}{(i-j)}^2$

(4)同质度：

$\mathrm{Homogeneity}=\underset{i,j=0}{\overset{L-1}{\mathrm{\Σ}}}\frac{p_{\mathrm{ij}}}{1+{(i-j)}^2}$

(5)逆差矩：

$\mathrm{Inverse}=\underset{i,j=0}{\overset{L-1}{\mathrm{\Σ}}}\frac{p_{i,j}}{1+|i-j|/L}$

(6)能量：

$\mathrm{Energy}=\sqrt{\underset{i,j=0}{\overset{L-1}{\mathrm{\Σ}}}{p_{\mathrm{ij}}}^2}$

(7)熵：

$\mathrm{Entropy}=-\underset{i,j=0}{\overset{L-1}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}\ln p_{\mathrm{ij}}$

(8)均值：

$\mathrm{Mean}=\underset{i,j}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}}$

(9)标准差：

$\mathrm{sd}=\sqrt{\underset{i,j=0}{\overset{L-1}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}{(i-\mathrm{Mean})}^2}$

(10)相关度：

$\mathrm{Correlation}=\underset{i,j=0}{\overset{L-1}{\mathrm{\Σ}}}\frac{p_{\mathrm{ij}}(i-\mathrm{Mean})(j-\mathrm{Mean})}{{\mathrm{Sd}}^2}$

以下统计量需要用到灰度共生矩阵的边缘概率分布，四种边缘概率分布的表达式分别是：

$p_x\left(i\right)=\underset{j}{\mathrm{\Σ}}{p}_{\mathrm{ij}}$

$p_y\left(j\right)=\underset{i}{\mathrm{\Σ}}{p}_{\mathrm{ij}}$

$p_{x+y}(k=i+j)=\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}{p}_{\mathrm{ij}}$

$p_{x-y}(k=|i-j\left|\right)=\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}{p}_{\mathrm{ij}}$

(11)和均值：

$\mathrm{SumMean}=\underset{i=0}{\overset{2L-2}{\mathrm{\Σ}}}{\mathrm{ip}}_{x+y}\left(i\right)$

(12)和标准差：

$\mathrm{SumSd}=\sqrt{\underset{i=0}{\overset{2L-2}{\mathrm{\Σ}}}{(i-\mathrm{SumMean})}^2{p}_{x+y}\left(i\right)}$

(13)和熵:

$\mathrm{SumEntropy}=-\underset{i=0}{\overset{2L-2}{\mathrm{\Σ}}}{p}_{x+y}\left(i\right){\log }_2{p}_{x+y}\left(i\right)$

(14)差均值：

$\mathrm{DifMean}=\underset{i=0}{\overset{2L-2}{\mathrm{\Σ}}}{\mathrm{ip}}_{x-y}\left(i\right)$

(15)差标准差：

$\mathrm{DifSd}=\sqrt{\underset{i=0}{\overset{2L-2}{\mathrm{\Σ}}}{(i-\mathrm{DifMean})}^2{p}_{x-y}\left(i\right)}$

(16)和熵:

$\mathrm{SumEntropy}=-\underset{i=0}{\overset{2L-2}{\mathrm{\Σ}}}{p}_{x-y}\left(i\right){\log }_2{p}_{x-y}\left(i\right)$

(17)边缘熵：

$\mathrm{HX}=\mathrm{HY}=-\underset{i}{\mathrm{\Σ}}{p}_x\left(i\right){\log }_2{p}_x\left(i\right)$

(18)联合边缘熵：

$\mathrm{HXY}1=-\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}{p}_{\mathrm{ij}}{\log }_2\left(p_x\left(i\right)p_y\left(j\right)\right)$

(19)独立边缘熵：

$\mathrm{HXY}2=-\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}{p}_x\left(i\right)p_y\left(i\right){\log }_2\left(p_x\left(i\right)p_y\left(j\right)\right)$

(20)信息相关度：

$\mathrm{INFH}1=\frac{\mathrm{Entropy}-\mathrm{HXY}1}{\max (\mathrm{HX},\mathrm{HY})}$

$\mathrm{INFH}2=\sqrt{1-\exp (-2(\mathrm{HXY}2-\mathrm{Entropy}))}$

梯度直方图：

(1)计算区域块中每个像素点的梯度幅度和梯度方向。为了计算梯度幅度，需要先求得垂直梯度和水平梯度：

垂直梯度：

f_h(x,y)＝f(x+1,y)-f(x-1,y)

水平梯度：

f_v(x,y)＝f(x,y+1)-f(x,y-1)

其中，f(x,y)表示在位置(x,y)的像素点的灰度值。在此基础上，梯度幅度和梯度方向则通过以下公式分别计算：

梯度幅度：

$m(x,y)=\sqrt{f_h{(x,y)}^2+f_v{(x,y)}^2}$

梯度方向：

$\mathrm{\θ}(x,y)=\arctan \frac{f_h(x,y)}{f_v(x,y)}$

(2)直方图生成。经过第一步的计算，每个像素点都可以计算出相应的梯度幅度和梯度方向，类似于灰度直方图，累加具有相同梯度方向的像素点的梯度幅度，从而生成基于梯度的直方图。为了降低计算的复杂度和提高运算速度，360度的方向平均量化为K个区间，落在相同区间的梯度方向被认为是相同的。累加完成后，该梯度直方图可以用一个K维向量表示。

(3)对K维向量做归一化。归一化后的K维向量就是最后用于描述梯度直方图的特征。

形状特征：在目前的研究中，较为常用而且有效的形状特征主要有以下四种：周长、面积、形状复杂度和似圆度。

形状复杂度和似圆度分别通过以下公式求得：

ShapeComplexity＝Perimeter²/Area

CircularDegree＝4π·Area/Perimenter²

其中Perimeter和Area分别表示归一化后的子区域周长和面积。

位置特征：通过下列表达式计算，

L_x＝1-2·|w/2-x|/w

L_y＝1-2·|h/2-y|/h

其中，w和h分别表示超声图像的宽度和高度，(x,y)则表示子区域的中心点位置。如果子区域所处位置越接近超声图像的中心点，则其位置特征计算出来的值越大。

步骤(4)特征选择：经过特征提取后，可以从乳腺肿瘤超声图像中提取出比较高维度的特征信息。然而，研究表明，以适当的方法降低特征的维度，去除冗余低效的特征，选择最具有表达能力的特征组合，不仅不会降低分类识别的准确率，而且还可以在一定程度上改善分类性能。特征选择方法有很多，如Variance Score、Laplacian Score和Fisher Score等，本发明选用了一种基于双聚类的特征排序和选择方法来选取最适合特定分类目的的特征子集，具体操作如下：

根据本发明所针对的特定的分类目的，采用提取双聚类的算法，该双聚类算法提取的是列元素相等的双聚类模式。该双聚类提取算法主要包括三个步骤：

(1)提取双聚类种子点：对每一单独列内的元素做聚类，每一个聚类的结果都被看做一个双聚类的种子点。对于给定的一个R行C列的数据矩阵，应用一种凝聚分等级聚类算法(Agglomerative Hierarchical Clustering)到目标矩阵的每一列。凝聚分等级聚类算法的具体描述如下：

1)计算所有原始数据的距离的最大值和最小值

2)将所有原始数据每个元素初始化为单独一类

3)建立距离数组，记录各类之间的距离，初始化为各个原始数据的距离

4)对所有数据做以下循环：找到距离数组的最小值，如果该值小于阈值τ，则将这两个类归为一类，重新计算距离数组。直到距离数组的最小值大于或等于τ，终止循环。

凝聚分等级聚类后，将得到目标矩阵的一系列双聚类种子点：

[C_s(i,j),N_cl(j)]＝HC(j,τ),j＝1,...,C

其中，HC代表的是凝聚分等级聚类算法，τ是聚类算法预设的阈值，C_s(i,j)表示第j列的第i个双聚类的种子点，N_cl(j)表示第j列双聚类种子点的总数目。

(2)启发式构造双聚类：利用一种启发式的方法对第一步中提取的双聚类种子点做扩展，根据预先设定的准则构造双聚类。扩展算法如下所述：

1)按照每个种子点包含的行数做升序排序。

2)从包含最小行数R_j的种子点开始，向其他列扩展相应的行，形成一个R_j行C列的矩阵M。

3)对于矩阵M，计算其均方残基得分(Mean Square Residue Score,MSRS)：

$\mathrm{MSRS}=\frac{1}{\left|R\right|\left|C\right|}\underset{i\∈R,j\∈C}{\mathrm{\Σ}}{(e_{\mathrm{ij}}-e_{\mathrm{iC}}-e_{\mathrm{Rj}}+e_{\mathrm{RC}})}^2$

$e_{\mathrm{iC}}=\frac{1}{\left|C\right|}\underset{j\∈C}{\mathrm{\Σ}}{e}_{\mathrm{ij}},e_{\mathrm{Rj}}=\frac{1}{\left|R\right|}\underset{i\∈R}{\mathrm{\Σ}}{e}_{\mathrm{ij}},e_{\mathrm{RC}}=\frac{1}{\left|R\right|\left|C\right|}\underset{i\∈R,j\∈C}{\mathrm{\Σ}}{e}_{\mathrm{ij}}$

其中，R是矩阵的行数，C是矩阵的列数，e_ij是该矩阵第i行第j列的元素值。如果计算得出的MSRS大于设定的阈值，则对于矩阵中的每一个点，计算除去该点所在的行或者列，计算新矩阵的MSRS,以新旧矩阵MSRS得分之差来衡量该点所在行或者列对MSRS的贡献值，删除MSRS贡献值最大的点所在的行或者列。不断迭代这一过程，直到矩阵的MSRS值少于预设定的阈值，则该矩阵为最终求得的一个双聚类。图3给出了双聚类种子点扩展示意图。

(3)移除冗余的双聚类：由于在第二步中生成了大量双聚类，而且存在一些比较小的冗余的双聚类被更大的双聚类所包含的情况。去除冗余的双聚类是十分必要的，因为重复的双聚类会造成后续特征排名方案中评估样本的分离性和特征的相关性时的重复计算。消除冗余的双聚类分两步进行：首先，对双聚类集合按照每个双聚类包含的列数目升序排序；然后，从包含最少列数目的双聚类开始迭代，检测它是否属于排在后面的双聚类的子集，如果是，则消除该双聚类。

对数据矩阵提取出双聚类之后，需要从这些双聚类中提取出能够衡量各个特征的信息。基于双聚类的特征排名方案考虑了特征的相关性和样本的分离性两个因素，其中特征的相关性用相关性得分表示，用于衡量特征子集中的特征之间联系的紧密程度，样本的分离性用分离性得分表示，用于衡量一个特征的独立程度。对于第k个特征，如果它出现在双聚类子集Z_k中的任一双聚类中，则该特征的相关性得分和分类性得分分别表示如下：

$\mathrm{correlation}\_\mathrm{score}=\underset{i=1}{\overset{n_{b,k}}{\mathrm{\Σ}}}\frac{n_{f,k}\left(i\right)}{n_c}$

$\mathrm{separability}\_\mathrm{score}=\frac{n_{s,k}}{n_r}\sqrt{\underset{i=1}{\overset{n_{b,k}}{\mathrm{\Σ}}}{({\mathrm{\μ}}_{i,k}-{\mathrm{\μ}}_{a,k})}^2/n_{b,k}}$

其中，n_b,k表示Z_k中双聚类的数目，n_f,k(i)表示Z_k中第i个双聚类包含的特征的数目，n_s,k表示Z_k中不重复的行的数目，μ_i,k表示Z_k中第i个双聚类的第k个特征的均值，μ_a,k表示μ_i,k的均值。

第k个特征的最终得分用双聚类分值表示：

$\mathrm{bicluster}\_\mathrm{score}=a\·\stackrel{\‾}{\mathrm{correlation}\_\mathrm{score}}+\stackrel{\‾}{\mathrm{separability}\_\mathrm{score}}$

其中，和分别表示correlation_score和separability_score归一化结果；a是一个常数，用于平衡correlation_score和separability_score对最终得分的贡献；

特征的bicluster_score分值越高，表明该特征越能表达该数据集，与其它特征的联系越紧密。

所述步骤(5)分类识别：选择出最有效的特征子集后，使用分类器对其进行学习分类，即可达到乳腺肿瘤区域自动识别的目的。常用的分类器有：K最近邻分类器(K-nearest Neighbor,KNN)，神经网络分类器(Neural Networks，NNs)和支持向量机分类器(Support Vector Machine，SVM)等。

为了验证本发明的正确性和有效性，做了以下实验：

实验中，使用VC++6.0对本发明进行实现。实验是基于46张乳腺肿瘤超声图像(恶性肿瘤和良性肿瘤各23张)来验证的，每张图像的大小均为400×300。这46张乳腺肿瘤超声图像是由中山大学附属肿瘤医院提供的，通过飞利浦医学超声图像仪器HDI5000SonoCT(换能器：L12-550mm宽频带线阵换能器；频率：7.1MHz)获取。实验测试主机参数为：***为Windows732位操作***，处理器为AMD5000+(2.60GHz)，内存为Kingston(2.00GB)。

乳腺肿瘤超声图像在经过图像分割后，图像中的每一子区域均被看做一个样本实例，每个样本实例被标记为肿瘤(正样本)或非肿瘤(负样本)，分别用1和-1表示；然后，对每一子区域提取出灰度直方图、灰度共生矩阵、梯度直方图、形状和位置等五类特征信息，每个样本实例一共可以提取出73维的特征，即每个样本实例可以用一个74维的向量表示(1维的类别信息和73维的特征信息)。对所有样本实例提取特征后，使用5折交叉验证进行算法评估：将所有样本集随机平均分为5组数据集，选取其中一组作为测试集，剩下的样本集则作为训练集；迭代五次，直到每一组数据集都被选过作为测试集为止。对于每次验证过程，先将基于双聚类的特征选择方法应用于训练集，依照特征的双聚类分值进行排序，然后分别对大小为1～73的特征子集在测试集上进行测试，通过其分类表现来验证本发明的正确性和有效性。

选取准确度(Accuracy)、敏感度(Sensibility)和特异度(Specificity)作为定量评价指标。其中，准确度表示在所有样本中被正确分类的样本占总样本的比例，敏感度表示正样本中被正确分类的样本占正样本的比例，特异度表示负样本中被正确分类的样本占负样本的比例。分别定义如下：

$\mathrm{Accuracy}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{FP}+\mathrm{TN}+\mathrm{FN}}$

$\mathrm{Sensibility}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$

$\mathrm{Specificity}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}$

其中，TP表示正样本被正确识别的数目，TN表示负样本被正确识别的数目，FN表示正样本被错误识别为负样本的数目，FP表示负样本被错误识别为正样本的数目。

图4是自动识别准确率随特征子集大小的变化图。从该图可以看出，准确率并没有随着特征数目的增加而提高，说明原始特征空间中存在一些冗余特征信息，这些特征不仅不能改善分类的效果，反而降低了自动识别的准确率；综合考虑准确度、敏感度和特异度等三个指标，从图4可以看出，当选择前25维特征(主要包括位置特征、形状特征和灰度共生矩阵)时，自动识别的准确率最高。

下表给出了选择前25维特征的条件下自动识别的准确度、敏感度和特异度，可以看出，三个指标都在97％以上，充分说明了本发明的正确性和有效性。

准确度	敏感度	特异度
			0.983±0.013	0.985±0.019	0.974±0.035

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (7)

1.一种基于超声图像的乳腺肿瘤区域自动识别方法，其特征在于，包含以下顺序的步骤：

S1.获取乳腺超声图像，并对其进行预处理；

S2.使用图像分割算法对预处理后的超声图像进行分割，得到多个分割子区域；

S3.提取超声图像的灰度直方图、纹理特征、梯度特征、形态特征，再加上二维位置信息，得到高维度的特征向量；

S4.对高维度的特征向量，通过基于双聚类的特征排序和选择方法来选取最有效的特征子集，具体包括以下步骤：

(1)提取双聚类种子点：对每一单独列内的元素做聚类，每一个聚类的结果都被看做一个双聚类的种子点；对于给定的一个R行C列的数据矩阵，应用一种凝聚分等级聚类算法到目标矩阵的每一列；

(2)启发式构造双聚类：利用一种启发式的方法对第一步中提取的双聚类种子点做扩展，根据预先设定的准则构造双聚类；

(3)消除冗余的双聚类：首先，对双聚类集合按照每个双聚类包含的列数目升序排序；然后，从包含最少列数目的双聚类开始迭代，检测它是否属于排在后面的双聚类的子集，如果是，则消除该双聚类；

(4)特征排序和选择：对数据矩阵提取出双聚类之后，需要从这些双聚类中提取出能够衡量各个特征的信息，即计算出每个特征相应的双聚类分值，分值越高，表明该特征越能表达该数据集，与其它特征的联系越紧密，然后根据所需要的特征数目选择出分值最高的特征子集，即为最有效的特征子集；

S5.选择出最有效的特征子集后，使用分类器对其进行学习分类，即可达到乳腺肿瘤区域自动识别的目的。

2.根据权利要求1所述的基于超声图像的乳腺肿瘤区域自动识别方法，其特征在于：步骤S1中，所述的预处理采用全变分模型，通过最小化全变分来达到滤波去噪的目的，全变分定义为

$\mathrm{TV}\left[u\right]=\underset{\mathrm{\Ω}}{\∫}|\▿u|\mathrm{dx}$

其中，Ω表示连续的信号域，表示图像梯度，dx则表示Ω的一个元素。

3.根据权利要求1所述的基于超声图像的乳腺肿瘤区域自动识别方法，其特征在于：步骤S2中，所述的图像分割算法阈值法、聚类法、最大后验法、马尔可夫随机场法、动态规划法、贝叶斯水平集法、基于图论的分割算法中的一种。

4.根据权利要求1所述的基于超声图像的乳腺肿瘤区域自动识别方法，其特征在于：步骤S4的分步骤(1)中，所述的一种凝聚分等级聚类算法具体如下：

a、计算所有原始数据的距离的最大值和最小值；

b、将所有原始数据每个元素初始化为单独一类；

c、建立距离数组，记录各类之间的距离，初始化为各个原始数据的距离；

d、对所有数据做以下循环：找到距离数组的最小值，如果该值小于阈值τ，则将这两个类归为一类，重新计算距离数组，直到距离数组的最小值大于或等于τ，终止循环；

凝聚分等级聚类后，将得到目标矩阵的一系列双聚类种子点：

[C_s(i,j),N_cl(j)]＝HC(j,τ),j＝1,...,C

其中，HC代表凝聚分等级聚类算法，τ是聚类算法预设的阈值，C_s(i,j)表示第j列的第i个双聚类的种子点，N_cl(j)表示第j列双聚类种子点的总数目。

5.根据权利要求1所述的基于超声图像的乳腺肿瘤区域自动识别方法，其特征在于：步骤S4的分步骤(2)中，所述的对双聚类种子点做扩展，其方法包含以下步骤：

a、按照每个种子点包含的行数做升序排序；

b、从包含最小行数R_j的种子点开始，向其他列扩展相应的行，形成一个R_j行C列的矩阵M；

c、对于矩阵M，计算其均方残基得分：

$\mathrm{MSRS}=\frac{1}{\left|R\right|\left|C\right|}\underset{i\∈R,j\∈C}{\mathrm{\Σ}}{(e_{\mathrm{ij}}-e_{\mathrm{iC}}-e_{\mathrm{Rj}}+e_{\mathrm{RC}})}^2$

$e_{\mathrm{iC}}=\frac{1}{\left|C\right|}\underset{j\∈C}{\mathrm{\Σ}}{e}_{\mathrm{ij}},e_{\mathrm{Rj}}=\frac{1}{\left|R\right|}\underset{i\∈R}{\mathrm{\Σ}}{e}_{\mathrm{ij}},e_{\mathrm{RC}}=\frac{1}{\left|R\right|\left|C\right|}\underset{i\∈R,j\∈C}{\mathrm{\Σ}}{e}_{\mathrm{ij}}$

其中，R是矩阵的行数，C是矩阵的列数，e_ij是该矩阵第i行第j列的元素值；如果计算得出的MSRS大于设定的阈值，则对于矩阵中的每一个点，计算除去该点所在的行或者列，计算新矩阵的MSRS,以新旧矩阵MSRS得分之差来衡量该点所在行或者列对MSRS的贡献值，删除MSRS贡献值最大的点所在的行或者列；不断迭代这一过程，直到矩阵的MSRS值少于预设定的阈值，则该矩阵为最终求得的一个双聚类。

6.根据权利要求1所述的基于超声图像的乳腺肿瘤区域自动识别方法，其特征在于：所述的步骤S4的分步骤(4)，具体如下：

A、对数据矩阵提取出双聚类之后，需要从这些双聚类中提取出能够衡量各个特征的信息：基于双聚类的特征排序方案考虑了特征的相关性和样本的分离性两个因素，其中特征的相关性用相关性得分表示，用于衡量特征子集中的特征之间联系的紧密程度，样本的分离性用分离性得分表示，用于衡量一个特征的独立程度；对于第k个特征，如果它出现在双聚类子集Z_k中的任一双聚类中，则该特征的相关性得分和分类性得分分别表示如下：

$\mathrm{correlation}\_\mathrm{score}=\underset{i=1}{\overset{n_{b,k}}{\mathrm{\Σ}}}\frac{n_{f,k}\left(i\right)}{n_c}$

$\mathrm{separability}\_\mathrm{score}=\frac{n_{s,k}}{n_r}\sqrt{\underset{i=1}{\overset{n_{b,k}}{\mathrm{\Σ}}}{({\mathrm{\μ}}_{i,k}-{\mathrm{\μ}}_{a,k})}^2/n_{b,k}}$

其中，n_b,k表示Z_k中双聚类的数目，n_f,k(i)表示Z_k中第i个双聚类包含的特征的数目，n_s,k表示Z_k中不重复的行的数目，μ_i,k表示Z_k中第i个双聚类的第k个特征的均值，μ_a,k表示μ_i,k的均值；

第k个特征的最终得分用双聚类分值表示：

$\mathrm{bicluster}\_\mathrm{score}=a\·\stackrel{\‾}{\mathrm{correlation}\_\mathrm{score}}+\stackrel{\‾}{\mathrm{separability}\_\mathrm{score}}$

其中，和分别表示correlation_score和separability_score归一化结果；a是一个常数，用于平衡correlation_score和separability_score对最终得分的贡献；

特征的bicluster_score分值越高，表明该特征越能表达该数据集，与其它特征的联系越紧密；

B、根据所需要的特征数目选择出分值最高的特征子集。

7.根据权利要求1所述的基于超声图像的乳腺肿瘤区域自动识别方法，其特征在于：步骤S5中，所述的分类器K最近邻分类器、神经网络分类器、支持向量机分类器中的一种。