CN101315699A - 一种增量式变分水平集快速医学图像分割方法 - Google Patents

Abstract

本发明提供了一种增量式变分水平集快速医学图像分割方法，首先选取初始边界，根据初始边界计算的子区域及其平均灰度，采用窄带法等快速算法求解水平集曲线演化过程，提取零水平集即新的边界，判断是否满足停止条件，若是则得到分割结果，若否则利用边界的运动导致区域的改变，根据增量在窄带范围内计算新区域的平均灰度，再进行窄带法等快速算法求解水平集曲线演化过程，最终的零水平集即为分割结果。本发明根据区域和区域中像素的动态变化，采用增量方式迭代求取区域的平均灰度，将其解析公式改进为递进迭代公式，从而可以采用窄带法等快速算法，较大地提高了分割效率，使得该模型更具实际意义。

Description

一种增量式变分水平集快速医学图像分割方法

(一)技术领域

本发明涉及医学图像分割领域，尤其涉及一种快速的医学图像分割方法。

(二)背景技术

在医学图像处理与分析应用中，图像分割技术起着关键的作用。医学图像分割的任务是从医学图像中提取包含重要诊断信息的感兴趣区域，为临床诊疗和病理学研究提供可靠的依据。由于医学图像的成像原理和人体组织结构本身的复杂性和差异性，医学图像与普通图像比较，不可避免地具有模糊、不均匀性等特点；同时，医学成像技术的快速发展使得获取各种复杂的海量医学图像数据成为可能，这些都对图像分割技术提出了更高的要求。

近年来，基于曲线演化几何流的变形模型成为医学图像分割领域的研究热点之一。基于几何流的变形模型又分为两类：一类是直接得到问题的几何流描述，以基于曲线演化的方法为代表，曲线演化方程即是一个几何流；另一类是来源于能量最小化的几何流，其思想是将所研究问题归结为一个泛函能量函数的极小问题，然后应用变分方法导出相应的几何流。近些年提出的水平集方法对基于几何流的图像分割方法的研究和应用具有深远的影响。

水平集的基本思想是：将运动的轮廓曲线(曲面)隐含地表达为高一维的曲面的零水平集，并通过曲面的运动来隐含地求解轮廓曲线(曲面)的演化。该方法是一种能自适应边界拓扑变化的稳定数值方法，不但提供了曲线演化的精确数值解法，而且很好地解决了非常棘手的拓扑变化问题，因此已发展为当今医学图像分割中最为活跃和成功的研究领域之一。

医学图像分割方法中，Chan-Vese模型是Chan和Vese在2001年提出的一种不依赖于梯度的分片常值简化Mumford-Shah模型，该模型根据区域的均值将图像仅分为两个区域(目标和背景)，即两相分割，并采用水平集方法进行数值求解。

基于变分水平集方法的Chan-Vese模型综合利用了整个图像区域的全部信息，因而其显著特点是全局优化性。初始轮廓曲线无需完全定位于目标的内部或者外部，而且仅需要一条初始轮廓曲线，就可以分割即使内部含有空洞的目标。另外，该模型不依赖于边缘梯度，即使图像中的边缘比较模糊，仍然能够得到比较理想的结果。

但该模型存在一个重要缺陷，即Chan-Vese模型的演化过程需要利用整个图像的全部信息，并且每次迭代都需要(且必须)在整个图像域上计算，因而无法采用窄带法或快速步进法等快速算法来提高分割效率。

(三)发明内容

本发明的目的在于提供一种根据区域和区域中像素的动态变化，采用增量方式迭代求取区域的平均灰度，将其解析公式改进为递进迭代公式，从而可以采用窄带法等快速算法，能够较大地提高分割效率的增量式变分水平集快速医学图像分割方法。

本发明的目的是这样实现的：首先选取初始边界，根据初始边界计算的子区域及其平均灰度，采用窄带法等快速算法求解水平集曲线演化过程，提取零水平集即新的边界，判断是否满足停止条件，若是则得到分割结果，若否则利用边界的运动导致区域的改变，根据增量在窄带范围内计算新区域的平均灰度，再进行窄带法等快速算法求解水平集曲线演化过程，最终的零水平集即为分割结果。

本发明的详细过程为：

首先根据实际分割问题选择合适的模型参数μ，v，λ₁，λ₂，进行如下分割过程：设ΔΩ₁，ΔΩ₂分别表示区域Ω₁，Ω₂增加或减少的部分，符号|Ω|表示区域Ω的面积，对于离散问题来说即区域上的点的个数；

设动态区域为有界连通开集，ΔΩ为其改变的区域，

(1)给定初始边界C₀，以此边界可得子区域Ω₁ ⁰，Ω₂ ⁰，分别计算其平均灰度ω₁ ⁰，ω₂ ⁰，并将ω₁ ⁰，ω₂ ⁰分别作为ω₁，ω₂；

(2)根据公式 $\left\{\begin{array}{c}{\∫\∫}_{\mathrm{\Δ}{\mathrm{\Ω}}_1}\tilde{H}\left({\mathrm{\ψ}}_1\right)\mathrm{dxdy}+{\∫\∫}_{{\mathrm{\Δ\Ω}}_2}\tilde{H}\left({\mathrm{\ψ}}_2\right)\mathrm{dxdy}=0\\ {\∫\∫}_{{\mathrm{\Δ\Ω}}_1}f\·\tilde{H}\left({\mathrm{\ψ}}_1\right)\mathrm{dxdy}+{\∫\∫}_{{\mathrm{\Δ\Ω}}_2}f\·\tilde{H}\left({\mathrm{\ψ}}_2\right)\mathrm{dxdy}=0\end{array}\right.$

对轮廓曲线C₀进行单步演化，实现上可以采用窄带法或快速步进法来快速求解；

(3)提取上一步演化结果的零水平集，并以此为新的边界，得到两个更新后的区域Ω₁，Ω₂，检查演化停止条件，满足则转到步骤(7)；

(4)分别计算区域Ω₁，Ω₂的变化区域ΔΩ₁，ΔΩ₂；

(5)采用递进方式，分别按照公式

${\mathrm{\ω}}_1=\frac{{\mathrm{\ω}}_1^0\left|{\mathrm{\Ω}}_1^0\right|+{\∫\∫}_{{\mathrm{\Δ\Ω}}_1}f\·\tilde{H}\left({\mathrm{\ψ}}_1\right)\mathrm{dxdy}}{\left|{\mathrm{\Ω}}_1^0\right|+{\∫\∫}_{{\mathrm{\Δ\Ω}}_1}\tilde{H}\left({\mathrm{\ψ}}_1\right)\mathrm{dxdy}}$ ${\mathrm{\ω}}_2=\frac{{\mathrm{\ω}}_2^0\left|{\mathrm{\Ω}}_2^0\right|+{\∫\∫}_{{\mathrm{\Δ\Ω}}_2}f\·\tilde{H}\left({\mathrm{\ψ}}_2\right)\mathrm{dxdy}}{\left|{\mathrm{\Ω}}_2^0\right|+{\∫\∫}_{{\mathrm{\Δ\Ω}}_2}\tilde{H}\left({\mathrm{\ψ}}_2\right)\mathrm{dxdy}}$

计算区域Ω₁，Ω₂的平均灰度ω₁，ω₂；

(6)转到步骤(2)；

(7)区域Ω₁，Ω₂即为分割结果。

本发明根据区域和区域中像素的动态变化，采用增量方式迭代求取区域的平均灰度，将其解析公式改进为递进迭代公式，从而可以采用窄带法等快速算法，较大地提高了分割效率。一种快速的医学图像分割方法，属于医学图像分割领域。本发明针对Chan-Vese模型在演化过程中每次迭代都必须在整个图像区域上计算而导致无法采用快速算法的缺陷，从基本理论公式的推导出发，提出一种增量式变分水平集快速医学图像分割方法。该方法根据区域和区域中像素的动态变化，采用增量方式迭代求取区域的平均灰度，将其解析公式改进为递进迭代公式，从而可以采用窄带法等快速算法，较大地提高了分割效率，使得该模型更具实际意义。本发明在医学图像分割领域有着广阔的应用前景。

本发明针对Chan-Vese模型无法采用快速算法的缺陷进行改进，提出一种增量式变分水平集快速医学图像分割方法。该方法给出了区域的平均灰度的增量式迭代求取公式，不再同时需要全部图像区域上的信息，因而能够利用窄带法或快速步进法来加快求解过程，解决了Chan-Vese模型因不能采用快速算法而导致分割过慢的重要不足。由于仍然利用图像的全局信息，从而在保持分割精度不变的同时较大地提高了分割速度，使得该模型在复杂的医学图像分割领域更具实际意义。

我们对区域(以及区域中像素)的动态变化进行考察。区域的变化是因轮廓曲线(区域边界)运动而产生的。设ΔΩ₁，ΔΩ₂分别表示区域Ω₁，Ω₂增加或减少的部分，符号|Ω|表示区域Ω的面积，对于离散问题来说即区域上的点的个数。

设动态区域

为有界连通开集，ΔΩ为其改变的区域，

若

(1)

则称

为ΔΩ关于Ω的符号隶属函数，简称ΔΩ的符号隶属函数。

符号隶属函数ψ的意义是，在区域Ω的变化区域ΔΩ中，对于进入Ω的像素，符号隶属值取1；对于离开Ω的像素，符号隶属值取-1。

令

为符号函数，即

$\tilde{H}\left(x\right)=\left\{\begin{array}{cc}1,& x>0\\ 0,& x=0\\ -1,& x<0\end{array}\right.$ (2)

设动态区域

为有界连通开集，ΔΩ为其改变的区域，若

$\left(3\right),A\left(\mathrm{\&Delta;\&Omega;}\right)={\&Integral;\&Integral;}_{\mathrm{\&Delta;\&Omega;}}\tilde{H}\left(\mathrm{\&psi;}\right)\mathrm{dxdy}$

则称

为ΔΩ的符号面积。

正的符号面积表示区域增大，负的符号面积表示区域缩小，零则说明区域保持大小(面积)不变的情况下改变了形状或位置。

利用式(1)和式(2)，ΔΩ₁，ΔΩ₂上的灰度变化总量可分别表示为

$\left\{\begin{array}{c}{S}_{{\mathrm{\&omega;}}_1}={\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_1}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}\\ S_{{\mathrm{\&omega;}}_2}={\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_2}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}\end{array}\right.---\left(4\right)$

其中ψ₁为ΔΩ₁关于Ω₁的符号隶属函数，ψ₂为ΔΩ₂关于Ω₂的符号隶属函数。

注意到，若区域Ω₁ ⁰经过变化区域ΔΩ₁后更新为Ω₁，则Ω₁上的平均灰度可按下式计算

${\mathrm{\&omega;}}_1=\frac{{\mathrm{\&omega;}}_1^0\left|{\mathrm{\&Omega;}}_1^0\right|+S_{{\mathrm{\&omega;}}_1}}{\left|{\mathrm{\&Omega;}}_1^0\right|+A\left({\mathrm{\&Delta;\&Omega;}}_1\right)}---\left(5\right)$

其中ω₁ ⁰为Ω₁ ⁰上的平均灰度。

将式(3)和(4)代入上式，得到

${\mathrm{\&omega;}}_1=\frac{{\mathrm{\&omega;}}_1^0\left|{\mathrm{\&Omega;}}_1^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_1}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}}{\left|{\mathrm{\&Omega;}}_1^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_1}\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}}---\left(6\right)$

同理可计算Ω₂ ⁰经过变化ΔΩ₂后的平均灰度

${\mathrm{\&omega;}}_2=\frac{{\mathrm{\&omega;}}_2^0\left|{\mathrm{\&Omega;}}_2^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_2}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}}{\left|{\mathrm{\&Omega;}}_2^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_2}\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}}---\left(7\right)$

综上分析，整理得

$\left\{\begin{array}{c}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=|\&dtri;\mathrm{\&phi;}|[\mathrm{\&mu;}\&dtri;\&CenterDot;\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}-v-{\mathrm{\&lambda;}}_1{(f-{\mathrm{\&omega;}}_1)}^2+{\mathrm{\&lambda;}}_2{(f-{\mathrm{\&omega;}}_2)}^2]\\ {\mathrm{\&omega;}}_1=\frac{{\mathrm{\&omega;}}_1^0\left|{\mathrm{\&Omega;}}_1^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_1}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}}{\left|{\mathrm{\&Omega;}}_1^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_1}\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}}\\ {\mathrm{\&omega;}}_2=\frac{{\mathrm{\&omega;}}_2^0\left|{\mathrm{\&Omega;}}_2^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_2}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}}{\left|{\mathrm{\&Omega;}}_2^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_2}\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}}\end{array}\right.---\left(8\right)$

我们将上式称为增量式Chan-Vese模型，并将利用该模型进行医学图像分割的方法称为增量式变分水平集快速医学图像分割方法。

事实上，在曲线演化过程中，区域的变化是由于边界曲线(即零水平集)C的单步运动产生的，因而变化区域ΔΩ₁，ΔΩ₂总是在零水平集邻近的小范围内，可以全部包含在一个窄带之中，因此可以采用窄带法或快速步进法来快速求解。

对于两相分割问题，进入Ω₁中的像素就是离开Ω₂的像素，因此有ΔΩ₁＝-ΔΩ₂，进一步，这两个变化区域上的符号面积及其像素值的和分别满足如下关系：

$\left\{\begin{array}{c}{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_1}\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_2}\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}=0\\ {\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_1}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_2}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}=0\end{array}\right.---\left(9\right)$

增量式Chan-Vese模型依然利用了图像全局信息，保留了全局优化的优点，且减少了计算量。一方面，迭代过程中平均灰度的求取大为简化。由式(8)和式(9)可知，仅仅需要在区域的变化的部分计算，而传统的计算公式必须对整个图像区域计算。另一方面，可以采用一般快速算法来加速水平集曲线演化过程，从而进一步降低计算量。

需要指出的是，从增量式Chan-Vese模型的推导过程可知，该模型与传统的Chan-Vese模型在理论上是完全等价的，因此分割精度不受影响。

本发明针对Chan-Vese模型在演化过程中每次迭代都必须在整个图像区域上计算而导致无法采用快速算法的缺陷，从基本理论公式的推导出发，提出一种增量式变分水平集快速医学图像分割方法。该方法根据区域和区域中像素的动态变化，采用增量方式迭代求取区域的平均灰度，将其解析公式改进为递进迭代公式，从而可以采用窄带法等快速算法，较大地提高了分割效率，使得该模型更具实际意义。本发明在医学图像分割领域有着广阔的应用前景。

(四)附图说明

图1为增量式变分水平集方法的算法流程图；

图2为脑MR(T₁)图像；

图3-图5为脑组织分割结果。

(五)具体实施方式

下面结合附图对本发明作进一步的说明：

本发明中提出的增量式变分水平集快速医学图像分割方法主要依据公式(9)进行数值计算，实施简洁，以下结合附图对本发明的具体实施方式作进一步的详细说明。

二维图像的两相分割问题为例，原图像如图2所示。首先根据实际分割问题选择合适的模型参数μ，v，λ₁，λ₂，本实施例中分别为μ＝1.0，v＝0，λ₁＝λ₂＝1.0，然后参照图1进行如下分割过程：

(1)给定初始边界C₀，以此边界可得区域Ω₁ ⁰，Ω₂ ⁰，分别计算其平均灰度ω₁ ⁰，ω₂ ⁰，并将ω₁ ⁰，ω₂ ⁰分别作为ω₁，ω₂；

(2)根据公式(9)对轮廓曲线C₀进行单步演化，实现上可以采用窄带法或快速步进法来快速求解；

(3)提取上一步演化结果的零水平集，并以此为新的边界，得到两个更新后的区域Ω₁，Ω₂，检查演化停止条件，满足则转到(7)；

(4)分别计算区域Ω₁，Ω₂的变化区域ΔΩ₁，ΔΩ₂；

(5)采用递进方式，分别按照公式(6)和公式(7)计算区域Ω₁，Ω₂的平均灰度ω₁，ω₂；

(6)转到(2)；

(7)区域Ω₁，Ω₂即为分割结果。

图3-图5给出了分割结果，从左至右分别为脑灰质、脑白质和脑脊液。

图1中：手工选取的初始边界101，根据初始边界计算的子区域及其平均灰度102，采用窄带法等快速算法求解水平集曲线演化过程103，提取零水平集即新的边界104，判断是否满足停止条件，若是则得到分割结果，若否则利用边界的运动导致区域的改变105，根据增量公式(6)和(7)在窄带范围内计算新区域的平均灰度106，再进行窄带法等快速算法求解水平集曲线演化过程103，最终的零水平集即为分割结果107。图3中分别为脑灰质301、脑白质302、脑脊液303。

表1给出了本发明方法与原Chan-Vese模型分割的时间比较，表中数据表明本发明方法减少了约22％的分割时间。虽然对于单幅图像来说分割时间以ms计，但实际医学图像通常都具有几十或几百个切片，因此节省的时间是非常可观的。

表1分割时间的比较

(单位：ms)

Claims (2)

1、一种增量式变分水平集快速医学图像分割方法，其特征在于：首先选取初始边界，根据初始边界计算的子区域及其平均灰度，采用窄带法等快速算法求解水平集曲线演化过程，提取零水平集即新的边界，判断是否满足停止条件，若是则得到分割结果，若否则利用边界的运动导致区域的改变，根据增量在窄带范围内计算新区域的平均灰度，再进行窄带法等快速算法求解水平集曲线演化过程，最终的零水平集即为分割结果。

2、根据权利要求1所述的一种增量式变分水平集快速医学图像分割方法，其特征在于：首先根据实际分割问题选择合适的模型参数μ，v，λ₁，λ₂，进行如下分割过程：设ΔΩ₁，ΔΩ₂分别表示区域Ω₁，Ω₂增加或减少的部分，符号|Ω|表示区域Ω的面积，对于离散问题来说即区域上的点的个数；设动态区域

为有界连通开集，ΔΩ为其改变的区域，

(1)给定初始边界C₀，以此边界可得子区域Ω₁ ⁰，Ω₂ ⁰，分别计算其平均灰度ω₁ ⁰，ω₂ ⁰，并将ω₁ ⁰，ω₂ ⁰分别作为ω₁，ω₂；

(2)根据公式 $\left\{\begin{array}{c}{\&Integral;\&Integral;}_{\mathrm{\&Delta;}{\mathrm{\&Omega;}}_1}\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}+{\&Integral;\&Integral;}_{\mathrm{\&Delta;}{\mathrm{\&Omega;}}_2}\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}=0\\ {\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_1}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_2}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}=0\end{array}\right.$

对轮廓曲线C₀进行单步演化，实现上可以采用窄带法或快速步进法来快速求解；

(3)提取上一步演化结果的零水平集，并以此为新的边界，得到两个更新后的区域Ω₁，Ω₂，检查演化停止条件，满足则转到步骤(7)；

(4)分别计算区域Ω₁，Ω₂的变化区域ΔΩ₁，ΔΩ₂；

(5)采用递进方式，分别按照公式

${\mathrm{\&omega;}}_1=\frac{{\mathrm{\&omega;}}_1^0\left|{\mathrm{\&Omega;}}_1^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_1}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}}{\left|{\mathrm{\&Omega;}}_1^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_1}\tilde{H}\left({\mathrm{\&psi;}}_1\right)\mathrm{dxdy}}$ ${\mathrm{\&omega;}}_2=\frac{{\mathrm{\&omega;}}_2^0\left|{\mathrm{\&Omega;}}_2^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_2}f\&CenterDot;\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}}{\left|{\mathrm{\&Omega;}}_2^0\right|+{\&Integral;\&Integral;}_{{\mathrm{\&Delta;\&Omega;}}_2}\tilde{H}\left({\mathrm{\&psi;}}_2\right)\mathrm{dxdy}}$

计算区域Ω₁，Ω₂的平均灰度ω₁，ω₂；

(6)转到步骤(2)；

(7)区域Ω₁，Ω₂即为分割结果。

Images