文章目錄

• 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算法的解就是將輪廓點移動到某個確定位置，使得以下等式成立：

$$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)+|?f{|}^2|\mathbf{v}??f{|}^{2}dxdy\text{\hspace{1em}(4)}$$

其中，$?f=[f_x,f_y]$，表示邊緣圖像（edge map）的梯度，任何圖像梯度算子都可以被用來計算這個變量。$u_x/u_y$、$v_x/v_y$為GVF得到的梯度場在$(x,y)$位置上$x$、$y$方向上的梯度值，由拉普拉斯方程計算得到，且其初始值定義為$?f$，即原圖像的二階導?？梢钥吹?#xff1a;

• 當$?f$很小時，能量由矢量場偏導數的平方和決定，從而產生一個緩慢變化的場；
• 當$?f$很大時，第二項占被積函數的主導地位，通過設置$\mathbf{v}=?f$可將其最小化，這可以再邊緣圖較大時保持幾乎等于邊緣圖的梯度，而在均勻區域中緩慢變化。

參數$\mu$是控制被積函數中第一項和第二項之間權衡的正則化參數，應根據圖像中存在的噪聲量來設置（噪聲越大$\mu$應越大）。可以預期由此最小化產生的矢量場既不是完全無旋場也不是完全無散場?？梢酝ㄟ^求解以下歐拉方程來求解GVF場（公式（4））：

$$\{\begin{array}{l}\mu {?}^{2}u?(u?f_x)(f_x^2+f_y^2)=0\\ \mu {?}^{2}v?(v?f_y)(f_x^2+f_y^2)=0\end{array}\text{\hspace{1em}(5)}$$

其中${?}^2$為拉普拉斯算子。

2.3 數值求解方法

將$u$和$v$看作時間的函數，可以將公示（5）（對應原文中公式（13a）、（13b））轉化為：

$$\{\begin{array}{l}{u}_t(x,y,t)=\mu {?}^{2}u(x,y,t)?[u(x,y,t)?f_x(x,y)]?[f_x(x,y{)}^2+f_y(x,y{)}^2]\\ v_t(x,y,t)=\mu {?}^{2}v(x,y,t)?[v(x,y,t)?f_y(x,y)]?[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 {?}^{2}u(x,y,t)?b(x,y)u(x,y,t)+c^1(x,y)\\ v_t(x,y,t)=\mu {?}^{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}{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)\\ {?}^{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})\\ {?}^{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 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实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。