CN102800100B - 基于距离势场和自适应气球力的图像分割方法 - Google Patents

Abstract

基于距离势场和自适应气球力的图像分割方法，它属于数字图像处理领域，本发明的目的是为了能够准确的分割长管状边缘，并且能够分割含椒盐噪声的图像。方法：一、构造距离势场算子；二、用构造出的距离势场算子求图像的距离势场，对势场求梯度获得距离力场，并对力场归一化；三、设置初始轮廓，计算自适应气球力；四、将从图像数据中计算的距离力场和自适应气球力带入到轮廓的力学平衡方程，采用有限差分法求解该方程，得到的稳态解组成的轮廓即为分割结果。本发明的优点是抗噪性能较高且可以分割长管状的图像目标。

Description

基于距离势场和自适应气球力的图像分割方法

技术领域

本发明涉及基于Snake模型的图像分割方法，属于数字图像处理领域。

背景技术

Snake模型，又称为活动轮廓模型、蛇模型。活动轮廓的含义是首先在待分割目标周围手动或自动设置初始轮廓，并赋予初始轮廓能量。轮廓在模型自身的拓扑内力和图像的灰度数据产生的外力的作用下发生形变，像蛇一样活动，达到轮廓能量最小时，轮廓的所在位置即是待分割目标的边缘。

Snake模型是将图像分割问题转化为求能量泛函最小化的过程。活动轮廓突出优点是一旦设置初始轮廓，后续的轮廓演化不需要人为参与，自动化程度高，而采用力学原理，对轮廓本身定义了弹性力和刚性力，这样，无论分割目标的结构如何复杂，分割的最终轮廓都是光滑封闭的，这符合一般情况下的自然图像，尤其是医学图像的客观规律。Snake模型的因计算效率高，适合建模，而被广泛应用于图像分割、目标跟踪、模式识别等应用。

Snake模型根据轮廓的表达方式的不同分为两种：参数Snake模型和几何Snake模型。参数Snake模型的轮廓由控制点连接而成，轮廓的变形特性由弹性和刚性的参数控制。Snake模型的轮廓能量定义如下：

E_snake＝E_int+E_image

Eint为内能，控制轮廓的拓扑特性，Eimage为外能，控制轮廓收敛特性，外能从图像数据获得，计算模型不唯一，是影响分割结果的重要因素。最早的参数Snake模型采用公式作为外能，该模型力场捕获范围小，且不能分割凹型边缘。

为了克服传统方法的不足，已经提出了许多方法。其中梯度矢量流(GVF)被引入用来作为轮廓的外力场。GVF的外力场V(x，y)＝(u(x，y)，v(x，y))通过最小化能量泛函获得：

$\mathrm{\ϵ}=\∫\∫\mathrm{\μ}(u_x^2+u_y^2+v_x^2+v_y^2)+{|\▿f|}^2{|V-\▿f|}^2\mathrm{dxdy}$

式中u_x，u_y，v_x，v_y分别为u(x，y)，v(x，y)关于x方向和y方向的导数，f为图像的边缘图，μ为常数。

GVF模型虽然具有较大的捕获范围，能够收敛凹型边缘，但还存在三点问题：1、GVF模型的力场范围不能通过参数控制。2、GVF模型的力场对噪声敏感。3、GVF模型的力场不能与气球力结合，难以分割长管状边缘。

为了克服孤立噪点以及提高分割效率，气球力模型(balloon)被引入作为轮廓的动态外力，气球力通过加在轮廓控制点法线方向的恒定大小的作用力来实现。气球力主要存在问题是力的膨胀或收缩方向一旦设定，则不能改变，在分割长管状边缘时，容易造成边界泄露。

发明内容

本发明的目的是为了能够准确分割长管状边缘的目标，同时能够分割含椒盐噪声的图像，进而提出了一种基于距离势场和自适应气球力的图像分割方法。

本发明基于距离势场和自适应气球力的图像分割方法，实现该方法的步骤如下：

一、构造距离势场算子；

二、用构造出的距离势场算子求图像的距离势场，对势场求梯度获得距离力场，并对力场归一化；

三、设置初始轮廓，计算自适应气球力；

四、将步骤二、步骤三得到距离力场和自适应气球力带入到轮廓的力学平衡方程，采用有限差分法求解该方程，得到的稳态解组成的轮廓即为分割结果。

其中步骤一所述构造距离势场算子的方法为：

a、构造一个边长为2×R+1的方形零值矩阵L_R，在矩阵中构造一个半径为R的圆形模板，半径R的选择一般不超过待处理图像的边长的一半；

b、将圆形模板内的元素(i，j)的值设置为l(i，j)，l(i，j)通过公式获得，其中h，p为常数，一般取0.5≤h≤2，1≤p≤3。得到的方形矩阵L_R即为距离势场算子；

其中步骤二所述用构造出的距离势场算子求图像的距离势场，对势场求梯度获得距离力场，并对力场归一化的方法为：

a、用步骤一构造出的距离势场算子L_R与图像I(x，y)做卷积，得到图像的距离势场E_DPF(x，y)。

b、对距离势场E_DPF(x，y)分别在x方向和y方向求偏导，得到距离力场 $f_{\mathrm{DPF}}(x,y)=(f_{\mathrm{DPFx}}(x,y),f_{\mathrm{DPFy}}(x,y))=(\frac{\∂E_{\mathrm{DPF}}(x,y)}{\∂x},$ $\frac{\∂E_{\mathrm{DPF}}(x,y)}{\∂y}),$ 其中f_DPFx(x，y)和f_DPFy(x，y)分别为图像中(x，y)点沿x方向和y方向的距离力场。

c、采用公式对距离力场归一化；

其中步骤三所述设置初始轮廓，计算自适应气球力的方法为：

首先手动在图像中待分割目标周围选取点{c₁，c₂...c_n}作为初始轮廓点，其中c₁＝c_n，

c_i＝(x_i，y_i)，i＝1，2，...，n为所选点的坐标，c_i点的自适应气球力采用式

$n_{\mathrm{balloon}}(x_i,y_i)=(n_{\mathrm{balloonx}}(x_i,y_i),n_{\mathrm{balloony}}(x_i,y_i))$

$=\mathrm{sign}\left(\mathrm{\θ}\right)\·(\frac{y_{i+1}-y_{i-1}}{\sqrt{{(x_{i+1}-x_{i-1})}^2+{(y_{i+1}-y_{i-1})}^2}}-\frac{x_{i+1}-x_{i-1}}{\sqrt{{(x_{i+1}-x_{i-1})}^2+{(y_{i+1}-y_{i-1})}^2}})$

获得，其中θ为c_i处轮廓法线方向与该处距离力方向的夹角， $\mathrm{sign}\left(\mathrm{\&theta;}\right)=\left\{\begin{array}{cc}1& \mathrm{\&theta;}<\mathrm{\&tau;}\\ -1& \mathrm{\&theta;}\&GreaterEqual;\mathrm{\&tau;}\end{array}\right.,$ τ为设定的阈值常数，一般取π/4≤τ≤3π/4。

其中步骤四所述将步骤二、步骤三得到距离力场和自适应气球力带入到轮廓的力学平衡方程，采用有限差分法求解该方程，得到稳态解的方法为：

将距离力和自适应气球力带入如下的轮廓力学平衡方程：

α(c_i-c_i-1)-α(c_i+1-c_i)

+β(c_i-2-2c_i-1+c_i)-2β(c_i-1-2c_i+c_i+1)+β(c_i-2c_i+1+c_i+2)

+λ(f_DPFx(c_i)，f_DPFy(c_i))+κ(n_balloonx(c_i)，n_balloony(c_i))＝0

其中λ和κ为常数，表示距离力和气球力的权重，一般取1≤λ≤3，1≤κ≤3，一般情况下κ要小于λ，上述方程写成矩阵形式为：

$\left\{\begin{array}{c}\mathrm{Ax}+\mathrm{\&lambda;}{f}_{\mathrm{DPFx}}(x,y)+\mathrm{\&kappa;}{n}_{\mathrm{balloonx}}(x,y)=0\\ \mathrm{Ay}+\mathrm{\&lambda;}{f}_{\mathrm{DPFy}}(x,y)+\mathrm{\&kappa;}{n}_{\mathrm{balloony}}(x,y)=0\end{array}\right.$

其中A为五对角带状矩阵：

$A=\left[\begin{array}{cccccccc}2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0& \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})\\ -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0& \mathrm{\&beta;}\\ \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0\\ 0& \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& ...& 0\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& .& .& .& .& .& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}\end{array}\right]$

引入时间辅助变量t，将静态偏微分方程转化为动态方程。将右式的0用曲线关于时间间隔Δt的一阶导数来替换，离散表达为

式中t——为迭代的次数；

Δt——第t次迭代和第t-1次迭代的时间间隔；

(x^t，y^t)——经过第t次迭代的轮廓曲线位置坐标；

(x^t-1，y^t-1)——经过第t-1迭次的轮廓曲线位置坐标。

假设在一次的迭代过程中距离力场和自适应气球力保持恒定，带入方程组中得到

Ax^t+λf_DPFx(x^t-1，y^t-1)+κn_balloonx(x^t-1，y^t-1)＝-γ(x^t-x^t-1)

Ay^t+λf_DPFy(x^t-1，y^t-1)+κn_balloony(x^t-1，y^t-1)＝-γ(y^t-y^t-1)

其中为常数，表示迭代步长。距离力场是分布在整个图像域中的，任意点的距离力f_DPF(x^t，y^t)可以通过对附近点的力线性插值获得。得到x，y的迭代方程：

x^t＝(A+γI)^-1(x^t-1-λf_dpfx(x^t-1，y^t-1)-κn_balloonx(x^t-1，y^t-1))

y^t＝(A+γI)^-1(y^t-1-λf_dpfy(x^t-1，y^t-1)-κn_balloony(x^t-1，y^t-1))

A+γI是带状矩阵，I为单位阵，它的逆可以通过LU分解获得；得到的稳态解组成的轮廓所包围的区域即为分割结果。

本发明使用归一化距离势场的力场作为轮廓外力场，具有可控的捕获范围，并且该力场对噪声不敏感，同时加入的自适应气球力可以自适应根据图像力场方向而改变方向，具有良好的边界特性，可以有效防止边界泄露，因此本发明具有能够分割长管状边缘的目标和抗噪性好的优点。

附图说明

图1为具体实施方式一中加噪声的U型图；

图2为具体实施方式一中加噪声的U型图的距离力场图；

图3为具体实施方式一中在加噪声U型图中设置的初始轮廓；

图4为具体实施方式一中基于距离势场和自适应气球力模型的分割结果；

图5为具体实施方式一中GVF模型对加噪声U型图的分割结果；

图6为具体实施方式二中长管状目标及初始轮廓；

图7为具体实施方式二中基于距离势场和自适应气球力模型的分割结果。

图8为具体实施方式二中传统气球力模型的分割结果。

图9为具体实施方式二中GVF模型的分割结果。

具体实施方式

具体实施方式一：本实施方式基于距离势场和自适应气球力的图像分割方法，以加有椒盐噪声的U型图作为仿真数据，实现该方法的步骤如下：

一、构造距离势场算子；构造方法为：

a、构造一个边长为2×R+1的方形零值矩阵L_R，在矩阵中构造一个半径为R的圆形模板，半径R的选择一般不超过待处理图像的边长的一半；

本实施例中选取R＝30；

b、将圆形模板内的元素(i，j)的值设置为l(i，j)，l(i，j)通过公式获得，其中h，p为常数，一般取0.5≤h≤2，1≤p≤3。得到的方形矩阵L_R即为距离势场算子；

本实例中选取h＝0.5，p＝2；

二、用构造出的距离势场算子求图像的距离势场，对势场求梯度获得距离力场，并对力场归一化的方法为：

a、用步骤一构造出的距离势场算子L_R与图像I(x，y)做卷积，得到图像的距离势场E_DPF(x，y)。

b、对距离势场E_DPF(x，y)分别在x方向和y方向求偏导，得到距离力场 $f_{\mathrm{DPF}}(x,y)=(f_{\mathrm{DPFx}}(x,y),f_{\mathrm{DPFy}}(x,y))=(\frac{\&PartialD;E_{\mathrm{DPF}}(x,y)}{\&PartialD;x},\frac{\&PartialD;E_{\mathrm{DPF}}(x,y)}{\&PartialD;y}),$ 其中f_DPFx(x，y)和f_DPFy(x，y)分别为图像中(x，y)点沿x方向和y方向的距离力场。

c、采用公式对距离力场归一化；

三、设置初始轮廓，计算自适应气球力的方法为：

首先手动在图像中待分割目标周围选取点{c₁，c₂...c_n}作为初始轮廓点，其中c₁＝c_n，

c_i＝(x_i，y_i)，i＝1，2，...，n为所选点的坐标，c_i点的自适应气球力采用式

$n_{\mathrm{balloon}}(x_i,y_i)=(n_{\mathrm{balloonx}}(x_i,y_i),n_{\mathrm{balloony}}(x_i,y_i))$

$=\mathrm{sign}\left(\mathrm{\&theta;}\right)\&CenterDot;(\frac{y_{i+1}-y_{i-1}}{\sqrt{{(x_{i+1}-x_{i-1})}^2+{(y_{i+1}-y_{i-1})}^2}}-\frac{x_{i+1}-x_{i-1}}{\sqrt{{(x_{i+1}-x_{i-1})}^2+{(y_{i+1}-y_{i-1})}^2}})$

获得，其中θ为c_i处轮廓法线方向与该处距离力方向的夹角， $\mathrm{sign}\left(\mathrm{\&theta;}\right)=\left\{\begin{array}{cc}1& \mathrm{\&theta;}<\mathrm{\&tau;}\\ -1& \mathrm{\&theta;}\&GreaterEqual;\mathrm{\&tau;}\end{array}\right.,$ τ为设定的阈值常数，一般取π/4≤τ≤3π/4。

本实施例中初始点的选取首先人工在周围选择12个点，在采用插值方法，每两个点之间插入5个点，总共72个点连接成初始轮廓。阈值τ取π/2。

四、将步骤二、步骤三得到距离力场和自适应气球力带入到轮廓的力学平衡方程，采用有限差分法求解该方程，得到稳态解的方法为：

将距离力和自适应气球力带入如下的轮廓力学平衡方程：

α(c_i-c_i-1)-α(c_i+1-c_i)

+β(c_i-2-2c_i-1+c_i)-2β(c_i-1-2c_i+c_i+1)+β(c_i-2c_i+1+c_i+2)

+λ(f_DPFx(c_i)，f_DPFy(c_i))+κ(n_balloonx(c_i)，n_balloony(c_i))＝0

其中λ和κ为常数，表示距离力和气球力的权重，一般取1≤λ≤3，1≤κ≤3，一般情况下κ要小于λ，上述方程写成矩阵形式为：

$\left\{\begin{array}{c}\mathrm{Ax}+\mathrm{\&lambda;}{f}_{\mathrm{DPFx}}(x,y)+\mathrm{\&kappa;}{n}_{\mathrm{balloonx}}(x,y)=0\\ \mathrm{Ay}+\mathrm{\&lambda;}{f}_{\mathrm{DPFy}}(x,y)+\mathrm{\&kappa;}{n}_{\mathrm{balloony}}(x,y)=0\end{array}\right.$

其中A为五对角带状矩阵：

$A=\left[\begin{array}{cccccccc}2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0& \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})\\ -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0& \mathrm{\&beta;}\\ \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0\\ 0& \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& ...& 0\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& .& .& .& .& .& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}\end{array}\right]$

引入时间辅助变量t，将静态偏微分方程转化为动态方程。将右式的0用曲线关于时间间隔Δt的一阶导数来替换，离散表达为

式中t——为迭代的次数；

Δt——第t次迭代和第t-1次迭代的时间间隔；

(x^t，y^t)——经过第t次迭代的轮廓曲线位置坐标；

(x^t-1，y^t-1)——经过第t-1迭次的轮廓曲线位置坐标。

假设在一次的迭代过程中距离力场和自适应气球力保持恒定，带入方程组中得到

Ax^t+λf_DPFX(x^t-1，y^t-1)+κn_balloonx(x^t-1，y^t-1)＝-γ(x^t-x^t-1)

Ay^t+λf_DPFy(x^t-1，y^t-1)+κn_balloony(x^t-1，y^t-1)＝-γ(y^t-y^t-1)

其中为常数，表示迭代步长。距离力场是分布在整个图像域中的，任意点的距离力f_DPF(x^t，y^t)可以通过对附近点的力线性插值获得。得到x，y的迭代方程：

x^t＝(A+γI)^-1(x^t-1-λf_dpfx(x^t-1，y^t-1)-κn_balloonx(x^t-1，y^t-1))

y^t＝(A+γI)^-1(y^t-1-λf_dpfy(x^t-1，y^t-1)-κn_balloony(x^t-1，y^t-1))

A+γI是带状矩阵，I为单位阵，它的逆可以通过LU分解获得。得到的稳态解组成的轮廓所包围的区域即为分割结果。

本实施例中α取0.4，β取0.3，γ取1，λ取2，κ取1.5。

本实施例是以加噪声的U型图为例(见图1)；距离势场的力场分布、初始轮廓见图2、3所示。

本实施例最终得到的分割结果如图4所示，图4中轮廓线即为U型边缘的分割结果，可以看出本实施例的方法可以准确的分割含椒盐噪声的图像，比采用GVF方法的分割结果(见图5)精度更高。

具体实施方式二：本实施方式基于距离势场和自适应气球力的图像分割方法，以真实的包含长管状目标的图像为例详述：

实现该方法的步骤如下：

一、构造距离势场算子；构造方法为：

a、构造一个边长为2×R+1的方形零值矩阵L_R，在矩阵中构造一个半径为R的圆形模板，半径R的选择一般不超过待处理图像的边长的一半；

本实施例中选取R＝100；

b、将圆形模板内的元素(i，j)的值设置为l(i，j)，l(i，j)通过公式获得，其中h，p为常数，一般取0.5≤h≤2，1≤p≤3。得到的方形矩阵L_R即为距离势场算子；

本实例中选取h＝1，p＝2；

二、用构造出的距离势场算子求图像的距离势场，对势场求梯度获得距离力场，并对力场归一化的方法为：

a、用步骤一构造出的距离势场算子L_R与图像I(x，y)做卷积，得到图像的距离势场E_DPF(x，y)。

b、对距离势场E_DPF(x，y)分别在x方向和y方向求偏导，得到距离力场 $f_{\mathrm{DPF}}(x,y)=(f_{\mathrm{DPFx}}(x,y),f_{\mathrm{DPFy}}(x,y))=(\frac{\&PartialD;E_{\mathrm{DPF}}(x,y)}{\&PartialD;x},\frac{\&PartialD;E_{\mathrm{DPF}}(x,y)}{\&PartialD;y}),$ 其中f_DPFx(x，y)和f_DPFy(x，y)分别为图像中(x，y)点沿x方向和y方向的距离力场。

c、采用公式对距离力场归一化；

三、设置初始轮廓，计算自适应气球力的方法为：

首先手动在图像中待分割目标周围选取点{c₁，c₂...c_n}作为初始轮廓点，其中c₁＝c_n，

c_i＝(x_i，y_i)，i＝1，2，...，n为所选点的坐标，c_i点的自适应气球力采用式

$n_{\mathrm{balloon}}(x_i,y_i)=(n_{\mathrm{balloonx}}(x_i,y_i),n_{\mathrm{balloony}}(x_i,y_i))$

$=\mathrm{sign}\left(\mathrm{\&theta;}\right)\&CenterDot;(\frac{y_{i+1}-y_{i-1}}{\sqrt{{(x_{i+1}-x_{i-1})}^2+{(y_{i+1}-y_{i-1})}^2}}-\frac{x_{i+1}-x_{i-1}}{\sqrt{{(x_{i+1}-x_{i-1})}^2+{(y_{i+1}-y_{i-1})}^2}})$

获得，其中θ为c_i处轮廓法线方向与该处距离力方向的夹角， $\mathrm{sign}\left(\mathrm{\&theta;}\right)=\left\{\begin{array}{cc}1& \mathrm{\&theta;}<\mathrm{\&tau;}\\ -1& \mathrm{\&theta;}\&GreaterEqual;\mathrm{\&tau;}\end{array}\right.,$ τ为设定的阈值常数，一般取π/4≤τ≤3π/4。

本实施例中初始点的选取首先人工在周围选择7个点，在采用插值方法，每两个点之间插入5个点，总共42个点连接成初始轮廓。阈值τ取5π/8。

四、将步骤二、步骤三得到距离力场和自适应气球力带入到轮廓的力学平衡方程，采用有限差分法求解该方程，得到稳态解的方法为：

将距离力和自适应气球力带入如下的轮廓力学平衡方程：

α(c_i-c_i-1)-α(c_i+1-c_i)

+β(c_i-2-2c_i-1+c_i)-2β(c_i-1-2c_i+c_i+1)+β(c_i-2c_i+1+c_i+2)

+λ(f_DPFx(c_i)，f_DPFy(c_i))+κ(n_balloonx(c_i)，n_balloony(c_i))＝0

其中λ和κ为常数，表示距离力和气球力的权重，一般取1≤λ≤3，1≤κ≤3，一般情况下κ要小于λ，上述方程写成矩阵形式为：

$\left\{\begin{array}{c}\mathrm{Ax}+\mathrm{\&lambda;}{f}_{\mathrm{DPFx}}(x,y)+\mathrm{\&kappa;}{n}_{\mathrm{balloonx}}(x,y)=0\\ \mathrm{Ay}+\mathrm{\&lambda;}{f}_{\mathrm{DPFy}}(x,y)+\mathrm{\&kappa;}{n}_{\mathrm{balloony}}(x,y)=0\end{array}\right.$

其中A为五对角带状矩阵：

$A=\left[\begin{array}{cccccccc}2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0& \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})\\ -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0& \mathrm{\&beta;}\\ \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0\\ 0& \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& ...& 0\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& .& .& .& .& .& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}\end{array}\right]$

引入时间辅助变量t，将静态偏微分方程转化为动态方程。将右式的0用曲线关于时间间隔Δt的一阶导数来替换，离散表达为

式中t——为迭代的次数；

Δt——第t次迭代和第t-1次迭代的时间间隔；

(x^t，y^t)——经过第t次迭代的轮廓曲线位置坐标；

(x^t-1，y^t-1)——经过第t-1迭次的轮廓曲线位置坐标。

假设在一次的迭代过程中距离力场和自适应气球力保持恒定，带入方程组中得到

Ax^t+λf_DPFx(x^t-1，y^t-1)+κn_balloonx(x^t-1，y^t-1)＝-γ(x^t-x^t-1)

Ay^t+λf_DPFy(x^t-1，y^t-1)+κn_balloony(x^t-1，y^t-1)＝-γ(y^t-y^t-1)

其中为常数，表示迭代步长。距离力场是分布在整个图像域中的，任意点的距离力f_DPF(x^t，y^t)可以通过对附近点的力线性插值获得。得到x，y的迭代方程：

x^t＝(A+γI)^-1(x^t-1-λf_dpfx(x^t-1，y^t-1)-κn_balloonx(x^t-1，y^t-1))

y^t＝(A+γI)^-1(y^t-1-λf_dpfy(x^t-1，y^t-1)-κn_balloony(x^t-1，y^t-1))

A+γI是带状矩阵，I为单位阵，它的逆可以通过LU分解获得。得到的稳态解组成的轮廓所包围的区域即为分割结果。

本实施例中α取0.4，β取0.3，γ取1，λ取2，κ取2.5。

本实施例是真实的长管状目标为例(见图6)。

本实施例最终得到的分割结果如图7所示。

可以看出本实施例的方法可以准确的分割长管状边缘，与传统气球力的分割结果(见图8)和GVF模型的分割结果(见图9)对比，准确度更高。

Claims (2)

1.基于距离势场和自适应气球力的图像分割方法，其特征在于，实现该方法的步骤如下：

一、构造距离势场算子；

二、用构造出的距离势场算子求图像的距离势场，对势场求梯度获得距离力场，并对力场归一化；

三、设置初始轮廓，计算自适应气球力；

四、将步骤二、步骤三得到的距离力场和自适应气球力带入到轮廓的力学平衡方程，采用有限差分法求解该方程，得到的稳态解组成的轮廓所包围的区域即为分割结果；

步骤一所述构造距离势场算子的方法为：

a、构造一个边长为2×R+1的方形零值矩阵L_R，在矩阵中构造一个半径为R的圆形模板，半径R的选择不超过待处理图像的边长的一半；

b、将圆形模板内的元素(i,j)的值设置为l(i,j)，l(i,j)通过公式获得，其中h，p为常数，0.5≤h≤2，1≤p≤3；得到的方形矩阵L_R即为距离势场算子；

步骤二所述用构造出的距离势场算子求图像的距离势场，对势场求梯度获得距离力场，并对力场归一化的方法为：

a、用步骤一构造出的距离势场算子L_R与图像I(x,y)做卷积，得到图像的距离势场E_DPF(x,y)；

b、对距离势场E_DPF(x,y)分别在x方向和y方向求偏导，得到距离力场 $f_{\mathrm{DPF}}(x,y)=(f_{\mathrm{DPF\; x}}(x,y),f_{\mathrm{DPF\; y}}(x,y))=(\frac{{\&PartialD;E}_{\mathrm{DPF}}(x,y)}{\&PartialD;x},\frac{{\&PartialD;E}_{\mathrm{DPF}}(x,y)}{\&PartialD;y}),$ 其中f_DPFx(x,y)和f_DPFy(x,y)分别为图像中(x,y)点沿x方向和y方向的距离力场；

c、采用公式对距离力场归一化；

步骤三所述设置初始轮廓，计算自适应气球力的方法为：

首先手动在图像中待分割目标周围选取点{c₁,c₂…c_n}作为初始轮廓点，其中c₁＝c_n，c_i＝(x_i,y_i),i＝1,2,...,n为所选点的坐标，c_i点的自适应气球力采用式

$\begin{array}{c}{n}_{\mathrm{balloon}}(x_i,y_i)=(n_{\mathrm{balloonx}}(x_i,y_i),n_{\mathrm{balloony}}(x_i,y_i))\\ =\mathrm{sign}\left(\mathrm{\&theta;}\right)\&CenterDot;(\frac{y_{i+1}-y_{i-1}}{\sqrt{{(x_{i+1}-x_{i-1})}^2+{(y_{i+1}+y_{i-1})}^2}}-\frac{x_{i+1}-x_{i-1}}{\sqrt{{(x_{i+1}-x_{i-1})}^2+{(y_{i+1}-y_{i-1})}^2}})\end{array}$

获得，其中θ为c_i处轮廓法线方向与该处距离力方向的夹角， $\mathrm{sign}\left(\mathrm{\&theta;}\right)=\left\{\begin{array}{cc}1& \mathrm{\&theta;}<\mathrm{\&tau;}\\ -1& \mathrm{\&theta;}\&GreaterEqual;\mathrm{\&tau;}\end{array}\right.,$ τ为设定的阈值常数，一般取π/4≤τ≤3π/4。

2.根据权利要求1所述的基于距离势场和自适应气球力的图像分割方法，其特征在于，步骤四所述的将步骤二、步骤三得到的距离力场和自适应气球力带入到轮廓的力学平衡方程，采用有限差分法求解该方程，得到稳态解的方法为：

将距离力和自适应气球力带入如下的轮廓力学平衡方程：

α(c_i-c_i-1)-α(c_i+1-c_i)

+β(c_i-2-2c_i-1+c_i)-2β(c_i-1-2c_i+c_i+1)+β(c_i-2c_i+1+c_i+2)

+λ(f_DPFx(c_i),f_DPFy(c_i))+κ(n_balloonx(c_i),n_balloony(c_i))＝0

其中λ和κ为常数，表示距离力和气球力的权重，一般取1≤λ≤3，1≤κ≤3，一般情况下κ要小于λ，上述方程写成矩阵形式为：

$\left\{\begin{array}{c}\mathrm{Ax}+{\mathrm{\&lambda;f}}_{\mathrm{DPF\; x}}(x,y)+{\mathrm{\&kappa;n}}_{\mathrm{balloonx}}(x,y)=0\\ \mathrm{Ay}+{\mathrm{\&lambda;f}}_{\mathrm{DPF\; y}}(x,y)+{\mathrm{\&kappa;n}}_{\mathrm{balloony}}(x,y)=0\end{array}\right.$

其中A为五对角带状矩阵：

$A=\left[\begin{array}{cccccccc}2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0& \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})\\ -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0& \mathrm{\&beta;}\\ \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& 0& ...& 0\\ 0& \mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}& -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& ...& 0\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ -(\mathrm{\&alpha;}+4\mathrm{\&beta;})& \mathrm{\&beta;}& .& .& .& .& .& 2\mathrm{\&alpha;}+6\mathrm{\&beta;}\end{array}\right]$

引入时间辅助变量t，将静态偏微分方程转化为动态方程；将右式的0用曲线关于时间间隔△t的一阶导数来替换，离散表达为

$0\&RightArrow;-\frac{x^t-x^{t-1}}{\mathrm{\&Delta;t}}$

$0\&RightArrow;-\frac{y^t-y^{t-1}}{\mathrm{\&Delta;t}}$

式中t——为迭代的次数；

△t——第t次迭代和第t-1次迭代的时间间隔；

(x^t,y^t)——经过第t次迭代的轮廓曲线位置坐标；

(x^t-¹,y^t-¹)——经过第t-1迭次的轮廓曲线位置坐标；

假设在一次的迭代过程中距离力场和自适应气球力保持恒定，带入方程组中得到

Ax^t+λf_DPFx(x^t-1,y^t-1)+κn_balloonx(x^t-1,y^t-1)＝-γ(x^t-x^t-1)

Ay^t+λf_DPFy(x^t-1,y^t-1)+κn_balloony(x^t-1,y^t-1)＝-γ(y^t-y^t-1)

其中为常数，表示迭代步长；距离力场是分布在整个图像域中的，任意点的距离力f_DPF(x^t,y^t)可以通过对附近点的力线性插值获得；得到x，y的迭代方程为：

x^t＝(A+γI)^-1(x^t-1-λf_dpfx(x^t-1,y^t-1)-κn_balloonx(x^t-1,y^t-1))

y^t＝(A+γI)^-1(y^t-1-λf_dpfy(x^t-1,y^t-1)-κn_balloony(x^t-1,y^t-1))

A+γI是带状矩阵，I为单位阵，它的逆可以通过LU分解获得；得到的稳态解组成的轮廓所包围的区域即为分割结果。