基于截断Dirichlet过程无限Student’s t 混合模型的脑部核磁共振图像分割方法

技术领域

本发明涉及计算机辅助图像处理技术，更具体地说，涉及一种基于截断Dirichlet过程无限Student’s t混合模型的脑部核磁共振图像分割方法。

背景技术

最近几十年，医学影像技术飞速发展，特别是核磁共振成像技术具有无创伤性、检查范围覆盖人体各个系统、成像数据丰富等优点，应用最为广泛。医学影像分析方法随之得到关注、研究和应用，医学图像分割是医学影像分析中的一个重要的研究内容。

医学图像分割是根据医学图像区域内的相似性和区域间的不同特性将图像分割成若干区域，主要包括研究解剖结构、识别感兴趣区域、观察肿瘤生长等。医学图像分割结果直接辅助诊治医师作为，对临床诊断、计算机辅助手术等都有重要意义。大脑又是人体最重要的部分，因此脑部核磁共振图像分割在医学图像分割领域尤为重要。

脑部核磁共振图像分割方法主要是利用图像的灰度特征，主要包括基于像素、基于区域、基于先验形状的方法。

基于像素的方法中主要用到的是基于统计学的方法，它单纯依据图像的灰度信息实现分割，方法灵活，能适应各种形状的医学图像的分割。

基于区域的方法利用图像相邻像素间的局部关系，因此适用于灰度均匀条件下的目标分割，不适用于高噪声的医学图像分割。

基于局部先验形状的方法的分割结果依赖于初始区域的选择，因此该方法对于交互分割比较理想；基于全局先验形状的方法只适合特定形状的图像分割。

现有技术中，关于基于统计模型的脑部核磁共振图像分割方法的代表性文献有：

文献(Sfikas G,Nikou C,Galatsanos N.Robust image segmentation with mixtures of student's t-distributions[C]//Image Processing,2007.ICIP 2007.IEEE International Conference on.IEEE,2007,1:I-273-I-276.)利用EM算法求解有限Student’s t混合模型，进行脑部核磁共振图像分割，然而这种方法明显的缺点是不能自动确定分割数。

文献(da Silva A R F.A Dirichlet process mixture model for brain MRI tissue classification[J].Medical Image Analysis,2007,11(2):169-182.)利用Dirichlet过程这种非参数贝叶斯模型实现了脑部核磁共振图像分割，但他用的是复杂的采样计算，而且他假设数据服从高斯分布，高斯分布对脑部核磁共振图像中的噪声敏感。

发明内容

本发明的目的在于克服现有技术的不足，提供一种针对基于统计模型，有良好抗噪性能的、能实现简单高效的自动分割的基于截断Dirichlet过程无限Student’s t混合模型的脑部核磁共振图像分割方法。

本发明的技术方案如下：

一种基于截断Dirichlet过程无限Student’s t混合模型的脑部核磁共振图像分割方法，基于Dirichlet过程无限Student’s t混合模型，将无限Student’s t混合模型中的分量数假设为预设的图像的分割数K，然后利用期望最大化算法对模型进行求解，再利用贝叶斯最大后验概率准则进行图像分割。

作为优选，包括如下在先步骤：

先在需要进行分割的待处理图像中去除非脑部组织的部分，得到预处理图像；

再利用预处理图像的灰度作为特征，将预处理图像的灰度矩阵拉直变为列向量，作为输入向量。

作为优选，假设输入向量符合Student’s t混合模型，利用期望最大化算法得到最大化模型的参数；再利用最大化模型，根据贝叶斯最大后验概率准则得到输入向量对应的标签向量；将标签向量转化为与待处理图像相同的矩阵形式，得到核磁共振图像的分割结果。

作为优选，假设输入向量符合Student’s t混合模型，以概率1离散的截断Dirichlet过程作为数据的先验分布，表示如下：

${\pi }_i=\left(\begin{array}{cc}{w}_i\underset{m=1}{\overset{i-1}{\Pi }}(1-w_m)& i=2:K\\ w_1& i=1\end{array}\right);$

w_j～Beta(1,α)；

其中，其中，α是模型中的超参数，K是预设的图像的分割数；

则，截断Dirichlet过程的无限Student’s t混合模型表示如下：

$L\left(\Psi \right)=\underset{j=1}{\overset{N}{\Sigma }}\underset{i=1}{\overset{K}{\Sigma }}{\mathrm{log\pi }}_{i}t(x;{\mu }_i,{\Sigma }_i,v_i);$

其中，Ψ表示模型中的所有参数，N是待处理图像中像素总数，x_j表示待处理图像中的第j个像素点，μ_i,Σ_i,ν_i是Student’s>i是先验概率，满足：

π_i≥0；

$\underset{i=1}{\overset{K}{\Sigma }}{\pi }_i=1.$

作为优选，向截断Dirichlet过程的无限Student’s t混合模型中引入两个丢失的信息z、u；

当x_j属于第i个类别时，z_ij＝1，否则z_ij＝0；

当z_ij＝1时，有x_j|u_j,zi_j＝1～Ν(μ_i,Σ_iu_j^-1)，且

则，截断Dirichlet过程的无限Student’s t混合模型的目标函数L(Ψ)为：

$L\left(\Psi \right)=\underset{j=1}{\overset{N}{\Sigma }}\underset{i=1}{\overset{K}{\Sigma }}{z}_{ij}{\mathrm{log\pi }}_{i}t(x_j;{\mu }_i,{\Sigma }_i,v_i)=\underset{j=1}{\overset{N}{\Sigma }}\underset{i=1}{\overset{K}{\Sigma }}{z}_{ij}{\mathrm{log\pi }}_{i}N(x_j;{\mu }_i,{\Sigma }_i{u_j}^{-1})\Gamma (u_j;\frac{v_i}{2},\frac{v_i}{2}).$

作为优选，利用期望最大化算法EM算法的E步对丢失的信息z和u进行估计：

第j个像素点属于类别i的后验概率表示为第j个像素点属于类别i概率占该像素点属于所有可能类别的概率的比例：

$z_{ij}=\frac{{\pi }_{ij}t(x_j;{\mu }_i,{\Sigma }_i,v_i)}{\underset{m=1}{\overset{K}{\Sigma }}{\pi }_{mj}t(x_j;{\mu }_m,{\Sigma }_m,v_m)};$

伽马分布作为u_ij的共轭先验，并利用伽马分布计算得到u_ij的后验，然后根据后验的均值计算得到

作为优选，利用期望最大化算法EM算法的M步对目标函数L(Ψ)求导，计算得到截断Dirichlet过程的无限Student’s t混合模型中的如下参数：

${\overrightarrow{\mu }}_i=\frac{{\Sigma }_{j=1}^N{z}_{ij}{u}_{ij}{\overrightarrow{x}}_j}{{\Sigma }_{j=1}^N{z}_{ij}{u}_{ij}};$

${\Sigma }_i=\frac{{\Sigma }_{j=1}^N{z}_{ij}{u}_{ij}({\overrightarrow{x}}_j-{\overrightarrow{\mu }}_i){({\overrightarrow{x}}_j-{\overrightarrow{\mu }}_i)}^T}{{\Sigma }_{j=1}^N{z}_{ij}{u}_{ij}};$

$w_i=\frac{{\Sigma }_{j=1}^N{z}_{ij}}{{\Sigma }_{j=1}^N{\Sigma }_{m=i+1}^K{z}_{mj}+\alpha -1+{\Sigma }_{j=1}^N{z}_{ij}},i=1:K-1;$

其中，是迭代期望最大化算法的E步中v_i的值。

作为优选，循环迭代期望最大化算法的E步和M步，直到算法收敛。

本发明的有益效果如下：

本发明所述的方法中，假设的Student’s t混合模型直接对应到脑部核磁共振图像不同的部分，而且Student’s t分布的高尾部特性决定了模型有较好的抗噪效果，从而实现脑部核磁共振图像分割。

求解Dirichlet过程无限Student’s t混合模型的过程中，现有技术常常需要用复杂的采样算法或变分推断方法。针对此问题，本发明用简单高效的期望最大化算法对基于截断Dirichlet过程无限Student’s t混合模型进行求解，使模型的求解变得更易实现。

本发明的实施，能够在PC端对一幅脑部核磁共振图像快速自动地进行分割。

附图说明

图1是本发明的方法流程图；

图2是脑部核磁共振图像的预处理图像；

图3是一幅预处理后的脑部核磁共振图像及其参考分割结果图；

图4是截断Dirichlet过程无限Student’s t混合模型示意图，图中，空心圆表示变量，阴影圆表示可观测量，圆角矩形表示参数，右下角表示循环次数，矩形框表示循环迭代；

图5是本发明的分割效果图；

图6是Dirichlet过程无限高斯混合模型方法的分割效果图；

图7是本发明方法与Dirichlet过程无限高斯混合模型方法分割结果对比图，其中给出了IBSRV01数据集中第30张和第50张图像的实验效果；

图8是本发明与截断Dirichlet过程无限高斯混合模型的分割结果对比图，其中给出了IBSRV01数据集中第30张和第50张图像的实验效果。

具体实施方式

以下结合附图及实施例对本发明进行进一步的详细说明。

本发明为了解决现有技术存在的不足，提供一种基于截断Dirichlet过程无限Student’s t混合模型的脑部核磁共振图像分割方法，使得用户可以简单快速地得到脑部核磁共振图像的分割结果。

本发明基于Dirichlet过程无限Student’s t混合模型，考虑到Dirichlet过程无限Student’s t混合模型实际上是可数无限的，因而，将无限Student’s t混合模型中的分量数假设为预设的图像的分割数K(适当的相当大的数K)，然后利用期望最大化算法对模型进行求解，再利用贝叶斯最大后验概率准则进行图像分割。

本发明所述的方法，如图1所示，步骤主要如下：

首先，在需要进行分割的待处理图像中去除非脑部组织的部分，利用简单分水岭算法来去除非脑部组织，得到预处理图像，如图2、图3所示；

其次，利用预处理图像的灰度作为特征，将预处理图像的灰度矩阵拉直变为列向量，作为输入向量；

再次，假设输入向量符合Student’s t混合模型，利用期望最大化算法得到最大化模型的参数；

然后，利用最大化模型，根据贝叶斯最大后验概率准则得到输入向量对应的标签向量；

最后，将标签向量转化为与待处理图像相同的矩阵形式，得到核磁共振图像的分割结果，如图5所示。

假设输入向量符合Student’s t混合模型，本发明以概率1离散的截断Dirichlet过程作为数据的先验分布，表示如下：

${\pi }_i=\left(\begin{array}{cc}{w}_i\underset{m=1}{\overset{i-1}{\Pi }}(1-w_m)& i=2:K\\ w_1& i=1\end{array}\right);$

w_j～Beta(1,α)；

其中，α是模型中的超参数，K是假定的适当的相当大的图像的分割数，例如假定的适当的相当大的图像的分割数；

则，截断Dirichlet过程的无限Student’s t混合模型表示如下：

$L\left(\Psi \right)=\underset{j=1}{\overset{N}{\Sigma }}\underset{i=1}{\overset{K}{\Sigma }}{\mathrm{log\pi }}_{i}t(x_j;{\mu }_i,{\Sigma }_i,v_i);$

其中，Ψ表示模型中的所有参数，N是待处理图像中像素总数，x_j表示待处理图像中的第j个像素点，μ_i,Σ_i,ν_i是Student’s>i是先验概率，满足：

π_i≥0；

$\underset{i=1}{\overset{K}{\Sigma }}{\pi }_i=1.$

为了用简单常用的期望最大化算法EM算法解得以上模型，则向截断Dirichlet过程的无限Student’s t混合模型中引入两个丢失的信息z、u；

当x_j属于第i个类别时，z_ij＝1，否则z_ij＝0；

当z_ij＝1时，有x_j|u_j,z_ij＝1～Ν(μ_i,Σ_iu_j^-1)，且

则，完整数据下的截断Dirichlet过程的无限Student’s t混合模型的目标函数L(Ψ)为：

$L\left(\Psi \right)=\underset{j=1}{\overset{N}{\Sigma }}\underset{i=1}{\overset{K}{\Sigma }}{z}_{ij}{\mathrm{log\pi }}_{i}t(x_j,·{\mu }_i,{\Sigma }_i,v_i)=\underset{j=1}{\overset{N}{\Sigma }}\underset{i=1}{\overset{K}{\Sigma }}{z}_{ij}{\mathrm{log\pi }}_{i}N(x_j;{\mu }_i,{\Sigma }_i{u_j}^{-1})\Gamma (u_j;\frac{v_i}{2},\frac{v_i}{2}).$

利用期望最大化算法EM算法的E步对丢失的信息z和u进行估计：

第j个像素点属于类别i的后验概率表示为第j个像素点属于类别i概率占该像素点属于所有可能类别的概率的比例：

$z_{ij}=\frac{{\pi }_{ij}t(x_j;{\mu }_i,{\Sigma }_i,v_i)}{\underset{m=1}{\overset{K}{\Sigma }}{\pi }_{mj}t(x_j;{\mu }_m,{\Sigma }_m,v_m)};$

伽马分布作为u_ij的共轭先验，并利用伽马分布计算得到u_ij的后验，然后根据后验的均值计算得到

因为要最大化截断Dirichlet过程的无限Student’s t混合模型的目标函数L(Ψ)，则利用期望最大化算法的M步对目标函数L(Ψ)求导，计算得到截断Dirichlet过程的无限Student’s t混合模型中的如下参数：

${\overrightarrow{\mu }}_i=\frac{{\Sigma }_{j=1}^N{z}_{ij}{u}_{ij}{\overrightarrow{x}}_j}{{\Sigma }_{j=1}^N{z}_{ij}{u}_{ij}};$

${\Sigma }_i=\frac{{\Sigma }_{j=1}^N{z}_{ij}{u}_{ij}({\overrightarrow{x}}_j-{\overrightarrow{\mu }}_i){({\overrightarrow{x}}_j-{\overrightarrow{\mu }}_i)}^T}{{\Sigma }_{j=1}^N{z}_{ij}{u}_{ij}};$

$w_i=\frac{{\Sigma }_{j=1}^N{z}_{ij}}{{\Sigma }_{j=1}^N{\Sigma }_{m=i+1}^K{z}_{mj}+\alpha -1+{\Sigma }_{j=1}^N{z}_{ij}},i=1:K-1;$

其中，是迭代期望最大化算法的E步中v_i的值。

进而，得到截断Dirichlet过程无限Student’s t混合模型，如图4所示。

循环迭代期望最大化算法的E步和M步，直到算法收敛。

经过以上步骤得到模型的参数后，根据贝叶斯最大后验概率准则得到输入向量数据对应的标签向量，第j个像素点根据得到它的标签i。

将得到的标签向量转化为待处理图像矩阵形式，得到核磁共振图像的分割结果。

如图6、图7所示，现有技术中的Dirichlet过程无限高斯混合模型方法，虽然能自动得正确得到图像分割的数目，但是由于高斯函数本身的轻尾部特性导致分割结果受图像中噪声的影响较大。而且这种方法由于采用了复杂的采样算法计算模型，耗时一般在十几秒。而本发明的方法在不到1秒的时间就能自动得到图像的分割结果，如表1所示。

表1：本发明与Dirichlet过程无限高斯混合模型方法分割准确率和分割时间的对比

在IBSRV01数据集的128张脑MRI图像上做实验，如图8所示，得到平均分割准确率和分割时间，如表2所示。

表2：本发明与截断无限高斯混合模型方法分割准确率的对比

上述实施例仅是用来说明本发明，而并非用作对本发明的限定。只要是依据本发明的技术实质，对上述实施例进行变化、变型等都将落在本发明的权利要求的范围内。