【主动轮廓模型 / Snake / Active Contour Model】算法原理与OpenCV实现_综合

文章目录

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 算法原理

首先需要人为地在图像上给出初始轮廓曲线，确切的说是一组用于控制曲线形状的控制点： $s\in[0,1]$ ，这些点收尾相连构成一个封闭的轮廓线。其中 $x (s)$ 和 $y (s)$ 分别表示每个控制点在图像中的坐标位置， $s$ 是以傅立叶变换形式描述边界的自变量，也可以理解为弧长。则Snake曲线的能量函数表示为：

$Esnake?=∫01Esnake(v(s))ds=∫01Eint(v(s))+Eimage(v(s))+Econ(v(s))ds=∫01Eint(v(s))+Eext(v(s))ds(1)\begin{aligned} 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{aligned} \tag{1}$

其中， $E_{int}$ 为内部能量， $E_{image}$ 为图像能量， $E_{con}$ 为外部约束能量。其中图像能量和外部约束能量统称为外部能量，即： $E_{ext}=E_{image}+E_{con}$ 。

2.1 内部能量 $E_{int}$

内部能量 $E_{int}$ 由保证曲线连续性的一阶导和保证曲线平滑的二阶导组成，表示为：

$Eint=12(α(s)∣vs(s)∣2+β(s)∣vss(s)∣2)(2)E_{int}=\frac{1}{2} (\alpha(s)|v_s(s)|^2+\beta(s)|v_{ss}(s)|^2) \tag{2}$

通过调整权值 $α(s)\alpha(s)$ 和 $β(s)\beta(s)$ 可以控制曲线的形状。例如将 $β(s)\beta(s)$ 置为0可以让曲线最终出现拐角，即曲线二阶不连续。

2.2 图像能量 $E_{image}$

在低层次的计算机视觉中，需要能够将轮廓吸引到特定图像特征的能量函数。原始的Snake轮廓模型中提出了3种不同的能量函数，分别将snake轮廓吸引到线、边和末端。完整的图像能量 $E_{image}$ 可以表示为这3个能量函数的权值组合，通过调整这3个权值，可以形成不同的轮廓形状。

$Eimage=ωlineEline+ωedgeEedge+ωtermEterm(3)E_{image}=\omega_{line}E_{line}+\omega_{edge}E_{edge}+\omega_{term}E_{term} \tag{3}$

（1）线函数 $E_{line}$

最简单直接且有用的图像能量函数是图像本身，即图像本身的灰度。令：

$Eline=I(x,y)(4)E_{line}=I(x,y) \tag{4}$

控制 $ωline\omega_{line}$ 的正负号可以控制轮廓被吸引到较暗的线或是较亮的线，也就是使轮廓试图靠近轮廓的最暗或最亮处。

然而，如果snake轮廓的一部分到达了一个低能量的图像特征位置，这个样条项将推动snake轮廓临近的部分朝着这个特征可能的延续方向移动，这会使其在一个最优的局部最小位置引入一个较大的能量。一种解决方案是允许snake轮廓和模糊能量函数平衡，然后慢慢降低模糊程度。

Marr和Hildreth在他们的论文里证明了“灰度的突然变化会在一阶导数中引起波峰或波谷，或在二阶导数中等效地引起零交点”。为了显示图像尺度空间连续性和Marr-Hildreth边缘检测理论的关系，snake模型中采用了模糊的边能量函数：

$Eline=?(Gσ??2I)2(6)E_{line}=-(G_{\sigma}*\nabla^2I)^2 \tag{6}$

其中 $GσG_{\sigma}$ 是标准差为 $σ\sigma$ 的高斯函数，该函数的最小值位于在Marr-Hildreth理论中被定义的 $Gσ??2IG_{\sigma}*\nabla^2I$ 的零交点处。在能量函数中加入这一项意味着snake轮廓在被吸引到零交点的同时，仍然受到它自己的平滑限制。

（2）边函数 $E_{edge}$

在图像上找边缘可以通过梯度来实现，令：

$Eedge=?∣?I(x,y)∣2(5)E_{edge}=-|\nabla I(x,y)|^2 \tag{5}$

在某个点上，梯度越大上式的能量越小，则snake轮廓将被吸引到梯度较大的区域。

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

为了找到轮廓的终止位置，作者将平滑过图像中等高线的曲率加入到能量函数中。令 $C(x,y)=Gσ(x,y)?I(x,y)C(x,y)=G_{\sigma}(x,y)*I(x,y)$ 是模糊后的图像， $θ=tan?1(CyCx)\theta=tan^{-1}(\frac{C_y}{C_x})$ 是梯度角， $n=(cosθ,sinθ)n=(cos\theta,sin\theta)$ 、 $n⊥=(?sinθ,cosθ)n_{\perp}=(-sin\theta,cos\theta)$ 是沿着和垂直梯度方向的单位向量。则 $C (x, y)$ 中等高线的曲率可以表示为：

$Eterm=?θ?n⊥=?C2/?n⊥2?C/?n=CyyCx2?2CxyCxCy+CxxCy2(Cx2+Cy2)3/2(6)\begin{aligned} E_{term}&=\frac{\partial \theta}{\partial n_{\perp}} \\ &=\frac{\partial C^2 / \partial n_{\perp}^2}{\partial C / \partial n} \\ &=\frac{C_{yy}C_x^2-2C_{xy}C_xC_y+C_{xx}C_y^2}{(C_x^2+C_y^2)^{3/2}} \end{aligned} \tag{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）的导数必须满足：

$αvss(s)?βvssss(s)??Eimage(v(s))??Econ(v(s))=0(7)\alpha v_{ss}(s) - \beta v_{ssss}(s) - \nabla E_{image}(v(s)) - \nabla E_{con}(v(s)) = 0 \tag{7}$

由 $v (s) = [x (s), y (s)]$ ，将上式改写成x和y两个方向有：

${αxss(s)?βxssss(s)??Eext?x=0αyss(s)?βyssss(s)??Eext?y=0(8)\begin{cases} \alpha x_{ss}(s) - \beta x_{ssss}(s) - \frac{\partial{E}_{ext}}{\partial x} = 0 \\ \alpha y_{ss}(s) - \beta y_{ssss}(s) - \frac{\partial{E}_{ext}}{\partial y} = 0 \end{cases} \tag{8}$

利用顺序连接的锚点的坐标差分来近似导数：

${xss(s)=x(s+1)+x(s?1)?2x(s)xssss(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))(9)\begin{cases} x_{ss}(s)=x(s+1)+x(s-1)-2x(s) \\ \begin{aligned} 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{aligned} \end{cases} \tag{9}$

将（9）带入（8），同时令 $fx=?Eext?x,fy=?Eext?yf_x=\frac{\partial E_{ext}}{\partial x}, f_y=\frac{\partial E_{ext}}{\partial y}$ 有：

${βx(s?2)?(α+4β)x(s?1)+(2α+6β)x(s)?(α+4β)x(s+1)+βx(s+2)+fx=0βy(s?2)?(α+4β)y(s?1)+(2α+6β)y(s)?(α+4β)y(s+1)+βy(s+2)+fy=0(10)\begin{cases} \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{cases} \tag{10}$

令：

${a=2α+6βb=?(α+4β)c=β(11)\begin{cases} a=2\alpha + 6\beta \\ b=-(\alpha+4\beta) \\ c=\beta \end{cases} \tag{11}$

则将（10）写成矩阵形式为：

${Ax+fx(x,y)=0Ay+fy(x,y)=0(12)\begin{cases} Ax+f_x(x,y)=0 \\ Ay+f_y(x,y)=0 \end{cases} \tag{12}$

其中A为五对角带状矩阵：

$\begin{bmatrix} a & b & c & \cdots & c & b \\ b & a & b & c & \cdots & c \\ c & b & a & b & c & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & c & b & a & b & c \\ c & \cdots & c & b & a & b \\ b & c & \cdots & c & b & a \end{bmatrix}, x= \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ \vdots \\ x_{n-1} \\ x_n \\ x_1 \end{bmatrix}, f_x= \begin{bmatrix} f_{x_1} \\ f_{x_2} \\ f_{x_3} \\ \vdots \\ f_{x_{n-1}} \\ f_{x_n} \\ f_{x_1} \end{bmatrix} \tag{13}$

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

${Axt+fx(xt?1,yt?1)=?γ(xt?xt?1)Ayt+fy(yt?1,yt?1)=?γ(yt?yt?1)(14)\begin{cases} Ax_t+f_x(x_{t-1},y_{t-1})=-\gamma(x_t-x_{t-1}) \\ Ay_t+f_y(y_{t-1},y_{t-1})=-\gamma(y_t-y_{t-1}) \end{cases} \tag{14}$

其中 $γ\gamma$ 是迭代步长，公式（14）的解为：

${xt=(A+γI)?1(xt?1?fx(xt?1,yt?1))yt=(A+γI)?1(yt?1?fy(xt?1,yt?1))(15)\begin{cases} 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{cases} \tag{15}$

由于 $A$ 为五对角带状矩阵，因此 $A+γIA+\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;
}

运行结果：

在这里插入图片描述

【主动轮廓模型 / Snake / Active Contour Model】算法原理与OpenCV实现

文章目录

1 概述

2 算法原理

2.1 内部能量 EintE_{int}Eint?

2.2 图像能量 EimageE_{image}Eimage?

（1）线函数ElineE_{line}Eline?

（2）边函数EedgeE_{edge}Eedge?

（3）末端函数EtermE_{term}Eterm?

2.3 外部能量 EconE_{con}Econ?