CN103793916A - 一种hifu治疗中子宫肌瘤超声图像分割方法 - Google Patents

Abstract

本发明公开了一种HIFU治疗中子宫肌瘤超声图像分割方法，其实现过程为：用户手动初始化肿瘤轮廓；构造图像高斯金字塔，得到粗尺度图像；构建形状约束能量，利用引入形状约束后的局域区域活动轮廓模型对粗尺度图像进行粗分割；将得到的粗尺度分割结果作为对原始图像进行分割的初始轮廓，再次利用引入形状约束后的局域区域活动轮廓模型对原始图像进行细分割；其细分割结果即是最后的肿瘤分割轮廓。该发明通过一种新的形状约束的引入提高了子宫肌瘤超声图像分割的准确性，避免了由于超声图像质量导致的边界泄露或过度收缩问题，而多尺度算法的实施极大地提高了分割效率。

Description

一种HIFU治疗中子宫肌瘤超声图像分割方法

技术领域

本发明涉及图像处理技术领域，涉及一种HIFU治疗中子宫肌瘤超声图像分割方法，一种HIFU治疗中基于形状约束的多尺度局域区域活动轮廓模型的子宫肌瘤超声图像分割方法。

背景技术

子宫肌瘤是困扰女性的一种常见的良性肿瘤。目前，高强度聚焦超声(HIFU)治疗作为一种新型的无创手术，由于其具有安全性和有效性，已被逐渐应用到子宫肌瘤的治疗中，减轻了传统手术中的患者痛苦。众所周知，超声图像由于其低信噪比，弱边界等问题，其准确分割一直是分割领域中的一个难题，至今仍未有完美的解决方案。而应用于HIFU治疗中的超声图像由于治疗过程中水介质干扰和角度恒定等问题，具有更加严重的噪声和边界的模糊性。然而，其图像中肿瘤靶区的分割是HIFU治疗中最为关键的一个步骤，其分割的准确性决定着手术的最终效果，而且，相比于其它HIFU系统，基于超声引导的HIFU系统的一个巨大优势就是实时性，而肿瘤区域分割的速度却极大影响着手术整体的进程。因此，一个应用于HIFU系统中高效而准确的分割方法亟待解决。

最近几年，活动轮廓模型广泛应用于图像分割。由于其平滑性和封闭的特征，被许多人应用到医学图像中并取得了较好的分割结果。活动轮廓模型首先由Kass等人提出，经过长时间的发展，按照所依据分割信息的不同，主要可以分为两类：基于边缘的模型和基于区域的模型。基于边缘的活动轮廓模型利用图像的梯度信息作为一种基于图像的力促使轮廓向目标边界移动，对于有清晰边界的目标有较好的分割效果。但由于梯度信息是一种高度局域化的图像信息，导致该模型有两个主要的缺点：对图像噪声的敏感性和对初始化轮廓的敏感性。将该模型应用于超声图像极易产生边界泄露的问题。基于区域的活动轮廓模型利用前景和背景区域的统计信息形成驱动力，适用于各区域灰度均匀的图像分割。其中最著名的方法是Chan and Vese提出的分段常数模型。该模型由于未使用图像的梯度信息，可以有效的分割含噪声和弱边界的图像，适用于超声图像等医学图像的分割。但由于该模型大都假定图像中各区域灰度是均匀的，并采用全局的统计信息，对于非均匀的图像容易产生错误的分割结果。近期，为了克服基于区域的活动轮廓模型难以分割不均匀目标的缺点，人们把更多注意力投向基于区域的局域化研究，提出了许多利用局域信息的新的基于区域活动轮廓模型。李纯明等人提出通过引入一个核函数来定义局域二元拟合能量在变分水平集框架中，从而嵌入局域的灰度信息到基于区域的活动轮廓模型中。之后，李纯明等人又在其基础上做了改进，深入研究了核函数的选择和基于区域范围大小的选择。S.Lankton等人则提出了一个局域化框架，允许利用全局信息的基于区域的能量式可以重新改写为局域化的形式，并深入分析了局域半径对于分割结果的影响。这些模型都引入了局域信息，更加适合于超声图像的分割难题，对灰度非均匀的图像已经有了较好的分割能力，但对于具有严重的低信噪比，低对比度，弱边界等问题的子宫肌瘤超声图像仍较易出现边界泄露或过度收缩的错误分割。

发明内容

针对上述子宫肌瘤超声图像分割较易出现的边界泄露或过度收缩等问题，本发明的目的是提供一种快速、准确的半自动分割方法应用于HIFU治疗中的子宫肌瘤超声图像。

本发明采用的技术方案为：一种HIFU治疗中子宫肌瘤超声图像分割方法，包括以下步骤：

步骤1：读取原始图像I₀，对所述的原始图像I₀进行粗尺度分割，其具体实现包括以下子步骤：

步骤1.1：初始化目标轮廓，在所述的原始图像I₀中划取一个椭圆作为原始图像I₀的初始轮廓C₀，使其能覆盖原始图像I₀中肿瘤的边缘轮廓；

步骤1.2：设所述的原始图像I₀大小为M×N，根据所述的原始图像I₀构造一个图像的高斯金字塔，得到图像大小为(M/4)×(N/4)的粗尺度图像I₂；

步骤1.3：将所述的原始图像I₀的初始轮廓C₀缩小到所述的粗尺度图像I₂大小，记为C₁，求解基于C₁的局域化区域活动轮廓模型，构建形状约束能量E_shape，并将所述的E_shape嵌入到所述的C₁的局域化区域活动轮廓模型中，得到基于C₁的总能量函数；

步骤1.4：通过变分法使所述的C₁的总能量函数最小化得到相应的水平集演化方程，通过该方程更新C₁的水平集函数，该水平集函数初始化为符号距离函数；

步骤1.5：判断所述的C₁的水平集函数是否收敛？

若否：则回转执行所述的步骤1.4，通过C₁的水平集演化方程不断更新C₁的水平集函数；

如是：则执行下述步骤1.6；

步骤1.6：判断C₁的水平集函数属于哪个尺度分割阶段？

若是细分割阶段，则得到细尺度分割轮廓，该细尺度分割轮廓作为最后的肿瘤分割轮廓，本方法执行完毕；

若是粗分割阶段，则由C₁的水平集函数的零水平集位置得到粗尺度图像I₂的粗尺度分割轮廓，转入细尺度分割阶段；

步骤2：细尺度分割，其具体实现包括以下子步骤：

步骤2.1：将所述的粗尺度图像I₂的分割轮廓的大小插值放大到对应原始图像I₀的大小，作为对原始图像I₀进行分割的初始轮廓C₂；

步骤2.2：求解基于C₂的局域化区域活动轮廓模型，利用所述的C₂构建相应的形状约束能量，将其嵌入到所述的C₂的局域化区域活动轮廓模型中，得到C₂的总能量函数；

步骤2.3：通过变分法得到所述的C₂的水平集演化方程，再通过该方程更新C₂的水平集函数；

步骤2.4：判断所述的C₂的水平集函数是否收敛？

若否：则回转执行所述的步骤2.3，通过C₂的水平集演化方程不断更新C₂的水平集函数；

如是：则执行下述步骤2.5；

步骤2.5：判断C₂的水平集函数属于哪个尺度分割阶段？

若是细分割阶段，则得到细尺度分割轮廓，该细尺度分割轮廓作为最后的肿瘤分割轮廓，本方法执行完毕；

若是粗分割阶段，则回转执行所述的步骤2.1；

作为优选，步骤1.3中所述的求解基于C₁的局域化区域活动轮廓模型，其具体实现过程为：对C₁曲线上的每一个点为单位独自计算其局域能量，为定义C₁曲线上每个点的局域区域，定义一个特征函数如下：

$B(x,y)=\left\{\begin{array}{cc}1,& \left|\right|x-y\left|\right|<r\\ 0,& \mathrm{otherwise}.\end{array}\right.$

其中x,y∈Ω，作为独立的空间变量各自代表一个点，r则代表半径参数，当点y在以x为中心，半径大小为r的圆内时，该特征函数Β(x,y)的值为1，否则为0；

根据特征函数Β(x,y)求出C₁上点x局域化下的区域内部的均值强度c_x1和区域外部的均值强度c_x2，得到局域化框架：

$c_{x1}\left(\mathrm{\&phi;}\right)=\frac{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;H\left(\mathrm{\&phi;}\left(y\right)\right)\&CenterDot;I\left(y\right)\mathrm{dy}}{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;H\left(\mathrm{\&phi;}\left(y\right)\right)\mathrm{dy}},$

$c_{x2}\left(\mathrm{\&phi;}\right)=\frac{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;(1-H\left(\mathrm{\&phi;}\left(y\right)\right))\&CenterDot;I\left(y\right)\mathrm{dy}}{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;(1-H\left(\mathrm{\&phi;}\left(y\right)\right))\mathrm{dy}},$

这里H(φ)是Heaviside函数，I(y)是局域范围内点y的灰度值；

将该局域化框架应用到C-V模型中，获得能量函数和C₁上点x相应的曲率流的局域化区域活动轮廓模型，简称LCV模型：

$E^{\mathrm{LCV}}(c_{x1},c_{x2},\mathrm{\&phi;})={\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right){\&Integral;}_{{\mathrm{\&Omega;}}_y}B(x,y)\&CenterDot;F_{\mathrm{region}}(I\left(y\right),\mathrm{\&phi;}\left(y\right))\mathrm{dydx}+\mathrm{\&mu;}{\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right)\left|\right|\&dtri;\mathrm{\&phi;}\left(x\right)\left|\right|\mathrm{dx}$

其中F_region=Hφ(y)(I(y)-c_x1)²+(1-Hφ(y))(I(y)-c_x2)²，δ(φ)是Dirac函数，也是H(φ)的导数，φ(x)是水平集函数，B(x,y)是表示局域范围的特征函数，F_region是基于面积的作用力，x表示整个图像中全局的一个像素点，y表示局域范围圆内的一个像素点，参数μ表示弧长项的权值，决定着曲线的平滑性。

作为优选，所述的特征函数Β(x,y)中，通过局域半径自适应选择函数R(x)，计算相应局域半径参数r，其中：

R(x)=10×arctan(0.28x-6)+24

x=k(||x_max-x_min||+||y_max-y_min||)

其中x是C₁在图像上x轴方向的最大值和最小值之差与y轴方向上最大值和最小值之差的和的一个比例值，k是控制比例值大小的系数，x_max,y_max分别是C₁在图像x,y轴方向上的最大值，x_min,y_min分别是C₁在图像x,y轴方向上的最小值。

作为优选，所述的k取值为0.25。

作为优选，步骤1.3中所述的构建形状约束能量E_shape，并将所述的E_shape嵌入到所述的C₁的局域化区域活动轮廓模型中，得到C₁的总能量函数，其总能量函数如下：

$E\left(\mathrm{\&phi;}\right)={\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right){\&Integral;}_{{\mathrm{\&Omega;}}_y}B(x,y)\&CenterDot;F_{\mathrm{region}}(I\left(y\right),\mathrm{\&phi;}\left(y\right))\mathrm{dydx}+\mathrm{\&mu;}{\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right)\left|\right|\&dtri;\mathrm{\&phi;}\left(x\right)\left|\right|\mathrm{dx}+E_{\mathrm{shape}}$

其中，

$E_{\mathrm{shape}}=\mathrm{\&beta;}{\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right)\&CenterDot;\mathrm{sign}({\mathrm{\&phi;}}_0,x){\left|\right|p_x-p_{\min }\left|\right|}^2\mathrm{dx}$

$\mathrm{sign}\left({\mathrm{\&phi;}}_{0}x\right)=\left\{\begin{array}{cc}1,& {\mathrm{\&phi;}}_0\left(x\right)>0\\ -1,& {\mathrm{\&phi;}}_0\left(x\right)<0\\ 0,& {\mathrm{\&phi;}}_0\left(x\right)=0.\end{array}\right.$

p_x是当前轮廓上点x在图像中的位置，p_min是初始轮廓C₁上到点x的最近的点的位置，φ(x)是水平集函数，φ₀(x)是初始轮廓C₁下的水平集函数表示，I(y)是图像在局域范围内点y的灰度值，参数β表示形状约束项的权值，决定着形状约束力的大小，该形状约束能力用来形成演化曲线在演化过程中朝向初始轮廓方向的作用力。

作为优选，步骤1.4中所述的通过变分法使所述的C₁的总能量函数最小化得到相应的水平集演化方程，其水平集演化方程为：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}\left(x\right)=\mathrm{\&delta;\&phi;}\left(x\right)[{\&Integral;}_{{\mathrm{\&Omega;}}_y}B(x,y)\mathrm{\&delta;\&phi;}\left(y\right)\&CenterDot;({(I\left(y\right)-c_{x1})}^2-{(I\left(y\right)-c_{x2})}^2)\mathrm{dy}+\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}\left(x\right)}{|\&dtri;\mathrm{\&phi;}\left(x\right)|}\right)+\mathrm{\&beta;sign}\left({\mathrm{\&phi;}}_{0}x\right){\left|\right|p_x-p_{\min }\left|\right|}^2$

其中，p_x是当前轮廓上点x在图像中的位置，p_min是初始轮廓C₁上到点x的最近的点的位置，φ(x)是水平集函数，φ₀(x)是初始轮廓C₁下的水平集函数表示，Β(x,y)是表示局域范围的特征函数，x表示整个图像中全局的一个像素点，y表示局域范围圆内的一个像素点，I(y)是图像在局域范围内点y的灰度值，参数μ表示弧长项的权值，决定着曲线的平滑性，参数β表示形状约束项的权值，决定着形状约束力的大小，δ(φ)是Dirac函数，c_x1，c_x2是根据特征函数Β(x,y)求出的轮廓上点x局域区域中轮廓内部和外部的均值强度。

作为优选，步骤2.2中所述的C₂的总能量函数的计算方法与步骤1.3中所述的C₁的总能量函数的计算方法相同。

作为优选，步骤2.3中所述的C₂的水平集演化方程的计算方法与步骤1.4中所述的C₁的水平集演化方程的计算方法相同。

与现有技术相比，本发明具有以下的有益效果：

本发明通过引入一种形状约束到局域化的基于区域的活动轮廓模型中，约束曲线演化获得一个更准确的分割结果，有效的避免了由于超声图像本身特性所导致的边界泄露或过度收缩等问题，为HIFU治疗中的子宫肌瘤超声图像提供了准确的分割结果。并通过多尺度分割算法的引入，极大的提高了分割效率，而且对局域化范围选择提出了有效的解决方案，实现了局域半径大小的自适应选取。

附图说明

图1：是本发明实施例的方法流程图。

图2：是本发明实施例中求解局域半径大小的自适应选择函数R(x)图像。

图3：是本发明实施例对合成图片的分割结果图。

图4-1：是本发明实施例未使用形状约束的LCV方法对HIFU中子宫肌瘤超声图像的分割结果图。

图4-2：是本发明实施例对HIFU中子宫肌瘤超声图像的分割结果图。

具体实施方式

以下结合附图和具体实施例对本发明做进一步的阐述。

请见图1，本发明所采用的技术方案是：一种HIFU治疗中子宫肌瘤超声图像分割方法，包括以下步骤：

步骤1：读取原始图像I₀，对原始图像I₀进行粗尺度分割，其具体实现包括以下子步骤：

步骤1.1：初始化目标轮廓，在原始图像I₀中划取一个椭圆作为原始图像I₀的初始轮廓C₀，使其能覆盖原始图像I₀中肿瘤的边缘轮廓；选择椭圆初始化是因为子宫肌瘤形状大都近似椭圆形状，更易获得较好的初始化效果；

步骤1.2：设原始图像I₀大小为M×N，根据原始图像I₀构造一个图像的高斯金字塔，得到图像大小为(M/4)×(N/4)的粗尺度图像I₂；

步骤1.3：将所述的原始图像I₀的初始轮廓C₀缩小到所述的粗尺度图像I₂大小，记为C₁，求解基于C₁的局域化区域活动轮廓模型，构建形状约束能量E_shape，并将所述的E_shape嵌入到所述的C₁的局域化区域活动轮廓模型中，得到基于C₁的总能量函数；

求解基于C₁的局域化区域活动轮廓模型，其具体实现过程为：对C₁曲线上的每一个点为单位独自计算其局域能量，为定义C₁曲线上每个点的局域区域，定义一个特征函数如下：

$B(x,y)=\left\{\begin{array}{cc}1,& \left|\right|x-y\left|\right|<r\\ 0,& \mathrm{otherwise}.\end{array}\right.$

其中x,y∈Ω，作为独立的空间变量各自代表一个点，r则代表半径参数，当点y在以x为中心，半径大小为r的圆内时，该特征函数Β(x,y)的值为1，否则为0；

根据特征函数Β(x,y)求出C₁上点x局域化下的区域内部的均值强度c_x1和区域外部的均值强度c_x2，得到局域化框架：

$c_{x1}\left(\mathrm{\&phi;}\right)=\frac{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;H\left(\mathrm{\&phi;}\left(y\right)\right)\&CenterDot;I\left(y\right)\mathrm{dy}}{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;H\left(\mathrm{\&phi;}\left(y\right)\right)\mathrm{dy}},$

$c_{x2}\left(\mathrm{\&phi;}\right)=\frac{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;(1-H\left(\mathrm{\&phi;}\left(y\right)\right))\&CenterDot;I\left(y\right)\mathrm{dy}}{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;(1-H\left(\mathrm{\&phi;}\left(y\right)\right))\mathrm{dy}},$

这里H(φ)是Heaviside函数，I(y)是局域范围内点y的灰度值；

将该局域化框架应用到C-V模型中，获得能量函数和C₁上点x相应的曲率流的局域化区域活动轮廓模型：

$E^{\mathrm{LCV}}(c_{x1},c_{x2},\mathrm{\&phi;})={\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right){\&Integral;}_{{\mathrm{\&Omega;}}_y}B(x,y)\&CenterDot;F_{\mathrm{region}}(I\left(y\right),\mathrm{\&phi;}\left(y\right))\mathrm{dydx}+\mathrm{\&mu;}{\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right)\left|\right|\&dtri;\mathrm{\&phi;}\left(x\right)\left|\right|\mathrm{dx}$

其中F_region=Hφ(y)(I(y)-c_x1)²+(1-Hφ(y))(I(y)-c_x2)²，δ(φ)是Dirac函数，也是H(φ)的导数，φ(x)是水平集函数，B(x,y)是表示局域范围的特征函数，F_region是基于面积的作用力，x表示整个图像中全局的一个像素点，y表示局域范围圆内的一个像素点，参数μ表示弧长项的权值，决定着曲线的平滑性。

特征函数Β(x,y)中，通过局域半径自适应选择函数R(x)，计算相应局域半径参数r，其中：

R(x)=10×arctan(0.28x-6)+24

x=k(||x_max-x_min||+||y_max-y_min||)

其中x是C₁在图像上x轴方向的最大值和最小值之差与y轴方向上最大值和最小值之差的和的一个比例值，k是控制比例值大小的系数，k取值为0.25，x_max,y_max分别是C₁在图像x,y轴方向上的最大值，x_min,y_min分别是C₁在图像x,y轴方向上的最小值。请见图2，根据子宫肌瘤图像中子宫肌瘤大小和实验结果，使由R(x)函数的求得的局域半径大小控制在10到40之间，有效地避免了当分割目标太大或太小时造成的局域半径过大或过小问题。

构建形状约束能量E_shape，并将所述的E_shape嵌入到所述的C₁的局域化区域活动轮廓模型中，得到C₁的总能量函数如下：

$E\left(\mathrm{\&phi;}\right)={\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right){\&Integral;}_{{\mathrm{\&Omega;}}_y}B(x,y)\&CenterDot;F_{\mathrm{region}}(I\left(y\right),\mathrm{\&phi;}\left(y\right))\mathrm{dydx}+\mathrm{\&mu;}{\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right)\left|\right|\&dtri;\mathrm{\&phi;}\left(x\right)\left|\right|\mathrm{dx}+E_{\mathrm{shape}}$

其中，

$E_{\mathrm{shape}}=\mathrm{\&beta;}{\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right)\&CenterDot;\mathrm{sign}({\mathrm{\&phi;}}_0,x){\left|\right|p_x-p_{\min }\left|\right|}^2\mathrm{dx}$

$\mathrm{sign}\left({\mathrm{\&phi;}}_{0}x\right)=\left\{\begin{array}{cc}1,& {\mathrm{\&phi;}}_0\left(x\right)>0\\ -1,& {\mathrm{\&phi;}}_0\left(x\right)<0\\ 0,& {\mathrm{\&phi;}}_0\left(x\right)=0.\end{array}\right.$

p_x是当前轮廓上点x在图像中的位置，p_min是初始轮廓C₁上到点x的最近的点的位置，φ(x)是水平集函数，φ₀(x)是初始轮廓C₁下的水平集函数表示，I(y)是图像在局域范围内点y的灰度值，参数β表示形状约束项的权值，决定着形状约束力的大小，该形状约束能力用来形成演化曲线在演化过程中朝向初始轮廓方向的作用力。

步骤1.4：通过变分法使C₁的总能量函数最小化得到相应的水平集演化方程，通过该方程更新C₁的水平集函数，该水平集函数初始化为符号距离函数；

其水平集演化方程为：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}\left(x\right)=\mathrm{\&delta;\&phi;}\left(x\right)[{\&Integral;}_{{\mathrm{\&Omega;}}_y}B(x,y)\mathrm{\&delta;\&phi;}\left(y\right)\&CenterDot;({(I\left(y\right)-c_{x1})}^2-{(I\left(y\right)-c_{x2})}^2)\mathrm{dy}+\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}\left(x\right)}{|\&dtri;\mathrm{\&phi;}\left(x\right)|}\right)+\mathrm{\&beta;sign}\left({\mathrm{\&phi;}}_{0}x\right){\left|\right|p_x-p_{\min }\left|\right|}^2$

其中，p_x是当前轮廓上点x在图像中的位置，p_min是初始轮廓C₁上到点x的最近的点的位置，φ(x)是水平集函数，φ0(x)是初始轮廓C₁下的水平集函数表示，Β(x,y)是表示局域范围的特征函数，x表示整个图像中全局的一个像素点，y表示局域范围圆内的一个像素点，I(y)是图像在局域范围内点y的灰度值，参数μ表示弧长项的权值，决定着曲线的平滑性，参数β表示形状约束项的权值，决定着形状约束力的大小，δ(φ)是Dirac函数，c_x1，c_x2是根据特征函数Β(x,y)求出的轮廓上点x局域区域中轮廓内部和外部的均值强度。

步骤1.5：判断所述的C₁的水平集函数是否收敛？

若否：则回转执行所述的步骤1.4，通过C₁的水平集演化方程不断更新C₁的水平集函数；

如是：则执行下述步骤1.6；

步骤1.6：判断C₁的水平集函数属于哪个尺度分割阶段？

若是细分割阶段，则得到细尺度分割轮廓，该细尺度分割轮廓作为最后的肿瘤分割轮廓，本方法执行完毕；

若是粗分割阶段，则由C₁的水平集函数的零水平集位置得到粗尺度图像I₂的粗尺度分割轮廓，转入细尺度分割阶段；

步骤2：细尺度分割，其具体实现包括以下子步骤：

步骤2.1：将粗尺度图像I₂的分割轮廓的大小插值放大到对应原始图像I₀的大小，作为对原始图像I₀进行分割的初始轮廓C₂；

步骤2.2：利用C₂构建相应的形状约束能量，求解基于C₂的局域化区域活动轮廓模型，将其嵌入到C₂的局域化区域活动轮廓模型中，得到C₂的总能量函数；C₂的总能量函数的计算方法与步骤1.3中C₁的总能量函数的计算方法相同。

步骤2.3：通过变分法得到C₂的水平集演化方程，再通过该方程更新C₂的水平集函数；C₂的水平集演化方程的计算方法与步骤1.4中C₁的水平集演化方程的计算方法相同。

步骤2.4：判断所述的C₂的水平集函数是否收敛？

若否：则回转执行所述的步骤2.3，通过C₂的水平集演化方程不断更新C₂的水平集函数；

如是：则执行下述步骤2.5；

步骤2.5：判断C₂的水平集函数属于哪个尺度分割阶段？

若是细分割阶段，则得到细尺度分割轮廓，该细尺度分割轮廓作为最后的肿瘤分割轮廓，本方法执行完毕；

若是粗分割阶段，则回转执行所述的步骤2.1；

请见图3，显示了本发明对合成图像的分割结果，表明所引入形状约束的作用。(a),(d)显示初始化轮廓，(b),(e)显示了未使用形状约束的分割结果；(c),(f)显示了使用形状约束后的分割结果。从图可以看到引入形状约束后使得最后分割结果的大致轮廓依然保持为初始轮廓椭圆。

请见图4-1、图4-2，显示了未使用形状约束的LCV方法同引入形状约束后的本发明算法对HIFU治疗中子宫肌瘤超声图像分割结果的比较，未引入形状约束前分割轮廓发生了边界溢出现象，引入形状约束后分割轮廓较好的分割出了肿瘤区域，避免了边界的溢出。

为定量分析，表1给出分割准确性和时间的比较。Dice Similarity Coefficient(DSC)值越接近于1，表明分割结果越好。本文发明1表示未使用多尺度分割算法，只使用形状约束。本文发明2表示在使用形状约束的基础上使用多尺度分割算法。

表1

方法	迭代次数	DSC	时间(s)
				LCV方法	400	0.89	140.32
本文发明1	400	0.95	114.69
				本文发明2	400	0.95	24.43

综上所述，本发明通过引入形状约束到局域区域活动轮廓模型中，克服了超声图像由于图像质量问题较易导致的边界泄露或过度收缩等问题，较大的提高了分割的准确性，并通过多尺度分割算法极大的提高了分割效率。

上述实例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明做出任何的修改和改变，都落入本发明的保护范围。

Claims (8)

1.一种HIFU治疗中子宫肌瘤超声图像分割方法，其特征在于，包括以下步骤：

步骤1：读取原始图像I₀，对所述的原始图像I₀进行粗尺度分割，其具体实现包括以下子步骤：

步骤1.1：初始化目标轮廓，在所述的原始图像I₀中划取一个椭圆作为原始图像I₀的初始轮廓C₀，使其能覆盖原始图像I₀中肿瘤的边缘轮廓；

步骤1.2：设所述的原始图像I₀大小为M×N，根据所述的原始图像I₀构造一个图像的高斯金字塔，得到图像大小为(M/4)×(N/4)的粗尺度图像I₂；

步骤1.3：将所述的原始图像I₀的初始轮廓C₀缩小到所述的粗尺度图像I₂大小，记为C₁，求解基于C₁的局域化区域活动轮廓模型，构建形状约束能量E_shape，并将所述的E_shape嵌入到所述的C₁的局域化区域活动轮廓模型中，得到基于C₁的总能量函数；

步骤1.4：通过变分法使所述的C₁的总能量函数最小化得到相应的水平集演化方程，通过该方程更新C₁的水平集函数，该水平集函数初始化为符号距离函数；

步骤1.5：判断所述的C₁的水平集函数是否收敛？

若否：则回转执行所述的步骤1.4，通过C₁的水平集演化方程不断更新C₁的水平集函数；

如是：则执行下述步骤1.6；

步骤1.6：判断C₁的水平集函数属于哪个尺度分割阶段？

若是细分割阶段，则得到细尺度分割轮廓，该细尺度分割轮廓作为最后的肿瘤分割轮廓，本方法执行完毕；

若是粗分割阶段，则由C₁的水平集函数的零水平集位置得到粗尺度图像I₂的粗尺度分割轮廓，转入细尺度分割阶段；

步骤2：细尺度分割，其具体实现包括以下子步骤：

步骤2.1：将所述的粗尺度图像I₂的分割轮廓的大小插值放大到对应原始图像I₀的大小，作为对原始图像I₀进行分割的初始轮廓C₂；

步骤2.2：求解基于C₂的局域化区域活动轮廓模型，利用所述的C₂构建相应的形状约束能量，将其嵌入到所述的C₂的局域化区域活动轮廓模型中，得到C₂的总能量函数；

步骤2.3：通过变分法得到所述的C₂的水平集演化方程，再通过该方程更新C₂的水平集函数；

步骤2.4：判断所述的C₂的水平集函数是否收敛？

若否：则回转执行所述的步骤2.3，通过C₂的水平集演化方程不断更新C₂的水平集函数；

如是：则执行下述步骤2.5；

步骤2.5：判断C₂的水平集函数属于哪个尺度分割阶段？

若是细分割阶段，则得到细尺度分割轮廓，该细尺度分割轮廓作为最后的肿瘤分割轮廓，本方法执行完毕；

若是粗分割阶段，则回转执行所述的步骤2.1。

2.根据权利要求1所述的HIFU治疗中子宫肌瘤超声图像分割方法，其特征在于：步骤1.3中所述的求解基于C₁的局域化区域活动轮廓模型，其具体实现过程为：对C₁曲线上的每一个点为单位独自计算其局域能量，为定义C₁曲线上每个点的局域区域，定义一个特征函数如下：

其中x,y∈Ω，作为独立的空间变量各自代表一个点，r则代表半径参数，当点y在以x为中心，半径大小为r的圆内时，该特征函数Β(x,y)的值为1，否则为0；

根据特征函数Β(x,y)求出C₁上点x局域化下的区域内部的均值强度c_x1和区域外部的均值强度c_x2，得到局域化框架：

$c_{x1}\left(\mathrm{\&phi;}\right)=\frac{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;H\left(\mathrm{\&phi;}\left(y\right)\right)\&CenterDot;I\left(y\right)\mathrm{dy}}{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;H\left(\mathrm{\&phi;}\left(y\right)\right)\mathrm{dy}},$

$c_{x2}\left(\mathrm{\&phi;}\right)=\frac{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;(1-H\left(\mathrm{\&phi;}\left(y\right)\right))\&CenterDot;I\left(y\right)\mathrm{dy}}{{\&Integral;}_{\mathrm{\&Omega;y}}B(x,y)\&CenterDot;(1-H\left(\mathrm{\&phi;}\left(y\right)\right))\mathrm{dy}},$

这里H(φ)是Heaviside函数，I(y)是局域范围内点y的灰度值；

将该局域化框架应用到C-V模型中，获得能量函数和C₁上点x相应的曲率流的局域化区域活动轮廓模型：

$E^{\mathrm{LCV}}(c_{x1},c_{x2},\mathrm{\&phi;})={\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right){\&Integral;}_{{\mathrm{\&Omega;}}_y}B(x,y)\&CenterDot;F_{\mathrm{region}}(I\left(y\right),\mathrm{\&phi;}\left(y\right))\mathrm{dydx}+\mathrm{\&mu;}{\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right)\left|\right|\&dtri;\mathrm{\&phi;}\left(x\right)\left|\right|\mathrm{dx}$ 其中F_region=Hφ(y)(I(y)-c_x1)²+(1-Hφ(y))(I(y)-c_x2)²，δ(φ)是Dirac函数，也是H(φ)的导数，φ(x)是水平集函数，B(x,y)是表示局域范围的特征函数，F_region是基于面积的作用力，x表示整个图像中全局的一个像素点，y表示局域范围圆内的一个像素点，参数μ表示弧长项的权值，决定着曲线的平滑性。

3.根据权利要求2所述的HIFU治疗中子宫肌瘤超声图像分割方法，其特征在于：所述的特征函数Β(x,y)中，通过局域半径自适应选择函数R(x)，计算相应局域半径参数r，其中：

R(x)=10×arctan(0.28x-6)+24

x=k(||x_max-x_min||+||y_max-y_min||)

其中x是C₁在图像上x轴方向的最大值和最小值之差与y轴方向上最大值和最小值之差的和的一个比例值，k是控制比例值大小的系数，x_max,y_max分别是C₁在图像x,y轴方向上的最大值，x_min,y_min分别是C₁在图像x,y轴方向上的最小值。

4.根据权利要求3所述的HIFU治疗中子宫肌瘤超声图像分割方法，其特征在于：所述的k取值为0.25。

5.根据权利要求2所述的HIFU治疗中子宫肌瘤超声图像分割方法，其特征在于：步骤1.3中所述的构建形状约束能量E_shape，并将所述的E_shape嵌入到所述的C₁的局域化区域活动轮廓模型中，得到C₁的总能量函数，其总能量函数如下：

$E\left(\mathrm{\&phi;}\right)={\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right){\&Integral;}_{{\mathrm{\&Omega;}}_y}B(x,y)\&CenterDot;F_{\mathrm{region}}(I\left(y\right),\mathrm{\&phi;}\left(y\right))\mathrm{dydx}+\mathrm{\&mu;}{\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right)\left|\right|\&dtri;\mathrm{\&phi;}\left(x\right)\left|\right|\mathrm{dx}+E_{\mathrm{shape}}$

其中，

$E_{\mathrm{shape}}=\mathrm{\&beta;}{\&Integral;}_{{\mathrm{\&Omega;}}_x}\mathrm{\&delta;\&phi;}\left(x\right)\&CenterDot;\mathrm{sign}({\mathrm{\&phi;}}_0,x){\left|\right|p_x-p_{\min }\left|\right|}^2\mathrm{dx}$

$\mathrm{sign}\left({\mathrm{\&phi;}}_{0}x\right)=\left\{\begin{array}{cc}1,& {\mathrm{\&phi;}}_0\left(x\right)>0\\ -1,& {\mathrm{\&phi;}}_0\left(x\right)<0\\ 0,& {\mathrm{\&phi;}}_0\left(x\right)=0\end{array}\right.$

p_x是当前轮廓上点x在图像中的位置，p_min是初始轮廓C₁上到点x的最近的点的位置，φ(x)是水平集函数，φ₀(x)是初始轮廓C₁下的水平集函数表_示，I(y)是图像在局域范围内点y的灰度值，参数β表示形状约束项的权值，决定着形状约束力的大小，该形状约束能力用来形成演化曲线在演化过程中朝向初始轮廓方向的作用力。

6.根据权利要求1所述的HIFU治疗中子宫肌瘤超声图像分割方法，其特征在于：步骤1.4中所述的通过变分法使所述的C₁的总能量函数最小化得到相应的水平集演化方程，其水平集演化方程为：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}\left(x\right)=\mathrm{\&delta;\&phi;}\left(x\right)[{\&Integral;}_{{\mathrm{\&Omega;}}_y}B(x,y)\mathrm{\&delta;\&phi;}\left(y\right)\&CenterDot;({(I\left(y\right)-c_{x1})}^2-{(I\left(y\right)-c_{x2})}^2)\mathrm{dy}+\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}\left(x\right)}{|\&dtri;\mathrm{\&phi;}\left(x\right)|}\right)+\mathrm{\&beta;sign}\left({\mathrm{\&phi;}}_{0}x\right){\left|\right|p_x-p_{\min }\left|\right|}^2$

其中，p_x是当前轮廓上点x在图像中的位置，p_min是初始轮廓C₁上到点x的最近的点的位置，φ(x)是水平集函数，φ₀(x)是初始轮廓C₁下的水平集函数表示，Β(x,y)是表示局域范围的特征函数，x表示整个图像中全局的一个像素点，y表示局域范围圆内的一个像素点，I(y)是图像在局域范围内点y的灰度值，参数μ表示弧长项的权值，决定着曲线的平滑性，参数β表示形状约束项的权值，决定着形状约束力的大小，δ(φ)是Dirac函数，c_x1，c_x2是根据特征函数Β(x,y)求出的轮廓上点x局域区域中轮廓内部和外部的均值强度。

7.根据权利要求1所述的HIFU治疗中子宫肌瘤超声图像分割方法，其特征在于：步骤2.2中所述的C₂的总能量函数的计算方法与步骤1.3中所述的C₁的总能量函数的计算方法相同。

8.根据权利要求1所述的HIFU治疗中子宫肌瘤超声图像分割方法，其特征在于：步骤2.3中所述的C₂的水平集演化方程的计算方法与步骤1.4中所述的C₁的水平集演化方程的计算方法相同。

Images