CN103593844B - 一种医学图像的多维多种特征提取方法 - Google Patents

Abstract

本发明提供一种颈动脉图像的多维多种特征提取方法，具体为：（1）依据颈动脉的医学三维超声体数据构建感兴趣区域二维序列图像和感兴趣区域三维图像；（2）对二维序列图像提取二维特征，对三维图像提取三维特征，所述二维特征包括二维纹理特征、二维形态特征和二维弹性特征，所述三维特征包括三维纹理特征；（3）验证各种特征是否对图像分类具有贡献，保留具有贡献的特征；以假发现率最小为目标，在具有贡献的特征中搜索得到全局最优特征组合。本发明针对图像进行了多个维度多类特征的提取，囊括了纹理、形态、弹性三大类特征，特征数目多，类型全，能够全面、客观、准确地反映图像特性，提高了分类的精确性，为临床应用提供重要的参考依据。

Description

一种医学图像的多维多种特征提取方法

技术领域

本发明属于医学图像信息处理领域，具体涉及一种医学图像的多维多种特征提取方法。

背景技术

图像特征提取正在成为计算机辅助诊断领域研究的热点之一，具有重要的意义。譬如，在颈动脉粥样硬化上，特征提取有助于易损斑块的自动识别，对医生的临床诊断具有有效的辅助作用。动脉粥样硬化斑块破裂引发血栓，极有可能导致急性心血管病事件的发生。然而，并非所有斑块都会破裂、形成血栓；斑块是否会破裂，取决于其易损性。易损斑块，即指容易破损、形成血栓、可能快速发展从而引发冠状动脉阻塞、死亡的斑块。因此，识别易损斑块是具有重要意义的前瞻性研究。提取易损斑块识别的关键技术在于显著性特征的提取，特征提取的好坏对斑块识别的准确性有很大的影响。目前，现有技术主要是基于形态、纹理、弹性上的特征提取，下面将详细介绍。

在纹理特征的提取上，纹理分析技术一直是计算机视觉、图像处理、图像分析、图像检索等的活跃研究领域。20世纪70年代以前出现了自相关函数法，功率谱方法和一些与灰度频率相关的方法。这些方法取得了一定的成功，但是没有具体的定义、描述或纹理模型，仅仅是某种数学变换。另外还有一些提取纹理特征的方法，也仅限于提取特定的图像属性，如纹理粗糙度，纹理直线性等。1973年，Haralick在利用陆地卫星图像研究美国加利福尼亚海岸带的土地利用问题时，开创性地提出了著名的灰度共生矩阵（Gray Level Co-occurrence Matrix，GLCM），它在纹理分析中是一个很好的方法，被广泛地应用于将灰度值转化为纹理信息。此外，这个阶段出现的方法主要还有灰度行程长度法（gray level runlength method，RLM）、灰度差分统计法（gray differential statistics），自回归模型法（Autoregressive model method，ARM）等，这些方法在纹理分类中有一定效果，但是这些方法的后续研究很少，在实际应用中也较少采用。在形态特征的提取上，以往的研究主要从斑块形态学的角度入手，发现易损斑块具有较大的外弹力膜面积、斑块面积、斑块负荷、偏心指数等。因此，提取这些形态学特征可用于识别易损斑块。但是，只提取斑块的形态学特征，并不能全面地衡量斑块的属性；而且在提取过程中，还需要较多人工介入，比如提取偏心指数时，医生需通过肉眼观察、手工描点测量斑块的最大厚度与最小厚度。这些缺点都是计算机自动识别易损斑块时需要着重克服的。

在弹性特征的提取上，血管搏动时，收缩压与舒张压的差异造成斑块的形变；不同斑块的弹性状况不同，因此形变也不同。可以从收缩期与舒张期的两幅超声（US）图像中检测斑块的形变，从而表征斑块的弹性特点。首先从US序列图像中自动抽取一对舒张期和收缩期图像，接着借助非刚性图像配准估计斑块的应变。最后从应变分布中提取表征斑块弹性状况的量化特征。

综上所述，以往的研究多从单个维度出发或只提取某一类型的图像特征，并不能全面客观地反映图像的特性。因此分类精确度较低，分类效果差，实用性不强。

发明内容

针对现有技术只能从单个维度出发或只提取某一类型的图像特征的缺陷，本发明提供一种图像的多维多种特征提取方法，全面、客观、准确地反映图像特性，提高了分类的精确性，为临床应用提供重要的参考依据。

一种颈动脉图像的多维多种特征提取方法，具体为：

（1）提取感兴趣区域：依据颈动脉的医学三维超声体数据构建感兴趣区域二维序列图像和感兴趣区域三维图像；

（2）提取图像特征：对感兴趣区域二维序列图像提取二维特征，以及对感兴趣区域三维图像提取三维特征，所述二维特征包括二维纹理特征、二维形态特征和二维弹性特征，所述三维特征包括三维纹理特征；

（3）最优特征选择：对各种图像特征进行数值归一化处理，处理后验证各种特征是否对图像分类具有贡献，保留具有贡献的特征；以假发现率最小为目标，在具有贡献的特征中搜索得到全局最优特征组合。

进一步地，

所述二维纹理特征包括灰度统计特征、频域统计特征和Laws能量特征；

所述二维形态特征包括周长面积类、径向特征类和椭圆拟合类；

所述二维弹性特征包括斑块横截面的顺从性和膨胀性；

所述三维纹理特征包括灰度共生矩阵中的角二阶矩、对比度和熵。

进一步地，

所述二维纹理特征的灰度统计特征包括一阶统计特征、灰度游程统计特征和灰度共生矩阵；

所述二维纹理特征的频域统计特征包括环状采样、楔状采样及两者的交叉采样；

所述二维纹理特征的Laws能量特征包括由灰度、边缘、点、波、涟漪和振荡滤波模板卷积产生的特征。

进一步地，

所述二维形态特征的周长面积类包括斑块面积、斑块负荷、内腔等效直径、血管等效直径、血管紧密度和内腔紧密度；

所述二维形态特征的径向特征类包括斑块的偏心指数和近似偏心指数；

所述二维形态特征的椭圆拟合类包括内腔长轴长、外膜长轴长、内腔短轴长、外膜短轴长、内腔偏心率和外膜偏心率。

进一步地，

所述血管紧密度表示为PE_E表示血管周长,EEMA表示外弹力膜面积；所述内腔紧密度表示为PE_L表示内膜周长，LA表示内腔面积。

进一步地，

所述步骤（3）验证各种特征是否对图像分类具有贡献的具体过程为：

提出假设H_o：特征F对图像分类具有贡献，

备择假设H_A：特征F不具有贡献；

根据统计学假设检验适用范围，选取H_o的概率标准即置信度α；

计算统计数据的u值：为分类结果为无班块的所有图像的特征F均值，为分类结果为有班块的所有图像的特征F均值，s₁ ²和s₁ ²分别为和对应的方差，n₁和n₂分别为分类结果为无班块和有班块的图像数目；

在正态分布Ν（0,1）的表中查询对应于1-α概率的接受区间D，如果u∈D，则特征F具有贡献，否则特征F不具有贡献。

进一步地，

所述步骤（3）去除冗余特征的具体过程为：以假发现率最小为目标，首先利用浮动搜索技术初始化候选特征组合，再利用遗传算法在候选特征组合中筛选以得到全局最优的特征组合。

进一步地，

所述步骤（1）具体为：

步骤（1.1）：读取三维超声体数据，依据三维超声体数据中的主颈动脉的血管走向，等间距切分三维超声体数据，得到二维横断面序列图像；

步骤（1.2）：在二维横断面序列图像的每一张图像中，分割得到各主颈动脉血管的内、外轮廓，将每一张二维血管的内、外轮廓之间的区域作为感兴趣区域二维序列图像；

步骤（1.3）：依据二维序列ROI图像及其二维序列的位置信息，按其空间位置关系三维重建，得到三维主颈动脉血管内、外轮廓，将三维血管的内、外轮廓之间的区域作为感兴趣区域三维ROI图像。

本发明与现有技术相比具有以下优点：

本发明针对劲动脉图像进行了多个维度多类特征的提取，囊括了纹理、形态、弹性三大类特征，特征数目多，类型全，能够更好地反映图像的特性。

进一步地，在二维形态特征中，本发明定义了血管紧密度和内腔紧密度，两特征反映了斑块的弹性，脂肪含量越高，弹性越强，斑块易损，两特征的提取能较好的反应图像特性。

进一步地，本发明使用了假设检验及浮动搜索的特征组合算法，对得到的多个特征组合进行组合，以便初步得到较优的特征组合集合，该较优的特征组合集作为后续的特征优化及分类的输入，保证了后续处理步骤的准确性，缩减了后续处理步骤的时间，以获得最优的特征组合。本发明使用了自组织、自适应、自学习性能都较强的遗传算法，只需根据应用，设定目标函数和适应度函数，该算法简单，运行时间较短，最终可以更少的特征，就能完整的描述病灶区域图像。

进一步地，本发明由水平集方法分割出颈动脉内外膜的轮廓，将内外膜轮廓的区域作为感兴趣区域ROI，以该非矩形ROI区域为特征提取原材料。该方法比现有研究中使用的手动截取的感兴趣矩形区域更为精确，提取的特征更能反映图像的特性。

附图说明

图1为本发明多维多类图像特征提取方法的流程图；

图2为内腔椭圆拟合绘图；

图3为外弹力膜椭圆拟合绘图；

图4为实施例的结果示意图；其中，图4（a）为原始三维颈动脉超声图像；图4（b）为图4（a）切分后的序列图像示意图（三张）；图4（c）为图4（b）上获得的颈动脉内外膜的轮廓序列图像示意图（三张）；图4（d）为图4（c）上选定ROI序列图像示意图（三张），即为自动勾勒内外轮廓的ROI序列图像；图4（e）为图4（d）重建的三维图像在不同传递函数下的显示效果示意图（三种传递函数）。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

为了行文简洁，约定下列名词简称：

P（i,j）：图像在位置(i,j)处的灰度值，这里i表示行，j表示列。

W：图像的宽度，

r_i：第i个灰度级，

N_G：可区分的灰度级数目，

N_R：游程长个数，

T_R：影像像素总数，

f(r_i)：图像灰度在各个灰度级上出现的频率，

M：灰度均值。

本发明特征提取方法步骤如图1所示，具体如下：

步骤（1）提取感兴趣区域：依据颈动脉的医学三维超声体数据构建感兴趣区域二维序列图像和感兴趣区域三维图像；

步骤（1.1）：读取三维超声体数据，依据三维超声体数据中的主颈动脉的血管走向，等间距切分三维超声体数据，得到二维横断面序列图像；

步骤（1.2）：在二维横断面序列图像的每一张图像中，分割得到各主颈动脉血管的内、外轮廓，将每一张二维血管的内、外轮廓之间的区域作为二维序列ROI图像；

步骤（1.3）：依据二维序列ROI图像及其二维序列的位置信息，按其空间位置关系三维重建，得到三维主颈动脉血管内、外轮廓，将三维血管的内、外轮廓之间的区域作为三维ROI图像；

进一步地，步骤（1.2）中，在二维横断面序列图像的每一张图像中，可采用全自动分割方法【如，参考：丁明跃,李鹤,杨鑫,等.超声颈动脉血管膜的自动分割方法[P].中国专利:申请号2012102161916】、半自动分割方法【如，参考：丁明跃，程洁玉，杨鑫，等.超声颈动脉斑块的自动分割方法[P].中国专利:申请号2012102163165.】以及手动分割方法【如，参考：丁明跃,杨鑫,金娇英,等.基于颈部超声图像的主颈动脉血管提取和厚度测量方法[P].中国专利:申请号2012102182221.；】。

下面以全自动分割方法为例详细说明如何分割得到各主颈动脉血管的内、外轮廓。在减少人工干预的前提下，利用连续帧颈动脉图像上灰度分布近似的特点，采用“分割—跟踪—分割”的策略，以当前帧图像上的血管壁分割结果自动指导下一帧图像上血管壁的分割过程，不断重复直到整个三维数据分割完成。具体为：

（1.2.1）实现对普通B超图片上目标区域的灰度分布估计。采用混合Gamma概率密度函数对感兴趣区域内的灰度直方图进行拟合，用拟合得到的混合Gamma概率密度函数参数来代替真实的灰度直方图用于指导后续的目标轮廓分割；

（1.2.2）利用新设计的水平集能量式，实现对每一帧图片上的血管外膜（MAB）与血管内腔（LIB）的分割。先利用人工置入的锚点或上一帧跟踪得到的轮廓分割MAB，再利用灰度分布特征和MAB与LIB之间位置关系的先验知识完成对LIB的分割；

（1.2.3）提出了一种连续帧上MAB轮廓的跟踪方法。根据连续帧图像之间灰度分布的相似性，利用当前帧图像上分割得到的MAB，预测下一帧图像上MAB轮廓的大致位置和形状，以实现连续二维图像上目标轮廓的自动分割；

（1.2.4）通过三维重建算法，将每一帧切片上分割得到的目标重建成完整的三维目标体，并计算出相应的临床特征值。最后，进行医生手动分割与之对比的研究。

步骤（2）提取图像特征：对感兴趣区域二维序列图像提取二维特征，以及对感兴趣区域三维图像提取三维特征，所述二维特征包括二维纹理特征、二维形态特征和二维弹性特征，所述三维特征包括三维纹理特征；

2.1二维纹理特征

二维纹理特征包括灰度统计特征、频域统计特征和Laws能量特征。其中，灰度统计特征包括：一阶统计、灰度游程统计和灰度共生矩阵等；频域统计特征包括：环状采样、楔状采样及两者的交叉采样；Laws能量特征包括：由灰度、边缘、点、波、涟漪、振荡滤波模板卷积产生的特征。

2.1.1灰度统计

对于图像的统计分析方法主要是基于图像像素的灰度值的分布与相互关系，找出反映这些关系的特征。基本原理是选择不同的统计量对纹理图像的统计特征进行提取。该方法主要适合医学图像中那些没有明显规则性的结构图像，特别适合随机的、非均匀性的结构。统计分析法中有灰度直方图、共生矩阵法、游程长度统计法和灰度差分统计法等。为了更详细地说明上述特征的提取方式，下面给出上述特征的通用计算方法。

1.一阶统计（FOS）

n阶矩：

$u_n=\underset{i=0}{\overset{N_G-1}{\mathrm{\Σ}}}{(r_i-m)}^{n}f\left(r_i\right)$

2.灰度游程统计（RL）

采用0°，45°，90°，135°四个方向，在每个方向上计算游程灰度和游程长度的联合概率密度P（i，j，θ）。

短游程优势：

${\mathrm{RF}}_1\left(\mathrm{\θ}\right)=\frac{1}{T_R}\underset{i=0}{\overset{N_G-1}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_R}{\mathrm{\Σ}}}\frac{P(i,j/\mathrm{\θ})}{j^2}---(2-1-1)$

长游程优势：

${\mathrm{RF}}_2\left(\mathrm{\θ}\right)=\frac{1}{T_R}\underset{i=0}{\overset{N_G-1}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_R}{\mathrm{\Σ}}}{j}^{2}p(i,j/\mathrm{\θ})---(2-1-2)$

灰度级非均匀性：

${\mathrm{RF}}_3\left(\mathrm{\θ}\right)=\frac{1}{T_R}\underset{j=0}{\overset{N_R}{\mathrm{\Σ}}}{\left[\underset{i=1}{\overset{N_G-1}{\mathrm{\Σ}}}p(i,j/\mathrm{\θ})\right]}^2---(2-1-3)$

游程长非均匀性：

${\mathrm{RF}}_4\left(\mathrm{\θ}\right)=\frac{1}{T_R}\underset{i=0}{\overset{N_G-1}{\mathrm{\Σ}}}{\left[\underset{j=1}{\overset{N_R}{\mathrm{\Σ}}}p(i,j/\mathrm{\θ})\right]}^2---(2-1-4)$

游程百分数：

${\mathrm{RF}}_5\left(\mathrm{\θ}\right)=\frac{1}{T_R}\underset{i=0}{\overset{N_G-1}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_R}{\mathrm{\Σ}}}p(i,j/\mathrm{\θ})---(2-1-5)$

其中：

$T_R=\underset{i=0}{\overset{N_G-1}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_R}{\mathrm{\Σ}}}p(i,j/\mathrm{\θ})---(2-1-6)$

3.灰度共生矩阵（GLCM）：

这些统计函数为:

（1）能量(Energy,or Uniformity,or Angular Second Moment)

$\mathrm{ASM}=\underset{i,j}{\mathrm{\Σ}}p{(i,j)}^2---(2-1-7)$

（2）对比度(Contrast,or Momentum)

$\mathrm{CON}=\underset{i,j}{\mathrm{\Σ}}{(i-j)}^{2}p(i,j)---(2-1-8)$

（3）相关性(Correlation)

$\mathrm{COR}=\frac{\underset{i,j}{\mathrm{\Σ}}\mathrm{ijp}(i,j)-{\mathrm{\μ}}_x{\mathrm{\μ}}_y}{{\mathrm{\σ}}_x{\mathrm{\σ}}_y}---(2-1-9)$

其中μ_x,μ_y,σ_x,σ_y分别是 $p_x\left(i\right)=\underset{x}{\mathrm{\Σ}}{p}_x(i,k)$ 和 $p_x\left(j\right)=\underset{k}{\mathrm{\Σ}}{p}_y(k,j)$ 的均值和方差.

（4）方差(Variance,or Sum of squares)

$\mathrm{VAR}=\underset{i,j}{\mathrm{\Σ}}{(i-\mathrm{\μ})}^{2}p(i,j)---(2-1-10)$

其中u是p（i,j）的均值

(5)逆差矩(Inverse Difference Moment)

$\mathrm{IDM}=\underset{i,j}{\mathrm{\Σ}}\frac{1}{1+{(i-j)}^2}p(i,j)---(2-1-11)$

(6)和平均(Sum Average)

$\mathrm{SA}=\underset{i=2}{\overset{2\mathrm{\ω}}{\mathrm{\Σ}}}{\mathrm{ip}}_{x+y}\left(i\right)---(2-1-12)$

其中 $p_{x+y}\left(i\right)=\stackrel{j+k=i}{\underset{j,k}{\mathrm{\Σ}}}p(j+k)$

(7)和方差(Sum Variance)

$\mathrm{SV}=\underset{i=2}{\overset{2\mathrm{\ω}}{\mathrm{\Σ}}}{(i+\underset{j=2}{\overset{\mathrm{\ω}}{\mathrm{\Σ}}}{p}_{x+y}\left(j\right)\lg \left(p_{x+y}\left(j\right)\right))}^2{p}_{x+y}\left(i\right)---(2-1-13)$

(8)和熵(Sum Entropy)

$\mathrm{SE}=-\underset{i=2}{\mathrm{\Σ}}\mathrm{\ω}{p}_{x+y}\left(i\right)\lg \left(p_{x+y}\left(i\right)\right)---(2-1-14)$

(9)熵(Entropy)

$\mathrm{ENT}=-\underset{i,j}{\mathrm{\Σ}}p(i,j)\lg \left(p(i,j)\right)---(2-1-15)$

(10)差方差(Difference Variance)

$\mathrm{DV}=\underset{i=2}{\overset{2\mathrm{\ω}}{\mathrm{\Σ}}}{(i+\underset{j=2}{\overset{\mathrm{\ω}}{\mathrm{\Σ}}}{p}_{x-y}\left(j\right)\lg \left(p_{x-y}\left(j\right)\right))}^2{p}_{x-y}\left(i\right)---(2-1-16)$

(11)差熵(Difference Entropy)

$\mathrm{DE}=-\underset{i=0}{\overset{\mathrm{\ω}-1}{\mathrm{\Σ}}}{p}_{x-y}\left(i\right)\lg \left(p_{x-y}\left(i\right)\right)---(2-1-17)$

(12)相关性信息度量(Information Measure of Correlation)

$\mathrm{IMC}=\frac{\mathrm{HXY}-\mathrm{HXY}1}{\max \{\mathrm{HX},\mathrm{HY}\}}---(2-1-18)$

其中HX和HY是p_x和p_y的熵

$\mathrm{HXY}=-\underset{i,j}{\mathrm{\Σ}}p(i,j)\lg \left(p(i,j)\right)---(2-1-19)$

$\mathrm{HXY}1=-\underset{i,j}{\mathrm{\Σ}}p(i,j)\lg \left(p_x\left(i\right)p_y\left(j\right)\right)---(2-1-20)$

$\mathrm{HXY}2=-\underset{i,j}{\mathrm{\Σ}}{p}_x\left(i\right)p_y\left(j\right)\lg \left(p_x\left(i\right)p_y\left(j\right)\right)---(2-1-21)$

(13)另一个相关性信息度量(Another Information Measure of Correlation)

$\mathrm{AIMC}=\sqrt{1-e^{2(\mathrm{HXY}2-\mathrm{HXY})}}---(2-1-22)$

(14)最大相关性系数(Maximal Correlation Coefficient)

$\mathrm{MCC}=\sqrt{\mathrm{\λ}\left(Q\right)}---(2-1-23)$

其中 $Q(i,j)=\underset{k}{\mathrm{\Σ}}\frac{p(i,k)p(j,k)}{p_x\left(i\right)p_y\left(k\right)}---(2-1-24)$

λ(Q)指Q的次特征值。

2.1.2基于傅里叶频谱分析法的纹理特征提取

提取纹理特征的功率谱方法主要包括环状采样、楔状采样及两者的交叉采样。通过极坐标下功率谱的环形区域积分的计算可以计算纹理的粗糙度，环形区域如式（2-1-25）所示；通过极坐标下楔型区域的积分可以计算纹理的方向性，楔形区域如式（2-1-26）所示。另外通过直角坐标下频率条内的积分可得到反应某频率特性的纹理特征。

$R_{r1,r2}=\left\{(u,v)\right|{r_1}^2\≤u_2+v_2\≤{r_2}^2,0\≤u,v\≤N-1\}---(2-1-25)$

$W_{{\mathrm{\θ}}_1,{\mathrm{\θ}}_2}=\left\{(u,v)\right|{{\mathrm{\θ}}_1}^2\≤\arctan (v/u)\≤{{\mathrm{\θ}}_2}^2,0\≤u,v\≤N-1\}---(2-1-26)$

2.1.3基于Laws能量度量的特征提取

Laws纹理测量的基本思想是设置两个窗口：

（1）一个是微窗口，可为3×3、5×5或7×7，常取5×5用来测量以像元为中心的小区域的灰度的不规则性，以形成属性，称为微窗口滤波；

（2）另一个为宏窗口，为15×15或32×32，用来在更大的窗口上求属性的一阶统计量（常为均值和标准偏差），称之为能量变换。

整个纹理分析过程为：

对f(x,y)进行微窗口滤波，得到F(x,y)，随即进行能量转换得到E(x,y)，再进行分量旋转得到C(x,y)，最终由分类得到特征量M(x,y).。

首先定义了一维滤波模板，然后通过卷积形成系列一维、二维滤波模板，用于检测和度量纹理的结构信息。

选定的一维滤波模板是：

其中，L代表灰度（Level），E代表边缘（Edge），S代表点（Spot），W代表波（Wave），R代表涟漪（Ripple），O代表振荡（Oscillation）。

L3=[121]

E3=[-101]

S3=[-12-1]

L5=[14641]

E5=[-1–2021]

S5=[-1020–1]

W5=[-120–21]

R5=[1–46–41]

L7=[1615201561

E7=[-1–4–50541]

S7=[-1-2141–2–1]

W7=[-1030–301]

R7=[1-2–14–1–21]

O7=[-16–1520–156–1]

1×3的矢量集是构成更大矢量的基础.每一个1×5的矢量可以由两个1×3矢量的卷积产生。1×7的矢量可以由1×3与1×5矢量卷积产生。垂直矢量和水平矢量可生成二维滤波模板。所得到的21个特征均由二维滤波模板产生分别为LL，EE，SS，WW，RR，OO，LE和EL的平均值，LS和SL的平均值，LW和WL的平均值，LR和RL的平均值，LO和OL的平均值，ES和SE的平均值，EW和WE的平均值，ER和RE的平均值，EO和OE的平均值，SW和WS的平均值，SR和RS的平均值，SO和OS的平均值，WR和RW的平均值，WO和OW的平均值，RO和OR的平均值。

2.2三维纹理特征

本发明中三维纹理特征主要是基于灰度共生矩阵的三维拓展。设f(x,y,z)为一幅三维数字图象，其大小为M×N×L，,则满足一定空间关系的灰度共生矩阵为P(i,j)=#{(x₁,y₁,z₁),(x₂,y₂,z₂)∈M×N×L|f(x₁,y₁,z₁)=i,f(x₂,y₂,z₂)=j}(3.2.1)其中i,j分别为空间中两点的灰度级，#(x)表示集合x中的元素个数，显然P为N_G×N_G的矩阵，若(x1,y1,z1)与(x2,y2,z2)间距离为d,两者与坐标横轴的夹角为θ，则可以得到各种间距及角度的灰度共生矩阵P（i，j，θ）。根据GLCM理论，图像中相距(Δx,Δy,Δz)的两个灰度像素同时出现的联合频率分布可以用灰度共生矩阵来表示，表示为M(Δx,Δy,Δz)(h,k)，其中位于(h,k)的元素mhk的值表示一个灰度为h而另一个灰度为k的两个相距为(Δx,Δy,Δz)的像素对出现的次数。这样就可以用灰度共生矩阵中的角二阶矩、对比度、熵等参数来描述三维图像的纹理特征。计算四个方向(θ=0°，45°,90°,135°)的灰度共生矩阵，并将各方向的灰度共生矩阵统计量取平均值作为三维的灰度共生矩阵，套用二维灰度共生矩阵纹理特征中的角二阶矩、对比度和熵公式。因此得到3个纹理特征值，分别用ASM、CON和ENT表示。

2.3二维形态特征

对二维序列ROI图像进行二维形态特征提取，包括周长面积类、径向特征类和椭圆拟合类。周长面积类包括斑块面积、斑块负荷、内腔等效直径、血管等效直径、血管紧密度和内腔紧密度等；径向特征类包括斑块的偏心指数和近似偏心指数等；椭圆拟合类包括内腔长轴长、外膜长轴长、内腔短轴长、外膜短轴长、内腔偏心率和外膜偏心率等。为了更详细地说明上述特征的提取方式，下面给出上述特征的通用计算方法。

2.3.1第一类：周长面积类

PEL即内腔周长；PEE即外弹力膜的周长：

$L=N_e+\sqrt{2}{N}_o---(2-3-1)$

其中，N_e和N_o分别是8方向链码中方向值为偶数和奇数的数目。

利用像素点数计算面积：

$A=\underset{x=1}{\overset{M}{\mathrm{\Σ}}}\underset{y=1}{\overset{N}{\mathrm{\Σ}}}f(x,y)---(2-3-2)$

斑块面积PA：PA=LA-EEMA (2-3-3)

斑块负荷PD：PD=PA/EEMA (2-3-4)

D_L内腔等效直径： $D_L=\sqrt{\mathrm{LA}}---(2-3-5)$

D_E血管等效直径： $D_E=\sqrt{\mathrm{EEMA}}---(2-3-6)$

CE血管紧密度： $\mathrm{CE}=\frac{{{\mathrm{PE}}_E}^2}{\mathrm{EEMA}}---(2-3-7)$

CL内腔紧密度： $\mathrm{CL}=\frac{{{\mathrm{PE}}_L}^2}{\mathrm{LA}}---(2-3-8)$

PE_E表示血管周长,EEMA表示外弹力膜面积，PE_L表示内膜周长，LA表示内腔面积。

特别说明的是，在二维形态特征中，本发明定义了血管紧密度和内腔紧密度，两特征反映了斑块的弹性，脂肪含量越高，弹性越强，斑块易损，两特征的提取能较好的反应图像特性。

2.3.2第二类：径向特征类

由内腔轮廓确定内腔中心。内腔中心即用求质心的方法求得，即为内腔轮廓坐标x，y的分别平均值，所确定的坐标。以内腔中心为起始点，沿着各角度的径线方向搜索两条轮廓与径线的交点，并计算两者间的距离(该角度的斑块厚度);找到最大距离与最小距离即为斑块的最大厚度WT_max与最小厚度WT_min。而WT_max所在径线反方向（相隔180度）的斑块厚度一般也非常薄，将其称为近似最小厚度，记为WT_min。从WT_max所在径线及其反方向径线上还可以计算内腔直径(LD)和血管直径(VD);与WT_max所在径线垂直的径线上则可以提取垂直厚度(pWT=pWT₁+pWT₂)、垂直内腔直径(pLD)和垂直血管直径(pVD)。

斑块的偏心指数(eccentricity index，EI):

EI=1-WT_min/WT_max (2-3-9)

近似偏心指数EI₂:

EI₂=1-WT_min2/WT_max (2-3-10)

极坐标系B(r,θ)转换：

$r=\sqrt{x^2+y^2}---(2-3-11)$

θ=arctan(y/x) (2-3-12)

2.3.3第三类：椭圆拟合类

1）最小二乘法原理：

把精度较高的观测量看作没有误差，并把这个观测量选作x，而把所有的误差只认为y的误差。其中c₁,c₂……c_m是m个要通过实验确定的参数。对于每组观测数据：{(x_i,y_i)|i=1,2,…,N}

y_i=f(x;c₁,c₂,…c_m) (2-3-13)

选取m组测量值代入式（2-3-13），便得到方程组。求m个方程的联立解即得m个参数的数值。

N<m时，参数不能确定。

N>m时，只能用曲线拟合的方法来处理。设测量中不存在系统误差，或者说已经修正，则y的观测值y_i围绕着期望值f(x;c₁,c₂,…,c_m)摆动，其分布为正态分布，则y_i的概率密度为

$p\left(y_i\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_i}\exp \{-\frac{{[y_i-<f(x_i;c_1,c_2,......,c_m)>]}^2}{2{\mathrm{\&sigma;}}_i^2}\}---(2-3-14)$

式(2-3-14)中σ_i是分布的标准误差。为简便起见，下面用C代表(c₁,c₂,…c_m)。考虑各次测量是相互独立的，故观测值(y₁,y₂,…y_N)的似然函数

$L=\frac{1}{{\left(\sqrt{2\mathrm{\&pi;}}\right)}^N{\mathrm{\&sigma;}}_1{\mathrm{\&sigma;}}_2...{\mathrm{\&sigma;}}_N}\exp \{-\frac{1}{2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{{[y_i-f(x;C)]}^2}{{\mathrm{\&sigma;}}_i^2}\}---(2-3-15)$

取似然函数L最大来估计参数C，应使

$\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{{\mathrm{\&sigma;}}_i^2}{[y_i-f(x_i;C)]}^2|=\min ---(2-3-16)$

$\frac{\&PartialD;}{{\&PartialD;c}_k}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{{\mathrm{\&sigma;}}_i^2}{[y_i-f(x_i;C)]}^2{|}_{c=\hat{c}}=0,(k=\mathrm{1,2},...,m)---(2-3-17)$

$\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{{\mathrm{\&sigma;}}_i^2}[y_i-f(x_i;C)]\frac{\&PartialD;f(x;C)}{{\&PartialD;C}_k}{|}_{c=\hat{c}}=0,(k=\mathrm{1,2},...,m)---(2-3-18)$

求得m个参数的估计值从而得到拟合的曲线方程 $f(x;{\hat{c}}_1,{\hat{c}}_2,...,{\hat{c}}_m).$

然而，对拟合的结果还应给予合理的评价。若y_i服从正态分布，可引入拟合的χ²量，

${\mathrm{\&chi;}}^2=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{{\mathrm{\&sigma;}}_i^2}{[y_i-f(x_i;C)]}^2---(2-3-19)$

把参数估计代入上式并比较式（3-25），便得到最小的χ²值

${\mathrm{\&chi;}}_{\min }^2=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{{\mathrm{\&sigma;}}_i^2}{[y_i-f(x_i;\hat{c})]}^2---(2-3-20)$

可以证明，服从自由度v＝N-m的χ²分布，由此可对拟合结果作χ²检验。

若x_min ²接近N-m（例如x_min ²≤N-m），则认为拟合结果是可接受的；

若则认为拟合结果与观测值有显著的矛盾。

2)非线性最小二乘法椭圆拟合：

目标函数：

$f(A,B,C,D,E)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{({{\mathrm{Ax}}_i}^2+{\mathrm{Bx}}_i{y}_i+{{\mathrm{Cy}}_i}^2+{\mathrm{Dx}}_i+{\mathrm{Ey}}_i+F)}^2---(2-3-21)$

极值原理：

$\frac{\&PartialD;f}{\&PartialD;B}=\frac{\&PartialD;f}{\&PartialD;C}=\frac{\&PartialD;f}{\&PartialD;D}=\frac{\&PartialD;f}{\&PartialD;E}=\frac{\&PartialD;f}{\&PartialD;F}=0---(2-3-22)$

求得方程系数A，B，C，D，E，F的值。

应用上述椭圆拟合的方式，分别对图像内轮廓，和外弹力膜进行椭圆拟合，分别得到如附图2所示结果。其中粗线表内腔轮廓，细线代表拟合的椭圆以及如附图3所示结果，其中粗线表外弹力膜轮廓，细线代表拟合的椭圆。进而可以得到椭圆拟合的相关参数MAALL内腔长轴长，MAALE外膜长轴长，MAILL内腔短轴长，MAILE外膜短轴长，ECCL内腔偏心率，ECCE外膜偏心率。

2.4弹性特征

二维弹性特征包括斑块横截面的顺从性（CC）和膨胀性（DC）两个弹性特征。

CC=πD(ΔD/2ΔP) (2-4-1)

DC=(2ΔD/D)/ΔP (2-4-2)

其中总颈动脉和右大腿总动脉的直径为D，P表示血管壁的压力。

步骤（3）最优特征选择：对各种图像特征进行数值归一化处理；验证各种特征是否对图像分类具有贡献，保留具有贡献的特征；以假发现率最小为目标，在具有贡献的特征中搜索得到全局最优特征组合。

3.1预处理

3.1.1剔除相差较远的点

为了减少噪声对特征有效性的影响，需要对提取的特征值进行剔除个别点的处理。

首先对统计同一特征的平均值。然后删除远离所该特征平均值的特征点，例如以特征值标准差的两倍为阈值，即在正态分布随机变量95%以外的点被剔除掉。

3.1.2特征值归一化

由于得到的特征值并不是均匀分布在均值周围，因此使用非线性方法：softmax比例，即下面所列公式，将特征值压缩到[0,1]之间。

$y=\frac{x_{\mathrm{ik}}-\stackrel{\&OverBar;}{x_k}}{r{\mathrm{\&sigma;}}_k},{\hat{x}}_{\mathrm{ik}}=\frac{1}{1+\exp (-y)}$

x_ik表示第i张二维图像的第k个特征，

表示第k个特征的平均值，

σ_k表示第k个特征的标准差，

表示归一化后第i张二维图像的第k个特征，

r是自己定义的常数，在本发明中r∈[0.8,1.2]。

3.2筛选具有贡献的特征

在对特征点的删除和归一化处理后，需要检验单个特征的有效性，如果该特征不具有贡献，即不同的图像特征值区分度不大，那么该特征就不列入接下来的特征选择步骤中。

对特征进行假设检验的理论基础是：根据总体的理论分布及小概率原理，对未知的总体提出两种彼此对立的假设，然后由样本的实际结果计算，作出在一定概率意义上应该接受哪种假设的推断。

具体过程为：

提出假设H_o：特征F对图像分类具有贡献；

备择假设H_A：特征F不具有贡献；

根据统计学假设检验适用范围，选取H_o的概率标准即置信度α；

计算统计数据的u值：为分类结果为无班块的所有图像的特征F均值，为分类结果为有班块的所有图像的特征F均值，s₁ ²和s₁ ²分别为和对应的方差，n₁和n₂分别为分类结果为无班块河有班块的图像数目；

在正态分布Ν（0,1）的表中查询对应于1-α概率的接受区间D，如果u∈D，则该特征具有贡献，否则该特征不具有贡献。

3.3寻找最优特征组合

本步骤的技术思路为：以假发现率最小为目标，首先利用浮动搜索技术确定初始特征组合，在初始特征组合的基础上再利用遗传算法进一步搜索得到全局最优特征组合。

3.3.1浮动搜索技术

此步骤选用了单个特征筛选与多个特征组合优化的混合方法，即首先使用假设检验对单个特征验证可分性，随即利用浮动搜索技术对特征进行组合，以得到较优组合。

该步骤的目的是：在由m个特征构成的特征集合中，寻找其中最好的k个特征构成特征子集，k=1,2，…，l≤m，使代价准则C最优。

选择特征子集的准则：将每个特征子集送入分类器中对应得到一个代价准则函数，即假发现率(False Discovery Rate,FDR)FDR=FP/(FP+TP)。（在双分类器中有4类可能的输出。如果输出的预测是p而真实的结果也是p，那么这就叫做真阳性（TP）；然而如果真实的结果是n，则这就叫做假阳性（FP）。相反的来说，一个真阴性发生在预测结果和实际结果都为n的时候，而假阴性是当预测输出是n而实际值是p的时候。）

令X_k={x₁,x₂,…x_k}是k个特征的最好组合集，Y_m-k是其余m-k个特征的集合。保留所有低维的最好子集，即分别对应于2,3，…k-1个特征向量的X_2，X_3，…X_k-1子集。方法核心是：在下一步中，通过从Y_m-k中借用一个特征来形成第k+1个最好子集X_k+1；然后，转向以前选择的低维子集，检验包含新特征后是否降低了代价准则值FDR。如果是，则使用新特征替换以前选择的特征。要使FDR值最小化，通过以下步骤实现算法：

（1）从中选择与Y_m-k最好C值对应的特征，它要与X_k组合，X_k+1={X_k，X_k+1}。

（2）检验：

1.即将它从X_k+1中剔除是，对代价的影响最小。

2.如果r=k+1，令k=k+1，并转向步骤（1）。

3.如果r≠k+1，且C(X_k+1-{x_r})<C(X_k)，则转向步骤（1）。即如果移走x_r并不能提高以前选择的k个最好组的代价，就不需要再进行后向搜索了。

4.如果k=2，令X_k=X_k+1-{x_r}，且C(X_k+1-{x_r})=C(X_k)；转到步骤（1）。

（3）排除：

1.X′k=X_k+1-{x_r}，即剔除x_r。

2.即在新组合中找到最不重要的特征。

3.如果C(X_k-{x_s})<C(X_k-1)，那么X_k=X′_k，并转向步骤（1），不再进行后向搜索。

4.令X′_k-1=X′_k-{x_s}，并且k=k-1。

5.如果k=2，令X_k=X′_k，并且C(X_k)=C(X′_k)，并且转到步骤（1）。

6.转到步骤（3）中的1.

通过执行循序前向算法建立X₂，可使算法初始化，该算法在找到l个特征以后结束。

3.3.2：对组合的特征使用优化算法对特征进行自动筛选，以得到全局最优的特征组合。

传统的优化方法和现代智能优化算法，其中传统方法包括穷举、贪婪、回朔、递归、分治，智能算法包括模拟退火、禁忌搜索、遗传算法、蚁群算法、粒子群算法等；本发明中，选用了自组织、自适应、自学习性能都较强的遗传算法。对上一步筛选出的k个特征作为初始组合，使用遗传算法在步骤2得到的特征集合的基础上进行挑选加以优化，以得到冗余最小化的最优组合。

遗传算法作为一种全局自适应随机搜索算法，通常使用二进制遗传编码，即等位基因Γ={0,1}，个体空间HL={0,1}L，且繁殖分为交叉与变异两个独立的步骤进行。其基本执行过程如下:

a)初始化算法参数。确定种群规模N、交叉概率Pc、变异概率Pm和终止进化准则，并置进化代数计数器t←0。

b)种群初始化。在上一步筛选出的k个特征基础上随机生成N个个体（这里将不同的特征组合（k=1,2，…，l≤m）编码成不同的二进制个体），作为初始种群X(0)。

c)种群进化。

(a)个体评价。计算X(t)种群中各个体的适应度。

(b)选择+交叉。运用锦标赛选择算子从X(t)中选出两个个体作为父个体，并依概率Pc进行交叉得到两个子个体。反复执行N/2次，直至得到N个子个体。这里规模为M的锦标赛选择算子是先从种群中随机选出M个个体，再通过对比M个个体的适应度值，选出其中最好的一个。通常情况下，M=2。

(c)变异。对上一步得到的N个子个体分别依概率Pm执行变异操作，得到规模大小为N的的新一代种群。

d)终止检验。如果满足终止准则，则输出X(t+1)中具有最大适应度的个体（特征组合）作为最优解，终止计算；否则置t←t+1，转c)，继续进化直至满足终止准则。

对优化的特征组合可使用分类器进行分类，评价其特征组合的有效性。验证的方法有多个分类器对比、留一法、随机法等等。上述步骤得到最优特征子集。

下面结合附图对本发明的一个具体实例进行介绍。

本发明实施例，包括以下五个步骤：

(1)载入待分解的原始超声图像图4（a）；依据三维超声体数据中的主颈动脉的血管走向，按垂直于血管走向方向，等间距切分三维超声体数据，得到二维横断面序列图像图4（b）；在二维横断面序列图像的每一张图像中，分别分割得到各主颈动脉血管的内、外轮廓图4（c），将每一张二维血管的内、外轮廓之间的区域作为二维序列ROI图像图4（d）；依据二维序列ROI图像，及其二维序列的位置信息，按其空间位置关系三维重建，得到三维主颈动脉血管内、外轮廓，将三维血管的内、外轮廓之间的区域作为三维ROI图像图4（e）；

(2)对二维序列ROI图像进行二维特征提取，对三维ROI图像进行三维特征提取。其中：二维纹理特征136个，三维纹理特征3个，二维形态特征26个，弹性特征2个。

(3)对提取的上述二维纹理、三维纹理、二维形态、二维弹性特征进行组合。选用了单个特征筛选与多个特征组合优化的混合方法，即首先使用假设检验对单个特征验证可分性，随即利用浮动搜索技术对特征组合进行组合，以得到较优组合。

(4)对组合的特征使用优化算法对特征进行自动筛选，以得到全局最优的特征组合。使用遗传算法对组合的特征进行挑选加以优化，以得到冗余最小化的最优组合。

(5)对得到的特征组合使用分类器进行留一法分类判决，已验证其特征的有效性，验证结果表明该方法对图像的特征提取、组合与优化具有良好的精确性。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种颈动脉图像的多维多种特征提取方法，具体为：

(1)提取感兴趣区域：依据颈动脉的医学三维超声体数据构建感兴趣区域二维序列图像和感兴趣区域三维图像；

(2)提取图像特征：对感兴趣区域二维序列图像提取二维特征，以及对感兴趣区域三维图像提取三维特征，所述二维特征包括二维纹理特征、二维形态特征和二维弹性特征，所述三维特征包括三维纹理特征；

(3)最优特征选择：对各种图像特征进行数值归一化处理，处理后验证各种特征是否对图像分类具有贡献，保留具有贡献的特征；以假发现率最小为目标，在具有贡献的特征中搜索得到全局最优特征组合；

所述二维纹理特征包括灰度统计特征、频域统计特征和Laws能量特征；

所述二维形态特征包括周长面积类、径向特征类和椭圆拟合类；

所述二维弹性特征包括斑块横截面的顺从性和膨胀性；

所述三维纹理特征包括灰度共生矩阵中的角二阶矩、对比度和熵；

所述二维纹理特征的灰度统计特征包括一阶统计特征、灰度游程统计特征和灰度共生矩阵；

所述二维纹理特征的频域统计特征包括环状采样、楔状采样及两者的交叉采样；

所述二维纹理特征的Laws能量特征包括由灰度、边缘、点、波、涟漪和振荡滤波模板卷积产生的特征；

所述二维形态特征的周长面积类包括斑块面积、斑块负荷、内腔等效直径、血管等效直径、血管紧密度和内腔紧密度；

所述二维形态特征的径向特征类包括斑块的偏心指数和近似偏心指数；

所述二维形态特征的椭圆拟合类包括内腔长轴长、外膜长轴长、内腔短轴长、外膜短轴长、内腔偏心率和外膜偏心率；

所述血管紧密度表示为PE_E表示血管周长,EEMA表示外弹力膜面积；所述内腔紧密度表示为PE_L表示内膜周长，LA表示内腔面积；

所述步骤(3)验证各种特征是否对图像分类具有贡献的具体过程为：

提出假设H_o：特征F对图像分类具有贡献，

备择假设H_A：特征F不具有贡献；

根据统计学假设检验适用范围，选取H_o的概率标准即置信度α；

计算统计数据的u值： 为分类结果为无班块的所有图像的特征F均值，为分类结果为有班块的所有图像的特征F均值，s₁ ²和s₂ ²分别为和对应的方差，n₁和n₂分别为分类结果为无班块和有班块的图像数目；

在正态分布Ν(0,1)的表中查询对应于1-α概率的接受区间D，如果u∈D，则特征F具有贡献，否则特征F不具有贡献；

所述步骤(3)去除冗余特征的具体过程为：以假发现率最小为目标，首先利用浮动搜索技术初始化候选特征组合，再利用遗传算法在候选特征组合中筛选以得到全局最优的特征组合。

2.根据权利要求1所述的颈动脉图像的多维多种特征提取方法，其特征在于，所述步骤(1)具体为：

步骤(1.1)：读取三维超声体数据，依据三维超声体数据中的主颈动脉的血管走向，等间距切分三维超声体数据，得到二维横断面序列图像；

步骤(1.2)：在二维横断面序列图像的每一张图像中，分割得到各主颈动脉血管的内、外轮廓，将每一张二维血管的内、外轮廓之间的区域作为感兴趣区域二维序列图像；

步骤(1.3)：依据二维序列ROI图像及其二维序列的位置信息，按其空间位置关系三维重建，得到三维主颈动脉血管内、外轮廓，将三维血管的内、外轮廓之间的区域作为感兴趣区域三维ROI图像。