CN109636816A - 一种超声图像分割方法 - Google Patents

Abstract

本发明公开了一种超声图像分割方法，实现步骤为：首先读取超声图像，其次手动初始化一条边缘轮廓曲线，通过使用多尺度滤波方法构建图像的边缘检测函数和演化模型的气球力因子，嵌入到水平集演化模型当中通过更新迭代水平集函数最终实现超声图像边缘轮廓的分割，本发明通过使用多尺度滤波方法和改进的气球力因子，避免了超声图像中由于弱边缘问题导致的边缘泄露现象，提高了超声图像的分割精度。

Description

一种超声图像分割方法

技术领域

本发明涉及图像处理技术领域，特别涉及一种超声图像分割方法。

背景技术

高强度聚焦超声(High Intensity Focused Ultrasound，HIFU)消融技术作为体外治疗肿瘤的医疗技术，具有无创性、治疗效果显著、副作用小、术后恢复快等优点，随着影像学的发展已经逐步应用在临床治疗肿瘤领域。HIFU治疗过程中需要使用手术引导系统为医生提供可视化的手术环境，由于超声成像实时性好，费用低，对人体无辐射伤害的特点，基于超声图像的HIFU引导系统得到了广泛应用。超声肿瘤图像的精确分割是HIFU手术成功与否的关键，然而受超声成像机制的影响，超声图像往往具有严重的斑点噪声和灰度低对比度，导致超声图像分辨率低，目标肿瘤边缘不清晰，给HIFU超声图像的精确分割增加了难度。

活动轮廓模型将图像的底层信息(灰度、梯度等)与高层信息(目标形状、位置)等有效的结合，因其平滑性和自适应的特点在超声图像分割领域得到了广泛应用。按照曲线构造方式的不同，活动轮廓模型可以分为参数活动轮廓模型和水平集活动轮廓模型，参数活动轮廓模型定义一条参数化曲线，结合曲线的内在特征和外部图像力将分割问题转变为能量泛函极小值问题，但是对曲线参数化的表征导致参数活动轮廓模型无法自主的处理曲线的拓扑演化问题。水平集活动轮廓模型使用高一维的水平集函数来隐含的表征初始曲线，曲线拓扑结构发生改变时，对高维水平集函数并无影响。Caselles等人和Malladi等人各自独立提出了基于水平集方法的几何活动轮廓模型，通过变分法和梯度下降流得到曲线稳态时的偏微分方程，但在水平集迭代过程需要反复重新初始化水平集函数，使水平集函数保持为符号距离函数，重复初始化过程增加了演化的计算复杂度，不利于图像的快速分割。Adalsteinsson等人提出了窄带法(Narrow Band)将计算区域限制在零水平及周围的一个小范围窄带内，缩小了计算量，Caselles等人提出了测地线活动轮廓模型，在几何活动轮廓模型的基础上增加了一个测地距离项，当曲线超出图像真实边缘时，该测地项能够牵引抑制曲线的延伸，一定程度上解决了图像边缘泄露问题。Li等人在水平集能量函数中添加了一个距离惩罚项，该惩罚项能够对偏离零水平集的情况进行修正，提出了距离正则项水平集模型，该模型使用势能函数定义正则项，具有更高的数值稳定性，避免了水平集演化过程中的重复初始化问题，降低了计算复杂度。学者们针对不同的医学图像应用实例，提出了不同的改进距离正则项模型，但对于具有严重斑点噪声，灰度低对比度、弱边缘的超声图像仍易出现边缘泄露问题，引发错误的分割结果。

发明内容

本发明旨在提供一种改进的多尺度滤波与弱边缘检测的超声图像分割方法，能够消除在HIFU超声图像分割中模型存在的边缘泄露问题。

本发明的具体技术实现方案如下：

一种用于HIFU治疗中的超声肿瘤图像分割方法，该方法包含以下步骤：

步骤1：读取原始图像I(x,y)。

步骤2：初始化轮廓曲线，手动在原始图像上选取一条闭合曲线，该曲线在目标边缘的附近，可与目标边缘相交。

步骤3：对图像进行预处理，构建多尺度滤波图像，在有效滤除图像斑点噪声的同时，更好的保留图像的边缘细节信息。根据多尺度滤波图像构建改进的边缘检测函数；

步骤4：由于二阶微分算子在图像的边缘两侧具有完全相反的符号，使用多尺度梯度下的图像拉普拉斯算子构造气球力因子，使气球力因子具有双向演化作用，提高模型弱边缘检测能力。

步骤5：使用步骤3与步骤4中获取的边缘检测函数与气球力因子，嵌入到距离正则项水平集活动轮廓模型中，得到改进的能量函数。

步骤6：基于变分法和梯度下降流得到能量函数最小化后的水平集演化方程，通过迭代更新水平集函数，当达到终止条件时，零水平集所表示的分割轮廓即是待分割目标边缘的真实轮廓。

步骤3包括：

步骤3-1，设置一组高斯标准差尺度σ₁,σ₂…σ_N，其中N表示尺度的个数，σ_N表示第N个高斯标准差尺度，初始化梯度矢量场x表示像素点的横坐标，y表示像素点的纵坐标，I表示原始图像，σ₁表示最大初始尺度；一般选取较大值，使之有效的滤出图像噪声，表示σ₁下的高斯滤波器。

步骤3-2，选取一定的尺度间隔η，在尺度由大到小的变化过程中σ_i+1＝σ_i-η，σ_i＞0，对每一图像像素点有：

公式中i＝1,2...N-1，V(x,y)表示初始梯度矢量场，V_i(x,y)表示σ₁下的梯度矢量场，表示V_i(x,y)的缩写。

该方法的判断依据是：随着i的提高，高斯标准差尺度σ_i逐渐减小，在尺度从大到小的过程中所产生的图像将渐渐受到噪声的影响(滤波不彻底)，而图像中噪声矢量方向是不规律的，因此在尺度从大到小的过程中，对任一像素点(x,y)，当两相邻尺度下的梯度矢量场之间的夹角小于45度时，认为该像素点的梯度是增强的，该梯度方向下的矢量场与图像的真实边缘与细节一致，从而对该点的矢量场进行叠加，反之，如果两个相邻尺度下的梯度矢量场之间的夹角大于45度时，说明低一级尺度下的该像素点的梯度方向与真实边缘不一致，判定是由噪声引起的，该像素点在低一级尺度下的梯度场被抑制。

步骤3-3，构建改进的边缘检测函数。

步骤3-3包括：构建如下改进的边缘检测函数g_k：

在原始图像的边缘区域，g_k近似为0，而在原始图像的均匀区域，g_k近似为1。

步骤4包括：构造如下气球力因子a_k：

a_k＝γ·sgn(div(V))，

其中，γ表示正常数，sgn表示符号函数，div表示散度算子。在原始图像的真实边缘附近，当演化轮廓线在真实边缘的外部时，a_k大于0，演化曲线保持收缩演化，而当演化轮廓线处于真实边缘的内部时，a_k小于0，演化曲线保持为膨胀演化，所以a_k能够使演化轮廓线进行双向演化，克服单向演化的缺陷。

步骤5中，改进的能量函数E(φ)的表示形式为：

其中，x表示像素点的横坐标，y表示像素点的纵坐标，φ表示构造的高维水平集函数，一般使用点(x,y)到初始曲线的最短距离构建，规定演化曲线内部的点(x,y)所对应的水平集函数φ为负值，演化曲线外部的点(x,y)对对应的水平集函数φ为正值，其零水平集可以隐含的表示给定的曲线，H(·)与δ(·)分别表示Heaviside函数和Dirac函数，α_k表示气球力因子，用于控制曲线演化的方向和速度；λ表示权重系数，μ表示距离惩罚项因子，一般取为1。R_p(φ)表示外部距离惩罚项，用来约束曲线的内部能量，避免演化过程中零水平的偏移，进而实现避免重复初始化的目的。

R_p(φ)计算公式如下：

其中p表示势能函数，由双势阱函数p(s)构造而成，表达式如下：

其中s表示内部参数。

步骤6中，使用变分法和梯度下降流将能量函数最小化得到相应的水平集演化方程为：

其中t表示时间。

作为优选，距离惩罚项因子μ取为1，应用中使用正则化的Dirac函数δ_ε(·)与正则化的Heaviside函数H_ε(·)代替广义的Dirac函数δ(·)和Heaviside函数H(·)，参数ε一般选为1，气球力因子a_k决定演化曲线的方向和速度，即γ的取值决定了演化速度与方向，初始尺度σ₁与尺度间隔η的选取根据图像信息而定。

本发明具有如下优点：通过使用多尺度滤波方法，在实现有效滤波的同时，对图像的边缘细节信息有一定的增强能力，避免了人为选择高斯标准差所带来的不足；使用多尺度滤波下的拉普拉斯算子来重构气球力因子，在实现双向演化的同时，提高了模型对弱边缘的检测能力，为HIFU治疗中超声肿瘤图像提供了精确的分割结果。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述或其他方面的优点将会变得更加清楚。

图1为本发明实施样例的流程图。

图2为实施样例中未实施本发明方法的超声图像分割结果图。

图3为本发明实施样例超声图像分割结果图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

本发明的核心思想是：提出一种基于多尺度滤波与弱边缘检测的超声图像分割方法，为HIFU治疗手术中医生提供精确的肿瘤边缘轮廓信息。本发明的优点如下所述：

使用多尺度滤波的方法，避免了人工选取高斯标准差带来的不足，在有效滤波的同时，对图像的边缘细节有一定的增强能力；使用多尺度滤波下的拉普拉斯算子来构造演化模型的气球力因子，消除了单向演化的缺陷，在实现双向演化的同时，增强了图像的弱边缘检测能力；相比于传统的距离正则项水平集演化模型，本发明避免了图像分割中边缘泄露现象，提高了超声图像的分割精度。

图1表示本发明的实施样例流程图，如图1所示，该实施样式包括以下步骤：

步骤1：读取原始待分割超声图像；

步骤2：初始化边缘轮廓线，该初始轮廓线一般设置在图像目标边缘轮廓附近。

步骤3：使用多尺度滤波方法，设置一组高斯标准差尺度σ₁,σ₂…σ_N，其中N表示尺度的个数，初始化梯度矢量场σ₁表示最大初始尺度，一般选取较大值，使之有效的滤出图像噪声，作为优选，本实施样例中σ₁取为10，选取一定的尺度间隔η，在尺度由大到小的变化过程中σ_i+1＝σ_i-η，σ_i＞0，对每一图像像素点有：

公式中i＝1,2...N-1，作为优选，本实施样例中η设置为2，N设置为4。使用该多尺度滤波方法，得到图像的边缘梯度图，构建边缘检测函数：

步骤4：使用多尺度滤波的图像拉普拉斯算子重新构造能量方程的气球力因子：

a_k＝γ·sgn(div(V))，γ表示正常数，用于控制演化曲线的演化方向，sgn表示符号函数，div表示散度算子，作为优选，本实施样例中γ设置为2。

步骤5：使用步骤3与步骤4中的边缘检测函数g_k和演化曲线的气球力因子a_k，

嵌入到距离正则项水平集活动轮廓模型中，得到本发明所构造的能量泛函。能量泛函的具体表示形式为：

这里φ表示构造的水平集函数，H(·)与δ(·)分别表示Heaviside函数和Dirac函数，λ表示权重系数，μ表示距离惩罚项因子，一般取为1。R_p(φ)表示外部距离惩罚项，用来约束曲线的内部能量，避免演化过程中零水平的偏移，进而实现避免重复初始化的目的。由构成，其中p表示势能函数，由双势阱函数构造而成，表达式如下：

步骤6：基于变分法以及梯度下降流得到能量泛函最小化后的水平集演化方程：

作为优选，距离惩罚项因子μ取为1，应用中使用正则化的Dirac函数δ_ε(·)与正则化的Heaviside函数H_ε(·)代替广义的Dirac函数δ(·)和Heaviside函数H(·)，参数ε选为1，气球力因子a_k决定演化曲线的方向和速度，即γ的取值决定了演化速度与方向。通过迭代更新水平集函数，当达到终止条件时，表示演化曲线已经收敛至图像的目标真实边缘轮廓。图2和图3分别显示了未用本发明的距离正则项水平集活动轮廓模型的分割结果和使用本发明的水平集活动轮廓模型的分割结果比较，由图2可知，未使用本发明的距离正则项水平集活动轮廓模型发生了边缘泄露现象，无法正确的分割出图像的弱边缘情形，而如图3所示，使用本发明的水平集活动轮廓模型则避免了演化模型的边缘泄露情形，综上所示，本发明通过使用多尺度滤波算法和改进的气球力因子，克服了超声图像中由于弱边缘导致的边缘泄露问题，提高了图像分割的准确率。

本发明提供了一种超声图像分割方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (7)

1.一种超声图像分割方法，其特征在于，包括以下步骤：

步骤1：读取原始图像；

步骤2：手动在原始图像上选取一条闭合曲线，作为初始化轮廓曲线；

步骤3：对图像进行预处理，构建多尺度滤波图像，根据多尺度滤波图像构建改进的边缘检测函数；

步骤4：构造气球力因子；

步骤5：使用步骤3与步骤4中获取的边缘检测函数与气球力因子，得到改进的能量函数；

步骤6：基于变分法和梯度下降流得到能量函数最小化后的水平集演化方程，通过迭代更新水平集函数，当达到终止条件时，零水平集所表示的分割轮廓即是待分割目标边缘的真实轮廓。

2.根据权利要求1所述的方法，其特征在于，步骤3包括：

步骤3-1，设置一组高斯标准差尺度σ₁,σ₂…σ_N，其中N表示尺度的个数，σ_N表示第N个高斯标准差尺度，初始化梯度矢量场x表示像素点的横坐标，y表示像素点的纵坐标，I表示原始图像，σ₁表示最大初始尺度；表示σ₁下的高斯滤波器；

步骤3-2，选取一定的尺度间隔η，在尺度由大到小的变化过程中σ_i+1＝σ_i-η，σ_i＞0，对每一图像像素点有：

公式中i＝1,2...N-1，V(x,y)表示初始梯度矢量场，V_i(x,y)表示σ₁下的梯度矢量场，表示V_i(x,y)的缩写；

步骤3-3，构建改进的边缘检测函数。

3.根据权利要求2所述的方法，其特征在于，步骤3-3包括：构建如下改进的边缘检测函数g_k：

在原始图像的边缘区域，g_k近似为0，而在原始图像的均匀区域，g_k近似为1。

4.根据权利要求3所述的方法，其特征在于，步骤4包括：构造如下气球力因子a_k：

a_k＝γ·sgn(div(V))，

其中，γ表示正常数，sgn表示符号函数，div表示散度算子。

5.根据权利要求4所述的方法，其特征在于，步骤5中，改进的能量函数E(φ)的表示形式为：

其中，x表示像素点的横坐标，y表示像素点的纵坐标，φ表示构造的高维水平集函数，H(·)与δ(·)分别表示Heaviside函数和Dirac函数，α_k表示气球力因子，用于控制曲线演化的方向和速度；λ表示权重系数，μ表示距离惩罚项因子，R_p(φ)表示外部距离惩罚项，用来约束曲线的内部能量。

6.根据权利要求5所述的方法，其特征在于，R_p(φ)计算公式如下：

其中p表示势能函数，由双势阱函数p(s)构造而成，表达式如下：

其中s表示内部参数。

7.根据权利要求6所述的方法，其特征在于，步骤6中，使用变分法和梯度下降流将能量函数最小化得到相应的水平集演化方程为：

其中t表示时间。