CN103729843A - 基于马尔可夫的医学图像分割方法 - Google Patents

Abstract

本发明提供一种基于马尔可夫的医学图像分割方法，包括步骤：步骤1、首先利用马尔可夫分割方法对原始图像进行初始分割；步骤2、利用水平集算法对图像进行处理，更新水平集函数、偏移场以及各区域灰度均值；步骤3、用原始图像除以当前更新的偏移场得到矫正后的图像；步骤4、利用马尔可夫分割算法对矫正后的图像进行分割；步骤5、判断当前水平集算法与马尔可夫分割算法的迭代次数之和是否大于预设值，如是，最终得到的矫正后的图像进行分割得到最终分割结果；如否，返回步骤2。实现了快速、稳健、准确地处理带有灰度不均匀的医学图像分割。

Description

基于马尔可夫的医学图像分割方法

技术领域

本发明属于图像处理技术领域，主要涉及马尔可夫及水平集技术中的图像分割方法。

背景技术

近些年来，医学影像技术飞速发展，其在临床中的应用也越来越广泛，计算机断层成像（CT）、B超、核磁共振成像（MRI）和电子内窥镜等多种现代成像技术已作为常规的辅助方法应用于临床诊断当中，这就使得医生可以更清晰，更直接地对病人组织器官观察，诊断，大大提高了诊断准确率，缩短了确诊时间。而且这些辅助手段的作用与影响也不断的变大，正带动着现代医学诊断的变革。目前更高维的，更高分辨率的的成像技术正在成为人们研究的热点。

由于医学成像技术的复杂性和人体组织器官的物质结构的特殊性使们医学图像不可避免的带有偏移场引起的灰度不均匀，器官运动产生的伪影与模糊以及人本内部其他物质带来的干扰影像等，其中影响最大的是偏移场。正是医学图像的这些特殊性使得大量的传统分割算法不再适用，因此，近些年来人们对于医学图像的分割算法的研究往往是在传统的经典分割算法的基础上充分考虑偏移场的影响来入手的。

医学图像最大的特点就是灰度不均匀。近年来，人们针对这一问题通过改进，融合，创新等方式在传统算法的基础上提出了大量非常实用的分割算法。在研究了医学图像中偏移场的特点后，人们提出的分割方法可以分为两大类：一类是先矫正灰度再分割的方式；另一类是灰度矫正和分割同时进行的方式。

对于前一类，通常是在矫正后再利用较为经典的分割方法，如模式识别中的各种聚类算法等，来进行分割，因此其关键在于灰度矫正算法。这类方法简单，容易实现，而且可以充分利用各种经典的传统方法。目前效果较好的，较实用的灰度矫正方法有N3算法，该方法就是对偏移场模型中的参数进行估计进而得到偏移场。详见文献：J.G.Sled,A.P.Zijdenbos,A.C.Evans.“Nonparametric method for automatic correction of intensity nonuniformity in MRIdata[J]”,IEEE Transactions on Medical Imaging,Vol.17,No.1,pp.[d]87-97,Feb.1998.

对于后一类，一般是在分割的过程中进行矫正，通常要在原有算法模型中引入对偏移场的估计部分。目前较为经典的方法有Li等人提出基于局部区域的灰度不均匀场水平集分割方法。该方法将图像看作是偏移场和真实图像的乘积，然后将该乘积引入到传统水平集C-V模型之中，然后在局部区域内应用该模型并扩展到全局构成新的能量函数，最后在最小化新的能量函数的过程中同时完成对偏移场的估计与分割。详见文献：C.M.Li,R.Huang,Z.H.Ding,etc.“A level set method for image segmentation in the presence of intensity inhomogeneities withapplication to MRI[J]”,IEEE Transactions on Image Processing,Vol.20,No.7,pp.2007-2016,July.2011.。水平集分割算法充分考虑到医学图像的灰度不均匀性，可以在分割的同时对医学图像进行灰度矫正，但其只利用了灰度信息并没有利用空间信息，而且处理分割区域差异很小（分割区域的灰度十分相近或对比度很低）的医学图像时对初始值的敏感性增加，分割准确性降低。另外，该算法是基于局部区域的，每个像素点要被计算多次，计算开销加大。

发明内容

本发明所要解决的技术问题是，提供一种快速、稳健、准确的医学图像分割方法。

本发明为解决上述技术问题所采用的技术方案是，基于马尔可夫的医学图像分割方法，包括以下步骤：

步骤1、首先利用马尔可夫分割方法对原始图像进行初始分割，利用初始分割结果对水平集算法进行初始化，得到初始化的水平集函数、偏移场以及各区域灰度均值；

步骤2、利用水平集算法对图像进行处理，更新水平集函数、偏移场以及各区域灰度均值；当第一次执行本步骤时，水平集算法处理的图像为原始图像，否则，水平集算法处理的图像为矫正图像；

步骤3、用原始图像除以当前更新的偏移场得到矫正后的图像；

步骤4、利用马尔可夫分割算法对矫正后的图像进行分割；

步骤5、判断当前水平集算法与马尔可夫分割算法的迭代次数之和是否大于预设值，如是，最终得到的矫正后的图像进行分割得到最终分割结果；如否，返回步骤2。

传统马尔可夫分割方法是一种常用的图像分割方法，其可以很好的利用空间统计信息。但由于它没有充分利用灰度信息，不适用于灰度不均匀的医学图像。在传统马尔可夫模型的基础上融入水平集模型，分割任务主要由马尔可夫模型来完成，水平集模型仅是对马尔可夫分割进行辅助修正。马尔可夫分割算法充分利用了空间信息，将水平集算法融入到马尔可夫算法中后，马尔可夫分割算法可提供一个接近准确结果的初值减少算法的迭代次数同时增强算法的鲁棒性，算法同时充分利用了空间统计信息和灰度信息，准确性也大大提高了。本发明将水平集算法的去偏移场矫正能力融入其中，使我们的算法在马尔可夫算法的基础上具有了分割灰度不均匀的医学图像的能力。

本发明的有益效果是，实现了快速、稳健、准确地处理带有灰度不均匀的医学图像分割；相对传统马尔可夫算法，能够对灰度不均匀的医学图像进行较准确地分割，相对水平集算法而言，对初值更不敏感，鲁棒性更好，准确性更高。

具体实施方式

为了方便地描述本发明内容，首先对以下一些现有技术进行简单介绍：

定义1：一幅M×N的图像记为Y={y_s|s∈S}，S={s=(i,j)|1≤i≤M,1≤j≤N}，y_s表示像素s的像素值，S为所有图像中像素点集合，图像的分割结果记为X={x_s|s∈S}，x_s表示图像像素点的所属类别，x_s取值范围记为L={1,2,...,N}，L代表图像像素点的所属类别。

定义2：8邻域系统及其势团能量。令δ(s)是像素点s的邻域，它是以位置s为中心，r为半径的一个圆形区域：δ(s)={s₁∈S|dist(s,s₁)≤r²,s≠s₁}，其中dist(s,s₁)代表s和s₁两点的欧式距离。依据此定义邻域系统：δ={δ(s)|s∈S}，它满足以下三个条件：1）

2）

箭头表示的是等价，s、t是一幅图像中的两个点3）故对于一个像素点s，其8邻域系统势团能量如下：

$\left\{\begin{array}{c}U\left(y_s\right)=\underset{s\∈S}{\mathrm{\Σ}}\underset{p\∈\mathrm{\δ}\left(s\right)}{\mathrm{\Σ}}V(y_s,y_p)\\ V(y_s,y_p)=\begin{array}{}\end{array}\left\{\begin{array}{c}-w,y_s=y_p\\ w,y_s\≠y_p\end{array}\right.\end{array}\right.;---\left(1\right)$

其中，y_s表示s点的像素值，y_p表示p点的像素值，p∈S；w为系数，可取为1，U(y_s)表示s点的8邻域系统势团能量。

定义3：Hamersley-Clifford定理：证明了Gibbs和马尔可夫模型MRF的等价关系，因此可以利用Gibbs分布计算出MRF中的先验概率P(X)。

定义4：水平集函数。设水平集函数φ(x)，x∈Ω，它满足：在封闭曲线C内部，φ(x)>0；在封闭曲线C外部，φ(x)<0；在曲线C上面，φ(x)=0。可用φ(x)的零水平集来表示曲线C。

定义5：Heaviside函数和狄拉克函数。Heaviside是一个单位阶跃函数，狄拉克函数是它的导数，即单位脉冲函数。由于阶跃函数不容易实现，因此在实际应用当中，我们使用式(2)和式(3)的平滑函数分别近似Heaviside函数H(x)和狄拉克函数δ(x)：

$H\left(x\right)=\frac{1}{2}[1+\frac{2}{\mathrm{\&pi;}}\arctan \left(\frac{x}{\mathrm{\&epsiv;}}\right)];---\left(2\right)$

$\mathrm{\&delta;}\left(x\right)=H^{\&prime;}\left(x\right)=\frac{1}{\mathrm{\&pi;}}\frac{\mathrm{\&epsiv;}}{{\mathrm{\&epsiv;}}^2+x^2};---\left(3\right)$

其中ε是参数，一般取1。

定义6：梯度下降法。就是利用负梯度方向来决定每次迭代的新的搜索方向，使得每次迭代能使待优化的目标函数逐步减小。

定义7：医学图像模型。由于受偏移场影响，医学图像灰度不均匀，因此，我们将医学图像I(x)看作是偏移场B(x)和真实图像（或称为本征图像）J(x)的乘积，如下：

I(x)=B(x)·J(x)；      (4)

定义8：Jaccard指标。国际上通用的可以用来衡量分割准确性的指标，其定义如下：

$\mathrm{Jaccard}(X,Y)=\frac{|X\&cap;Y|}{|x\&cup;|};---\left(5\right)$

其中X和Y分别表示真实结果和分割结果。Jaccard指标应该在[0 1]区间内，其值越大分割结果就越准确。

（1）传统马尔可夫模型

其实质就是将图像分割问题转化为求像素点的最大后验概率P(X|Y)，即

$\hat{x}=\arg \max \left\{P\left(X\right|Y)\right\};---\left(6\right)$

其中

表示最优分割结果。再根据贝叶斯公式：P(X|Y)∝P(Y|X)P(X)可将上式转化为如下：

$\hat{x}=\arg \max \left\{P\left(Y\right|X)P\left(X\right)\right\};---\left(7\right)$

P(X)是MRF的先验概率，可根据Hamersley-Clifford定理可得：

$P\left(X{=x}_s\right)=\frac{1}{Z}\exp [-U(X=x_s)];---\left(8\right)$

其中 $Z=\underset{x_s}{\mathrm{\&Sigma;}}\exp [-U\left(x_s\right)],$ 是归一化系数。

P(Y|X)是条件概率，其满足高斯分布：

$P\left(y_s\right|x_s)=\frac{1}{\sqrt{2\mathrm{\&pi;}{{\mathrm{\&sigma;}}_l}^2}}\exp [-\frac{{(y_s-{\mathrm{\&mu;}}_l)}^2}{2{{\mathrm{\&sigma;}}_l}^2}];---\left(9\right)$

其中μ_l，σ_l分别是第l类的均值和标准差，表达式如下：

$\left\{\begin{array}{c}{\mathrm{\&mu;}}_l=\frac{1}{\left|S_l\right|}\underset{s\&Element;S_l}{\mathrm{\&Sigma;}}{y}_s\\ {{\mathrm{\&sigma;}}_l}^2=\frac{1}{\left|S_l\right|}\underset{s\&Element;S_l}{\mathrm{\&Sigma;}}{(y_s-{\mathrm{\&mu;}}_l)}^2\end{array}\right.---\left(10\right)$

其中S_l是属于第l类的像素点所构成的区域，它是一个归一化值。

我们假设随机变量y_s间相互独立，可得：

$\begin{array}{c}\hat{X}=\arg \max \left[\log (P\left(Y\right|X)\&CenterDot;P\left(x\right))\right]\\ =\arg \min [U\left(Y\right|X)+U\left(X\right)]\end{array};---\left(11\right)$

$U\left(Y\right|X)=\underset{S}{\mathrm{\&Sigma;}}[\ln \left(\sqrt{2\mathrm{\&pi;}{{\mathrm{\&sigma;}}_l}^2}\right)+\frac{{(y_s-{\mathrm{\&mu;}}_l)}^2}{2{{\mathrm{\&sigma;}}_l}^2}];---\left(12\right)$

显然，根据最大后验概率MAP准则可知，我们的最优分割结果为

这便是马尔可夫分割模型，我们将采用条件迭代算数（ICM）来求得最分割结果。对每个像素点y_s执行如下三步，直到所得结果变化很小为止。

1）根据式(10)更新参数μ_l和σ_l；

2）根据式(1)和式(12)计算U(x_s)和U(y_s|x_s)；

3）更新第l类的像素点以满足式(11)最小。

传统马尔可夫分割方法可以很好的利用空间统计信息，但没有充分利用灰度信息，同时也不适用于灰度不均匀的医学图像。因此可以融入灰度矫正来解这一问题。

（2）水平集模型

这里使用的是一种基于局部区域的灰度不均匀场水平集模型。其主要思想是最小化如下能量函数：

$\left\{\begin{array}{c}E(\mathrm{\&Phi;},B,C)={\&Integral;}_{\mathrm{\&Omega;}}(\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}{e}_i\&CenterDot;M_i\left(\mathrm{\&Phi;}\left(x\right)\right))\mathrm{dx}+\mathrm{\&mu;L}\left(\mathrm{\&Phi;}\left(x\right)\right)+\mathrm{vP}\left(\mathrm{\&Phi;}\left(x\right)\right)\\ e_i={\&Integral;}_{\mathrm{\&Omega;}}K(y-x){|I\left(x\right)-B\left(y\right)c_i|}^2\mathrm{dy}\end{array}\right.;---\left(13\right)$

其中I(x)，Ω→R是原始图；Φ(x)=(φ₁(x),φ₂(x))是两个水平集函数的向量形式；C=(c₁,c₂,c₃)是三个区域的灰度均值。

是长度约束项；是距离符号函数约束项。K(x)是高斯核函数，如下：

$K\left(x\right)=\left\{\begin{array}{cc}\frac{1}{a}\&CenterDot;e^{-{\left|x\right|}^2/2{\mathrm{\&sigma;}}^2},& \left|x\right|\&le;r\\ 0,& \left|x\right|>r\end{array}\right.;---\left(14\right)$

其中σ是高斯函数中的标准差，r是某一像素点的领域的半径。M_i(Φ(x))是像素点归属函数，如下：

$\left\{\begin{array}{c}{M}_1\left(\mathrm{\&Phi;}\left(x\right)\right)=H\left({\mathrm{\&phi;}}_1\left(x\right)\right)\&CenterDot;H\left({\mathrm{\&phi;}}_2\left(x\right)\right)\\ M_2\left(\mathrm{\&Phi;}\left(x\right)\right)=H\left({\mathrm{\&phi;}}_1\left(x\right)\right)\&CenterDot;(1-H\left({\mathrm{\&phi;}}_2\left(x\right)\right))\\ M_3\left(\mathrm{\&Phi;}\left(x\right)\right)=1-H\left({\mathrm{\&phi;}}_1\left(x\right)\right)\end{array}\right.;---\left(15\right)$

显然，当该能量函数最小时，既能驱使轮廓线向目标边界移动，又能保持曲线的光滑，我们可以得到一个分割结果和偏移场的估计，进而我们使用原始图像除以偏移场就可以得到矫正后的图像。

能量函数中存在三个变量：Φ，B，C。我们的最小化能量函数的方法是：每个迭代过程都固定两个变量，然后对第三个变量求偏导，最后得到这个变量的表达式，这样不停的交替进行，直到能量函数达到最小值或不再变化为止。C，B的表达式和Φ的演化方程如下：

$c_i=\frac{{\&Integral;}_{\mathrm{\&Omega;}}(B\left(x\right)*K\left(x\right))\&CenterDot;I\left(x\right)\&CenterDot;M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}}{{\&Integral;}_{\mathrm{\&Omega;}}(B^2\left(x\right)*K\left(x\right))\&CenterDot;M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}};---\left(16\right)$

$\left\{\begin{array}{c}B\left(x\right)=\frac{(I\left(x\right)\&CenterDot;J\left(x\right))*K\left(x\right)}{J^2\left(x\right)*K\left(x\right)}\\ J\left(x\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{c}_i\&CenterDot;M_i\left(\mathrm{\&phi;}\left(x\right)\right)\end{array}\right.;---\left(17\right)$

$\left\{\begin{array}{c}\frac{\&PartialD;{\mathrm{\&phi;}}_m\left(x\right)}{\&PartialD;t}=-\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}\frac{\&PartialD;M_i\left({\mathrm{\&phi;}}_m\left(x\right)\right)}{\&PartialD;{\mathrm{\&phi;}}_m\left(x\right)}\&CenterDot;e_i+\mathrm{\&mu;}\&CenterDot;\mathrm{\&delta;}\left({\mathrm{\&phi;}}_m\left(x\right)\right)\&CenterDot;\mathrm{div}\left(\frac{\&dtri;{\mathrm{\&phi;}}_m\left(x\right)}{|\&dtri;{\mathrm{\&phi;}}_m\left(x\right)|}\right)\\ +v\&CenterDot;\mathrm{div}(\mathrm{dp}\left(\right|\&dtri;{\mathrm{\&phi;}}_m\left(x\right)\left|\right)\&CenterDot;\&dtri;{\mathrm{\&phi;}}_m\left(x\right))\\ \mathrm{dp}\left(z\right)=\frac{d}{\mathrm{dz}}P\left(z\right)\end{array}\right.;---\left(18\right)$

其中“*”表示卷积运算，i=1，2，3。

上面的水平集分割算法充分考虑到医学图像的灰度不均匀性，可以在分割的同时对医学图像进行灰度矫正，但其只利用了灰度信息并没有利用空间信息，而且处理分割区域差异很小（分割区域的灰度十分相近或对比度很低）的医学图像时对初始值的敏感性增加，分割准确性降低。另外，该算法是基于局部区域的，每个像素点要被计算多次，计算开销加大。

本发明的分割算法包括以下步骤：

步骤1.初始化参数Φ（水平集函数），C（各区灰度均值），B（偏移场）。

首先利用马尔可夫分割方法对图像进行分割，用它的分割结果对水平集进行初始化，由于此结果已接近正确分割结果，这就大大减少水平集的迭代次数。初始化只需要对马尔可夫分割的三个区域赋三个不同值即可，这样便可得到水平集的初始值Φ₀，然后利用统计的方法求得C的初值C₀。最后我们利用式(17)计算得到B的初始值B₀。设置算法总的迭代次数为K，并将其初始化为1。

步骤2.利用水平集算法对图像进行处理，更新水平集函数Φ，偏移场B，各区域灰度均值C；当第一次执行本步骤时，水平集算法处理的图像为原始图像，否则，水平集算法处理的图像为矫正图像；

这里将该步骤中水平集的迭代次数设置为10，同时这里的迭代次数计入算法总的迭代次数K中，也即执行10次水平集迭代就结束本步骤并跳至步骤3。

步骤3.用原始图像I(x)除以步骤2中求得的偏移场B，得到矫正后的图像，然后对该图像作归一化处理，使其像素值在[0255]间；

步骤4.利用马尔可夫分割算法对矫正后的图像进行分割；

将本步骤中的马尔可夫算法的迭代次数设置为8，同时这里的迭代次数也计入算法总的迭代次数K中，即执行8次马尔可夫算法的迭代就结束本步骤并跳至步骤5；

步骤5.判断总的迭代次数K大于60，如是，最终得到的矫正后的图像进行分割得到最终分割结果；如否，返回步骤2。

步骤2中在更新各个参数时，我们选择最后更新B，这样可以更好的保证B的准确性，使得我们对马尔可夫结果的修正更好。

在整个算法中，由于我们使用的是两种算法相互修正，因此为了减少开销，我们将每种方法的迭代次数设置为一个较小的值。尤其上面的水平集算法采用了局部区域的思想，每个像素点会被重复计算多次，计算开销加大，因此，将水平集算法的部分的迭代次数仅设置为10，这样可以大大提高算法的速度。

Claims (1)

1.基于马尔可夫的医学图像分割方法，其特征在于，包括以下步骤：

步骤1、首先利用马尔可夫分割方法对原始图像进行初始分割，利用初始分割结果对水平集算法进行初始化，得到初始化的水平集函数、偏移场以及各区域灰度均值；

步骤2、利用水平集算法对图像进行处理，更新水平集函数、偏移场以及各区域灰度均值；当第一次执行本步骤时，水平集算法处理的图像为原始图像，否则，水平集算法处理的图像为矫正图像；

步骤3、用原始图像除以当前更新的偏移场得到矫正后的图像；

步骤4、利用马尔可夫分割算法对矫正后的图像进行分割；

步骤5、判断当前水平集算法与马尔可夫分割算法的迭代次数之和是否大于预设值，如是，最终得到的矫正后的图像进行分割得到最终分割结果；如否，返回步骤2。