CN1471034A - 基于水平集和分水岭方法的医学图像分割方法 - Google Patents

Abstract

一种基于水平集和Watershed方法的医学图像分割方法，包括步骤：各向异性扩散滤波，用于去除噪声；采用Watershed方法，对图像进行过度分割；堆数据结构的建立，用于定位窄带中具有最小时间T的网格点；采用Fast Marching方法，对图像进行最后的分割。本发明利用Watershed方法与改进的Fast Marching方法进行医学图像分割，不仅能够大大地提高分割的速度，同时还具有广泛的适应性，无论是CT还是MR影像。在计算机辅助诊断和治疗等领域有着重要的应用价值。

Description

基于水平集和分水岭方法的医学图像分割方法

技术领域

本发明涉及图像处理，特别涉及结合Watershed(分水岭)方法和水平集来进行医学图像分割的方法。

背景技术

在对图像的研究和应用中，人们往往只对图像中的某些部分感兴趣。图像分割就是利用图像信息中部分特征将图像中若干感兴趣的区域提取出来，感兴趣的区域通常称为目标或前景(其余部分称为背景)，在此基础上才能对目标进行特征提取和测量，从而使更高层的分析和理解成为可能。图像分割可用的特征包括图像灰度、颜色、纹理、局部统计特征或频谱特征等。

图像分割是由图像处理到图像分析的关键步骤，也是计算机视觉领域低层次视觉中的主要问题，同时它又是一个经典难题。由于此问题的重要性和困难性，多年来图像分割问题一直吸引着国内外研究人员为之努力，至今已提出了上千个各种类型的分割算法。但到由于图像分割只能利用图像信息中的部分特征分割区域，因此各种方法必然带有局限性和针对性，只能针对各种实际应用领域的需求来选择各种方法，到目前为止还不存在一个通用的分割方法。

水平集方法主要是从界面传播等研究领域中逐步发展起来的，它是处理封闭运动界面随时间演化过程中几何拓扑变化的有效的计算工具。Osher和Sethian首先提出依赖时间的运动界面的水平集描述。水平集方法的基本思想是将平面闭合曲线隐含的表达为二维曲面函数的水平集，即具有相同函数值的点集，通过水平集函数曲面的进化隐含的求解曲线的运动。尽管这种转化使得问题在形式上变得复杂，但在问题的求解上带来很多优点，其最大的优点在于曲线的拓扑变化能够得到很自然的处理，而且可以获得唯一的满足熵条件的弱解。

水平集函数的演化满足如下的基本方程：

                  Φ_t+F|Φ|＝0               (1)

其中Φ为水平集函数，其零水平集表示目标轮廓曲线，即：

                  Γ(t)＝{x|Φ(x，t)＝0}        (2)

|Φ|表示水平集函数的梯度范数；F为曲面法线方向上的速度函数，控制曲线的运动，一般F包括两项：与图像有关的项(如梯度信息)以及与曲线的几何形状有关的项(如曲线的曲率)。

假定Γ是一个超平面的初始位置，F是在法线方向上的速度。在Level Set算法中，我们视Γ是一个高维函数Φ(x，y，z)的零水平集。然后，根据链式规则，可以产生如式(1)的运动的超平面的演化方程。

水平集算法中最普通的一种就是所谓的Fast Marching(快速匹配)方法。现在考虑一种界面运动的特殊情况，那就是界面的运动速度F＞0。假定T是界面经过一个指定点(x，y)的时间，这样T就满足如下的方程：

                  |T|F＝1                     (3)

这个公式简单的说明了到达时间的梯度和界面的运动速度成反比。从广义上说，有两种方法可以用来近似运动的界面随时间变化的位置：一种是通过迭代和数字近似公式(1)中微商来解决；另一种是构建公式(3)中到达时间T的解决方案。而Fast Marching方法依赖于后一种方法。

公式(3)是著名的Eikonal方程的一种形式，Sethian首先提出，要得到公式(3)中的到达时间T，等价于求解下面的二次方程。 $[{\max {(D_{\mathrm{ij}}^{-x}T,0)}^2+\min {(D_{\mathrm{ij}}^{+x}T,0)}^2+\max {(D_{\mathrm{ij}}^{-y}T,0)}^2+{\min (D_{\mathrm{ij}}^{+y}T,0)}^2]}^{1/2}=1/F_{i,j}----\left(4\right)$

这里D^-和D⁺分别是后向差分和前向差分算子。

下面我们给出Sethian首先提出的Fast Marching方法的具体步骤：

初始化：

1)活动点：活动点就是所有网格点中时间T固定的点。在我们算法中，也就是用户指定的种子点，时间T(x，y)＝0，见图2，以二维情况为例；

2)窄带点：所有在窄带中的点叫做窄带点。在我们算法中，也就是所有种子点的4-邻接的点，时间T(x，y)＝1/F(x，y)，见图2，以二维情况为例；

3)远离点：除了活动点和窄带点外，所有其它的网格点为远离点，T(x，y)＝TIME_MAX，见图2，以二维情况为例；

图2中的红色格子表示活动点，绿色格子表示窄带点，而白色的格子表示远离点。图2(a)表示Fast Marching方法的开始，图2(b)和图2(c)表示Fast Marching方法的过程。

循环：

1)开始循环：设点(i，j)是窄带中具有最小时间T的点；

2)标识点(i，j)为活动点，并从窄带中删除；

3)标识点(i，j)的4-邻接点：如果点(i，j)的邻接点为活动点，则不改变时间；如果其邻接点为窄带点，则按照公式(4)更新邻接点的时间；如果其邻接点为远离点，则标识该邻接点为窄带点，同时按照公式(4)更新该邻接点的时间；

4)如果一点的到达时间超过指定的阈值，则结束循环，否则跳到1)；

从Fast Marching方法的具体步骤我们很容易知道，在最坏的情况下，Fast Marching方法需要遍历一幅图像的所有像素。R.Malladi首先定义了如下的速度函数： $F(x,y)=e^{-\mathrm{\α}|\▿G_{\mathrm{\σ}}*I(x,y)|}----\left(5\right)$

这里I(x，y)为原始图像，是梯度算子，G_σ ^*I(x，y)表示对原始图像进行高斯平滑，而0＜α＜1是权重系数。

发明内容

本发明的目的是提供一种快速、准确的医学图像分割算法，能够对CT和MR等模态的图像进行分割。

为实现上述目的，基于水平集和Watershed方法的医学图像分割方法包括步骤：

各向异性扩散滤波，用于去除噪声；

采用Watershed方法，对图像进行过度分割；

堆数据结构的建立，用于定位窄带中具有最小时间T的网格点；

采用Fast Marching方法，对图像进行最后的分割。

本发明利用Watershed方法与改进的Fast Marching方法进行医学图像分割，不仅能够大大地提高分割的速度，同时还具有广泛的适应性，无论是CT还是MR影像。在计算机辅助诊断和治疗等领域有着重要的应用价值。

附图说明

图1是利用Watershed方法与改进Fast marching方法进行图像分割的结构框图；

图2是活动、窄带以及远离点示意图；

图3是Fast Marching方法跨越图像边界示意图；

图4是Watershed变换后相邻区域示意图；

图5是膝关节CT图像分割；

图6是脑部MR图像：脑回分割；

图7是脑部MR图像：脑白质分割；

图8是脑部MR图像：脑白质分割。

具体实施方式

本发明的核心思想是将Watershed算法与改进的Fast Marching方法结合起来进行图像分割。由于引入了Watershed算法对图像进行过度分割，该方法只需计算种子点到其所在的区域边界的到达时间，对于区域内部的其它像素点，种子点的到达时间则无需计算，这样算法的运算速度就得到大大的改善。而且，我们根据区域之间的统计特性的相似度重新定义了FastMarching方法的速度函数。

下面结合附图详细描述本发明的医学图像分割算法。作为一种具体的实现方案，结构框图见图1。主要包括四个步骤：各向异性扩散滤波，Watershed方法进行过度分割(分割的区域数目超过图像中包含的实际对象数)，堆数据结构的建立，改进的Fast Marching方法进行最后的分割。下面对其逐一介绍。步骤一：各向异性扩散滤波

医学图像通常含有随机噪声，噪声有可能造成边缘检测错误，为减少随机噪声对分割算法的影响，我们可以先对图像进行滤波。但是，利用低通滤波算子平滑图像噪声的同时通常也会平滑掉边缘，这样将导致边界信息被破坏，从而给图像的后处理带来困难。为避免此情况的发生，在我们的算法中将采用各向异性扩散滤波来去除噪声，该方法在去除噪声的同时还能保持边缘位置和增强边缘。

各向异性扩散滤波最早是由Perona和Malik提出的，他们在传统的各向同性扩散方程基础上提出了下面的各向异性扩散方程，并将其用于图像的边缘增强。 $\frac{\∂I(x,y,t)}{\∂t}=\mathrm{div}[g\left(\right||\▿I|\left|\right)\▿I]----\left(6\right)$

其中|I|表示图像梯度幅值，扩散系数是一个与梯度有关的函数g(·)。g(·)是一个用来控制平滑程度的单调递减函数，满足g(0)＝1，并且在x→∞时g(x)→0。这意味着扩散系数随图像梯度的增大而减小，这样可以保证在区域内部(I小)以较快速度扩散，而在边缘点(I大)则不再扩散，从而起到边缘增强的作用，因此g(·)也称为边缘停止函数。步骤二：Watershed方法进行过度分割

Watershed(分水岭)算法是一种基于数学形态学的方法，最早由Digabel和Lantuejou提出。Watershed算法也可以分类为基于区域的分割方法。Watershed方法的思想来源于地理学：它将梯度幅值图像看成一幅地形图，而梯度幅值对应海拔高度，图像中不同梯度值的区域就对应于山峰和山谷间盆地。设想在各个局部极小值点的位置打一个洞，然后将地形图逐渐浸入一个湖中，全局极小值点的盆地先进水。水位逐渐升高漫过盆地，当相邻两个盆地的水即将合并时，这时在两个盆地间建坝拦截。此过程将图像划分为许多个山谷盆地，分水岭就是分隔这些盆地的堤坝。Watershed算法的基本步骤如下：1)非线性扩散滤波平滑图像；2)计算平滑后图像的梯度图；3)计算梯度图的最大值MAX；4)选择合适的水淹阈值DT＝MAX*l/n；5)对每一个像素，如果当前的像素和其8-邻居的值都小于DT，则合并当前像素和其8-邻居6)合并过小区域(如：区域包含的点数＜10).

分水岭算法将产生原始图像的过度分割图(分割的区域数目超过图像中包含的实际对象数)，过度分割的区域数目取决于淹没水平面高度的选择。步骤三：堆数据结构的建立

Fast Marching方法的关键之处在于如何快速的定位窄带中具有最小时间T的网格点。这里我们使用下面的堆数据结构：

struct

{

     double time；种子点到网格点的到达时间

     int x；    网格点的X方向坐标

     int y；    网格点的Y方向坐标

}HeapNode；

我们存储种子点到其它网格点的到达时间T，同时保存该网格点在网格结构中的位置(X方向的坐标和Y方向的坐标)。下面是一些对堆数据结构的操作：

heap_init：初始化堆数据结构，将种子点作为活动点，计算种子点到边缘点的到达时间，并按从小到大地顺序存入堆数据结构中；

heap_push：将远离点转化为窄带点，并按从小到大地顺序存入堆数据结构中；

heap_pop：从堆中删除具有最小时间的网格点，并将其改为活动点；

heap_update：更新种子点到窄带中网格点的到达时间。步骤四：改进的Fast Marching方法进行最后的分割

从公式(5)可知，速度函数F(x，y)仅仅依靠图像中的边界信息，也就是图像的梯度信息而没有充分的利用图像区域的全局信息。因此，对于图像中的边界模糊或者离散状，就很难获得理想的分割效果。如图3所示：

图3(a)表示膝盖图像原始图，绿色箭头所指的分别是种子点选择和边界模糊的地方；图3(b)表示在图像边界模糊的地方分割结果跨越了边界，这是和速度函数的定义以及参数α的选择息息相关的，Fast Marching方法的关键在于速度函数的选取。由于其速度函数选取了图像的梯度信息，也就是图像的边界信息，而没有充分的利用图像区域的全局信息，必然导致了在图像灰度比较接近，边界模糊的地方分割结果发生错误。

我们知道，一般来说，图像是由一块块的小区域构成，每块小区域都具有相同的性质，像灰度值比较接近，纹理结构相似等。如何有效的利用这些相邻的小区域的信息对最后的分割结果至关重要。为此，我们引入Watershed变换，先把图像过度分割(分割的区域数目超过图像中包含的实际对象数)成许多这样的小区域。

对图像进行过度分割有三个方面的好处。一是对Fast Marching方法而言，假设种子点(活动点)位于一个小区域的内部，我们只需计算种子点到其边界像素点的到达时间，对于这个小区域内部的其它像素点，由于其具有相同的性质(在我们算法中，指的是像素点的灰度值比较接近)，所以种子点到这些像素点的到达时间无需计算，这样Fast Marching方法的速度大大提高；二是由于过度分割后这些小区域的边界有可能是潜在的正确的分割结果，所以我们约定最后的分割结果的边界一定位于由过度分割得到的区域的边界上，从而减少了分割的盲目性，提高分割的准确性；三是相邻的这些小区域的局部统计特征的相似度对于Fast Marching方法的速度函数具有很好的参考价值。

记Watershed变换对图像进行过度分割后的区域为R₁，R₂，R₃，...，R_n，n为过度分割后的区域数。每个区域R_i的统计特征记为：{SF_i，j|l i n，l j m}，其中m表示区域统计特征的个数。SF_i，1是该区域内象素的灰度均值，SF_i，2是该区域内象素的灰度方差。记录的特征的个数可以根据图像的特性来决定，例如图像中有明显的纹理的话，可以在区域特征集中增加区分纹理的特征。在我们的实验中，只用到SF_i，1和SF_i，2两个特征。

图4所示的是Watershed变换后两个相邻的区域R₁和R₂，种子点(活动点)s₁设在区域R₁中。图中红线表示的是区域R₁和R₂共同的边界，p₁，p₂，p₃...是边界上的像素点，而q₁，q₂，q₃...是区域R₁其它边界上的像素点。

当种子点s₁到达区域边界时，我们可以计算出区域R₁和R₂的统计特征SF_1，j和SF_2，j，并按如下公式计算出它们之间的相似度： ${\mathrm{Sm}}_{i,j}=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}|{\mathrm{SF}}_{i,m}-{\mathrm{SF}}_{j,m}|----\left(7\right)$

这里Sm_i，j表示区域R_i和区域R_j之间的相似度，m表示区域R_i和区域R_j的统计特征个数，SF_i，m表示区域R_i的统计特征，SF_j，m表示区域R_j的统计特征。得到区域R_i和区域R_j之间的相似度，我们就可以定义区域R_i和区域R_j共同的边界上像素点的速度函数： $F_{i,j}=e^{-\mathrm{\β}\left|{\mathrm{Sm}}_{i,j}\right|}----\left(8\right)$

这里是0β1的正数。

综上所述，改进后的Fast Marching方法的步骤如下：

1.初始化：

1)活动点：标识用户指定的种子点(可以有多个)及种子点所在的区域的内部像素点为活动点，时间T(x，y)＝0；

2)窄带点：标识所有种子点所在区域的边界像素点为窄带点。时间T(x，y)＝1/F(x，y)；

3)远离点：除了活动点和窄带点外，所有其它的网格点为远离点，T(x，y)＝TIME_MAX；

2.循环：

1)开始循环：设点(i，j)是窄带中具有最小时间T的点，不妨设其所在的区域为R_i，R_j；

2)计算区域R_i和R_j的相似度，标识点(i，j)为活动点，并从窄带中删除；

3)标识点(i，j)的所在区域的所有边界点：如果点(i，j)的所在区域的边界点为活动点，则不改变时间；如果其所在区域的边界点为窄带点，则按照公式(4)更新边界点的时间；如果其所在区域的边界点为远离点，则标识该边界点为窄带点，同时按照公式(4)更新该边界点的时间；

4)如果一点的到达时间超过指定的阈值，则结束循环，否则跳到1)；运行结果

我们用C++语言实现了本发明所描述的算法，并且在几个不同的数据集上作了实验。所有的实验都是在一台PIII 800，128MB内存，操作系统为Windows 2000的PC机上完成的。

图5、6、7、8为几组实例。图5(a)为原始的CT膝关节图像，图5(b)是利用Watershed方法的过渡分割图，图5(c)是利用改进的FastMarching方法的最后分割图，和图3想比，改进后的Fast Marching方法能够在边界模糊的地方获得满意的分割结果。图6(a)是脑部MR原图像，图6(b)是Watershed过度分割图，图6(c)是最后的脑回分割结果。图7(a)是脑部MR原图像，图7(b)是Watershed过度分割图，图7(c)是最后的脑白质分割图。图8如同图7。

这几组实例充分证明了我们算法的有效性。

Claims (5)

1.一种基于水平集和Watershed方法的医学图像分割方法，包括步骤：

各向异性扩散滤波，用于去除噪声；

采用Watershed方法，对图像进行过度分割；

堆数据结构的建立，用于定位窄带中具有最小时间T的网格点；

采用改进的Fast Marching方法，对图像进行最后的分割。

2.按权利要求1所述的方法，其特征在于所述的改进的Fast Marching方法包括如下步骤：

1.初始化：

1)活动点：标识用户指定的种子点及种子点所在的区域的内部像素点为活动点；

2)窄带点：标识所有种子点所在区域的边界像素点为窄带点；

3)远离点：除了活动点和窄带点外，所有其它的网格点为远离点；

2.循环：

1)开始循环：设点(i，j)是窄带中具有最小时间T的点；

2)计算区域R_i和R_j的相似度，标识点(i，j)为活动点，并从窄带中删除；

3)标识点(i，j)所在区域的所有边界点：如果点(i，j)的所在区域的边界点为活动点，则不改变时间；如果其所在区域的边界点为窄带点，更新边界点的时间；如果其所在区域的边界点为远离点，则标识该边界点为窄带点，更新该边界点的时间；

4)如果种子点到某一点的到达时间超过指定的阈值，则结束循环，否则重新开始循环；

3.按权利要求2所述的方法，其特征在于所述的更新边界点的时间按下述公式进行： $[{\max {(D_{\mathrm{ij}}^{-x}T,0)}^2+\min {(D_{\mathrm{ij}}^{+x}T,0)}^2+\max {(D_{\mathrm{ij}}^{-y}T,0)}^2+{\min (D_{\mathrm{ij}}^{+y}T,0)}^2]}^{1/2}=1/F_{i,j}$

4.按权利要求2所述的方法，其特征在于所述的计算区域Ri和Rj的相似度按下述公式进行： ${\mathrm{Sm}}_{i,j}=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}|{\mathrm{SF}}_{i,m}-{\mathrm{SF}}_{j,m}|$

5.按权利要求2所述的方法，其特征在于所述的改进的Fast Marching方法的初始速度按照下述公式进行： $\underset{\‾}{F_{i,j}=e^{-\mathrm{\β}\left|{\mathrm{Sm}}_{i,j}\right|}}$

Images