文章目录

1 问题引入
- 1.1 传统Snake模型的缺陷
- 1.2 亥姆霍兹定理（Helmholtz theorem）
2 GVF Snake
- 2.1 边缘图（Edge Map）
- 2.2 梯度矢量流（Gradient Vector Flow，GVF）
- 2.3 数值求解方法
3 OpenCV实现

1 问题引入

1.1 传统Snake模型的缺陷

传统的Snake算法（参考博客）存在以下两个主要缺陷，在《Snakes, Shapes, and Gradient Vector Flow》这篇文章中提出的GVF Snake模型解决了这两个问题：

初始轮廓必须在真实轮廓附近，或对初始位置敏感；
很难进入凹陷的边界。

如图1所示，文章中作者用一个U形轮廓作为示例来直观解释传统Snake的缺陷。在图（c）中可以看到，U形轮廓凹陷处的力达到了左右平衡即方向相反大小相同，无法使轮廓再继续向凹陷处移动。

我们知道，Snake算法的解就是将轮廓点移动到某个确定位置，使得以下等式成立：

Fint+Fext(g)=0(1)F_{int} + F_{ext}^{(g)} = 0 \tag{1} Fint+Fext(g)=0(1)

其中外部力Fext(g)F_{ext}^{(g)}Fext(g)的选取对于Snake模型的表现有巨大影响，外部力又可以被分为两种：静态（static）和动态（dynamic）。其中静态力从图像数据计算并且不随模型的更新改变，动态力是随着模型更新而改变的力。传统的Snake模型外部力就是一种静态力，GVF Snake中所描述的外部力也是一种静态力，它不会随时间变化或依赖于Snake模型轮廓的位置。

1.2 亥姆霍兹定理（Helmholtz theorem）

GVF Snake中所描述外部力的数学原理来自于亥姆霍兹定理，即最一般的静态向量场可以被分解为两个分量：无旋（无旋度）分量和无散（无散度）分量。

由于传统Snake模型的变分公式长生的外力是势函数的梯度，因此它必须作为静态无旋场加入力平衡方程（公式(1)）。因此，可以通过允许包括无旋分量和无散分量来获得更一般的静态场。文章中，设计了一个外力场，使其具有大的捕获范围的同时又具有指向边界凹陷的力。

2 GVF Snake

GVF Snake中定义Fext(g)=v(x,y)F_{ext}^{(g)}=\mathbf{v}(x,y)Fext(g)=v(x,y)，仅这里与传统的Snake不同，称为梯度矢量流（Gradient Vector Flow，GVF）。

xt(s,t)=αx′′(s,t)−βx′′′′(s,t)+v(2)x_t(s,t)=\alpha x^{''}(s,t)-\beta x^{''''}(s,t)+\mathbf{v} \tag{2} xt(s,t)=αx′′(s,t)−βx′′′′(s,t)+v(2)

求解上述方程的参数曲线即为GVF Snake，求解方式与传统的Snake相同。由于v(x,y)\mathbf{v}(x,y)v(x,y)通常不是一个无旋场，因此该方程通常不表示（1）中的能量最小化问题的欧拉方程。

2.1 边缘图（Edge Map）

方法很多，比如可以用：

f(x,y)=−Eext(i)(x,y)(3)f(x,y)=-E_{ext}^{(i)}(x,y) \tag{3} f(x,y)=−Eext(i)(x,y)(3)

2.2 梯度矢量流（Gradient Vector Flow，GVF）

定义GVF为使以下能量函数最小化的矢量场v(x,y)=[u(x,y),v(x,y)]\mathbf{v}(x,y)=[u(x,y),v(x,y)]v(x,y)=[u(x,y),v(x,y)]：

ε=∫∫μ(ux2+uy2+vx2+vy2)+∣∇f∣2∣v−∇f∣2dxdy(4)\varepsilon = \int \int \mu (u_x^2+u_y^2+v_x^2+v_y^2)+| \nabla f |^2| \mathbf{v} - \nabla f |^2 dxdy \tag{4} ε=∫∫μ(ux2+uy2+vx2+vy2)+∣∇f∣2∣v−∇f∣2dxdy(4)

其中，∇f=[fx,fy]\nabla f=[f_x,f_y]∇f=[fx,fy]，表示边缘图像（edge map）的梯度，任何图像梯度算子都可以被用来计算这个变量。ux/uyu_x/u_yux/uy、vx/vyv_x/v_yvx/vy为GVF得到的梯度场在(x,y)(x,y)(x,y)位置上xxx、yyy方向上的梯度值，由拉普拉斯方程计算得到，且其初始值定义为∇f\nabla f∇f，即原图像的二阶导。可以看到：

当∇f\nabla f∇f很小时，能量由矢量场偏导数的平方和决定，从而产生一个缓慢变化的场；
当∇f\nabla f∇f很大时，第二项占被积函数的主导地位，通过设置v=∇f\mathbf{v}=\nabla fv=∇f可将其最小化，这可以再边缘图较大时保持几乎等于边缘图的梯度，而在均匀区域中缓慢变化。

参数μ\muμ是控制被积函数中第一项和第二项之间权衡的正则化参数，应根据图像中存在的噪声量来设置（噪声越大μ\muμ应越大）。可以预期由此最小化产生的矢量场既不是完全无旋场也不是完全无散场。可以通过求解以下欧拉方程来求解GVF场（公式（4））：

{μ∇2u−(u−fx)(fx2+fy2)=0μ∇2v−(v−fy)(fx2+fy2)=0(5)\begin{cases} \mu \nabla ^2 u - (u-f_x)(f_x^2+f_y^2) = 0 \\ \mu \nabla ^2 v - (v-f_y)(f_x^2+f_y^2) = 0 \end{cases} \tag{5} {μ∇2u−(u−fx)(fx2+fy2)=0μ∇2v−(v−fy)(fx2+fy2)=0(5)

其中∇2\nabla ^2∇2为拉普拉斯算子。

2.3 数值求解方法

将uuu和vvv看作时间的函数，可以将公示（5）（对应原文中公式（13a）、（13b））转化为：

{ut(x,y,t)=μ∇2u(x,y,t)−[u(x,y,t)−fx(x,y)]⋅[fx(x,y)2+fy(x,y)2]vt(x,y,t)=μ∇2v(x,y,t)−[v(x,y,t)−fy(x,y)]⋅[fx(x,y)2+fy(x,y)2](6)\begin{cases} u_t(x,y,t)=\mu \nabla^2 u(x,y,t)-[u(x,y,t)-f_x(x,y)]\cdot[f_x(x,y)^2+f_y(x,y)^2] \\ v_t(x,y,t)=\mu \nabla^2 v(x,y,t)-[v(x,y,t)-f_y(x,y)]\cdot[f_x(x,y)^2+f_y(x,y)^2] \end{cases} \tag{6} {ut(x,y,t)=μ∇2u(x,y,t)−[u(x,y,t)−fx(x,y)]⋅[fx(x,y)2+fy(x,y)2]vt(x,y,t)=μ∇2v(x,y,t)−[v(x,y,t)−fy(x,y)]⋅[fx(x,y)2+fy(x,y)2](6)

上式进一步简化：

{ut(x,y,t)=μ∇2u(x,y,t)−b(x,y)u(x,y,t)+c1(x,y)vt(x,y,t)=μ∇2v(x,y,t)−b(x,y)v(x,y,t)+c2(x,y)(7)\begin{cases} u_t(x,y,t)=\mu \nabla^2 u(x,y,t)-b(x,y)u(x,y,t)+c^1(x,y) \\ v_t(x,y,t)=\mu \nabla^2 v(x,y,t)-b(x,y)v(x,y,t)+c^2(x,y) \end{cases} \tag{7} {ut(x,y,t)=μ∇2u(x,y,t)−b(x,y)u(x,y,t)+c1(x,y)vt(x,y,t)=μ∇2v(x,y,t)−b(x,y)v(x,y,t)+c2(x,y)(7)

其中，

{b(x,y)=fx(x,y)2+fy(x,y)2c1(x,y)=b(x,y)fx(x,y)c2(x,y)=b(x,y)fy(x,y)(8)\begin{cases} b(x,y)=f_x(x,y)^2+f_y(x,y)^2 \\ c^1(x,y)=b(x,y)f_x(x,y) \\ c^2(x,y)=b(x,y)f_y(x,y) \end{cases} \tag{8} ⎩⎨⎧b(x,y)=fx(x,y)2+fy(x,y)2c1(x,y)=b(x,y)fx(x,y)c2(x,y)=b(x,y)fy(x,y)(8)

公式（7）中的系数b(x,y)b(x,y)b(x,y)、c1(x,y)c^1(x,y)c1(x,y)和c2(x,y)c^2(x,y)c2(x,y)在整个迭代过程中是固定的，只需要计算一次。为方便计算，作以下假设：

{i,j,n:对应x、y、tΔx、Δy:像素间隔spacingΔt:每次迭代的时间间隔\begin{cases} i,j,n: 对应x、y、t \\ \Delta x、\Delta y: 像素间隔spacing \\ \Delta t: 每次迭代的时间间隔 \end{cases} ⎩⎨⎧i,j,n:对应x、y、tΔx、Δy:像素间隔spacingΔt:每次迭代的时间间隔

则每次迭代的偏微分方程可以写为：

{ut=1Δt(ui,jn+1−ui,jn)vt=1Δt(vi,jn+1−vi,jn)∇2u=1ΔxΔy(ui+1,j+ui,j+1+ui−1,j+ui,j−1−4ui,j)∇2v=1ΔxΔy(vi+1,j+vi,j+1+vi−1,j+vi,j−1−4vi,j)(9)\begin{cases} u_t=\frac{1}{\Delta t} (u_{i,j}^{n+1} - u_{i,j}^n) \\ v_t=\frac{1}{\Delta t} (v_{i,j}^{n+1} - v_{i,j}^n) \\ \nabla^2u = \frac{1}{\Delta x \Delta y}(u_{i+1,j}+u_{i,j+1}+u_{i-1,j}+u_{i,j-1}-4u_{i,j}) \\ \nabla^2v = \frac{1}{\Delta x \Delta y}(v_{i+1,j}+v_{i,j+1}+v_{i-1,j}+v_{i,j-1}-4v_{i,j}) \end{cases} \tag{9} ⎩⎨⎧ut=Δt1(ui,jn+1−ui,jn)vt=Δt1(vi,jn+1−vi,jn)∇2u=ΔxΔy1(ui+1,j+ui,j+1+ui−1,j+ui,j−1−4ui,j)∇2v=ΔxΔy1(vi+1,j+vi,j+1+vi−1,j+vi,j−1−4vi,j)(9)

将这些近似值带入公式（7）可以得出GVF的迭代解：

{ui,jn+1=(1−bi,jΔt)ui,jn+r(ui+1n+ui,j+1n+ui−1,jn+ui,j−1n−4ui,jn)+ci,j1Δtvi,jn+1=(1−bi,jΔt)vi,jn+r(vi+1n+vi,j+1n+vi−1,jn+vi,j−1n−4vi,jn)+ci,j2Δt(10)\begin{cases} u_{i,j}^{n+1}=(1-b_{i,j}\Delta t)u_{i,j}^n+r(u_{i+1}^n+u_{i,j+1}^n+u_{i-1,j}^n+u_{i,j-1}^n-4u_{i,j}^n)+c_{i,j}^1 \Delta t \\ v_{i,j}^{n+1}=(1-b_{i,j}\Delta t)v_{i,j}^n+r(v_{i+1}^n+v_{i,j+1}^n+v_{i-1,j}^n+v_{i,j-1}^n-4v_{i,j}^n)+c_{i,j}^2 \Delta t \end{cases} \tag{10} {ui,jn+1=(1−bi,jΔt)ui,jn+r(ui+1n+ui,j+1n+ui−1,jn+ui,j−1n−4ui,jn)+ci,j1Δtvi,jn+1=(1−bi,jΔt)vi,jn+r(vi+1n+vi,j+1n+vi−1,jn+vi,j−1n−4vi,jn)+ci,j2Δt(10)

其中：

r=μΔtΔxΔy(11)r=\frac{\mu \Delta t}{\Delta x \Delta y} \tag{11} r=ΔxΔyμΔt(11)

为了保证收敛，必须确保Δt≤ΔxΔy4μ\Delta t \leq \frac{\Delta x \Delta y}{4\mu}Δt≤4μΔxΔy。

3 OpenCV实现

这里给出作为外力的梯度矢量流（GVF）的OpenCV实现。

std::vector<cv::Mat> GVFSnake::GVF(cv::Mat edgeImg, double mu, int iterNum)
{cv::Mat fImg = edgeImg.clone();fImg.convertTo(fImg, CV_32FC1);// a.Calculate the gradient of the edge mapqDebug() << "Calculate the gradient of the edge map";cv::Mat fx, fy;cv::Sobel(fImg, fx, fImg.depth(), 1, 0, 1, 0.5, 0, cv::BORDER_CONSTANT); // fImg must be CV_32FC1 formatcv::Sobel(fImg, fy, fImg.depth(), 0, 1, 1, 0.5, 0, cv::BORDER_CONSTANT);// b.Squared magnitude of the gradient fieldqDebug() << "Squared magnitude of the gradient field";cv::Mat SqrMagf = fx.mul(fx) + fy.mul(fy);// c.Initialize GVF to the gradientqDebug() << "Initialize GVF to the gradient";cv::Mat u = fx.clone();cv::Mat v = fy.clone();// d.Iteratively solve for the GVF u,vqDebug() << "Iteratively solve for the GVF u,v";for (int i = 0; i < iterNum; i++){// equal to del2 in matlabcv::Mat mKernal = (cv::Mat_<float>(3, 3) << 0, 1 / 4, 0,1 / 4, -1, 1 / 4,0, 1 / 4, 0);cv::Mat uLap, vLap;cv::filter2D(u, uLap, u.depth(), mKernal);cv::filter2D(v, vLap, v.depth(), mKernal);// why *4 ? ref: https://blog.csdn.net/qq_41634276/article/details/80659218u += mu * 4 * uLap - SqrMagf.mul(u - fx);v += mu * 4 * vLap - SqrMagf.mul(v - fy);}std::vector<cv::Mat> uvMatVector;uvMatVector.push_back(u);uvMatVector.push_back(v);return uvMatVector;
}

【主动轮廓模型（二）】《GVF Snake》算法原理与OpenCV实现相关推荐

离散化-利用计算机求解y=x,基于边缘的主动轮廓模型——从零到一用python实现snake...
从零到一实现snake算法 1.Snake算法原理 Kass等人最早于1988年提出了主动轮廓模型,该方法通过构造能量函数,从而将图像分割转化为求解能量泛函极值的问题,其核心思想在于,作者认为和之前的 ...
Active Contour Models 主动轮廓模型（snake模型）
主动轮廓模型主要用于解决图像中目标物体的分割操作.理论上是可以解决二维乃至多维的情况,不过最初的模型是在二维图像上建立的. 主动轮廓模型(Active Contour Model),又被称为Snake ...
水平集方法引入主动轮廓模型算法中的两种方法
水平集方法引入主动轮廓模型算法中的两种方法 1.传统的基于主动轮廓模型和水平集理论的方法 2.变分水平集方法在讲解水平集理论在主动轮廓模型中的应用前,我们先用流程图说明一下常见的处理主动轮廓模型的流 ...
Active Contour Models 主动轮廓模型
参考博客: https://www.mathworks.com/matlabcentral/fileexchange/19567-active-contour-segmentation 数字图像处理- ...
《Matlab图像处理》part1 Snakes：Active Contour Models 主动轮廓模型
<Matlab图像处理>part1 Snakes:Active Contour Models 主动轮廓模型参考博客: 数字图像处理-图像分割:Snake主动轮廓模型 Matlab代码及运 ...
基于带有信息熵和联合矢量的LBF主动轮廓模型的PET-CT成像中对静脉血管肺结节分割（笔记四）
-----------------------------------------------------------------SUV 标准吸收值-------------------------- ...
主动轮廓模型 matlab,主动轮廓模型的功能.ppt
主动轮廓模型的功能 Amelioration d'images ultrasonores via alignement de vasculature - Julien Jomier - UNC - 2 ...
OpenCV学习笔记（Python）———— 主动轮廓模型
本文包含主动轮廓模型代码以及实例分割代码原图: 效果图: 主动轮廓模型: morphsnakes.py # -*- coding: utf-8 -*- """ ==== ...
十三种基于直方图的图像全局二值化算法原理、实现、代码及效果（转）
源:十三种基于直方图的图像全局二值化算法原理.实现.代码及效果.
十三种基于直方图的图像全局二值化算法原理、实现、代码及效果
十三种基于直方图的图像全局二值化算法实现 1. 什么是基于直方图的图像全局二值化算法 2. 灰度平均值 3. 百分比阈值(P-Tile法) 3. 基于双峰的阈值 3.1 基于平均值的阈值 3.2 基于 ...

【主动轮廓模型（二）】《GVF Snake》算法原理与OpenCV实现