CN103700104A - 一种利用先验知识的人类外侧膝状体自动分割方法 - Google Patents

Abstract

本发明公开了一种利用先验知识的人类外侧膝状体LGN自动分割方法。该方法首先将人脑的结构像数据进行偏场校正，然后将经过校正的图像进行大脑子结构的分割从而得到腹侧间脑区的模板。接着将校正过的结构像数据和腹侧间脑区的模板配准到MNI标准空间。根据外侧膝状体解剖结构的先验知识得到LGN在MNI空间内的区域界限，接着在这个区域内分别用区域增长法、k均值法、大津算法和图割法在分割出LGN。然后用将几种分割结果融合得到LGN区域的估计值作为分割结果。最后将分割结果变换到结构像的原始空间，即为LGN的最终分割结果。

Description

一种利用先验知识的人类外侧膝状体自动分割方法

技术领域

本发明涉及计算机图像处理领域，特别涉及一种利用外侧膝状体（LGN）的解剖先验知识来在人脑结构像中自动分割出LGN的医学图像分割方法。

背景技术

图像分割是指利用图像信息将图像中“感兴趣的目标物体”从复杂的景象中提取出来。感兴趣的目标物体一般被称为前景，图像中其余部分被称为背景。整体而言，图像分割方法可以分为两种体系：人机交互式分割和计算机全自动分割。前者过于依赖专家经验且耗时太长，后者自主分割的结果往往不能不能令人满意。在实际应用领域中，人们针对不同需求来选择适应的方法，到目前为止还不存在一种通用的分割方法。

外侧膝状体是前段视路传入神经纤维的终止点，交叉与未交叉的神经纤维最终在该处换神经元后向后形成视放射。作为一个视路的第四个神经元（前三个分别是视网膜内的椎体细胞、杆体细胞，视网膜内的双极细胞和视网膜内的节细胞），外侧膝状体在视觉通路中扮演着重要的角色。但由于作为丘脑一部分的LGN被颞叶和海马回遮挡且LGN自身体积很小，LGN很难被观察到。通过解剖人脑，人们发现LGN必然在腹侧间脑区（VDC）内。1997年Andrews等人发现LGN大致在一个5mm*5mm*5mm的立方体区域内,而2004年Kastner等人发现LGN中心点在Talairach空间中的位置大概在左侧（-23.33，-21，-4.66）和右侧（22.88，-21.3,4.63）附近。

区域增长法、k均值法、大津算法和图割算法都是图像分割领域中常用的算法。区域增长法利用图像的灰度信息将位置接近且灰度相似的像素聚集成一个区域来得到分割结果。k均值法是一种典型的聚类算法，它利用图像特征将图像的体素聚成指定的类别数。大津算法是图像分割中阈值选取的最佳算法，它遍历图像的灰度值来寻找一个最佳阈值将图像分成背景和前景两部分从而使得背景和前景之间的类间方差最大。图割算法将图论引入图像分割领域，图像的体素作为图的节点，体素之间的相关性作为连接节点边的权重，寻找一个切割边的方式将节点分成两部分，切割损失能量最小的方式即为最佳切割方式。

在医学图像领域中，人工分割目标物体非常耗时，但是自动分割方法准确度不高，为了解决低准确度的问题，我们利用若干不同方式的分割结果来计算出目标物体真实分割的概率估计值，把这个估计值作为目标物体的分割结果。

发明内容

（一）要解决的技术问题

本发明所要解决的技术问题是提供一种自动的LGN分割方法，以对人脑结构像中的LGN进行准确、自动的分割。

（二）技术方案

为解决上述的技术问题，本发明提出一种利用先验知识的人类外侧膝状体自动分割方法，用于将人脑结构像中的LGN分割出来。该方法包括如下步骤：步骤S1：将人脑的结构像数据进行不均匀场校正；步骤S2：将校正过的结构像配准到MNI标准空间；步骤S3：利用LGN解剖结构的先验知识得到LGN在MNI空间内的大致区域；步骤S4：在这个区域内分别用区域增长法、k均值法、大津算法和图割法进行分割；步骤S5：将四种分割结果融合得到对LGN的真实分割的估计；步骤S6：将估计出的结果变换到结构像的原始空间，得到LGN的最终分割结果。

（三）有益效果

本发明利用LGN解剖结构的先验知识制作LGN先验知识模板，将分割区域限定在这个模板中。这一措施将解剖结构先验知识融入LGN分割中，既避免了人机交互分割存在的耗时问题，也提高了传统自动分割的准确性。本发明还将几种分割方法融合来得到LGN真实分割的估计值，进一步提高自动分割的准确性。

附图说明

图1是本发明的基于先验知识的人类外侧膝状体自动分割方法的流程图;

图2是LGN分割过程中图像的变化图，其中（a）是人脑的原始MR结构像；（b）是经过不均匀场校正的MR结构像；（c）是配准到MNI标准空间的MR结构像；（d）是MNI标准空间下LGN的先验知识模板；（e）是融合四种分割方法得到的对LGN分割的估计值；（f）是原始空间下LGN的最终分割结果；

图3是LGN分割结果的三维显示图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

本发明的核心思想是将LGN解剖结构的先验知识融入LGN的自动分割中，从而解决了人机交互的耗时问题和提高自动分割的准确性。本发明还融合几种分割方法来估计LGN的真实分割，从而进一步提高自动分割的准确性。为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明做进一步的详细说明。图1显示了本发明的外侧膝状体自动分割方法的流程图。如图1所示，本发明的方法包括如下步骤：

步骤S1：将人脑的结构像数据进行不均匀场校正。

利用核磁共振技术得到人脑的结构像数据，然后利用N3（nonparametric nonuniform intensity normalization，非参数不均匀灰度校正）算法对不均匀场进行校正。具体方法参见Sled,J.G.,A.P.Zijdenbos,et al.(1998)."A nonparametric method for automatic correction of intensitynonuniformity in MRI data."IEEE Trans Med Imaging17(1):87-97。

步骤S2：将校正过的结构像配准到MNI标准空间。

将图像配准到MNI标准空间是指将图像配准到MNI模板上。MNI模板是由蒙特利尔神经学研究机构(Montreal Neurological Institute)制作的。该机构利用305个正常人的大脑图像通过手工标记关键点和线性配准等方式得到最后的平均模板，这个模板被称为MNI模板。具体MNI模板制作方法参见Evans,A.C.,D.L.Collins,et al.(1993)."3d StatisticalNeuroanatomical Models from305Mri Volumes."Nuclear ScienceSymposium&Medical Imaging Conference,Vols1-3:1813-1817。

首先利用12个参数的仿射变换将校正后的图像配准到MNI模板，这个变换过程解决了图像和MNI模板大脑尺寸和位置的主要差异问题。然后用非线性配准方法处理图像和MNI模板的一些小尺寸结构差异。具体方法参见Friston,K.J.,J.Ashburner,et al.(1995)."Spatial registration andnormalization of images."Hum Brain Mapp3(3):165-189。

步骤S3：利用LGN解剖结构的先验知识得到LGN在MNI空间内的大致区域。

由于传统的医学分割方法中，全自动分割的准确性较低，人机交互分割因依赖专家经验太耗时，本发明为了解决这两个问题，将专家经验即LGN解剖结构的先验知识转化成计算机可以操作的步骤，从而制作MNI空间下的LGN先验知识模板。具体操作如下：

步骤S3a：LGN必然在腹侧间脑区（VDC）内，所以将经过校正的图像进行大脑子结构的分割从而得到VDC模板，并将VDC模板配准到MNI标准空间。

本发明利用Bruce Fischl等人制作的脑图谱将大脑分成45个子结构，然后挑出左右VDC区域制作VDC模板。接着利用步骤S2提到的方法，将VDC模板配准到MNI标准空间。利用脑图谱分割大脑的方法参见Fischl,B.,D.H.Salat,et al.(2002)."Whole brain segmentation:automated labelingof neuroanatomical structures in the human brain."Neuron33(3):341-355。

步骤S3b：在Talairach空间中LGN的中心点位置大概在左侧（-23.33，-21，-4.66）和右侧（22.88，-21.3,4.63）附近。将这两个坐标转换成MNI空间坐标为左侧P_l（-23，-22，-7）和右侧P_r（26，-22，-8）；

步骤S3c：左右两侧LGN大致在一个5mm*5mm*5mm的立方体空间内，考虑到转换带来的误差，分别构建以为P_l和P_r为中心，7mm为边长的立方体模板；

步骤S3d：将VDC模板和立方体模板求交集，构建LGN先验知识模板。

步骤S4：在这个区域内分别用区域增长法、k均值法、大津算法和图割法进行分割。

步骤S4a：在区域增长法中，以P_l和P_r为起始点，半径为5mm来进行区域增长，将区域增长结果和LGN先验知识模板求交集得到分割结果。

把LGN先验知识模板中的第i个体素P_i(i=1,2,…,N)的灰度值表示为I_i(i=1,2,…,N)，其坐标向量表示为V_i，区域增长的约束条件是候选点与起始点P_s的距离不超过半径R₀，且灰度值的差异不超过I₀。如下：

$\begin{array}{c}\mathrm{Region}={\∪}_{i=1}^N{P}_i,(i=\mathrm{1,2},\·\·\·,N)\\ \mathrm{if}\left|\right|V_i-V_s\left|\right|\≤R_0,|I_i-I_s|\≤I_0\end{array}$

这里Region是区域增长法的结果，N是LGN先验知识模板中体素的个数，V_s和I_s是起始点的坐标向量和灰度值。本发明中R₀取5mm，I₀为30。

步骤S4b：在k均值法中，体素特征由体素的灰度值和其在MNI空间中的坐标构成，利用k均值法和体素特征将LGN先验知识模板中的体素聚成两类，体素数目多的一类为分割结果。

步骤S4c：在大津算法中，遍历LGN先验知识模板中的灰度值来寻找分割LGN的最佳阈值从而使得分割得到结果的类间方差最大，体素数目多的一类为分割结果。

大津算法是图像分割中阈值选取的最佳算法。它是按图像的灰度特性,将图像分成背景和前景两部分。背景和前景之间的类间方差越大,说明构成图像的两部分的差别越大,当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。

设LGN先验知识模板限定区域中图像的灰度级是L，则灰度范围为[0,L-1]，若阈值为t，则前景和背景的类间方差定义为：

f(t)=w₀(t)*[u₀(t)-u]²+w₁(t)*[u₁(t)-u]²,t∈[0,L-1]

其中w₀为背景所占比例，u₀为背景的灰度均值，w₁为前景所占比例，u₁为前景的灰度均值，u为LGN先验知识模板限定区域中所有体素的均值。遍历所有阈值t，使得使以上表达式值最大的t，即为分割图像的最佳阈值。

利用最佳阈值将图像分成两类，体素数目多的一类为分割结果。

步骤S4d：在图割算法中，LGN先验知识模板中的体素作为图的节点，利用体素坐标信息和体素灰度信息来计算体素间的相似性，把相似性作为节点之间的边的权重，利用图割算法把图分成两部分，使得分割后能量损失最小，体素数目多的一类为分割结果。

本发明定义体素i和体素j间相似性为

$\mathrm{similarity}(i,j)=w_c\·\exp (-\frac{{\left|\right|f_{\mathrm{ci}}-f_{\mathrm{cj}}\left|\right|}^2}{2{\mathrm{\σ}}_c^2})+w_g\·\exp (-\frac{{\left|\right|f_{\mathrm{gi}}-f_{\mathrm{gi}}\left|\right|}^2}{2{\mathrm{\σ}}_g^2})$

其中f_ci=[x_i,y_i,z_i]和fc_j=[x_j,y_j,z_j]分别是体素i和体素j坐标信息构成的特征，f_gi=[grayvalue_i,6-neighborvalue_i,meanvalue_i,medvalue_i,stdvalue_i]和f_gj=[grayvalue_j,6-neighborvalue_j,meanvalue_j,medvalue_j,stdvalue_j]分别是体素i和体素j附近灰度信息构成的特征，grayvalue是该体素对应的灰度值，6-neighborvalue是该体素上下左右前后六个体素的灰度值，meanvalue是该体素和周围六个体素的灰度平均值，medvalue是该体素和周围六个体素的灰度中值，stdvalue是该体素和周围六个体素的灰度标准差。w_c,w_g分别是坐标特征和灰度特征对应的权重，σ_c,σ_g分别是坐标特征和灰度特征对应的核宽；

定义能量损失函数E为

$E=\frac{\mathrm{cut}(A,B)}{\mathrm{assov}(A,V)}+\frac{\mathrm{cut}(A,B)}{\mathrm{assov}(B,V)}$

其中V是原始连通图，A和B是图切后得到的两部分，cut(A,B)是图V中连接A和B中节点的所有边的权重之和，assov(A,V)是有一个节点在A中的所有边的权重之和。

将分割结果向量记为S=[s₁,s₂,…,s_N]^T，每个分量s_i代表V图中每个节点的分割结果：如果s_i值是1，则对应节点属于一个子图，如果s_i值是

则对应节点属于另一个子图。将V图的相似性矩阵记为

，用

表述次数矩阵（每个对角线元素是相应节点权重之和，非对角线元素为0），则

而能量损失函数E转换成

为了计算出使得E最小的二值向量S，先计算出使得E最小的实数向量S，然后通过设置阈值决定每个节点属于哪一边。经过变换可以知道E的最小值应该是矩阵

的第二小特征值（因为最小特征值是0，而E不会是0），求出对应的特征向量即为实数向量S。然后遍历S中的值作为阈值，将S转成只有1和-b的二值向量，计算对应的损失函数E，从中选择一个阈值使得E最小，此时对应的S则将图像分成两类，体素数目多的一类为图割算法的分割结果。

步骤S5：融合四种分割方法得到LGN区域的估计值作为分割结果。

为了进一步提到自动分割算法的准确性，这里把若干不同方式的分割结果融合来估计LGN的真实分割。

步骤S5a：构造LGN分割的似然估计函数。

步骤S4将LGN先验知识模板区域中的所有体素(N个)用四种分割算法独立赋予类标志Label,Label∈{0,1}，对于某体素i(i=1,2,…,N)，T(i)为其真实类标志，记D_m(i)为第m种分割算法赋予第i个体素的类标志，将先验知识模板区域中的N个体素的四种分割结果组成N×4维的决策矩阵D，记p=[p₁,p₂,p₃,p₄]^T和q=[q₁,q₂,q₃,q₄]^T，p_j和q_j代表第j(j=1,2,3,4)种分割算法的敏感度和灵敏度，分别用(D,T)和f(D,T|p,q)来表示图像分割的完全数据和其概率分布函数，而似然估计函数则为ln[f(D,T|p,q)]。

步骤S5b：利用EM算法进行最大化估计得到LGN分割的估计值。

EM算法中第k次迭代的E步骤计算公式为

${W_i}^{\left(k\right)}=\frac{{a_i}^{\left(k\right)}}{{a_i}^{\left(k\right)}+{b_i}^{\left(k\right)}}$

这个式子用来估算第i个体素被标记为1的概率W_i ^(k)，其中

${a_i}^{\left(k\right)}=f(T_i=1)\underset{j:D_{\mathrm{ij}}=1}{\mathrm{\Π}}{p_j}^{\left(k\right)}\underset{j:D_{\mathrm{ij}}=0}{\mathrm{\Π}}(1-{p_j}^{\left(k\right)}),$

${b_i}^{\left(k\right)}=f(T_i=0)\underset{j:D_{\mathrm{ij}}=0}{\mathrm{\Π}}{q_j}^{\left(k\right)}\underset{j:D_{\mathrm{ij}}=1}{\mathrm{\Π}}(1-{q_i}^{\left(k\right)}),$

这里f(T_i=1)和f(T_i=0)分别表示第i(i=1,2,…,N)个体素真实标志是1和0的概率，D_ij(i=1,2,…,N;j=1,2,3,4)是N×4维决策矩阵D中的元素，k是迭代次数，循环迭代，直至W_i ^(k)，p^(k)和q^(k)收敛。

EM算法中第k次迭代的M步骤计算公式为

$\begin{array}{cc}{p_j}^{(k+1)}=\frac{{\mathrm{\Σ}}_{i:D_{\mathrm{ij}}=1}{W_i}^{\left(k\right)}}{{\mathrm{\Σ}}_i{W_i}^{\left(k\right)}}& {q_j}^{(k+1)}=\frac{{\mathrm{\Σ}}_{i:D_{\mathrm{ij}}=0}(1-{W_i}^{\left(k\right)})}{{\mathrm{\Σ}}_i(1-{W_i}^{\left(k\right)})}\end{array}$

这两个式子用来估计第k+1次迭代中使用的敏感度p^(k+1)和灵敏度q^(k+1)，其中W_i ^(k)是第i个体素被标记为1的概率，D_ij(i=1,2,…,N;j=1,2,3,4)是N×4维决策矩阵D中的元素，k是迭代次数，循环迭代，直至W_i ^(k)，p^(k)和q^(k)收敛。

在本步骤中先将W_i ^(k)，p^(k)和q^(k)初始化，然后利用E步骤和M步骤的算式计算，循环迭代，直至收敛。最后利用阈值0.5将W_i ^(k)二值化，若W_i ^(k)>0.5则第i个体素在LGN区域内，反之，则第i个体素不在LGN区域内。

更多细节参见Warfield,S.K.,K.H.Zou,et al.(2004)."Simultaneoustruth and performance level estimation(STAPLE):an algorithm for thevalidation of image segmentation."IEEE Trans Med Imaging23(7):903-921.

步骤S6：将估计出来的分割结果变换到结构像的原始空间，得到LGN的最终分割结果。

因为本发明为了利用MNI空间内的解剖先验知识，对LGN的分割过程均是在MNI空间下进行，最终要利用步骤S2得到的变换信息将MNI空间下的LGN分割结果反变换到结构像的原始空间，完成整个LGN的分割。

下面通过一个具体应用来进一步说明本发明的方法。在一台配置有Intel Xeon E5-2687W八核3.10GHz的处理器、64Gb的内存、操作系统为Linux CentOS6.3的工作站上使用FreeSurfer和MATLAB软件完成的。

对一个正常人的大脑结构像进行LGN分割，大脑结构像的原始图像如图2（a）所示。

步骤S1：利用N3算法对人脑的结构像数据进行不均匀场校正，校正后图像如图2（b）所示。

步骤S2：利用12个参数的放射变换和非线性配准方法将校正后的图像配准到MNI标准空间，配准后图像如图2（c）所示。

步骤S3：对原始结构像分割出VDC区域并将VDC模板配准到MNI空间。然后利用LGN中心点位置构造立方体模板，将VDC模板和立方体模板求交集得到LGN的先验知识模板，如图2（d）所示。

步骤S4：在LGN的先验知识模板内分别用区域增长法、k均值法（kmeans）、大津算法（otsu）和图割法（Ncut）进行分割。

步骤S5：融合四种分割结果中得到对LGN的真实分割的估计，这个估计值即为MNI空间下LGN的分割结果。图2（e）是融合四种分割方法得到的对LGN分割的估计值。

步骤S6：将MNI空间下估计出的LGN分割结果变换到结构像的原始空间，得到LGN的最终分割结果，如图2（f）和图（3）所示。

由上述对于本发明的具体实施方式，并辅以实施例的说明可见，本发明利用LGN解剖结构的先验知识制作LGN先验知识模板，将先验知识融入LGN分割中，解决了人机交互分割存在的耗时问题并提高传统自动分割的准确性。然后将多种分割方法融合得到LGN真实分割的估计值，进一步提高自动分割的准确性。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种利用先验知识的人类外侧膝状体LGN自动分割方法，该方法包括如下步骤：

步骤S1：将人脑的结构像数据进行不均匀场校正；

步骤S2：将校正过的结构像配准到MNI标准空间；

步骤S3：利用LGN解剖结构的先验知识得到LGN在MNI空间内的大致区域；

步骤S4：在所述区域内分别用区域增长法、k均值法、大津算法和图割法进行分割；

步骤S5：将几种分割算法融合得到LGN区域的估计值作为分割结果；

步骤S6：将估计的LGN分割结果变换到结构像的原始空间，得到LGN的最终分割结果。

2.如权利要求1所述的方法，其特征在于，所述步骤S3中的利用LGN解剖结构的先验知识得到LGN在MNI空间内的大致区域包括：

步骤S3a：LGN必然在腹侧间脑区VDC内，所以将经过校正的图像进行大脑子结构的分割从而得到VDC模板，并将VDC模板配准到MNI标准空间；

步骤S3b：在Talairach空间中LGN的中心点位置在左侧（-23.33，-21，-4.66）和右侧（22.88，-21.3,4.63）附近，将这两个坐标转换成MNI空间坐标为左侧P_l（-23，-22，-7）和右侧P_r（26，-22，-8）；

步骤S3c：左右两侧LGN在一个5mm*5mm*5mm的立方体空间内，分别构建以为P_l和P_r为中心，7mm为边长的立方体模板；

步骤S3d：将VDC模板和立方体模板求交集，构建LGN先验知识模板。

3.如权利要求1所述的方法，其特征在于，所述步骤S4中的步骤包括：

步骤S4a：在区域增长法中，以P_l和P_r为起始点，半径为5mm来进行区域增长，将区域增长结果和LGN先验知识模板求交集得到分割结果；

步骤S4b：在k均值法中，体素特征由体素的灰度值和其在MNI空间中的坐标构成，利用k均值法和体素特征将LGN先验知识模板中的体素聚成两类，体素数目多的一类为分割结果；

步骤S4c：在大津算法中，遍历LGN先验知识模板中的灰度值来寻找分割LGN的最佳阈值从而使得两类的类间方差最大，体素数目多的一类为分割结果；

步骤S4d：在图割算法中，LGN先验知识模板中的体素作为图的节点，利用体素坐标信息和体素灰度信息来计算体素间的相似性，把相似性作为节点之间的边的权重，利用图割算法把图分成两部分，使得分割后能量损失最小，体素数目多的一类为分割结果。

4.如权利要求3所述的方法，其特征在于，在步骤S4d中体素i和体素j间相似性的定义是

$\mathrm{similarity}(i,j)=w_c\·\exp (-\frac{{\left|\right|f_{\mathrm{ci}}-f_{\mathrm{cj}}\left|\right|}^2}{2{\mathrm{\σ}}_c^2})+w_g\·\exp (-\frac{{\left|\right|f_{\mathrm{gi}}-f_{\mathrm{gi}}\left|\right|}^2}{2{\mathrm{\σ}}_g^2})$

其中f_ci=[x_i,y_i,z_i]和fc_j=[x_j,y_j,z_j]分别是体素i和体素j坐标信息构成的特征，f_gi=[grayvalue_i,6-neighborvalue_i,meanvalue_i,medvalue_i,stdvalue_i]和f_gj=[grayvalue_j,6-neighborvalue_j,meanvalue_j,medvalue_j,stdvalue_j]分别是体素i和体素j附近灰度信息构成的特征，grayvalue是该体素对应的灰度值，6-neighborvalue是该体素上下左右前后六个体素的灰度值，meanvalue是该体素和周围六个体素的灰度平均值，medvalue是该体素和周围六个体素的灰度中值，stdvalue是该体素和周围六个体素的灰度标准差,w_c,w_g分别是坐标特征和灰度特征对应的权重，σ_c,σ_g分别是坐标特征和灰度特征对应的核宽；

在步骤S4d中能量损失函数E的定义是

$E=\frac{\mathrm{cut}(A,B)}{\mathrm{assov}(A,V)}+\frac{\mathrm{cut}(A,B)}{\mathrm{assov}(B,V)}$

其中V是原始连通图，A和B是图切后得到的两部分，cut(A,B)是图V中连接A和B中节点的所有边的权重之和，assov(A,V)是有一个节点在A中的所有边的权重之和。

5.如权利要求1所述的方法，其特征在于，所述步骤S5中将四种分割算法融合来估计出LGN的真实分割，步骤包括：

步骤S5a：构造LGN分割的似然估计函数；

步骤S5b：利用EM算法进行最大化估计得到LGN分割的估计值。

6.如权利要求5中所述的方法，其特征在于，在步骤S5a中构造LGN分割的似然估计函数过程是：

步骤S4将LGN先验知识模板区域中的所有体素(N个)用四种分割算法独立赋予类标志Label,Label∈{0,1}，对于某体素i(i=1,2,…,N)，T(i)为其真实类标志，记D_m(i)为第m种分割算法赋予第i个体素的类标志，将先验知识模板区域中的N个体素的四种分割结果组成N×4维的决策矩阵D，记p=[p₁,p₂,p₃,p₄]^T和q=[q₁,q₂,q₃,q₄]^T，p_j和q_j代表第j(j=1,2,3,4)种分割算法的敏感度和灵敏度，分别用(D,T)和f(D,T|p,q)来表示图像分割的完全数据和其概率分布函数，而似然估计函数则为ln[f(D,T|p,q)]。

7.如权利要求5所述的方法，其特征在于，在步骤S5b中利用EM算法中第k次迭代的E步骤计算公式为

${W_i}^{\left(k\right)}=\frac{{a_i}^{\left(k\right)}}{{a_i}^{\left(k\right)}+{b_i}^{\left(k\right)}}$

这个式子用来估算第i个体素被标记为1的概率W_i ^(k)，其中

${a_i}^{\left(k\right)}=f(T_i=1)\underset{j:D_{\mathrm{ij}}=1}{\mathrm{\Π}}{p_j}^{\left(k\right)}\underset{j:D_{\mathrm{ij}}=0}{\mathrm{\Π}}(1-{p_j}^{\left(k\right)}),$

${b_i}^{\left(k\right)}=f(T_i=0)\underset{j:D_{\mathrm{ij}}=0}{\mathrm{\Π}}{q_j}^{\left(k\right)}\underset{j:D_{\mathrm{ij}}=1}{\mathrm{\Π}}(1-{q_i}^{\left(k\right)}),$

这里f(T_i=1)和f(T_i=0)分别表示第i(i=1,2,…,N)个体素真实标志是1和0的概率，D_ij(i=1,2,…,N;j=1,2,3,4)是N×4维决策矩阵D中的元素，k是迭代次数，循环迭代，直至W_i ^(k)，p^(k)和q^(k)收敛,

步骤S5b中利用EM算法中第k次迭代的M步骤计算公式为

$\begin{array}{cc}{p_j}^{(k+1)}=\frac{{\mathrm{\Σ}}_{i:D_{\mathrm{ij}}=1}{W_i}^{\left(k\right)}}{{\mathrm{\Σ}}_i{W_i}^{\left(k\right)}}& {q_j}^{(k+1)}=\frac{{\mathrm{\Σ}}_{i:D_{\mathrm{ij}}=0}(1-{W_i}^{\left(k\right)})}{{\mathrm{\Σ}}_i(1-{W_i}^{\left(k\right)})}\end{array}$

这两个式子用来估计第k+1次迭代中使用的敏感度p^(k+1)和灵敏度q^(k+1)，其中W_i ^(k)是第i个体素被标记为1的概率，Di_j(i=1,2,…,N;j=1,2,3,4)是N×4维决策矩阵D中的元素，k是迭代次数，循环迭代，直至W_i ^(k)，p^(k)和q^(k)收敛,

步骤S5b中先将W_i ^(k)，p^(k)和q^(k)初始化，然后利用E步骤和M步骤的算式计算，循环迭代，直至收敛，最后利用阈值0.5将W_i ^(k)二值化，若W_i ^(k)>0.5则第i个体素在LGN区域内，反之，则第i个体素不在LGN区域内。

Images