文章目錄

• 1 概述
• 2 算法原理
• • 2.1 內部能量 $E_{int}$
  • 2.2 圖像能量 $E_{image}$
  • • （1）線函數$E_{line}$
    • （2）邊函數$E_{edge}$
    • （3）末端函數$E_{term}$
  • 2.3 外部能量 $E_{con}$
• 3 模型求解
• 4 算法實現（OpenCV3）

---

1 概述

主動輪廓模型（也稱Active Contour Model、Snake）是Kass等人在1988年提出的，該算法將圖像分割問題轉換為求解能量泛函最小值的問題。主要思路是通過構造能量泛函，經過算法迭代，輪廓曲線由初始位置逐漸向使能量函數最小（或局部極小）的圖像邊緣逼近，最終分割出目標。

2 算法原理

首先需要人為地在圖像上給出初始輪廓曲線，確切的說是一組用于控制曲線形狀的控制點：$v(s)=[x(s),y(s)]s\in [0,1]$，這些點收尾相連構成一個封閉的輪廓線。其中$x(s)$和$y(s)$分別表示每個控制點在圖像中的坐標位置，$s$是以傅立葉變換形式描述邊界的自變量，也可以理解為弧長。則Snake曲線的能量函數表示為：

$$\begin{array}{rl}{E}_{snake}^{?}& ={\int }_0^1{E}_{snake}(v(s))ds\\ & ={\int }_0^1{E}_{int}(v(s))+E_{image}(v(s))+E_{con}(v(s))ds\\ & ={\int }_0^1{E}_{int}(v(s))+E_{ext}(v(s))ds\end{array}\text{\hspace{1em}(1)}$$

其中，$E_{int}$為內部能量，$E_{image}$為圖像能量，$E_{con}$為外部約束能量。其中圖像能量和外部約束能量統稱為外部能量，即：$E_{ext}=E_{image}+E_{con}$。

2.1 內部能量 $E_{int}$

內部能量$E_{int}$由保證曲線連續性的一階導和保證曲線平滑的二階導組成，表示為：

$$E_{int}=\frac{1}{2}(\alpha (s)|v_s(s){|}^2+\beta (s)|v_{ss}(s){|}^2)\text{\hspace{1em}(2)}$$

通過調整權值$\alpha (s)$和$\beta (s)$可以控制曲線的形狀。例如將$\beta (s)$置為0可以讓曲線最終出現拐角，即曲線二階不連續。

2.2 圖像能量 $E_{image}$

在低層次的計算機視覺中，需要能夠將輪廓吸引到特定圖像特征的能量函數。原始的Snake輪廓模型中提出了3種不同的能量函數，分別將snake輪廓吸引到線、邊和末端。完整的圖像能量 $E_{image}$可以表示為這3個能量函數的權值組合，通過調整這3個權值，可以形成不同的輪廓形狀。

$$E_{image}={\omega }_{line}{E}_{line}+{\omega }_{edge}{E}_{edge}+{\omega }_{term}{E}_{term}\text{\hspace{1em}(3)}$$

（1）線函數$E_{line}$

最簡單直接且有用的圖像能量函數是圖像本身，即圖像本身的灰度。令：

$$E_{line}=I(x,y)\text{\hspace{1em}(4)}$$

控制${\omega }_{line}$的正負號可以控制輪廓被吸引到較暗的線或是較亮的線，也就是使輪廓試圖靠近輪廓的最暗或最亮處。

然而，如果snake輪廓的一部分到達了一個低能量的圖像特征位置，這個樣條項將推動snake輪廓臨近的部分朝著這個特征可能的延續方向移動，這會使其在一個最優的局部最小位置引入一個較大的能量。一種解決方案是允許snake輪廓和模糊能量函數平衡，然后慢慢降低模糊程度。

Marr和Hildreth在他們的論文里證明了“灰度的突然變化會在一階導數中引起波峰或波谷，或在二階導數中等效地引起零交點”。為了顯示圖像尺度空間連續性和Marr-Hildreth邊緣檢測理論的關系，snake模型中采用了模糊的邊能量函數：

$$E_{line}=?(G_{\sigma }?{?}^{2}I{)}^2\text{\hspace{1em}(6)}$$

其中$G_{\sigma }$是標準差為$\sigma$的高斯函數，該函數的最小值位于在Marr-Hildreth理論中被定義的$G_{\sigma }?{?}^{2}I$的零交點處。在能量函數中加入這一項意味著snake輪廓在被吸引到零交點的同時，仍然受到它自己的平滑限制。

（2）邊函數$E_{edge}$

在圖像上找邊緣可以通過梯度來實現，令：

$$E_{edge}=?|?I(x,y){|}^2\text{\hspace{1em}(5)}$$

在某個點上，梯度越大上式的能量越小，則snake輪廓將被吸引到梯度較大的區域。

（3）末端函數$E_{term}$

為了找到輪廓的終止位置，作者將平滑過圖像中等高線的曲率加入到能量函數中。令$C(x,y)=G_{\sigma }(x,y)?I(x,y)$是模糊后的圖像，$\theta =ta{n}^{?1}(\frac{C_y}{C_x})$是梯度角，$n=(cos\theta ,sin\theta )$、$n_{\perp }=(?sin\theta ,cos\theta )$是沿著和垂直梯度方向的單位向量。則$C(x,y)$中等高線的曲率可以表示為：

$$\begin{array}{rl}{E}_{term}& =\frac{?\theta }{?n_{\perp }}\\ & =\frac{?C^2/?n_{\perp }^2}{?C/?n}\\ & =\frac{C_{yy}{C}_x^2?2C_{xy}{C}_x{C}_y+C_{xx}{C}_y^2}{(C_x^2+C_y^2{)}^{3/2}}\end{array}\text{\hspace{1em}(6)}$$

通過組合$E_{edge}$和$E_{term}$，我們能夠創建一個被邊和末端吸引的snake曲線。

2.3 外部能量 $E_{con}$

外部能量$E_{con}$來自于外部的約束力。在論文原文中，作者給了一個外部約束的例子， 包括外部的固定點、連接兩條snake曲線的錨點或鼠標拖動的點。例如，為了在$x_1$和$x_2$點之間創建一個連接的彈性力，就可以將$?k(x_1?x_2{)}^2$添加到外部能量$E_{con}$中。

3 模型求解

由歐拉方程，求解能量$E_{snake}^{?}$的最小值（局部極小值），公式（1）的導數必須滿足：

$$\alpha v_{ss}(s)?\beta v_{ssss}(s)??E_{image}(v(s))??E_{con}(v(s))=0\text{\hspace{1em}(7)}$$

由$v(s)=[x(s),y(s)]$，將上式改寫成x和y兩個方向有：

$$\{\begin{array}{l}\alpha x_{ss}(s)?\beta x_{ssss}(s)?\frac{?E_{ext}}{?x}=0\\ \alpha y_{ss}(s)?\beta y_{ssss}(s)?\frac{?E_{ext}}{?y}=0\end{array}\text{\hspace{1em}(8)}$$

利用順序連接的錨點的坐標差分來近似導數：

$$?\begin{array}{l}{x}_{ss}(s)=x(s+1)+x(s?1)?2x(s)\\ \begin{array}{rl}{x}_{ssss}(s)=& (x(s+2)+x(s)?2x(s+1))\\ & +(x(s)+x(s?2)?2x(s?1))\\ & ?2(x(s+1)+x(s?1)?2x(s))\end{array}\end{array}\text{\hspace{1em}(9)}$$

將（9）帶入（8），同時令$f_x=\frac{?E_{ext}}{?x},f_y=\frac{?E_{ext}}{?y}$有：

$$\{\begin{array}{l}\beta x(s?2)?(\alpha +4\beta )x(s?1)+(2\alpha +6\beta )x(s)?(\alpha +4\beta )x(s+1)+\beta x(s+2)+f_x=0\\ \beta y(s?2)?(\alpha +4\beta )y(s?1)+(2\alpha +6\beta )y(s)?(\alpha +4\beta )y(s+1)+\beta y(s+2)+f_y=0\end{array}\text{\hspace{1em}(10)}$$

令：

$$?\begin{array}{l}a=2\alpha +6\beta \\ b=?(\alpha +4\beta )\\ c=\beta \end{array}\text{\hspace{1em}(11)}$$

則將（10）寫成矩陣形式為：

$$\{\begin{array}{l}Ax+f_x(x,y)=0\\ Ay+f_y(x,y)=0\end{array}\text{\hspace{1em}(12)}$$

其中A為五對角帶狀矩陣：

$$A=?\begin{array}{cccccc}a& b& c& ?& c& b\\ b& a& b& c& ?& c\\ c& b& a& b& c& ?\\ ?& ?& ?& ?& ?& ?\\ ?& c& b& a& b& c\\ c& ?& c& b& a& b\\ b& c& ?& c& b& a\end{array}?,x=?\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ?\\ x_{n?1}\\ x_n\\ x_1\end{array}?,f_x=?\begin{array}{c}{f}_{x_1}\\ f_{x_2}\\ f_{x_3}\\ ?\\ f_{x_{n?1}}\\ f_{x_n}\\ f_{x_1}\end{array}?\text{\hspace{1em}(13)}$$

利用梯度下降法，在第t次迭代有：

$$\{\begin{array}{l}A{x}_t+f_x(x_{t?1},y_{t?1})=?\gamma (x_t?x_{t?1})\\ A{y}_t+f_y(y_{t?1},y_{t?1})=?\gamma (y_t?y_{t?1})\end{array}\text{\hspace{1em}(14)}$$

其中$\gamma$是迭代步長，公式（14）的解為：

$$\{\begin{array}{l}{x}_t=(A+\gamma I{)}^{?1}(x_{t?1}?f_x(x_{t?1},y_{t?1}))\\ y_t=(A+\gamma I{)}^{?1}(y_{t?1}?f_y(x_{t?1},y_{t?1}))\end{array}\text{\hspace{1em}(15)}$$

由于$A$為五對角帶狀矩陣，因此$A+\gamma I$也是一個五對角帶狀矩陣，原文中作者用LU分解來求其逆矩陣。

4 算法實現（OpenCV3）

網上關于snake算法的實現多是matlab代碼。OpenCV2中可用cvSnakeImage來實現，但這個函數在OpenCV3中已經被刪除。本文將在OpenCV3.14中實現snake算法代碼。

cv::Mat Interate(cv::Mat image,cv::Mat xs,cv::Mat ys,double alpha,double beta,double gamma,double kappa,double wl,double we,double wt,int iterations ) {// 相關參數int N = iterations;cv::Mat smth = image.clone();// 圖像大小qDebug() << "Calculating size of image";cv::Size size = image.size();int row = size.height;int col = size.width;// 計算外部力（圖像力）qDebug() << "Computing external forces";cv::Mat E_line = smth.clone(); // E_line is simply the image intensitiescv::Mat gradx, grady;cv::Sobel(smth, gradx, smth.depth(), 1, 0, 1, 1, 0, cv::BORDER_CONSTANT);cv::Sobel(smth, grady, smth.depth(), 0, 1, 1, 1, 0, cv::BORDER_CONSTANT);qDebug() << "Computing gradx and grady";cv::Mat E_edge(row, col, CV_32FC1);for (int i = 0; i < gradx.rows; i++){for (int j = 0; j < gradx.cols; j++){float v_gradx = gradx.at<float>(i, j);float v_grady = grady.at<float>(i, j);E_edge.at<float>(i, j) = -1 * std::sqrt(v_gradx * v_gradx + v_grady * v_grady); // E_edge is measured by gradient in the image}}// 導數maskqDebug() << "masks for taking various derivatives";cv::Mat m1 = (cv::Mat_<float>(1, 2) << -1, 1);cv::Mat m2 = (cv::Mat_<float>(2, 1) << -1, 1);cv::Mat m3 = (cv::Mat_<float>(1, 3) << 1, -2, 1);cv::Mat m4 = (cv::Mat_<float>(3, 1) << 1, -2, 1);cv::Mat m5 = (cv::Mat_<float>(2, 2) << 1, -1, -1, 1);cv::Mat cx, cy, cxx, cyy, cxy;filter2D(smth, cx, -1, m1);filter2D(smth, cy, -1, m2);filter2D(smth, cxx, -1, m3);filter2D(smth, cyy, -1, m4);filter2D(smth, cxy, -1, m5);// 計算 E_termcv::Mat E_term(row, col, CV_32FC1);for (int i = 0; i < row; i++){for (int j = 0; j < col; j++){int v_cx = cx.at<float>(i, j);int v_cy = cy.at<float>(i, j);int v_cxx = cxx.at<float>(i, j);int v_cyy = cyy.at<float>(i, j);int v_cxy = cxy.at<float>(i, j);E_term.at<float>(i, j) = (v_cyy*v_cx*v_cx - 2 * v_cxy*v_cx*v_cy + v_cxx * v_cy*v_cy) / (std::pow((1 + v_cx * v_cx + v_cy * v_cy), 1.5));}}// 計算E_extcv::Mat E_ext = (wl*E_line + we * E_edge - wt * E_term);// 計算梯度cv::Mat fx, fy;cv::Sobel(E_ext, fx, E_ext.depth(), 1, 0, 1, 0.5, 0, cv::BORDER_CONSTANT);cv::Sobel(E_ext, fy, E_ext.depth(), 0, 1, 1, 0.5, 0, cv::BORDER_CONSTANT);cv::transpose(xs, xs);cv::transpose(ys, ys);int m = xs.rows;int n = 1;int mm = fx.cols;int nn = fx.rows;// 計算五對角狀矩陣，b(i)表示vi系數(i = i - 2 到 i + 2)double b[5];b[0] = beta;b[1] = -(alpha + 4 * beta);b[2] = 2 * alpha + 6 * beta;b[3] = b[1];b[4] = b[0];cv::Mat A = cv::Mat::eye(m, m, CV_32FC1);cv::Mat eyeMat0 = cv::Mat::eye(m, m, CV_32FC1);circRowShift(eyeMat0, 2);eyeMat0.convertTo(eyeMat0, CV_32FC1);A = b[0] * eyeMat0;cv::Mat eyeMat1 = cv::Mat::eye(m, m, CV_32FC1);circRowShift(eyeMat1, 1);eyeMat1.convertTo(eyeMat1, CV_32FC1);A = A + b[1] * eyeMat1;cv::Mat eyeMat2 = cv::Mat::eye(m, m, CV_32FC1);circRowShift(eyeMat2, 0);eyeMat2.convertTo(eyeMat2, CV_32FC1);A = A + b[2] * eyeMat2;cv::Mat eyeMat3 = cv::Mat::eye(m, m, CV_32FC1);circRowShift(eyeMat3, -1);eyeMat3.convertTo(eyeMat3, CV_32FC1);A = A + b[3] * eyeMat3;cv::Mat eyeMat4 = cv::Mat::eye(m, m, CV_32FC1);circRowShift(eyeMat4, -2);eyeMat4.convertTo(eyeMat4, CV_32FC1);A = A + b[4] * eyeMat4;// 計算矩陣的逆cv::Mat Ainv(A.size(), CV_32FC1);A = A + gamma * cv::Mat::eye(m, m, CV_32FC1);cv::invert(A, Ainv); // Computing Ainv cv::Mat srcImg = cv::imread("D:/endo_image.jpg");// 迭代更新曲線for (int i = 0; i < N; i++){cv::Mat intFx(fx.size(), CV_32FC1);cv::Mat intFy(fy.size(), CV_32FC1);cv::remap(fx, intFx, xs, ys, cv::INTER_LINEAR, cv::BORDER_CONSTANT);cv::remap(fy, intFy, xs, ys, cv::INTER_LINEAR, cv::BORDER_CONSTANT);cv::Mat ssx(xs.size(), CV_32FC1);cv::Mat ssy(ys.size(), CV_32FC1);for (int k = 0; k < xs.rows; k++){for (int l = 0; l < xs.cols; l++){ssx.at<float>(k, l) = gamma * xs.at<float>(k, l) - kappa * intFx.at<float>(k, l);ssy.at<float>(k, l) = gamma * ys.at<float>(k, l) - kappa * intFy.at<float>(k, l);}}// 更新曲線位置xs = Ainv * ssx;ys = Ainv * ssy;cv::Mat resultImg = srcImg.clone();for (int j = 0; j < xs.rows; j++){cv::Point center = cv::Point(xs.at<float>(j, 0), ys.at<float>(j, 0));cv::circle(resultImg, center, 4, cv::Scalar(0, 255, 0), 2);}// 顯示cv::imshow("result", resultImg);cv::waitKey(30);}return image; }

運行結果：

總結

以上是生活随笔為你收集整理的【主动轮廓模型（一）】《Snakes: Active Contour Models》算法原理与OpenCV实现的全部內容，希望文章能夠幫你解決所遇到的問題。

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