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

---

1 问题引入

1.1 传统Snake模型的缺陷

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

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

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

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

$$F_{int}+F_{ext}^{(g)}=0\text{\hspace{1em}(1)}$$

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

1.2 亥姆霍兹定理（Helmholtz theorem）

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

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

2 GVF Snake

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

$$x_t(s,t)=\alpha x^{{}^{\prime \prime }}(s,t)-\beta x^{{}^{\prime \prime \prime \prime }}(s,t)+\mathbf{v}\text{\hspace{1em}(2)}$$

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

2.1 边缘图（Edge Map）

方法很多，比如可以用：

$$f(x,y)=-E_{ext}^{(i)}(x,y)\text{\hspace{1em}(3)}$$

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

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

$$\epsilon =\int \int \mu (u_x^2+u_y^2+v_x^2+v_y^2)+|\nabla f{|}^2|\mathbf{v}-\nabla f{|}^{2}dxdy\text{\hspace{1em}(4)}$$

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

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

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

$$\{\begin{array}{l}\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{array}\text{\hspace{1em}(5)}$$

其中${\nabla }^2$为拉普拉斯算子。

2.3 数值求解方法

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

$$\{\begin{array}{l}{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{array}\text{\hspace{1em}(6)}$$

上式进一步简化：

$$\{\begin{array}{l}{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{array}\text{\hspace{1em}(7)}$$

其中，

$$\{\begin{array}{l}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{array}\text{\hspace{1em}(8)}$$

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

$$\{\begin{array}{l}i,j,n:对应x、y、t\\ \Delta x、\Delta y:像素间隔spacing\\ \Delta t:每次迭代的时间间隔\end{array}$$

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

$$\{\begin{array}{l}{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 }^{2}u=\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 }^{2}v=\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{array}\text{\hspace{1em}(9)}$$

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

$$\{\begin{array}{l}{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{array}\text{\hspace{1em}(10)}$$

其中：

$$r=\frac{\mu \Delta t}{\Delta x\Delta y}\text{\hspace{1em}(11)}$$

为了保证收敛，必须确保$\Delta t\le \frac{\Delta x\Delta y}{4\mu }$。

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 map
	qDebug() << "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 format
	cv::Sobel(fImg, fy, fImg.depth(), 0, 1, 1, 0.5, 0, cv::BORDER_CONSTANT);

	// b.Squared magnitude of the gradient field
	qDebug() << "Squared magnitude of the gradient field";
	cv::Mat SqrMagf = fx.mul(fx) + fy.mul(fy);

	// c.Initialize GVF to the gradient
	qDebug() << "Initialize GVF to the gradient";
	cv::Mat u = fx.clone();
	cv::Mat v = fy.clone();

	// d.Iteratively solve for the GVF u,v
	qDebug() << "Iteratively solve for the GVF u,v";
	for (int i = 0; i < iterNum; i++)
	{
		// equal to del2 in matlab
		cv::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/80659218
		u += 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;
}