CN1595451A - 核磁共振多谱图像分割方法 - Google Patents

Abstract

一种用于图像处理领域的核磁共振多谱图像分割方法。具体如下：对于核磁共振图像的多个加权谱建立三维矢量灰度空间，运用矢量梯度算子求得图像中任意一点的矢量梯度值；运用灰度矢量和矢量梯度值生成训练用的特征向量并输入神经网络，运用动态步幅函数的方法动态调节神经网络的学习率；采用层次规约方法对神经网络的神经元节点进行规约；以矢量梯度大小作为筛选标准进一步进行部分有监督的神经网络训练；对训练完成的神经网络用于图像的分割和分类，输出分割完成的图像结果。本发明充分利用核磁共振影像序列中的多维梯度信息进行神经网络训练，同时加入有监督神经网络训练方式，从而实现对多谱核磁共振影像准确分割的目的。

Description

核磁共振多谱图像分割方法

技术领域

本发明涉及一种用于图像处理技术领域的方法，特别是一种核磁共振多谱图像分割方法。

技术背景

核磁共振是一种探测组织氢原子核磁化性能的非电离损害技术，核磁共振图像(MRI)在诊断各种病变及准确定位方面远优于包括CT在内的各种检查技术。MRI成像具有多谱特性，其灰度值主要取决于身体内部固有的组织参数：组织内可以移动的质子密度PD、纵向驰豫时间T1、横向驰豫时间T2。通过改变输入脉冲的时间间隔和组成方式，可获得这三种固有参数的加权图像序列。多谱MRI图像联合分析能充分利用MR各序列的数据，达到提高分割效果的目的。另一方面神经网络技术具有高速并行处理能力、通过实例学习能力和最优搜索能力，并对随机噪声有很强的鲁棒性，将其运用在磁共振图像中可以达到图像定量分析和图像可视化技术前期处理的目的。

经文献检索发现，Suchendra M.Bhandarka等人在IEEE Computer-BasedMedical Systems(《计算机医疗系统》)2001：294-299上发表的“Segmentationof Multispectral MR Images Using a Hierarchical Self-Organizing Map”(《利用多层自组织神经网络的多谱核磁共振图像分割》)，该文中提出利用多层自组织神经网络进行分割的方法，即对原有自组织神经网络的特征向量运用了多谱化的方法，使磁共振多谱图像的信息得到充分的利用。但是该项技术有以下三个缺点：(1)该技术在特征向量中融入了图像的空间位置信息，而对应于实际的图像，有各种可能的成像形式，一概将空间位置信息加入会扰乱图像的正常分割和识别。(2)该技术未能有效利用图像的边界信息，图像边界一般都有较高的梯度值，同时也是图像信息熵较低的部分，如充分发掘这一特性可以加快自组织神经网络的训练过程。(3)该技术为纯粹的无监督分割方法，无法控制神经网络最后训练的精度，无法进行人工交互并融入专家的经验信息，造成训练成功率的下降。

发明内容

本发明针对背景技术的上述不足之处，提出一种核磁共振多谱图像分割方法。使其利用多谱图像中提取出的矢量梯度信息，动态调整自组织神经网络的步幅函数，达到高效无监督训练的目的，而为了进一步提高神经网络的分割能力，还可对已经训练完成的自组织神经网络再进行部分有监督的分割训练。

本发明是通过以下技术方案实现的，方法步骤如下：

步骤一：对于核磁共振图像的多个加权谱建立三维矢量灰度空间；

步骤二：在该三维矢量灰度空间中，运用矢量梯度算子求得图像中任意一点在三维空间中的矢量梯度值；

步骤三：由步骤一和步骤二所生成的灰度矢量和矢量梯度值生成自组织神经网络输入训练用的特征向量；

步骤四：将特征向量输入神经网络，运用动态步幅函数的方法动态调节神经网络的学习率，以提高神经网络的学习效果；

步骤五：由于图像分割得到的实际区域数目远大于目标数目，需要进行合并，采用层次规约的方法对神经网络的神经元节点进行规约。

步骤六：在已形成的神经网络基础上需进一步进行部分有监督的神经元训练，以矢量梯度大小作为筛选标准将训练的像素点集中于图像的边界区，减少有监督训练的训练量，

步骤七：将训练完成的神经网络用于图像的分割和分类，输出分割完成的图像结果。

以下对本发明作进一步的说明，具体内容如下：

所述的步骤二，具体实现如下：

在三维灰度空间中，运用下列矢量梯度算子求得图像任意点的矢量梯度值：

$\mathrm{\ΔH}=\left[\begin{array}{ccc}{H}_{-}& 0& H_{+}\end{array}\right],\mathrm{\ΔV}=\left[\begin{array}{c}{V}_{-}\\ 0\\ V_{+}\end{array}\right]$

其中矢量元素H+(x₀，y₀)、H-(x₀，y₀)、V+(x₀，y₀)、V-(x₀，y₀)代表像素点(x₀，y₀)在(2ω+1)×(2ω+1)滑窗内的平均灰度三维矢量值，可由下式计算；参数ω为正整数，代表了滑窗的大小。

$H_{+}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0-\mathrm{\ω}}{\overset{y=y_0+\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0+1}{\overset{x=x_0+\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

$H_{-}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0-\mathrm{\ω}}{\overset{y=y_0+\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0-1}{\overset{x=x_0-\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

$V_{+}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0+1}{\overset{y=y_0+\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0-\mathrm{\ω}}{\overset{x=x_0+\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

$V_{-}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0-1}{\overset{y=y_0-\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0-\mathrm{\ω}}{\overset{x=x_0+\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

定义(x₀，y₀)点的三维空间梯度值为

$G(x_0,y_0)=\sqrt{{\left|\right|\mathrm{\ΔH}(x_0,y_0)\left|\right|}^2+{\left|\right|\mathrm{\ΔV}(x_0,y_0)\left|\right|}^2}$

    ΔH(x₀，y₀)＝H₊(x₀，y₀)-H_(x₀，y₀)

其中                                   。

    ΔV(x₀，y₀)＝V₊(x₀，y₀)-V_-(x₀，y₀)

所述的步骤四，具体实现如下：

将特征向量输入自组织神经网络进行训练，其中算法采用改进的步幅函数，形式如下：

$\mathrm{\α}(k,G(x,y))=\frac{\mathrm{\α}\left(0\right)}{k}\frac{2}{1+\exp \left(\frac{G(x,y)-\mathrm{\λ}}{\mathrm{\η}}\right)}$

其中λ为梯度阈值，η为倾斜坡度参数，用于调整步幅函数在阈值点附近变化的快慢。

所述的步骤六，具体实现如下：

(1)算法按迭代节拍k＝1，2，…进行，对于每个节拍k，依次或随机地从训练集样本中取出一个训练样本并标明所属的理想类别；

(2)训练样本送入神经元阵列训练，得到一个获胜端和一个次胜端；

(3)对获胜端和次胜端所连神经元权向量进行调整，调整的具体形式如下，其中C_opt(k)是该点所属的理想类别，各神经元权为W_j(k)，获胜端为j*(k)，次胜端为i*(k)。

①若C_f(k)＝C_opt(k)，则

其中θ是一个小正常数。

②若C_f(k)≠C_opt(k)且C_s(k)＝C_opt(k)，

③若C_f(k)≠C_opt(k)且C_s(k)≠C_opt(k)，则

其中l(k)是与类别C_opt(k)相应的神经元编号。

本发明方法主要具有以下特点和作用：(1)采用神经网络分割方式，具有高速并行处理能力、通过实例学习能力和最优搜索能力，并对随机噪声有很强的鲁棒性。(2)充分利用了多谱图像中提取出的矢量梯度信息，动态调整自组织神经网络的步幅函数，达到高效无监督训练的目的。(3)可对已经训练完成的自组织神经网络再进行部分有监督的分割训练，进一步提高了神经网络的分割能力。

具体实施方式

结合本发明方法的具体内容提供以下实施例：

核磁共振设备采集的多谱加权影像序列是医学影像序列的一种。人体不同器官的正常和病理组织的脉冲纵向驰豫时间(T1)和横向驰豫时间(T2)相对恒定，互有差别。组织间弛豫时间和自旋质子密度的差别是核磁共振的成像基础。通过施加不同的脉冲激发序列(如长TR、短TE)，可得到不同加权组合的核磁共振成像。

本例多谱核磁共振影像序列为某某医院利用马可尼医疗设备公司1.5T核磁共振影像设备(Eclipse 1.5T，Marconi Medical Systems，Inc，)采集的17岁男孩的脑部序列。

运行本发明方法的设备为普通个人计算机，其具体配置为P41.8G，256M内存，80G硬盘。发明方法的具体实现平台为Windows2000下的Matlab计算平台。

实施情况如下：

步骤1：对于核磁共振图像的多个加权谱建立三维矢量灰度空间，对于本例实际提取的为T1、T2和质子加权(PD)成像，大小为256×256像素，对应图像中某点坐标(115，189)的点，其灰度三维矢量值M(115，189)＝(138，250，21)。

步骤2：同时在该三维灰度空间中，运用矢量梯度算子求得该点(x₀，y₀)在其矢量梯度值。在三维灰度空间中，运用下列矢量梯度算子求得图像任意点的矢量梯度值：

$\mathrm{\ΔH}=\left[\begin{array}{ccc}{H}_{-}& 0& H_{+}\end{array}\right],\mathrm{\ΔV}=\left[\begin{array}{c}{V}_{-}\\ 0\\ V_{+}\end{array}\right]$

其中矢量元素H+(x₀，y₀)、H-(x₀，y₀)、V+(x₀，y₀)、V-(x₀，y₀)代表像素点(x₀，y₀)在(2ω+1)×(2ω+1)滑窗内的平均灰度三维矢量值，可由下式计算；参数ω为正整数，代表了滑窗的大小，本例中采用ω＝1，划窗大小为3×3。

$H_{+}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0-\mathrm{\ω}}{\overset{y=y_0+\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0+1}{\overset{x=x_0+\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

$H_{-}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0-\mathrm{\ω}}{\overset{y=y_0+\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0-1}{\overset{x=x_0-\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

$V_{+}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0+1}{\overset{y=y_0+\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0-\mathrm{\ω}}{\overset{x=x_0+\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

$V_{-}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0-1}{\overset{y=y_0-\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0-\mathrm{\ω}}{\overset{x=x_{0+}\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

定义(x₀，y₀)点的三维空间梯度值为

$G(x_0,y_0)=\sqrt{{\left|\right|\mathrm{\ΔH}(x_0,y_0)\left|\right|}^2+{\left|\right|\mathrm{\ΔV}(x_0,y_0)\left|\right|}^2}$

    ΔH(x₀，y₀)＝H+(x₀，y₀)-H_(x₀，y₀)

其中                                           ，本例中实际求得G(115，189)＝2143.2

    ΔV(x₀，y₀)＝V+(x₀，y₀)-V_(x₀，y₀)

步骤3：由步骤1和步骤2所生成的灰度矢量和矢量梯度值生成自组织神经网络输入训练用的特征向量。当得到全部样本像素点的灰度三维矢量值和矢量梯度值后，按顺序对每个样本点(x，y)将M(x，y)和G(x，y)输入自组织神经网络。本例中将前面求得的(115，189)点获得的三维矢量值M(115，189)＝(138，250，21)和G(115，189)＝2143.2输入神经网络进行训练。

步骤4：由步骤3所得到的特征向量输入神经网络，运用动态步幅函数的方法动态调节神经网络的学习率，以达到提高神经网络学习效果的目的。具体如下：

(1)随机设置初始权向量W_j(1)，j＝1~L。L为神经元的个数，本例中为300个。设置最大迭代计算次数K(K≥P)，本例中为500次。

(2)按照迭代节拍k＝1，2，…，K，进行下列迭代计算：

对于每一节拍k，由训练集中依次或随机的取出一个样本输入向量并且表示为X(k)，然后用下列迭代公式由已知的Wj(k)求Wj(k+1)：

W_j(k+1)＝W_j(k)+α(k)Λ(j，j*(k)，k)[X(k)-W_j(k)]，j＝1～L

(3)当k＝K，迭代结束并输出Wj(K)，j＝1~L作为学习所得的诸神经元权向量。

算法中α(k)为步幅函数。Λ(j，j*(k)，k)为领域函数，Λ(j，j*(k)，k)按下式计算：

$\mathrm{\Λ}(j,j^{*}\left(k\right),k)=\exp \{-\frac{d_{{\mathrm{jj}}^{*}\left(k\right)}^2\left(k\right)}{{2\mathrm{\σ}}^2\left(k\right)}\}$

称为宽度系数，随着迭代节拍k的增加而下降，影响的领域也越小。

算法同时采用改进的步幅函数，形如下式：

$\mathrm{\α}(k,G(x,y))=\frac{\mathrm{\α}\left(0\right)}{k}\frac{2}{1+\exp \left(\frac{G(x,y)-\mathrm{\λ}}{\mathrm{\η}}\right)}$

本例中，当网络训练到第43步，即k＝43时，σ(43)＝0.113，λ＝1000，η＝13，求得α＝1.27。由此对应各神经元求出领域函数Λ，进一步求出训练后神经元权值Wj(44)

自组织神经网络输出层经过训练后将分成不同的区域，而各区域对应输入有不同的响应特征，将原有特征空间中的聚类中心映射到输出层而保持拓扑结构不变，达到对目标固有特征作出客观的划分的目的。

步骤5：该步采用层次规约的方法步骤四得到的神经网络的神经元节点进行规约。

1)对于自组织图中的神经元统计其最佳命中次数直方图；

2)初始化构造L个类，各包含一个神经元；

3)计算类中心间的距离矩阵，找出类间距离最小的两类将其合并，累计两类的最佳命中次数之和作为新类的命中次数，同时更新类中心的位置；

4)通过随机取样剔除孤立点。如果一个类命中次数增长过慢，就去掉该孤立点；

5)若类已归约到指定数目则停止，否则转3)。

步骤6：为了得到更为精确的分割结果同时为了减少有监督训练的训练量，在已形成的神经网络基础上进一步进行部分有监督的神经元训练。该步采用以矢量梯度大小作为筛选标准将训练的像素点集中于图像的边界区的方法。

算法仍按迭代节拍k＝1，2，…进行。对于每个节拍k，依次或随机地从训练集样本中取出一矢量梯度大于特定值的(x，y)，其中Copt(k)是该点所属的理想类别。将(x，y)送入各神经元权为Wj(k)的阵列，可以有一个获胜端j*(k)，其相应的类别为Cf(k)，还可以有一个次胜端i*(k)，其相应类别为Cs(k)。获胜端的权向量Wj*_(k)和次胜端Wi*_(k)(k)满足下列关系：

‖W_j*(k)(k)-M(k)‖＜‖W_i*(k)(k)-M(k)‖＜‖W_j(k)-M(k)‖

j≠j*(k)，j≠i^*(k)

按下式进行Wj(k+1)的迭代运算，并更新自组织神经网络：

①若C_f(k)＝C_opt(k)，则

其中θ是一个小正常数。

②若C_f(k)≠C_opt(k)且C_s(k)＝C_opt(k)，

③若C_f(k)≠C_opt(k)且C_s(k)≠C_opt(k)，则

对于本例中，在第15节拍取矢量梯度的阈值为1500，取一新点(211，34)进行训练，得到获胜端为j*(15)和次胜端i*(15)，对应类别分别为Cf(15)：白质和Cs(k)：脑脊液，结果Copt(15)正是白质，所以对应上述情况为①，神经元权值得到加强，求得所有Wj(16)的值。

对Wj(k)调整的结果是使那些分布在交叠区域内的矢量，向着正确分类的那一面靠近，这样通过有监督的学习，使交叠区域内能在竞争中获胜的神经元数目越来越少，其分界面也越来越明确。

步骤7：对训练完成的神经网络用于图像的分割和分类，输出分割完成的图像结果。

本发明方法取得了较好的实施效果：最终分割结果显示出了对脑部几种重要的生物组织(灰质(Graymatter)、白质(Whitematter)、皮层(Cortex)、脑脊液(CSF))的多谱序列分割效果良好，在不同组织之间的敏感边缘处都得到了较好的反映。

下表给出了分类结果与参考分类的分类结果的定量比较：

	动态无监督分割  部分有监督分割  参考分割
		灰质白质皮层脑脊液背景	22.21％         19.09％         22.65％14.66％         17.92％         17.97％4.08％          3.58％          4.64％3.05％          2.73％          2.87％56.03％         56.64％         51.90％

Claims (5)

1、一种核磁共振多谱图像分割方法，其特征在于，方法步骤如下：

步骤一：对于核磁共振图像的多个加权谱建立三维矢量灰度空间；

步骤二：在该三维矢量灰度空间中，运用矢量梯度算子求得图像中任意一点在三维空间中的矢量梯度值；

步骤三：由步骤一和步骤二所生成的灰度矢量和矢量梯度值生成自组织神经网络输入训练用的特征向量；

步骤四：将特征向量输入神经网络，运用动态步幅函数的方法动态调节神经网络的学习率，以提高神经网络的学习效果；

步骤五：由于图像分割得到的实际区域数目远大于目标数目，需要进行合并，采用层次规约的方法对神经网络的神经元节点进行规约。

步骤六：在已形成的神经网络基础上需进一步进行部分有监督的神经元训练，以矢量梯度大小作为筛选标准将训练的像素点集中于图像的边界区，减少有监督训练的训练量，

步骤七：将训练完成的神经网络用于图像的分割和分类，输出分割完成的图像结果。

2、根据权利要求1所述的核磁共振多谱图像分割方法，其特征是，所述的步骤二，具体实现如下：

在三维灰度空间中，运用下列矢量梯度算子求得图像任意点的矢量梯度值：

$\mathrm{\ΔH}=[H\_0H_{+}],\mathrm{\ΔV}=\left[\begin{array}{c}{V}_{-}\\ 0\\ V_{+}\end{array}\right]$

其中矢量元素H+(x₀，y₀)、H-(x₀，y₀)、V+(x₀，y₀)、V-(x₀，y₀)代表像素点(x₀，y₀)在(2ω+1)×(2ω+1)滑窗内的平均灰度三维矢量值，参数ω为正整数，代表了滑窗的大小。

3、根据权利要求2所述的核磁共振多谱图像分割方法，其特征是，像素点(x₀，y₀)在(2ω+1)×(2ω+1)滑窗内的平均灰度三维矢量值，通过以下公式获得：

$H_{+}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0-\mathrm{\ω}}{\overset{y=y_0+\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0+1}{\overset{x=x_0+\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

$H_{-}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0-\mathrm{\ω}}{\overset{y=y_0+\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0-1}{\overset{x=x_0-\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

$V_{+}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{y=y_0+1}{\overset{y=y_0+\mathrm{\ω}}{\mathrm{\Σ}}}\underset{x=x_0-\mathrm{\ω}}{\overset{x=x_0+\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

$V_{-}(x_0,y_0)=\frac{1}{\mathrm{\ω}(2\mathrm{\ω}+1)}\underset{{y=y}_0-1}{\overset{y=y_0-\mathrm{\ω}}{\mathrm{\Σ}}}\underset{{x=x}_0-\mathrm{\ω}}{\overset{x=x_0+\mathrm{\ω}}{\mathrm{\Σ}}}M(x,y)$

定义(x₀，y₀)点的三维空间梯度值为

$G(x_0,y_0)=\sqrt{{\left|\right|\mathrm{\ΔH}(x_0,y_0)\left|\right|}^2+{\left|\right|\mathrm{\ΔV}(x_0,y_0)\left|\right|}^2}$

        ΔH(x₀，y₀)＝H+(x₀，y₀)-H-(x₀，y₀)。

其中

        ΔV(x₀，y₀)＝V+(x₀，y₀)-V-(x₀，y₀)

4、根据权利要求1所述的核磁共振多谱图像分割方法，其特征是，所述的步骤四，具体实现如下：

将特征向量输入自组织神经网络进行训练，其中算法采用改进的步幅函数，形式如下：

$\mathrm{\α}(k,G(x,y))=\frac{\mathrm{\α}\left(0\right)}{k}\frac{2}{1+\exp \left(\frac{G(x,y)-\mathrm{\λ}}{\mathrm{\η}}\right)}$

其中λ为梯度阈值，η为倾斜坡度参数，用于调整步幅函数在阈值点附近变化的快慢。

5、根据权利要求1所述的核磁共振多谱图像分割方法，其特征是，所述的步骤六，具体实现如下：

(1)算法按迭代节拍k＝1，2，…进行，对于每个节拍k，依次或随机地从训练集样本中取出一个训练样本并标明所属的类别；

(2)训练样本送入神经元阵列训练，得到一个获胜端和一个次胜端；

(3)对获胜端和次胜端所连神经元权向量进行调整，调整的具体形式如下，其中C_opt(k)是该点所属的理想类别，各神经元权为W_j(k)，获胜端为j*(k)，次胜端为i*(k)：

①若C_f(k)＝C_opt(k)，则

其中θ是一个小正常数；

②若C_f(k)≠C_opt(k)且C_s(k)＝C_opt(k)，

③若C_f(k)≠C_opt(k)且C_s(k)≠C_opt(k)，则

其中l(k)是与类别C_opt(k)相应的神经元编号。