CN102737378A - 测地线活动轮廓分割图像的方法和系统 - Google Patents

Abstract

一种测地线活动轮廓分割图像的方法，包括：获取图像，设置图像的初始轮廓；根据所述图像计算所述初始轮廓的距离场，根据所述距离场初始化二维格子玻尔兹曼模型演化方程的演化参数，所述演化参数包括平衡分布函数；根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程；根据所述推演结果计算所述目标轮廓的距离场，根据所述距离场判断所述目标轮廓的误差是否小于阈值，若是，则输出所述目标轮廓；否则，根据所述推演结果更新所述演化参数，继续执行所述根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程的步骤。此外，还提供了一种测地线活动轮廓分割图像的系统。上述测地线活动轮廓分割图像的方法和系统能提高执行效率。

Description

测地线活动轮廓分割图像的方法和系统

技术领域

本发明涉及图像分割领域，特别是涉及一种测地线活动轮廓分割图像的方法和系统。

背景技术

使用测地线活动轮廓模型能够得到图像中分割目标的封闭轮廓，非常利于图像的后期分析。测地线活动轮廓模型通常可用水平集方法表示。水平集方法的优点在于可适应曲线演化过程中物体表面的拓扑变化，并且作为一种全局方法，它可以有效对付由图像噪音或重复特征引起的错误匹配，具有较强的鲁棒性。

然而，利用水平集方法表示后，接着需要求解得到相应的偏微分方程。由于数字图像固有的不连续性，偏微分方程的非线性，解析解很难得到或者是根本不存在的。这时有必要借助于数值计算以获取该问题的近似解。由于计算区域是零水平集的高纬空间，水平集方法的数值计算量非常大。水平集方法中常用的数值计算方法包括：1）显式差分方法，其优点在于实现简单且直观。缺点在于受算法稳定性的限制，只能使用很小的迭代步长。整个计算过程需要多次迭代导致计算效率非常低；2）隐式差分方法，其优点在于稳定性好，能实现大步长的迭代计算以提高计算效率。缺点在于解隐式方程需求系数矩阵的逆，例如大小为M×N的图像，其系数矩阵将是(MN)²维，此计算量仍然较大。而且这种方法以求解线性方程组为基础，每个未知数的解都和全局有关，很难进行分块并行计算。

因此，传统技术中，分割图像的方法执行效率较低，活动轮廓分割图像的速度较慢。

发明内容

基于此，有必要提供一种能提高执行效率的测地线活动轮廓分割图像的方法。

一种测地线活动轮廓分割图像的方法，包括：

获取图像，设置图像的初始轮廓；

根据所述图像计算所述初始轮廓的距离场，根据所述距离场初始化二维格子玻尔兹曼模型演化方程的演化参数，所述演化参数包括平衡分布函数；

根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程；

根据所述推演结果计算所述目标轮廓的距离场，根据所述距离场判断所述目标轮廓的误差是否小于阈值，若是，则输出所述目标轮廓；否则，根据所述推演结果更新所述演化参数，继续执行所述根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程的步骤。

在其中一个实施例中，所述演化参数还包括松弛因子和外力项；

所述二维格子玻尔兹曼模型演化方程为：

$I_i(x+c_i,n+1)-I_i(x,n)=\mathrm{\ω}[I_i^{\mathrm{eq}}(x,n)-I_i(x,n)]+F_i$

其中，i为预设的演化方向，x为节点像素坐标，所述c_i为速度矢量，ω为松弛因子，F_i为外力项,为平衡分布函数，n为推演过程中的迭代次数，I_i(x,n)为迭代次数为n时，节点x的演化方向i上的分布函数。

在其中一个实施例中，所述根据所述图像计算所述初始轮廓的距离场，根据所述距离场初始化二维格子玻尔兹曼模型演化方程的演化参数的步骤为：

根据公式：

$D(p(x,y),S)=\underset{s\∈S}{\min }\left\{\right|p-s\left|\right\}=\underset{s(i,j)\∈S}{\min }\left\{\sqrt{{(x-i)}^2+{(y-j)}^2}\right\},\∀p(x,y)\∈M$

计算所述初始轮廓的距离场；其中，M为获取到的图像，p(x,y)为图像M上的像素点，S为初始轮廓，s(i,j)为初始轮廓S上的某个像素点；

根据预设的演化方向和所述初始轮廓的距离场计算演化初值、平衡分布函数；

计算所述图像的边缘截止函数，根据所述边缘截止函数和所述预设的演化方向计算所述松弛因子和外力项。

在其中一个实施例中，所述根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程的步骤为：

根据所述二维格子玻尔兹曼模型演化方程的碰撞过程和迁移过程推演所述二维格子玻尔兹曼模型演化方程。

在其中一个实施例中，所述根据所述推演结果更新所述演化参数的步骤为：

根据公式：

$\left\{\begin{array}{cc}{I}_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=\mathrm{1,2,3,4})\\ I_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/36& (i=\mathrm{5,6,7,8})\\ I_i^{\mathrm{eq}}(x,n+1)=4\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=0)\end{array}\right.$

更新所述平衡分布函数；其中，n为推演过程中的迭代次数，

为更新后的迭代次数为n+1时演化方向为i的平衡分布函数，I_i(x,n+1)为推演所得的迭代次数为n+1时的分布函数。

在其中一个实施例中，所述根据所述二维格子玻尔兹曼模型演化方程输出所述目标轮廓的步骤为：

根据所述目标轮廓的距离场提取零水平集；

根据所述零水平集输出目标轮廓。

在其中一个实施例中，所述方法还包括：

每迭代推演预设次数时，根据所述目标轮廓的距离场重置所述分布函数和演化参数，继续执行所述根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程的步骤；

所述演化参数包括演化初值、平衡分布函数。

此外，还有必要提供一种能提高执行效率的测地线活动轮廓分割图像的系统。

一种测地线活动轮廓分割图像的系统，包括：

图像获取模块，用于获取图像，设置图像的初始轮廓；

初始化模块，用于根据所述图像计算所述初始轮廓的距离场，根据所述距离场初始化二维格子玻尔兹曼模型演化方程的演化参数，所述演化参数包括平衡分布函数；

演化方程推演模块，用于根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程；

目标轮廓输出模块，用于根据所述推演结果计算所述目标轮廓的距离场，在所述目标轮廓的误差小于阈值时，输出所述目标轮廓；在所述目标轮廓的误差大于等于阈值时，根据所述推演结果更新所述演化参数，调用所述演化方程推演模块。

在其中一个实施例中，所述演化参数还包括松弛因子和外力项；

所述二维格子玻尔兹曼模型演化方程为：

$I_i(x+c_i,n+1)-I_i(x,n)=\mathrm{\ω}[I_i^{\mathrm{eq}}(x,n)-I_i(x,n)]+F_i$

其中，i为预设的演化方向，x为节点像素坐标，所述c_i为速度矢量，ω为松弛因子，F_i为外力项,

为平衡分布函数，n为推演过程中的迭代次数，I_i(x,n)为迭代次数为n时，节点x的演化方向i上的分布函数。

在其中一个实施例中，所述初始化模块还用于根据公式：

$D(p(x,y),S)=\underset{s\∈S}{\min }\left\{\right|p-s\left|\right\}=\underset{s(i,j)\∈S}{\min }\left\{\sqrt{{(x-i)}^2+{(y-j)}^2}\right\},\∀p(x,y)\∈M$

计算所述初始轮廓的距离场；其中，M为获取到的图像，p(x,y)为图像M上的像素点，S为初始轮廓，s(i,j)为初始轮廓S上的某个像素点，根据预设的演化方向和所述初始轮廓的距离场计算演化初值、平衡分布函数，计算所述图像的边缘截止函数，根据所述边缘截止函数和所述预设的演化方向计算所述松弛因子和外力项。

在其中一个实施例中，所述演化方程推演模块用于根据所述二维格子玻尔兹曼模型演化方程的碰撞过程和迁移过程推演所述二维格子玻尔兹曼模型演化方程。

在其中一个实施例中，所述目标轮廓输出模块还用于根据公式：

$\left\{\begin{array}{cc}{I}_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=\mathrm{1,2,3,4})\\ I_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/36& (i=\mathrm{5,6,7,8})\\ I_i^{\mathrm{eq}}(x,n+1)=4\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=0)\end{array}\right.$

更新所述平衡分布函数；其中，n为推演过程中的迭代次数，

为更新后的迭代次数为n+1时演化方向为i的平衡分布函数，I_i(x,n+1)为推演所得的迭代次数为n+1时的分布函数。

在其中一个实施例中，所述目标轮廓输出模块还用于根据所述目标轮廓的距离场提取零水平集，根据所述零水平集输出目标轮廓。

在其中一个实施例中，还包括演化方程重置模块，用于每迭代推演预设次数时，根据所述目标轮廓的距离场重置所述分布函数和演化参数，并继续调用所述演化方程推演模块；

所述演化参数包括平衡分布函数。

上述测地线活动轮廓分割图像的方法和系统，通过迭代的方式推演二维格子玻尔兹曼模型演化方程来输出目标轮廓。迭代推演过程中，每个像素节点的分布函数的更新只需要相邻节点的数据（8邻域的像素点），使得可通过并行计算来完成推演过程，例如，可使用GPU平台进行并行计算，从而提高了活动轮廓分割图像的执行效率。

附图说明

图1为一个实施例中测地线活动轮廓分割图像的方法的流程图；

图2为一个实施例中图像的初始轮廓和初始轮廓的距离场；

图3为另一个实施例中根据图像计算初始轮廓的距离场，根据距离场初始化二维格子玻尔兹曼模型演化方程的演化参数的流程图；

图4为一个实施例中迭代次数阈值为50次时输出的图像的目标轮廓和目标轮廓的距离场；

图5为另一个实施例中迭代次数阈值为100次时输出的图像的目标轮廓；

图6为另一个实施例中迭代次数阈值为150次时输出的图像的目标轮廓；

图7为另一个实施例中测地线活动轮廓分割图像的系统的结构示意图。

具体实施方式

在一个实施例中，如图1所示，一种测地线活动轮廓分割图像的方法，包括：

步骤S102，获取图像，设置图像的初始轮廓。

可将获取到的图像的边缘设置为初始轮廓。如图2所示，方形的边缘即为设置的初始轮廓。优选的，初始轮廓为圆形或方形，可以减少分割图像时的计算量。

步骤S104，根据图像计算初始轮廓的距离场，根据距离场初始化二维格子玻尔兹曼模型演化方程的演化参数，演化参数包括演化初值、平衡分布函数。

在一个实施例中，演化参数还包括松弛因子和外力项。二维格子玻尔兹曼模型演化方程为：

$I_i(x+c_i,n+1)-I_i(x,n)=\mathrm{\ω}[I_i^{\mathrm{eq}}(x,n)-I_i(x,n)]+F_i$

其中，i为预设的演化方向，x为节点像素坐标，c_i为速度矢量，ω为松弛因子，F_i为外力项,

为平衡分布函数，n为推演过程中的迭代次数，I_i(x,n)为代次数为n时，节点x的演化方向i上的分布函数。

在本实施例中，演化方向有9个，包括节点x的8个邻域方向和其自身不变的方向。演化方向i的取值范围为0至8。获取到的图像被划分为多个离散化的网格，网格中的每个节点与图像相应的像素对应。c_i为演化方向变化的速度矢量，可通过公式：

$\left\{\begin{array}{c}{c}_0=0\\ c_{\mathrm{1,2,3,4}}=(\cos \frac{(i-1)\mathrm{\π}}{2},\sin \frac{(i-1)\mathrm{\π}}{2})\\ c_{\mathrm{5,6,7,8}}=\sqrt{2}(\cos (\frac{(i-5)\mathrm{\π}}{2}+\frac{\mathrm{\π}}{4}),\sin (\frac{(i-5)\mathrm{\π}}{2}+\frac{\mathrm{\π}}{4}))\end{array}\right.$

计算演化方向变化的速度矢量c_i。其中，c₀为节点x自身不变方向的变化的速度矢量，c_{1，2，3，4}为节点x上下左右4个邻域方向的速度矢量，c_5，6，7,8为节点x剩余的4个邻域方向上速度矢量。

本实施例中，如图3所示，根据图像计算初始轮廓的距离场，根据距离场初始化二维格子玻尔兹曼模型演化方程的演化参数的步骤可具体为：

步骤S202，根据公式：

$D(p(x,y),S)=\underset{s\∈S}{\min }\left\{\right|p-s\left|\right\}=\underset{s(i,j)\∈S}{\min }\left\{\sqrt{{(x-i)}^2+{(y-j)}^2}\right\},\∀p(x,y)\∈M$

计算初始轮廓的距离场；其中，M为获取到的图像，p(x,y)为图像M上的像素点，S为初始轮廓，s(i,j)为初始轮廓S上的某个像素点。

步骤S204，根据预设的演化方向和初始轮廓的距离场计算演化初值、平衡分布函数。

演化初值为二维格子玻尔兹曼模型演化方程在迭代次数为0时的分布函数，即I_i(x,0)。本实施例中，可根据公式：

$\left\{\begin{array}{cc}{I}_i(x,0)=I_i^{\mathrm{eq}}(x,0)=D/9& (i=\mathrm{1,2,3,4})\\ I_i(x,0)=I_i^{\mathrm{eq}}(x,0)=D/36& (i=\mathrm{5,6,7,8})\\ I_i(x,0)=I_0^{\mathrm{eq}}(x,0)=4D/9& (i=0)\end{array}\right.$

计算演化初值和平衡分布函数。其中，I_i(x,0)为演化方向i上演化初值，为演化方向为i的平衡分布函数，D为由步骤S202计算得到的初始轮廓的距离场。

步骤S206，计算图像的边缘截止函数，根据边缘截止函数和预设的演化方向计算松弛因子和外力项。

本实施例中，可根据公式：

$g=\frac{1}{1+{\left(\right|\▿G_{\mathrm{\σ}}*M|/K)}^2}$

计算图像的边缘截止函数。其中，G_σ是方差为σ的高斯核（σ为预设值），M为图像的灰度值，K为阈值。例如，σ可优选为1，K可优选为5。

本实施例中，可根据公式：

$\mathrm{\ω}=\frac{2}{1+3C\·g}$

计算松弛因子。其中，ω为演化参数中的松弛因子，g为已获取的图像的边缘截止函数，C为预设的系数。

本实施例中，可根据公式：

$\left\{\begin{array}{cc}{F}_i=g/9,& (i=\mathrm{1,2,3,4})\\ F_i=g/36,& (i=\mathrm{5,6,7,8})\end{array}\right.$

计算演化参数中的外力项。其中，i为演化方向，g为已获取的图像的边缘截止函数。（i=0时的演化方向为非移动项，不需要迭代更新。故只有8个方向演化方程，所以不需要演化方向为i=0这个外力项）。

步骤S106，根据演化参数推演二维格子玻尔兹曼模型演化方程。

本实施例中，可根据二维格子玻尔兹曼模型演化方程的碰撞过程和迁移过程推演该二维格子玻尔兹曼模型演化方程。

本实施例中，在二维格子玻尔兹曼模型演化方程的碰撞过程中，可通过公式：

$I_i(x,n+1)=I_i(x,n)+\mathrm{\ω}[I_i^{\mathrm{eq}}(x,n)-I_i(x,n)]+F_i$

推演迭代次数加1后的分布函数。其中，I_i(x,n)为迭代次数为n+1时的分布函数，n为迭代次数，I_i(x,n)为迭代次数为n时的分布函数，ω为松弛因子，F_i为外力项，

为迭代次数为n时的平衡分布函数。

本实施例中，在二维格子玻尔兹曼模型演化方程的迁移过程中，若节点x处于计算网格的边缘（图像边缘），则可根据公式：

I_i(x+c_i,n+1)＝I_i(x+c_-i,n+1)

计算经过碰撞过程得出的迭代次数为n+1时的结合演化方向的速度矢量的分布函数。

若节点x处于计算网格的非边缘区域，则可根据公式：

I_i(x+c_i,n+1)＝I_i(x,n+1)

计算经过碰撞过程得出的迭代次数为n+1时的结合演化方向的速度矢量的分布函数。

需要说明的是，c_i只是的代表了迁移的方向，因此，x+c_i也为一个向量，包含在x中。

步骤S108，根据推演结果计算目标轮廓的距离场。

步骤S110，根据距离场判断目标轮廓的误差是否小于阈值，若是，则执行步骤S112，输出目标轮廓；否则，执行步骤S114，根据推演结果更新演化参数，并继续执行步骤S106，根据演化参数推演二维格子玻尔兹曼模型演化方程。

本实施例中，可根据公式：

$D_{\mathrm{old}}(n+1)=\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)$

计算目标轮廓的距离场。其中，D_old(n+1)为目标轮廓的距离场，i为演化方向，n为迭代次数，I_i(x,n+1)为迭代次数为n+1时的分布函数。

若D_old(n+1)与D_old(n)的差小于阈值，则目标轮廓的误差小于阈值，则本实施例中，可根据公式：

S_new＝contour_zero(D_old(n+1))

输出目标轮廓。其中，D_old为距离场，contour_zero为零水平集提取函数，S_new为提取出的目标轮廓。

若D_old(n+1)与D_old(n)的差大于阈值，则表示对二维格子玻尔兹曼模型方程需要继续推演。在推演之前，需要更新演化参数，具体为更新平衡分布函数。本实施例中，可根据公式：

$\left\{\begin{array}{cc}{I}_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=\mathrm{1,2,3,4})\\ I_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/36& (i=\mathrm{5,6,7,8})\\ I_i^{\mathrm{eq}}(x,n+1)=4\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=0)\end{array}\right.$

计算迭代次数加1时的平衡分布函数。其中，n为迭代次数，

为迭代次数加1（即迭代次数为n+1）时演化方向为i的平衡分布函数，I_i(x,n+1)为求解所得的迭代次数加1时的分布函数。

在另一个实施例中，还可在每迭代推演预设的次数时，根据距离场获取目标轮廓，可在获取到的目标轮廓的误差小于阈值时，输出目标轮廓。例如，可在每迭代推演5次时，根据计算得到的距离场获取目标轮廓的坐标序列，然后将该坐标序列与前次获取的目标轮廓的坐标序列相减，若方差小于阈值，则目标轮廓收敛且误差小于阈值，即可输出目标轮廓。

更新平衡分布函数之后，继续执行步骤S106，此时迭代次数由n变为n+1。也就是说，可迭代执行步骤S106、步骤S108、步骤S110、步骤S114，直至误差小于阈值，然后执行步骤S112，输出目标轮廓。

在另一个实施例中，每迭代推演预设次数时，可根据目标轮廓的距离场重置分布函数和演化参数，继续执行步骤S106。演化参数包括平衡分布函数。

在本实施例中，每迭代推演预设次数时，可根据前述的初始化演化方程的方法根据目标轮廓的距离场重置分布函数和平衡分布函数。例如，可每迭代推演10次，重新计算目标轮廓的距离场，然后根据距离场重新计算I_i(x,10)和

而不根据推演结果获取，并根据重新计算的I_i(x,10)和

继续推演过程，计算I_i(x,11)和

重新计算的I_i(x,10)和

的方法与前述的根据初始轮廓的距离场计算演化初值和平衡分布函数的公式和方法相同，仅仅在于迭代次数n由0变为10，在此不再赘述。

也就是说，每迭代执行步骤S106、步骤S108、步骤S110、步骤S114达到预设的次数时，可根据当前的目标轮廓的距离场重新设置二维格子玻尔兹曼模型演化方程的分布函数和演化参数，然后再根据重置后的分布函数和演化参数继续迭代执行步骤S106、步骤S108、步骤S110、步骤S114，直至误差小于阈值，然后执行步骤S112，输出目标轮廓。

由于演化过程中得到的目标轮廓的距离场相对于随机划定的初始轮廓的距离场更加接近图像真实的轮廓，因此，每迭代预设的次数重置演化参数，可以使二维格子玻尔兹曼模型演化方程可以随着演化的进程进行调整，从而使得迭代得到的目标轮廓更加接近图像中真实的轮廓。

本实施例中，请同时参考图4、图5、图6，分别为迭代了50次、100次、150次后输出的目标轮廓。可以看出，目标轮廓随着迭代推演的次数的增加逐步趋近图像中的真实轮廓。

在另一个实施例中，如图7所示，一种测地线活动轮廓分割图像的系统，包括：

图像获取模块102，用于获取图像，设置图像的初始轮廓。

可将获取到的图像的边缘设置为初始轮廓。如图2所示，方形的边缘即为设置的初始轮廓。优选的，初始轮廓为圆形或方形，可以减少分割图像时的计算量。

初始化模块104，用于根据图像计算初始轮廓的距离场，根据距离场初始化二维格子玻尔兹曼模型演化方程的演化参数，演化参数包括演化初值、平衡分布函数。

在一个实施例中，演化参数还包括松弛因子和外力项。二维格子玻尔兹曼模型演化方程为：

$I_i(x+c_i,n+1)-I_i(x,n)=\mathrm{\ω}[I_i^{\mathrm{eq}}(x,n)-I_i(x,n)]+F_i$

其中，i为预设的演化方向，x为节点像素坐标，c_i为速度矢量，ω为松弛因子，F_i为外力项,

为平衡分布函数，n为推演过程中的迭代次数，I_i(x,n)为迭代次数为n时，节点x的演化方向i上的分布函数。

在本实施例中，演化方向有9个，包括节点x的8个邻域方向和其自身不变的方向。演化方向i的取值范围为0至8。获取到的图像被划分为多个离散化的网格，网格中的每个节点与图像相应的像素对应。c_i为演化方向变化的速度矢量，初始化模块104可用于通过公式：

$\left\{\begin{array}{c}{c}_0=0\\ c_{\mathrm{1,2,3,4}}=(\cos \frac{(i-1)\mathrm{\π}}{2},\sin \frac{(i-1)\mathrm{\π}}{2})\\ c_{\mathrm{5,6,7,8}}=\sqrt{2}(\cos (\frac{(i-5)\mathrm{\π}}{2}+\frac{\mathrm{\π}}{4}),\sin (\frac{(i-5)\mathrm{\π}}{2}+\frac{\mathrm{\π}}{4}))\end{array}\right.$

计算演化方向变化的速度矢量c_i。其中，c₀为节点x自身不变方向的变化的速度矢量，c_{1，2，3，4}为节点x上下左右4个邻域方向的速度矢量，c_5，6，7,8为节点x剩余的4个邻域方向上速度矢量。

本实施例中，根据图像计算初始轮廓的距离场，初始化模块104可用于根据公式：

$D(p(x,y),S)=\underset{s\∈S}{\min }\left\{\right|p-s\left|\right\}=\underset{s(i,j)\∈S}{\min }\left\{\sqrt{{(x-i)}^2+{(y-j)}^2}\right\},\∀p(x,y)\∈M$

计算初始轮廓的距离场；其中，M为获取到的图像，p(x,y)为图像M上的像素点，S为初始轮廓，s(i,j)为初始轮廓S上的某个像素点。

本实施例中，初始化模块104还可用于根据预设的演化方向和初始轮廓的距离场计算演化初值、平衡分布函数。

演化初值为二维格子玻尔兹曼模型演化方程在迭代次数为0时的分布函数，即I_i(x,0)。本实施例中，初始化模块104还可用于根据公式：

$\left\{\begin{array}{cc}{I}_i(x,0)=I_i^{\mathrm{eq}}(x,0)=D/9& (i=\mathrm{1,2,3,4})\\ I_i(x,0)=I_i^{\mathrm{eq}}(x,0)=D/36& (i=\mathrm{5,6,7,8})\\ I_i(x,0)=I_0^{\mathrm{eq}}(x,0)=4D/9& (i=0)\end{array}\right.$

计算演化初值和平衡分布函数。其中，I_i(x,0)为演化方向i上演化初值，

为演化方向为i的平衡分布函数，D为计算得到的初始轮廓的距离场。

本实施例中，初始化模块104还可用于计算图像的边缘截止函数，根据边缘截止函数和预设的演化方向计算松弛因子和外力项。

本实施例中，初始化模块104可用于根据公式：

$g=\frac{1}{1+{\left(\right|\▿G_{\mathrm{\σ}}*M|/K)}^2}$

计算图像的边缘截止函数。其中，G_σ是方差为σ的高斯核（σ为预设值），M为图像的灰度值，K为阈值。例如，σ可优选为1，K可优选为5。

本实施例中，初始化模块104可用于根据公式：

$\mathrm{\ω}=\frac{2}{1+3C\·g}$

计算松弛因子。其中，ω为演化参数中的松弛因子，g为已获取的图像的边缘截止函数，C为预设的系数。

本实施例中，初始化模块104可用于根据公式：

$\left\{\begin{array}{cc}{F}_i=g/9,& (i=\mathrm{1,2,3,4})\\ F_i=g/36,& (i=\mathrm{5,6,7,8})\end{array}\right.$

计算演化参数中的外力项。其中，i为演化方向，g为已获取的图像的边缘截止函数。（i=0时的演化方向为非移动项，不需要迭代更新。故只有8个方向演化方程，所以不需要演化方向为i=0这个外力项）。

演化方程推演模块106，用于根据演化参数推演二维格子玻尔兹曼模型演化方程。

本实施例中，演化方程推演模块106可用于根据二维格子玻尔兹曼模型演化方程的碰撞过程和迁移过程推演该二维格子玻尔兹曼模型演化方程。

本实施例中，在二维格子玻尔兹曼模型演化方程的碰撞过程中，演化方程推演模块106可用于通过公式：

$I_i(x,n+1)=I_i(x,n)+\mathrm{\ω}[I_i^{\mathrm{eq}}(x,n)-I_i(x,n)]+F_i$

推演迭代次数加1后的分布函数。其中，I_i(x,n)为迭代次数为n+1时的分布函数，n为迭代次数，I_i(x,n)为迭代次数为n时的分布函数，ω为松弛因子，F_i为外力项，

为迭代次数为n时的平衡分布函数。

本实施例中，在二维格子玻尔兹曼模型演化方程的迁移过程中，若节点x处于计算网格的边缘（图像边缘），则演化方程推演模块106可用于根据公式：

I_i(x+c_i,n+1)＝I_i(x+c_-i,n+1)

计算经过碰撞过程得出的迭代次数为n+1时的结合演化方向的速度矢量的分布函数。

若节点x处于计算网格的非边缘区域，则演化方程推演模块106可用于根据公式：

I_i(x+c_i,n+1)＝I_i(x,n+1)

计算经过碰撞过程得出的迭代次数为n+1时的结合演化方向的速度矢量的分布函数。

需要说明的是，c_i只是的代表了迁移的方向，因此，x+c_i也为一个向量，包含在x中。

目标轮廓输出模块108，用于根据推演结果计算目标轮廓的距离场，在目标轮廓的误差小于阈值时，输出目标轮廓；在目标轮廓的误差大于等于阈值时，根据推演结果更新演化参数，调用演化方程推演模块。

本实施例中，目标轮廓输出模块108可用于根据公式：

$D_{\mathrm{old}}(n+1)=\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)$

计算目标轮廓的距离场。其中，D_old(n+1)为目标轮廓的距离场，i为演化方向，n为迭代次数，I_i(x,n+1)为迭代次数为n+1时的分布函数。

若D_old(n+1)与D_old(n)的差小于阈值，则目标轮廓的误差小于阈值，则本实施例中，可根据公式：

S_new＝contour_zero(D_old(n+1))

输出目标轮廓。其中，D_old为距离场，contour_zero为零水平集提取函数，S_new为提取出的目标轮廓。

若D_old(n+1)与D_old(n)的差大于阈值，则表示对二维格子玻尔兹曼模型方程需要继续推演。在推演之前，需要更新演化参数，具体为更新平衡分布函数。本实施例中，可根据公式：

$\left\{\begin{array}{cc}{I}_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=\mathrm{1,2,3,4})\\ I_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/36& (i=\mathrm{5,6,7,8})\\ I_i^{\mathrm{eq}}(x,n+1)=4\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=0)\end{array}\right.$

计算迭代次数加1时的平衡分布函数。其中，n为迭代次数，

为迭代次数加1（即迭代次数为n+1）时演化方向为i的平衡分布函数，I_i(x,n+1)为求解所得的迭代次数加1时的分布函数。

在另一个实施例中，目标轮廓输出模块108还可用于在每迭代推演预设的次数时，根据距离场获取目标轮廓，可在获取到的目标轮廓的误差小于阈值时，输出目标轮廓。例如，可在每迭代推演5次时，根据计算得到的距离场获取目标轮廓的坐标序列，然后将该坐标序列与前次获取的目标轮廓的坐标序列相减，若方差小于阈值，则目标轮廓收敛且误差小于阈值，即可输出目标轮廓。

更新平衡分布函数之后，目标轮廓输出模块108调用演化方程推演模块106继续推演，此时迭代次数由n变为n+1。也就是说，演化方程推演模块106和目标轮廓输出模块108以迭代的方式推演二维格子玻尔兹曼演化方程，直至输出误差小于阈值的目标轮廓。

在另一个实施例中，如图7所示，测地线活动轮廓分割图像的系统还包括演化方程重置模块110，可用于在每迭代推演预设次数时，可根据目标轮廓的距离场重置分布函数和演化参数，并调用演化方程推演模块。演化参数包括平衡分布函数。

在本实施例中，每迭代推演预设次数时，演化方程重置模块110可用于根据前述的初始化演化方程的方法根据目标轮廓的距离场重置分布函数和平衡分布函数。例如，可每迭代推演10次，重新计算目标轮廓的距离场，然后根据距离场重新计算I_i(x,10)和

而不根据推演结果获取，并根据重新计算的I_i(x,10)和

继续推演过程，计算I_i(x,11)和

重新计算的I_i(x,10)和

的方法与前述的根据初始轮廓的距离场计算演化初值和平衡分布函数的公式和方法相同，仅仅在于迭代次数n由0变为10，在此不再赘述。

也就是说，演化方程推演模块106和目标轮廓输出模块108迭代推演预设的次数时，演化方程重置模块110可用于根据当前的目标轮廓的距离场重新设置二维格子玻尔兹曼模型演化方程的分布函数和演化参数，然后演化方程推演模块106再根据重置后的分布函数和演化参数继续迭代推演，直至输出误差小于阈值的目标轮廓。

由于演化过程中得到的目标轮廓的距离场相对于随机划定的初始轮廓的距离场更加接近图像真实的轮廓，因此，每迭代预设的次数重置演化参数，可以使二维格子玻尔兹曼模型演化方程可以随着演化的进程进行调整，从而使得迭代得到的目标轮廓更加接近图像中真实的轮廓。

本实施例中，请同时参考图4、图5、图6，分别为迭代了50次、100次、150次后输出的目标轮廓。可以看出，目标轮廓随着迭代推演的次数的增加逐步趋近图像中的真实轮廓。

上述测地线活动轮廓分割图像的方法和系统，通过迭代的方式推演二维格子玻尔兹曼模型演化方程来输出目标轮廓。迭代推演过程中，每个像素节点的分布函数的更新只需要相邻节点的数据（8邻域的像素点），使得可通过并行计算来完成推演过程，例如，可使用GPU平台进行并行计算，从而提高了活动轮廓分割图像的执行效率。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (14)

1.一种测地线活动轮廓分割图像的方法，包括：

获取图像，设置图像的初始轮廓；

根据所述图像计算所述初始轮廓的距离场，根据所述距离场初始化二维格子玻尔兹曼模型演化方程的演化参数，所述演化参数包括平衡分布函数；

根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程；

根据所述推演结果计算所述目标轮廓的距离场，根据所述距离场判断所述目标轮廓的误差是否小于阈值，若是，则输出所述目标轮廓；否则，根据所述推演结果更新所述演化参数，继续执行所述根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程的步骤。

2.根据权利要求1所述的测地线活动轮廓分割图像的方法，其特征在于，所述演化参数还包括松弛因子和外力项；

所述二维格子玻尔兹曼模型演化方程为：

$I_i(x+c_i,n+1)-I_i(x,n)=\mathrm{\ω}[I_i^{\mathrm{eq}}(x,n)-I_i(x,n)]+F_i$

其中，i为预设的演化方向，x为节点像素坐标，所述c_i为速度矢量，ω为松弛因子，F_i为外力项,

为平衡分布函数，n为推演过程中的迭代次数，I_i(x,n)为迭代次数为n时，节点x的演化方向i上的分布函数。

3.根据权利要求2所述的测地线活动轮廓分割图像的方法，其特征在于，所述根据所述图像计算所述初始轮廓的距离场，根据所述距离场初始化二维格子玻尔兹曼模型演化方程的演化参数的步骤为：

根据公式：

$D(p(x,y),S)=\underset{s\∈S}{\min }\left\{\right|p-s\left|\right\}=\underset{s(i,j)\∈S}{\min }\left\{\sqrt{{(x-i)}^2+{(y-j)}^2}\right\},\∀p(x,y)\∈M$

计算所述初始轮廓的距离场；其中，M为获取到的图像，p(x,y)为图像M上的像素点，S为初始轮廓，s(i,j)为初始轮廓S上的某个像素点；

根据预设的演化方向和所述初始轮廓的距离场计算演化初值、平衡分布函数；

计算所述图像的边缘截止函数，根据所述边缘截止函数和所述预设的演化方向计算所述松弛因子和外力项。

4.根据权利要求2所述的测地线活动轮廓分割图像的方法，其特征在于，所述根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程的步骤为：

根据所述二维格子玻尔兹曼模型演化方程的碰撞过程和迁移过程推演所述二维格子玻尔兹曼模型演化方程。

5.根据权利要求2所述的测地线活动轮廓分割图像的方法，其特征在于，所述根据所述推演结果更新所述演化参数的步骤为：

根据公式：

$\left\{\begin{array}{cc}{I}_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=\mathrm{1,2,3,4})\\ I_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/36& (i=\mathrm{5,6,7,8})\\ I_i^{\mathrm{eq}}(x,n+1)=4\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=0)\end{array}\right.$

更新所述平衡分布函数；其中，n为推演过程中的迭代次数，

为更新后的迭代次数为n+1时演化方向为i的平衡分布函数，I_i(x,n+1)为推演所得的迭代次数为n+1时的分布函数。

6.根据权利要求2所述的测地线活动轮廓分割图像的方法，其特征在于，所述根据所述二维格子玻尔兹曼模型演化方程输出所述目标轮廓的步骤为：

根据所述目标轮廓的距离场提取零水平集；

根据所述零水平集输出目标轮廓。

7.根据权利要求2所述的测地线活动轮廓分割图像的方法，其特征在于，所述方法还包括：

每迭代推演预设次数时，根据所述目标轮廓的距离场重置所述分布函数和演化参数，继续执行所述根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程的步骤；

所述演化参数包括演化初值、平衡分布函数。

8.一种测地线活动轮廓分割图像的系统，其特征在于，包括：

图像获取模块，用于获取图像，设置图像的初始轮廓；

初始化模块，用于根据所述图像计算所述初始轮廓的距离场，根据所述距离场初始化二维格子玻尔兹曼模型演化方程的演化参数，所述演化参数包括平衡分布函数；

演化方程推演模块，用于根据所述演化参数推演所述二维格子玻尔兹曼模型演化方程；

目标轮廓输出模块，用于根据所述推演结果计算所述目标轮廓的距离场，在所述目标轮廓的误差小于阈值时，输出所述目标轮廓；在所述目标轮廓的误差大于等于阈值时，根据所述推演结果更新所述演化参数，调用所述演化方程推演模块。

9.根据权利要求8所述的测地线活动轮廓分割图像的系统，其特征在于，所述演化参数还包括松弛因子和外力项；

所述二维格子玻尔兹曼模型演化方程为：

$I_i(x+c_i,n+1)-I_i(x,n)=\mathrm{\ω}[I_i^{\mathrm{eq}}(x,n)-I_i(x,n)]+F_i$

其中，i为预设的演化方向，x为节点像素坐标，所述c_i为速度矢量，ω为松弛因子，F_i为外力项,为平衡分布函数，n为推演过程中的迭代次数，I_i(x,n)为迭代次数为n时，节点x的演化方向i上的分布函数。

10.根据权利要求9所述的测地线活动轮廓分割图像的系统，其特征在于，所述初始化模块还用于根据公式：

$D(p(x,y),S)=\underset{s\∈S}{\min }\left\{\right|p-s\left|\right\}=\underset{s(i,j)\∈S}{\min }\left\{\sqrt{{(x-i)}^2+{(y-j)}^2}\right\},\∀p(x,y)\∈M$

计算所述初始轮廓的距离场；其中，M为获取到的图像，p(x,y)为图像M上的像素点，S为初始轮廓，s(i,j)为初始轮廓S上的某个像素点，根据预设的演化方向和所述初始轮廓的距离场计算演化初值、平衡分布函数，计算所述图像的边缘截止函数，根据所述边缘截止函数和所述预设的演化方向计算所述松弛因子和外力项。

11.根据权利要求9所述的测地线活动轮廓分割图像的系统，其特征在于，所述演化方程推演模块用于根据所述二维格子玻尔兹曼模型演化方程的碰撞过程和迁移过程推演所述二维格子玻尔兹曼模型演化方程。

12.根据权利要求9所述的测地线活动轮廓分割图像的系统，其特征在于，所述目标轮廓输出模块还用于根据公式：

$\left\{\begin{array}{cc}{I}_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=\mathrm{1,2,3,4})\\ I_i^{\mathrm{eq}}(x,n+1)=\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/36& (i=\mathrm{5,6,7,8})\\ I_i^{\mathrm{eq}}(x,n+1)=4\left(\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{I}_i(x,n+1)\right)/9& (i=0)\end{array}\right.$

更新所述平衡分布函数；其中，n为推演过程中的迭代次数，

为更新后的迭代次数为n+1时演化方向为i的平衡分布函数，I_i(x,n+1)为推演所得的迭代次数为n+1时的分布函数。

13.根据权利要求9所述的测地线活动轮廓分割图像的系统，其特征在于，所述目标轮廓输出模块还用于根据所述目标轮廓的距离场提取零水平集，根据所述零水平集输出目标轮廓。

14.根据权利要求9所述的测地线活动轮廓分割图像的系统，其特征在于，还包括演化方程重置模块，用于每迭代推演预设次数时，根据所述目标轮廓的距离场重置所述分布函数和演化参数，并继续调用所述演化方程推演模块；所述演化参数包括平衡分布函数。

Images