CN102509113B - 一种脑瘤mib-1指数范围检测方法 - Google Patents

Abstract

本发明涉及一种通过分析脑瘤患者的磁共振图像得到脑瘤MIB-1指数范围的检测方法，包括步骤：采集脑瘤患者的磁共振图像，构造病变区域图像训练样本；提取病变区域图像训练样本的图像特征，根据提取的病变区域图像训练样本的图像特征，训练得到支持向量机模型；构造病变区域图像检测样本，提取病变区域图像检测样本的图像特征；根据病变区域图像检测样本的图像特征，用支持向量机模型检测并获取该病变区域图像检测样本的MIB-1指数所处的范围。本发明克服了现有技术中只能通过手术获得脑瘤病理组织后利用免疫组织化学检测MIB-1指数的缺陷，避免了在免疫组织化学检测中引入的检测者主观思维及标准化不足的问题。

Description

一种脑瘤MIB-1指数范围检测方法

技术领域

本发明涉及图像处理及识别技术，特别涉及一种通过分析脑瘤患者的磁共振图像得到脑瘤MIB-1指数范围的检测方法。

背景技术

脑胶质瘤是中枢神经系统最常见的肿瘤，约占原发性脑肿瘤的42％。尽管随着医疗技术的发展，脑胶质瘤的治疗方法逐步发展为以手术治疗为主，结合放疗，化疗，免疫治疗等的综合疗法，但很多患者的预后并没有明显的改善。临床通常检测phosphatase and tensin homolog deleted onchromosome 10(PTEN)、epidermal growth factor receptor(EGFR)、theO6-methylguanine-DNA methyltransferase(MGMT)、tumor protein 53(P53)、monoclonal antibody of cell proliferation associated nuclear antigen(MIB-1)等蛋白表达状况，并综合上述蛋白表达的信息来评价胶质瘤的恶性程度和患者的预后状况。其中，单克隆抗体MIB-1的表达状况即MIB-1指数的值能反映胶质瘤的增殖活性，其作为一种中间结果信息，对于评价胶质瘤的恶性程度和患者预后具有一定的指导意义。

现有技术中，临床使用最为广泛的MIB-1指数的检测方法为免疫组织化学技术，该技术需要手术获得患者的胶质瘤病理切片后才能进行检测，对患者造成较大创伤，且无法指导制定术前治疗方案。2003年于《诊断病理学杂志》的232-235页刊登的一篇《免疫组化在病理诊断中的正确应用》中表明，免疫组织化学技术在标准化和结果量化方面存在不足，检测结果容易受到检测人员的主观影响。目前尚没有基于图像处理和模式识别技术进行MIB-1指数范围检测的方法。

2008年于《科技广场》的9-11页刊登的一篇《基于支持向量机的诊断方法研究》中指出：所谓支持向量机的诊断就是根据已知的故障样本，对未知工作状态样本数据进行分类。文中还提到：支持向量机的核函数形式主要有三类：多项式形式的核函数、径向基形式的核函数和S型核函数；其模型分别对应：q阶多项式分类器、径向基函数分类器和两层的感知器神经网络；其中，采用径向基函数的支持向量仅需确定惩罚因子和核宽度两个模型参数。如文中所述，确定模型参量一般有两种有效的方法：一是交叉有效；二是网络搜索。

发明内容

本发明提供了一种能够对采集得到的脑瘤患者的磁共振图像进行处理及识别的脑瘤MIB-1指数范围检测方法；克服了现有技术中只能通过手术获得脑瘤病理组织后利用免疫组织化学检测MIB-1指数的缺陷，避免了在免疫组织化学检测中引入的检测者主观思维及标准化不足的问题。

一种脑瘤MIB-1指数范围检测方法，包括以下步骤：

(1)采集脑瘤患者的磁共振图像，构造病变区域图像训练样本；

所述的脑瘤患者的磁共振图像的体数据分辨率可以为512x512x16体素；

所述的构造病变区域图像训练样本，包括以下步骤：

a.截取磁共振图像中的病变区域图像；

b.对所述的病变区域图像进行分类，对每个类别选取一定数量的包含该类别的病变区域图像，得到病变区域图像训练样本；

对所述的病变区域图像进行分类，包括以下步骤：

i、对病变区域图像的MIB-1指数进行检测，确定MIB-1指数的分界点；

ii、根据所述的MIB-1指数的分界点，对病变区域图像进行分类；

所述的MIB-1指数的分界点的选择可以根据临床的具体需求来确定；一般情况下，MIB-1≤5％或MIB＞5％是判断胶质瘤增殖活性强弱的有效分界点，因此可以优选5％为MIB-1指数的分界点；

(2)提取病变区域图像训练样本的图像特征，根据提取的病变区域图像训练样本的图像特征，训练得到支持向量机模型；

所述的提取病变区域图像训练样本的图像特征，可以为提取病变区域图像训练样本的纹理特征，具体包括步骤：

使用基本灰度统计信息提取病变区域图像训练样本的图像特征，包括：均值、方差、偏度、能量、绝对梯度均值以及绝对梯度方差6个图像特征；

使用灰度共生矩阵提取病变区域图像训练样本的图像特征，包括：角二阶矩、对比度、相关系数、方差、逆差矩、和平均、和方差、和熵、熵、差平均、惯性、差方差、差熵13个图像特征；

使用灰度-梯度共生矩阵提取病变区域图像训练样本的图像特征，包括：灰度平均值、梯度平均值、灰度方差、梯度方差、小梯度优势、大梯度优势、灰度分布的不均匀性、梯度分布的不均匀性、能量、相关系数、灰度熵、梯度熵、混合熵、惯性、逆差矩15个图像特征；

使用游程长度矩阵提取病变区域图像训练样本的图像特征，包括：短游程优势度量、长游程优势度量、灰度分布、游程分布、游程百分比5个图像特征；

使用闵可夫斯基泛函提取病变区域图像训练样本的图像特征，包括：周长最大值、周长最小值、欧拉数最大值、欧拉数最小值、周长起始点对应阈值、周长终止点对应阈值、欧拉数起始点对应阈值、欧拉数终止点对应阈值8个图像特征；

根据提取的病变区域图像训练样本的图像特征，训练得到支持向量机模型，包括以下步骤：

a、根据提取的病变区域图像训练样本的图像特征构造病变区域图像特征样本集；

b、根据病变区域图像特征样本集确定支持向量机模型的参数；

由于病变区域图像特征样本集所含的维数较高，可能含有一些冗余的图像特征，这些冗余的图像特征一方面可能降低分类的精度，另一方面会大大增加支持向量机的计算开销；因此作为优选，可基于离散粒子群算法对病变区域图像特征样本集进行优化，得到病变区域图像优化特征集；根据该病变区域图像优化特征集确定支持向量机模型的参数，有效地降低图像特征的复杂程度；

基于离散粒子群算法对病变区域图像特征样本集进行优化，得到病变区域图像优化特征集，包括以下步骤：

i、确定离散粒子群算法的粒子群规模，每个粒子对应一个随机选取的病变区域图像特征样本集子集；

ii、将每个粒子对应的病变区域图像特征样本集子集分别放入支持向量机中用留一交叉检验法进行分类，计算每个病变区域图像特征样本集子集的分类准确率，并将该准确率作为对应粒子的适应度；基于粒子的适应度迭代更新病变区域图像特征样本集子集达一定次数，将适应度最大的粒子对应的特征样本集子集作为病变区域图像优化特征集；

(3)构造病变区域图像检测样本，提取病变区域图像检测样本的图像特征；根据病变区域图像检测样本的图像特征，用支持向量机模型识别并获取该病变区域图像检测样本的MIB-1指数所处的范围。

本发明的有益效果为：

一、本发明基于脑瘤患者的磁共振图像，能够获得MIB-1指数的中间结果信息，具有快速及时的优点；

二、通过采集和分析脑瘤患者的磁共振图像来检测MIB-1指数范围，可以在获得肿瘤组织切片前无创地获取关于脑瘤恶性程度和患者预后的中间结果信息，效率高且操作简易；

三、客观地获得脑瘤MIB-1指数所处的范围，避免了在免疫组织化学技术中引入检测者的主观影响，也避免了免疫组织化学检测方法标准化不足的问题；

四、仅需通过图像分析获得MIB-1指数的范围，不需要消耗化学试剂等，具有成本低的优点。

附图说明

图1为本发明实施例1的脑瘤MIB-1指数范围检测方法流程示意图图；

图2为实施例1的S1中构造病变区域图像训练样本的流程示意图；

图3为实施例1中由T1加权序列的磁共振图像中截取的病变区域图像的示意图；

图4为实施例1中使用的SVM分类器在线性可分情况下的最优分类面的示意图；

图5为实施例1中用闵可夫斯基泛函提取图像特征的示意图；

图6为实施例2中对病变区域图像特征样本集进行优化的流程示意图；

图7为使用未优化的病变区域图像特征样本集训练支持向量机，然后对验证样本集进行分类得到的ROC曲线。

图8为对病变区域图像特征样本集进行特征优选后，使用优选后的训练样本集训练支持向量机，然后对验证样本集进行分类得到的ROC曲线；

具体实施方式

下面结合附图详细介绍本发明的具体实施过程。

实施例1：

如图1所示的一种脑瘤MIB-1指数范围检测方法，包括以下步骤：

S1采集脑瘤患者的磁共振图像，构造病变区域图像训练样本；

在本步骤中，磁共振图像包括T1加权序列、T1增强序列、FLAIR序列中的任一种或任几种，具体的采集方法如下：

使用磁共振扫描仪(例如GE Healthcare，1.5T)采集脑胶质瘤患者的横断位、冠状位或矢状位的磁共振图像，该磁共振图像包括T1加权序列、T1增强序列和FLAIR序列。其中，T1加权序列的成像参数优选为Repetition Time＝1966.1ms，Echo Time＝21.088ms，Inversion Time＝750ms；T1增强序列的成像参数优选为Repetition Time＝1967.25ms，EchoTime＝7.264ms，Inversion Time＝750ms；FLAIR序列图像的成像参数优选为Repetition Time＝8002ms，Echo Time＝122.904ms，Inversion Time＝2000ms。优选采集脑胶质瘤患者的横断位的磁共振图像。对患者进行扫描时，每个序列的磁共振图像的体数据分辨率为512x512x16体素，即每个二维切面图像的分辨率为512x512，共有16个二维切面图像。磁共振图像的格式一般为DICOM。

在本步骤中，构造病变区域图像训练样本，如图2所示，包括以下步骤：

S101截取磁共振图像中的病变区域图像；

在T1加权序列、T1增强序列和FLAIR序列的任一种或任几种序列中截取磁共振图像中的病变区域图像，具体方法如下：

在T1加权序列、T1增强序列和FLAIR序列的任一种或任几种序列的横断位、冠状位或矢状位的磁共振图像的体数据中的每个二维切面图像上至多截取一个病变区域图像，即图3中的小方框内的区域；病变区域图像的格式优选为：尺寸为16x16像素、灰度级为256级、图像格式为tif。由于T1加权序列、T1增强序列和FLAIR序列中的部分序列上，磁共振图像的伪影比较严重，容易影响图像分析和分类的结果，因此不同序列的病变区域图像的截取个数不尽相同。

S102利用免疫组织化学技术对病变区域图像的MIB-1指数进行检测，确定MIB-1指数的分界点；MIB-1指数的分界点的选择可以根据临床的具体需求来确定；一般情况下，MIB-1≤5％或MIB＞5％是判断胶质瘤增殖活性强弱的有效分界点，因此优选5％为MIB-1指数的分界点；

S103根据所述的MIB-1指数的分界点，基于S102所得的检测结果对病变区域图像进行分类；分类结果是将病变区域图像标记为MIB-1≤5％或MIB＞5％；

S104提取MIB-1≤5％和MIB＞5％的两个类别中一定数量的包含该类别的病变区域图像构成病变区域图像训练样本；

步骤S104中的每个类别选取为训练样本的数量满足如下公式：

$\frac{1}{2}m\≤n\≤\frac{3}{4}m$

其中，m为每个类别中病变区域图像的个数；n为每个类别中选取为病变区域图像训练样本的个数。

S2提取病变区域图像训练样本的图像特征；

在本实施例中，对病变区域图像训练样本提取的图像特征为病变区域图像训练样本的纹理特征，包括对病变区域图像训练样本提取基本灰度统计信息、灰度共生矩阵、灰度-梯度共生矩阵、游程长度矩阵、以及闵可夫斯基泛函。以下对上述提取方式作进一步说明：

a、用基本灰度统计信息提取病变区域图像训练样本的图像特征：

病变区域通常为二维数字图像，以下以一幅数字图像为例进行进一步说明。此处以f(x，y)表示一幅二维数字图像，假设其大小为M×N。基本的灰度统计信息主要有均值、方差、偏度、能量、绝对梯度均值以及绝对梯度方差。上述每个特征的公式如下所示：

1)均值：

$\mathrm{\μ}=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}f(i,j)$

2)方差：

${\mathrm{\σ}}^2=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{[f(i,j)-\mathrm{\μ}]}^2$

3)偏度：

${\mathrm{\μ}}_3=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left[\frac{f(i,j)-\mathrm{\μ}}{\mathrm{\σ}}\right]}^3$

4)能量：

$E=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left[f(i,j)\right]}^2$

对于给定的邻域窗口大小w×w，绝对梯度矩阵的元素定义为：

得到大小为M×N的绝对梯度矩阵AG，通过绝对梯度矩阵计算绝对梯度均值和绝对梯度方差，公式如下：

5)绝对梯度均值：

${\mathrm{\μ}}_{\mathrm{ag}}=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{AG}(i,j)$

6)绝对梯度方差：

${\mathrm{\σ}}_{\mathrm{ag}}^2=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{[\mathrm{AG}(i,j)-{\mathrm{\μ}}_{\mathrm{ag}}]}^2$

b、用灰度共生矩阵提取病变区域图像训练样本的图像特征：

仍以一幅二维数字图像为例进行说明，f(x，y)表示一幅二维数字图像，假设其大小为M×N，最高灰度级为第N_g级。在二维数字图像中，在某个方向上相隔一定距离的一对像素点的灰度出现的统计规律，从一定程度上可以反映这个二维数字图像的图像特征。这个统计规律可以用一个矩阵描述，即灰度共生矩阵。

在二维数字图像中，任意取一像素点(x，y)以及偏离它的另一像素点(x+a，y+b)形成一个点对。设该点对的灰度值为(i，j)，即像素点(x，y)的灰度值为i，像素点(x+a，y+b)的灰度值为j。固定a和b，令点(x，y)在整幅二维数字图像上移动，则会得到各种(i，j)值。假如二维数字图像的灰度级别为N_g，则i与j的组合共有种。在整幅二维数字图像中，统计每一种组合出现的频率为P(i，j，d，θ)，则构成大小为N_g×N_g的灰度共生矩阵，其中d是点对之间的距离

θ为点对构成的向量与坐标横轴之间的夹角，即点对的方向。灰度共生矩阵本质上就是两个像素点的联合直方图，距离差分值(a，b)取不同的数值组合，都可以得到二维数字图像沿一定方向θ、相隔一定距离

的灰度共生矩阵。灰度共生矩阵的数学表达式为：

P(i，j，d，θ)＝#{(x，y)，(x+a，y+b)∈M×N|f(x，y)＝i，f(x+a，y+b)＝j}

其中#{x}表示集合x中元素的个数。P为N_g×N_g的矩阵。若(x，y)与(x+a，y+b)之间距离为d，点对构成的向量与坐标横轴之间的夹角为θ，则可以得到各种间距及角度的灰度共生矩阵P(i，j，d，θ)。

通常可以设置d＝1，分别计算0°、45°、90°、135°四个方向的图像特征，对这四个方向的特征求取均值，得到与方向无关的特征。利用灰度共生矩阵提取的13个图像特征，包括角二阶矩、对比度、相关系数、方差、逆差矩、和平均、和方差、和熵、熵、差平均、惯性、差方差、差熵。在计算上述每个特征时，设置d＝1，分别计算0°、45°、90°、135°四个方向的图像特征，对这四个方向的特征求取均值，得到与方向无关的特征，则每个特征公式如下所示：

1)角二阶矩：

$f_1=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}{\left\{P(i,j)\right\}}^2$

2)对比度：

$f_2=\underset{n=0}{\overset{N_g}{\mathrm{\Σ}}}{n}^2\left\{\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j)\right\},|i-j|=n$

3)相关系数：

$f_3=\{\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j)-{\mathrm{\μ}}_x{\mathrm{\μ}}_y\}/{\mathrm{\σ}}_x{\mathrm{\σ}}_y$

式中，μ_x，σ_x分别是{P_x(i)；i＝1，2，…，N_g}的均值和标准差，μ_y，σ_y分别是{P_y(j)；j＝1，2，…，N_g}的均值和标准差。其中， $P_x\left(i\right)=\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j),i=\mathrm{1,2},...,N_g,$ $P_y\left(j\right)=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j),j=\mathrm{1,2},...,N_g.$

4)方差：

$f_4=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}{(i-\mathrm{\μ})}^{2}P(i,j)$

式中，μ是P(i，j)的均值。

5)逆差矩：

$f_5=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}\frac{1}{1+{(i-j)}^2}P(i,j)$

6)和平均：

$f_6=\underset{i=2}{\overset{2N_g}{\mathrm{\Σ}}}{\mathrm{iP}}_{x+y}\left(i\right)$

其中，

$P_{x+y}\left(k\right)=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j),k=\mathrm{2,3},...,2N_g,|i+j|=k$

7)和方差：

$f_7=\underset{i=2}{\overset{{2N}_g}{\mathrm{\Σ}}}{(1-f_6)}^2{P}_{x+y}\left(i\right)$

其中，

$P_{x+y}\left(k\right)=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j),k=\mathrm{2,3},...,2N_g,|i+j|=k$

8)和熵：

$f_8=-\underset{i=2}{\overset{2N_g}{\mathrm{\Σ}}}{P}_{x+y}\left(i\right){\log }_2\left[P_{x+y}\left(i\right)\right]$

其中，

$P_{x+y}\left(k\right)=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j),k=\mathrm{2,3},...,2N_g,|i+j|=k$

9)熵：

$f_9=-\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j){\log }_2\left[P(i,j)\right]$

10)差平均：

$f_{10}=\underset{i=0}{\overset{N_g}{\mathrm{\Σ}}}{\mathrm{iP}}_{x-y}\left(i\right)$

其中，

$P_{x-y}\left(k\right)=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j),k=\mathrm{0,1},...,N_g-1,|i-j|=k$

11)惯性：

$f_{11}=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}{(i,j)}^{2}P(i,j)$

12)差方差：

$f_{12}=\underset{i=0}{\overset{N_g-1}{\mathrm{\Σ}}}{(i-f_{10})}^2{P}_{x-y}\left(i\right)$

其中，

$P_{x-y}\left(k\right)=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j),k=\mathrm{0,1},...,N_g-1,|i-j|=k$

13)差熵：

$f_{13}=-\underset{i=0}{\overset{N_g-1}{\mathrm{\Σ}}}{P}_{x-y}\left(i\right){\log }_2\left[P_{x-y}\left(i\right)\right]$

其中，

$P_{x-y}\left(k\right)=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}P(i,j),k=\mathrm{0,1},...,N_g-1,|i-j|=k$

c、用灰度-梯度共生矩阵提取病变区域图像训练样本的图像特征：

仍以一幅二维数字图像为例进行说明，f(x，y)表示一幅二维数字图像，假设其大小为M×N，最高灰度级为第N_g级，对二维数字图像的灰度矩阵进行正规化变换：

F(x，y)＝[f(x，y)×(N_g-1)/f_max]+1

式中，[x]表示x的整数部分，f_max为二维数字图像的最大灰度值。

二维数字图像各像素的梯度计算使用Sobel算子，设二维数字图像的梯度矩阵为g(x，y)，其中x＝1，2，…，N，y＝1，2，…，M。为了使g(x，y)梯度值分布在更大的离散间隔N_s个等级中，我们对二维数字图像的梯度矩阵进行正规化变换：

G(x，y)＝[g(x，y)×(N_s-1)/g_max]+1

式中，[x]表示x的整数部分，g_max为二维数字图像的梯度矩阵的最大值，N_s是正规化变换后梯度矩阵的最大值。

通常，将正规化变换后的二维数字图像的灰度矩阵F(x，y)简称为正规化灰度矩阵，将正规化变换后的二维数字图像的梯度矩阵G(x，y)简称为正规化梯度矩阵，将正规化灰度矩阵和正规化梯度矩阵进行结合就可以得到灰度-梯度共生矩阵：

{H(i，j)；i＝1，2，…，N_g，j＝1，2，…，N_s}

其中，(i，j)表示灰度-梯度共生矩阵的第i行第j列元素，H(i，j)表示正规化灰度矩阵中灰度值为i，并且正规化梯度矩阵中梯度值为j的像素点的个数。N_s是正规化梯度矩阵的最大值，N_g是正规化灰度矩阵的最大值。

将灰度-梯度共生矩阵H(i，j)进行正规化变换，得到：

p(i，j)＝H(i，j)/(N_g×N_s)，i＝1，2，…，N_g，j＝1，2，…，N_s

利用灰度-梯度共生矩阵提取的图像特征主要有15个，分别是灰度平均值、梯度平均值、灰度方差、梯度方差、小梯度优势、大梯度优势、灰度分布的不均匀性、梯度分布的不均匀性、能量、相关系数、灰度熵、梯度熵、混合熵、惯性、逆差矩。上述每个特征的公式如下所示：

1)灰度平均值：

$T_1={\mathrm{\μ}}_G=\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}i\left[\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}p(i,j)\right]$

2)梯度平均值：

$T_2={\mathrm{\μ}}_F=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\left[\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}p(i,j)\right]$

3)灰度方差：

$T_3={\mathrm{\σ}}_G={\left\{\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_G)}^2\right[\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}p(i,j)\left]\right\}}^{1/2}$

4)梯度方差：

$T_4={\mathrm{\σ}}_F={\left\{\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_F)}^2\right[\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}p(i,j)\left]\right\}}^{1/2}$

5)小梯度优势：

$T_5=\left[\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}\frac{H(i,j)}{j^2}\right]/H$

其中，

$H=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}H(i,j)$

6)大梯度优势：

$T_6=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}H(i,j)$

7)灰度分布的不均匀性：

$T_7=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}{\left[\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}H(i,j)\right]}^2/H$

8)梯度分布的不均匀性：

$T_8=\underset{j=1}{\overset{N_g}{\mathrm{\Σ}}}\left[\underset{i=1}{\overset{N_s}{\mathrm{\Σ}}}H(i,j)\right]/H$

9)能量：

$T_9=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}{\left[p(i,j)\right]}^2$

10)相关系数：

$T_{10}=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}(i-{\mathrm{\σ}}_F)(j-{\mathrm{\σ}}_G)p(i,j)$

11)灰度熵：

$T_{11}=-\left\{\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\right[\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}p(i,j)\left]{\log }_2\right[\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}p(i,j)\left]\right\}$

12)梯度熵：

$T_{12}=-\left\{\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}\right[\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}p(i,j)\left]{\log }_2\right[\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}p(i,j)\left]\right\}$

13)混合熵：

$T_{13}=-\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}p(i,j){\log }_{2}p(i,j)$

14)惯性：

$T_{14}=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}{(i,j)}^{2}p(i,j)$

15)逆差矩：

$T_{15}=\underset{i=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_s}{\mathrm{\Σ}}}\frac{1}{1+{(i-j)}^2}p(i,j)$

d、用游程长度矩阵提取病变区域图像训练样本的图像特征：

仍以一幅二维数字图像为例进行说明，f(x，y)表示一幅二维数字图像，假设其大小为M×N。在数字图像中，称沿一确定方向θ并具有相同灰度的一批邻接像素集合为一游程，集合中点的个数为游程长度。

对于某一方向θ，其灰度游程长度矩阵R(θ)的元素定义为：

r(i，j |θ)＝#{(n₁，m₁)，(n₂，m₂)，...(n_j，m_j)|f(n_k，m_k)＝i}

其中，(n₁，m₁)，(n₂，m₂)，...(n_j，m_j)是指沿方向θ，灰度值为i，且长度为j的任一游程，显然r(i，j|θ)反映了图像中灰度值为i，沿θ方向，游程为j的像素点个数。分别计算0°、45°、90°、135°四个方向的灰度游程长度矩阵，进而统计：这样就构成了灰度游程长度矩阵R，其大小为m×n，然后再从这个灰度游程长度矩阵R中提取一些纹理特征进行分析。提取的特征如下所示：

1)短游程优势度量：

$L_1=\left[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{R(i,j)}{j^2}\right]/\left[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}R(i,j)\right]$

2)长游程优势度量：

$L_2=\left[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{j}^{2}R(i,j)\right]/\left[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}R(i,j)\right]$

3)灰度分布：

$L_3=\left\{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\left[\underset{j=1}{\overset{n}{\mathrm{\Σ}}}R(i,j)\right]}^2\right\}/\left[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}R(i,j)\right]$

4)游程分布：

$L_4=\left\{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}R(i,j)\right]}^2\right\}/\left[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}R(i,j)\right]$

5)游程百分比：

$L_5=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{R(i,j)}{M\×N}$

e、用闵可夫斯基泛函提取病变区域图像训练样本的图像特征：

仍以一幅二维数字图像为例进行说明，f(x，y)表示一幅二维数字图像，假设其大小为M×N，灰度级范围为[N_s，N_g]。

闵可夫斯基泛函有3个泛函公式，分别表示图像中目标点的面积(area)，周长(perimeter)以及欧拉数(Euler number)即连通性的信息。闵可夫斯基泛函是针对黑白二值图像进行计算的，因此对于二维灰度图像，首先应采用阈值变换法将图像变换成二值图像。假设q是N_s与N_g之间的数，以q为阈值变换得到的二值图像f_binary(x，y|q)定义为：

$f_{\mathrm{binary}}(x,y|q)=\left\{\begin{array}{cc}1,& f(x,y)\&GreaterEqual;q\\ 0,& f(x,y)<q\end{array}\right.$

其中，f_binary(x，y|q)＝0的点为背景点，f_binary(x，y|q)＝1的点为目标点。

通常在二维数字图像的灰度级范围[N_s，N_g]内取100个分位点的值作为阈值，通过阈值变换可以得到相应的100幅二值图像。

对于每一幅二值图像，闵可夫斯基泛函的公式定义如下：

mf.area＝n_s，

mf.perim＝-4n_s+2n_e，

mf.Euler＝n_s-n_e+n_v，

假设二值图像中每个像素点为有1个面，4条边，4个顶点的正方形，则公式中n_s为目标点的数目，n_e是目标点的边的条数，n_v是目标点的顶点数。例如图5中所示的二值图像中，n_s＝3，n_e＝10，n_v＝8。由上述计算可知，任意一个灰度阈值t对应一幅二值图像的mf.area、mf.perim、mf.Euler共3个闵可夫斯基泛函值。

因此对一幅二维数字图像进行闵可夫斯基泛函运算可以得到3个特征函数：MF1(t，mf.area)，MF2(t，mf.perim)，MF3(t，mf.Euler)，分别提取其中MF2，MF3函数曲线的最大值，最小值，曲线起始点对应的阈值，曲线终止点对应的阈值共8维特征。

基于上述方式提取的病变区域图像的图像特征归类如表1所示，并依次标明了图像的特征序号。

表1.病变区域图像的图像特征

S3根据提取的病变区域图像训练样本的图像特征，训练得到支持向量机模型；

通过S2中的a、b、c、d和e的五种图像特征提取方式，每个病变区域图像训练样本上可以提取图像特征有47个，其中包括使用基本灰度统计信息提取的6个、使用灰度共生矩阵提取的13个、使用灰度-梯度共生矩阵提取的15个、使用游程长度矩阵提取的5个、使用闵可夫斯基泛函提取的8个。对每一个病变区域图像训练样本提取47个图像特征，并且对其进行依据分类MIB-1≤5％或MIB＞5％对其进行标记，具体为：如果病变区域图像训练样本标记为MIB-1≤5％，则将该病变区域图像训练样本的类别记为-1；如果病变区域图像训练样本标记为MIB＞5％，则将该病变区域图像训练样本记为+1。通过上述处理，所有病变区域图像训练样本的图像特征和类别标记结合，构成病变区域图像特征样本集。病变区域图像特征样本集可以表示为S＝{(f_i，l_i)|i＝1，2，…，n}，其中f_i是第i个病变区域图像训练样本的47个的图像特征，表示为f_i＝[f_i1，f_i2，…，f_ij…，f_i47]，f_ij是第i个病变区域图像训练样本的第j个图像特征的特征值；l_i是第i个病变区域图像的类别标记；如果第i个病变区域图像标记为MIB-1≤5％，则l_i＝-1；如果第i个病变区域图像标记为MIB＞5％，则l_i＝+1；n是病变区域图像特征样本集中的训练样本的个数。

支持向量机(Support Vector Machines，SVM)作为样本的二分类器，是统计学习理论的重要成果，它从线性可分的最优分类面发展而来，目标是为了产生一个能够对不同类的例子进行有效分类的分类器。

病变区域图像特征样本集作为一给定的训练集，表示为S＝{(f_i，l_i)|i＝1，2，…，n}。假设病变区域图像特征样本集可被一个超平面线性划分，该超平面记为(w·x)+b＝0。如果训练集中的所有向量均能被某超平面正确划分，并且距离超平面最近的异类向量之间的距离最大(即边缘最大化)，则该超平面为最优超平面。两维情况可如图4所示。

图4中，实心点和空心点代表两不同类别的样本，H为分类面，H₁和H₂分别为过各类中离分类线最近的样本且平行于分类线的直线，它们之间的距离叫做分类间隔。所谓最优分类面就是要求分类面不但能够将两类正确分开(训练错误率为0)，而且使分类间隔最大。其中在H₁和H₂这两条直线上的向量被称为支持向量。一组支持向量可以唯一地确定一个超平面。

对于线性可分问题，不失一般性，其分类方程满足：

y_i[(w·x_i)+b]-1≥0，i＝1，...，n

此时分类间隔等于2/||w ||，使分类间隔最大等价于使||w ||²最小。满足分类方程和||w||²最小的分类面就是最优分类面，所以构造最优分类面的问题转化为在分类方程的约束条件下求||w ||²最小。使用Lagrange乘子方法可以归结为一个二次规划问题，即在约束条件

i＝1，...n下面对ai求解下列函数的最大值：

$Q\left({\mathrm{\&alpha;}}_i\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i-\frac{1}{2}\underset{i,j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\mathrm{\&alpha;}}_j{y}_i{y}_j(x_i\&CenterDot;x_j)$

a_i为原问题中与每个分类方程对应的Lagrange乘子。这是一个不等式约束下二次函数寻优问题，存在唯一解。求解上述问题后得到的最优分类函数是：

$f\left(x\right)=\mathrm{sgn}\{(w\&CenterDot;x)+b\}=\mathrm{sgn}\{\underset{i,j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{y}_i(x_i,x)+b\}$

上式中的求和实际上只对支持向量进行。b是分类阈值，可以用任意一个支持向量求得，或通过两类中任意一对支持向量取中值求得。

本发明中对应的问题为非线性，此时可以通过非线性变换将其转换为某个高维空间中的线性问题，在变换空间中求解最优分类面。在线性SVM中，向量之间只进行点积运算。在高维特征空间中，若能找到一个核函数K，使得K(x_i，x_j)＝Φ(x_i)·Φ(x_j)，则只需要进行内积运算。这种内积运算可以用原空间中的函数实现，甚至没有必要获取变换Φ的形式。根据泛函的有关理论，只要一种核函数K(x_i，x_j)满足Mercer条件，它就对应某一变换空间中的内积。因此，在最优分类面中采用适当的内积函数K(x_i，x_j)后，目标函数就变为：

$Q\left({\mathrm{\&alpha;}}_i\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i-\frac{1}{2}\underset{i,j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\mathrm{\&alpha;}}_j{y}_i{y}_{j}K(x_i,x_j)$

而相应的分类函数也变为：

$f\left(x\right)=\mathrm{sgn}\{\underset{i,j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{y}_{i}K(x_i,x)+b\}$

上面的方法在保证训练样本全部被正确分类，即在经验风险为0的前提下，通过最大化分类间隔来获得最好的推广性能。如果希望在经验风险和推广性能之间求得某种平衡，可以通过引入正的松驰因子ξ_i来允许错分样本的存在，这时，分类方程可表示为：

y_i[(w·x_i)+b]-1+ξ_i≥0，i＝1，...，n

而在目标函数里，也就是最小化

里加入惩罚项

这里ξ_i看作训练样本关于分类超平面的偏差，上式折衷考虑样本偏差(即最大分类正确率)和最大分类间隔(即机器泛化能力)，C＞0是个常数，控制对错分样本惩罚的程度。这样SVM最一般的表述可表示成：

Maximize： $Q\left(\mathrm{\&alpha;}\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i-\frac{1}{2}\underset{i,j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\mathrm{\&alpha;}}_j{y}_i{y}_{j}K(x_i,x_j)$

S.t. $\underset{i,j=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i{\mathrm{\&alpha;}}_i=\mathrm{0,0}\&le;{\mathrm{\&alpha;}}_i\&le;C,i=1,...,n$

由此，分类器决策函数如下：

$f\left(x\right)=\mathrm{agn}\{\underset{x_i\&Element;\mathrm{SVs}}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^{*}{y}_{i}K(x_i,x)+b^{*}\}$

其中，

$b^{*}=-\frac{1}{2}\underset{x_i\&Element;\mathrm{SVs}}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^{*}{y}_i[K(x_r,x_i)+K(x_s,x_i)]$

式中x_r与x_s分属不同类的支持向量。SVM的最终输出只由少数的支持向量所确定。

在非线性的SVM算法中，不同的内积函数将形成不同的算法，此处支持向量机的核函数优选为径向基函数，惩罚因子C和核宽度σ两个参数采用网格搜索的方法进行选取。网格搜索的方法具体为，将惩罚因子C和核宽度σ分别取N个值和M个值，对N×M个C和σ的组合分别训练不同的支持向量机，采用留一法评估支持向量机的推广能力，选择最高分类准确率所对应的C和σ的组合，作为支持向量机的最优参数。其中，留一法，即留一法交叉检验(leave-one-out cross validation strategy，LOOCV)，通常对样本个数较少的样本集进行分类时使用留一法，以下以一个样本集为例进行进一步说明。只使用样本集中的一个样本作为验证样本，而剩余样本作为训练样本，重复训练支持向量机和对验证样本进行验证的过程，一直持续到每个样本都被当作一次验证样本；根据每个样本的验证结果计算得到留一法的分类准确率。留一法的分类准确率可用于评估支持向量机的推广能力。参数选择过程使用的训练样本集为病变区域图像特征样本集。

实验证明，C与σ的取值按指数增长的方式排列时，能更有效的找到C和σ的最优组合(例如，C＝2^-5，2^-3，…，2¹⁵，σ＝2^-15，2^-13，…，2³)。最初可以进行粗网格搜索，即C和σ的取值间隔较大(如，C＝2^-5，2，2⁷…，2¹⁹，σ＝2^-15，2^-9，2^-3…2⁹)，得到较优的C和σ取值范围。然后在较优的取值范围中，对C和σ的值进行细化的网格搜索。此处优选采用N＝5，M＝6，惩罚因子的取值范围优选为C∈{2^-1，2⁰，2¹，2²，2³}，核宽度的取值范围优选为σ∈{2^-4，2^-3，2^-2，2^-1，2⁰，2¹}。

根据上述最优参数设定支持向量机的参数，使用病变区域图像特征样本集对支持向量机进行训练，得到训练后的支持向量机。

SVM程序选用的是台湾林智仁(Chih-Jen Lin)教授于2001年开发的支持向量机库(LIBSVM)。支持向量机模型训练可通过下列描述语句实现：

Model＝SVM_Train(S)

其中，输入的S为得到的病变区域图像特征样本集，输出的为训练好的支持向量机模型，具体过程此处不再赘述。

S4构造验证样本集，根据验证样本集对训练后的支持向量机模型进行验证；

将病变区域图像中除提取为病变区域图像训练样本以外的其他图像构造为验证样本集，使用训练后的支持向量机模型对所述验证样本集进行分类，获得各验证样本的脑瘤MIB-1指数的状况。

S5构造病变区域图像检测样本，具体步骤同S1；

S6提取病变区域图像检测样本的图像特征，具体步骤同S2；

S7支持向量机模型根据输入的病变区域图像检测样本的图像特征，识别并获取该病变区域图像检测样本的类别，即其MIB-1指数所处的范围。

基于上述方法进行的具体实施过程如下：

临床采集24个胶质瘤患者的磁共振图像，其中包括T1加权序列、T1增强序列和FLAIR序列。使用免疫组织化学技术检测并分类：24个胶质瘤患者中14位患者MIB-1≤5％，10位患者MIB-1＞5％。由于T1加权序列、T1增强序列和FLAIR序列中的部分序列上，磁共振图像的伪影比较严重，容易影响图像分析和分类的结果，因此不同序列的病变区域图像个数不尽相同。从上述磁共振图像中获得的病变区域图像中，T1加权序列上，MIB-1≤5％和MIB＞5％的病变区域图像个数分别为76和62个；T1增强序列上，MIB-1≤5％和MIB＞5％的病变区域图像个数分别为75个和63个；FLAIR序列上，MIB-1≤5％和MIB＞5％的病变区域图像个数分别为81个和64个。

从每个类别中任选的一定数量包含该类别的病变区域图像构成病变区域图像训练样本，训练样本的样本个数选取如下：T1加权序列上，MIB-1≤5％和MIB＞5％的训练样本个数分别为38个和31个；T1增强序列上，MIB-1≤5％和MIB＞5％的训练样本个数分别为38个和32个；FLAIR序列上，MIB-1≤5％和MIB＞5％的训练样本个数分别为41个和32个。

用S2的方法对上述训练样本提取图像特征，并基于该图像特征，构造病变区域图像特征样本集，并利用该变区域图像特征样本集训练得到支持向量机模型。

S4中的验证样本集的选取如下：T1加权序列上，验证样本集中MIB-1≤5％和MIB＞5％的样本个数分别为38个和31个；T1增强序列上，验证样本集中MIB-1≤5％和MIB＞5％的样本个数分别为37个和31个；FLAIR序列上，验证样本集中MIB-1≤5％和MIB＞5％的样本个数分别为40个和32个。

使用受试者工作特征(receiver operating characteristic，ROC)评价该分类结果的准确性，使用ROC曲线下面积A_z来反映分类的有效性。当A_z取值等于0.5时，分类没有意义，相当于随机分类；当A_z取值等于1是，分类完善，对所有样本都能准确分类。A_z值越大，则分类效果越好。

经过计算得到，对来自T1加权序列的验证样本集进行分类，A_z达到0.6238；对来自T1增强序列的验证样本集进行分类，A_z达到0.6580；对来自FLAIR序列的验证样本集进行分类，A_z达到0.6716，详见图7，支持向量机对验证样本的分类较为准确，有一定的识别率。

实施例2：

由于病变区域图像特征样本集所含的维数较高，可能含有一些冗余的图像特征，这些冗余的图像特征一方面可能降低分类的精度，另一方面会大大增加支持向量机的计算开销；因此作为优选，可基于离散粒子群算法对病变区域图像特征样本集进行优化，得到病变区域图像优化特征集；根据该病变区域图像优化特征集确定支持向量机模型的参数，有效的降低图像特征的复杂程度。

在粒子群算法中，优化问题的每个潜在解都可以想象成搜索空间中的一个点，称为粒子。粒子当前位置的好坏由目标函数进行评估，目标函数根据粒子的位置计算出相应的适应值Ap。每个粒子都知道自己目前为止所发现的最好位置(即历史飞行中适应值Ap最大时对应的位置)，这个可以看作是粒子自己的飞行经验。除此之外，每个粒子还知道目前为止整个群体中所有粒子发现的最好位置，这个可以看作粒子同伴的飞行经验。粒子在搜索空间中以一定的速度飞行，这个速度根据它本身的飞行经验和同伴的飞行经验来动态调整，继而被用来计算粒子的新位置，向空间中最优的位置靠近。优化搜索正是由一群随机初始化形成的粒子所组成的种群中，以迭代的方式进行的，直到满足一定的终止条件，例如达到指定的迭代次数。最后得到的整个群体中所有粒子发现的最好位置即为优化的最优解。在本发明中，粒子是在特征空间中进行搜索的，目的是搜索到使分类结果最优的特征子集，因此适应度是由支持向量机的分类准确率来决定的。

对病变区域图像训练样本集进行特征优化时，病变区域图像训练样本集表示为S_m＝{(f_i，l_i)|i＝1，2，…，n}，其中f_i是第i个训练样本的m维的特征向量，表示为f_t＝[f_i1，f_i2，…，f_ij，…，f_im]，f_ij是第i个训练样本的第j个特征的特征值；l_i是第i个训练样本的类别标记；对于MIB＞5％的训练样本，l_i＝+1；对于MIB-1≤5％的训练样本，l_i＝-1；n是训练样本数量。

据上所述，在实施例1的S2与S3之间，如图6所示，对病变区域图像训练样本的特征优化，包括以下步骤：

S201确定初始的病变区域图像特征样本集子集；

令病变区域图像特征样本集为S_p，S_p中的特征有p维，令p＝m；在离散粒子群算法中，每个粒子即为p维特征搜索空间中的1个点，可以表示为长度为p的二进制位串，每一位的取值为0或为1，1表示相应的特征被选中，0表示相应的特征未被选中。例如，总特征有5维，则粒子(10010)表示第1维和第4维特征被选中构成特征子集。令粒子群中粒子的数目为N，初始化粒子群，即对每个粒子随机赋以一个长度为p的二进制位串。综合实验结果及计算量，本发明优选粒子群的规模为N＝30，即初始化随机产生30个长度为p的二进制位串，得到30个病变区域图像特征样本集子集。

S202将病变区域图像特征样本集子集分别放入支持向量机中用留一交叉检验法进行分类，计算每个病变区域图像特征样本集子集的分类准确率，并将该准确率作为对应粒子的适应度；基于粒子的适应度迭代更新病变区域图像特征样本集子集达一定次数，将适应度最大的粒子对应的特征样本集子集作为病变区域图像优化特征集。

对于每个粒子，将其对应选出的特征子集分别放入支持向量机用留一法进行分类，将得到的分类准确率作为该粒子的适应度Ap。对于N个粒子在第一次迭代计算中得到的N个适应度Ap，选出Ap的最大值Ag_max，并将Ag_max记为群体当前最优值Ag_best，将Ag_best对应粒子的位置记为群体最优位置g_best。对于每个粒子的Ap值，记为每个粒子各自的个体当前最优值Ap_best，将各自Ap_best对应的位置记为个体最优位置p_best。

本步骤中，每一次迭代结束后，更新粒子的速度和位置：

其中，v_ij(k+1)为粒子i在第(k+1)次迭代中第j维的速度，是[0，1]区间的随机数；a是[0.6，0.3]区间上随迭代次数递减的数；x_ij(k)为粒子i在第k次迭代中第j维的当前位置；p_bestij为粒子i在第j维的个体最优位置；g_bestj为整个粒子群在第j维的全局最优位置；并且，对每次迭代完的粒子随机选取10％进行变异(随机赋予二进制值)，以增强粒子跳出局部最优点的能力。

在本实施例中，指定的迭代次数为100次，停止计算时的g_best所对应的病变区域图像优化特征集即为优化得到的病变区域图像优化特征集。

通过对病变区域图像训练样本集的特征优化，特征维数大幅下降，特征的复杂程度得到有效地降低。例如，本发明中T1加权序列、T1增强序列、FLAIR序列上的训练样本集的特征维数为47。经过特征优选，T1加权序列上优选后的训练样本集的特征维数为16，T1增强序列上优选后的训练样本集的特征维数为15，FLAIR序列上优选后的训练样本集的特征维数为22。对于不同序列，具体优选出的特征序号如表2所示，表2中的特征标号基于表1的特征标号。

表2经优化得到的病变区域图像优化特征集

特征优化后，在S3中，使用病变区域图像优化特征集对支持向量机进行训练。首先，选择支持向量机的参数。支持向量机的核函数优选为径向基函数，惩罚因子C和核宽度σ两个参数采用网格搜索的方法进行选取。网格搜索的方法具体为，将惩罚因子C和核宽度σ分别取N个值和M个值，对N×M个C和σ的组合分别训练不同的支持向量机，采用留一法评估支持向量机的推广能力，选择最高分类准确率所对应的C和σ的组合，作为支持向量机的最优参数。需要注意的是，参数选择过程中仅使用优化后的训练样本集，本发明优选采用N＝5，M＝6，惩罚因子的取值优选为C∈{2^-1，2⁰，2¹，2²，2³}，核宽度的取值优选为σ∈{2^-4，2^-3，2^-2，2^-1，2⁰，2¹}。

然后，根据最优参数，设定支持向量机的参数，使用T1加权序列、T1增强序列或FLAIR序列中任一种序列的优化后的训练样本集对支持向量机进行训练，得到优化训练后的支持向量机。

采用未经特征优选的病变区域图像训练样本集对支持向量机的分类计算带来了很大的负担。例如，T1加权序列上，验证样本集中MIB-1≤5％和MIB＞5％的样本分别为38个和31个，使用训练后的支持向量机对验证样本进行分类，需要0.96秒，而使用优选训练后的支持向量机对新的验证样本进行分类，需要0.73秒；T1增强序列上，验证样本集中MIB-1≤5％和MIB＞5％的样本分别为37个和31个，那么使用训练后的支持向量机对验证样本进行分类，需要0.91秒，而使用优选训练后的支持向量机对新的验证样本进行分类，只需要0.68秒；FLAIR序列上，验证样本集中MIB-1≤5％和MIB＞5％的样本分别为40个和32个，那么使用训练后的支持向量机对验证样本进行分类，需要1.15秒，而使用优选训练后的支持向量机对新的验证样本进行分类，只需要0.75秒。因此，为了提高方法的效率，本实施例使用特征优选后的训练样本集对支持向量机进行训练，并使用优选训练后的支持向量机对新的验证样本进行分类。

使用受试者工作特征评价该分类结果的准确性，经过计算得到，对来自T1加权序列的优选后的验证样本集进行分类，A_z达到0.8291；对来自T1增强序列的优选后的验证样本集进行分类，A_z达到0.7977；对来自FLAIR序列的优选后的验证样本集进行分类，A_z达到0.8066。基于特征优化后得到的ROC曲线如图8所示，与未经特征优化得到的ROC曲线(如图7所示)比较，可知：支持向量机对优选后的验证样本能做到更为准确的分类，尤其在T1序列上具有比较高的识别率。

Claims (9)

1.一种脑瘤MIB-1指数范围检测方法，包括以下步骤：

（1）采集脑瘤患者的磁共振图像，构造病变区域图像训练样本；

（2）提取病变区域图像训练样本的图像特征，根据提取的病变区域图像训练样本的图像特征，训练得到支持向量机模型；

（3）构造验证样本集，根据验证样本集对训练后的支持向量机模型进行验证；

将病变区域图像中除提取为病变区域图像训练样本以外的其他图像构造为验证样本集，使用训练后的支持向量机模型对所述验证样本集进行分类，获得各验证样本的脑瘤MIB-1指数的状况；

（4）构造病变区域图像检测样本，提取病变区域图像检测样本的图像特征；根据病变区域图像检测样本的图像特征，用训练得到的支持向量机模型识别并获取该病变区域图像检测样本的MIB-1指数所处的范围。

2.根据权利要求1所述的脑瘤MIB-1指数范围检测方法，其特征在于：在步骤（1）中，所述脑瘤患者的磁共振图像的体数据分辨率为512x512x16体素。

3.根据权利要求1所述的脑瘤MIB-1指数范围检测方法，其特征在于：在步骤（1）中，所述的构造病变区域图像训练样本，包括以下步骤：

a.截取磁共振图像中的病变区域图像；

b.对所述的病变区域图像进行分类，对每个类别选取一定数量的病变区域图像，得到病变区域图像训练样本。

4.根据权利要求3所述的脑瘤MIB-1指数范围检测方法，其特征在于：在步骤b中，所述的对所述的病变区域图像进行分类，包括以下步骤：

i、对病变区域图像的MIB-1指数进行检测，确定MIB-1指数的分界点；

ii、根据所述的MIB-1指数的分界点，对病变区域图像进行分类。

5.根据权利要求1所述的脑瘤MIB-1指数范围检测方法，其特征在于：在步骤（2）中，所述的提取病变区域图像训练样本的图像特征，包括提取病变区域图像训练样本的纹理特征。

6.根据权利要求5所述的脑瘤MIB-1指数范围检测方法，其特征在于：所述的提取病变区域图像训练样本的纹理特征，包括步骤：

使用基本灰度统计信息提取病变区域图像训练样本的图像特征；

使用灰度共生矩阵提取病变区域图像训练样本的图像特征；

使用灰度-梯度共生矩阵提取病变区域图像训练样本的图像特征；

使用游程长度矩阵提取病变区域图像训练样本的图像特征；

使用闵可夫斯基泛函提取病变区域图像训练样本的图像特征。

7.根据权利要求1所述的脑瘤MIB-1指数范围检测方法，其特征在于：在步骤（2）中，根据提取的病变区域图像训练样本的图像特征，训练得到支持向量机模型，包括以下步骤：

a、根据提取的病变区域图像训练样本的图像特征构造病变区域图像特征样本集；

b、根据病变区域图像特征样本集确定支持向量机模型的参数。

8.根据权利要求7所述的脑瘤MIB-1指数范围检测方法，其特征在于：在步骤b中，所述的根据病变区域图像特征样本集确定支持向量机模型的参数，包括：基于离散粒子群算法对病变区域图像特征样本集进行优化，得到病变区域图像优化特征集。

9.根据权利要求8所述的脑瘤MIB-1指数范围检测方法，其特征在于：基于离散粒子群算法对病变区域图像特征样本集进行优化，得到病变区域图像优化特征集，包括以下步骤：

i、确定离散粒子群算法的粒子群规模，每个粒子对应一个随机选取的病变区域图像特征样本集子集；

ii、将每个粒子对应的病变区域图像特征样本集子集分别放入支持向量机中用留一交叉检验法进行分类，计算每个病变区域图像特征样本集子集的分类准确率，并将该准确率作为对应粒子的适应度；基于粒子的适应度迭代更新病变区域图像特征样本集子集达一定次数，将适应度最大的粒子对应的特征样本集子集作为病变区域图像优化特征集。

Images