文章目录
1 概述2 算法原理2.1 内部能量 EintE_{int}Eint2.2 图像能量 EimageE_{image}Eimage（1）线函数ElineE_{line}Eline（2）边函数EedgeE_{edge}Eedge（3）末端函数EtermE_{term}Eterm2.3 外部能量 EconE_{con}Econ3 模型求解4 算法实现（OpenCV3）1 概述
主动轮廓模型（也称Active Contour Model、Snake）是Kass等人在1988年提出的，该算法将图像分割问题转换为求解能量泛函最小值的问题。主要思路是通过构造能量泛函，经过算法迭代，轮廓曲线由初始位置逐渐向使能量函数最小（或局部极小）的图像边缘逼近，最终分割出目标。
2 算法原理
首先需要人为地在图像上给出初始轮廓曲线，确切的说是一组用于控制曲线形状的控制点：v(s)=[x(s),y(s)]s∈[0,1]v(s)=[x(s),y(s)] s\in[0,1]v(s)=[x(s),y(s)]s∈[0,1]，这些点收尾相连构成一个封闭的轮廓线。其中x(s)x(s)x(s)和y(s)y(s)y(s)分别表示每个控制点在图像中的坐标位置，sss是以傅立叶变换形式描述边界的自变量，也可以理解为弧长。则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} Esnake∗=∫01Esnake(v(s))ds=∫01Eint(v(s))+Eimage(v(s))+Econ(v(s))ds=∫01Eint(v(s))+Eext(v(s))ds(1)
其中，EintE_{int}Eint为内部能量，EimageE_{image}Eimage为图像能量，EconE_{con}Econ为外部约束能量。其中图像能量和外部约束能量统称为外部能量，即：Eext=Eimage+EconE_{ext}=E_{image}+E_{con}Eext=Eimage+Econ。
2.1 内部能量 EintE_{int}Eint
内部能量EintE_{int}Eint由保证曲线连续性的一阶导和保证曲线平滑的二阶导组成，表示为：
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} Eint=21(α(s)∣vs(s)∣2+β(s)∣vss(s)∣2)(2)
通过调整权值α(s)\alpha(s)α(s)和β(s)\beta(s)β(s)可以控制曲线的形状。例如将β(s)\beta(s)β(s)置为0可以让曲线最终出现拐角，即曲线二阶不连续。
2.2 图像能量 EimageE_{image}Eimage
在低层次的计算机视觉中，需要能够将轮廓吸引到特定图像特征的能量函数。原始的Snake轮廓模型中提出了3种不同的能量函数，分别将snake轮廓吸引到线、边和末端。完整的图像能量 EimageE_{image}Eimage可以表示为这3个能量函数的权值组合，通过调整这3个权值，可以形成不同的轮廓形状。
Eimage=ωlineEline+ωedgeEedge+ωtermEterm(3)E_{image}=\omega_{line}E_{line}+\omega_{edge}E_{edge}+\omega_{term}E_{term} \tag{3} Eimage=ωlineEline+ωedgeEedge+ωtermEterm(3)
（1）线函数ElineE_{line}Eline
最简单直接且有用的图像能量函数是图像本身，即图像本身的灰度。令：
Eline=I(x,y)(4)E_{line}=I(x,y) \tag{4} Eline=I(x,y)(4)
控制ωline\omega_{line}ωline的正负号可以控制轮廓被吸引到较暗的线或是较亮的线，也就是使轮廓试图靠近轮廓的最暗或最亮处。
然而，如果snake轮廓的一部分到达了一个低能量的图像特征位置，这个样条项将推动snake轮廓临近的部分朝着这个特征可能的延续方向移动，这会使其在一个最优的局部最小位置引入一个较大的能量。一种解决方案是允许snake轮廓和模糊能量函数平衡，然后慢慢降低模糊程度。
Marr和Hildreth在他们的论文里证明了“灰度的突然变化会在一阶导数中引起波峰或波谷，或在二阶导数中等效地引起零交点”。为了显示图像尺度空间连续性和Marr-Hildreth边缘检测理论的关系，snake模型中采用了模糊的边能量函数：
Eline=−(Gσ∗∇2I)2(6)E_{line}=-(G_{\sigma}*\nabla^2I)^2 \tag{6} Eline=−(Gσ∗∇2I)2(6)
其中GσG_{\sigma}Gσ是标准差为σ\sigmaσ的高斯函数，该函数的最小值位于在Marr-Hildreth理论中被定义的Gσ∗∇2IG_{\sigma}*\nabla^2IGσ∗∇2I的零交点处。在能量函数中加入这一项意味着snake轮廓在被吸引到零交点的同时，仍然受到它自己的平滑限制。
（2）边函数EedgeE_{edge}Eedge
在图像上找边缘可以通过梯度来实现，令：
Eedge=−∣∇I(x,y)∣2(5)E_{edge}=-|\nabla I(x,y)|^2 \tag{5} Eedge=−∣∇I(x,y)∣2(5)
在某个点上，梯度越大上式的能量越小，则snake轮廓将被吸引到梯度较大的区域。
（3）末端函数EtermE_{term}Eterm
为了找到轮廓的终止位置，作者将平滑过图像中等高线的曲率加入到能量函数中。令C(x,y)=Gσ(x,y)∗I(x,y)C(x,y)=G_{\sigma}(x,y)*I(x,y)C(x,y)=Gσ(x,y)∗I(x,y)是模糊后的图像，θ=tan−1(CyCx)\theta=tan^{-1}(\frac{C_y}{C_x})θ=tan−1(CxCy)是梯度角，n=(cosθ,sinθ)n=(cos\theta,sin\theta)n=(cosθ,sinθ)、n⊥=(−sinθ,cosθ)n_{\perp}=(-sin\theta,cos\theta)n⊥=(−sinθ,cosθ)是沿着和垂直梯度方向的单位向量。则C(x,y)C(x,y)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} Eterm=∂n⊥∂θ=∂C/∂n∂C2/∂n⊥2=(Cx2+Cy2)3/2CyyCx2−2CxyCxCy+CxxCy2(6)
通过组合EedgeE_{edge}Eedge和EtermE_{term}Eterm，我们能够创建一个被边和末端吸引的snake曲线。
2.3 外部能量 EconE_{con}Econ
外部能量EconE_{con}Econ来自于外部的约束力。在论文原文中，作者给了一个外部约束的例子， 包括外部的固定点、连接两条snake曲线的锚点或鼠标拖动的点。例如，为了在x1x_1x1和x2x_2x2点之间创建一个连接的弹性力，就可以将−k(x1−x2)2-k(x_1-x_2)^2−k(x1−x2)2添加到外部能量EconE_{con}Econ中。
3 模型求解
由欧拉方程，求解能量Esnake∗E_{snake}^*Esnake∗的最小值（局部极小值），公式（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} αvss(s)−βvssss(s)−∇Eimage(v(s))−∇Econ(v(s))=0(7)
由v(s)=[x(s),y(s)]v(s)=[x(s),y(s)]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)−βxssss(s)−∂x∂Eext=0αyss(s)−βyssss(s)−∂y∂Eext=0(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} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧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)
将（9）带入（8），同时令fx=∂Eext∂x,fy=∂Eext∂yf_x=\frac{\partial E_{ext}}{\partial x}, f_y=\frac{\partial E_{ext}}{\partial y}fx=∂x∂Eext,fy=∂y∂Eext有：
{β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} {β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)
令：
{a=2α+6βb=−(α+4β)c=β(11)\begin{cases} a=2\alpha + 6\beta \\ b=-(\alpha+4\beta) \\ c=\beta \end{cases} \tag{11} ⎩⎪⎨⎪⎧a=2α+6βb=−(α+4β)c=β(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} {Ax+fx(x,y)=0Ay+fy(x,y)=0(12)
其中A为五对角带状矩阵：
A=[abc⋯cbbabc⋯ccbabc⋯⋮⋮⋮⋮⋮⋮⋯cbabcc⋯cbabbc⋯cba],x=[x1x2x3⋮xn−1xnx1],fx=[fx1fx2fx3⋮fxn−1fxnfx1](13)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} A=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡abc⋮⋯cbbab⋮c⋯ccba⋮bc⋯⋯cb⋮abcc⋯c⋮babbc⋯⋮cba⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤,x=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡x1x2x3⋮xn−1xnx1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤,fx=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡fx1fx2fx3⋮fxn−1fxnfx1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤(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} {Axt+fx(xt−1,yt−1)=−γ(xt−xt−1)Ayt+fy(yt−1,yt−1)=−γ(yt−yt−1)(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} {xt=(A+γI)−1(xt−1−fx(xt−1,yt−1))yt=(A+γI)−1(yt−1−fy(xt−1,yt−1))(15)
由于AAA为五对角带状矩阵，因此A+γIA+\gamma IA+γ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;}
运行结果：
