【主动轮廓模型（一）】《Snakes: Active Contour Models》算法原理与OpenCV实现

---

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}^{\ast }& ={\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 }\ast {\nabla }^{2}I{)}^2\text{\hspace{1em}(6)}$$

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

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

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

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

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

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

为了找到轮廓的终止位置，作者将平滑过图像中等高线的曲率加入到能量函数中。令$C(x,y)=G_{\sigma }(x,y)\ast 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{\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}_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}^{\ast }$的最小值（局部极小值），公式（1）的导数必须满足：

$$\alpha v_{ss}(s)-\beta v_{ssss}(s)-\nabla E_{image}(v(s))-\nabla 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{\partial E_{ext}}{\partial x}=0\\ \alpha y_{ss}(s)-\beta y_{ssss}(s)-\frac{\partial E_{ext}}{\partial 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{\partial E_{ext}}{\partial x},f_y=\frac{\partial E_{ext}}{\partial 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=\left[\begin{array}{cccccc}a& b& c& \cdots & c& b\\ b& a& b& c& \cdots & c\\ c& b& a& b& c& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ \cdots & c& b& a& b& c\\ c& \cdots & c& b& a& b\\ b& c& \cdots & c& b& a\end{array}\right],x=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ⋮\\ x_{n-1}\\ x_n\\ x_1\end{array}\right],f_x=\left[\begin{array}{c}{f}_{x_1}\\ f_{x_2}\\ f_{x_3}\\ ⋮\\ f_{x_{n-1}}\\ f_{x_n}\\ f_{x_1}\end{array}\right]\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 intensities

	cv::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
		}
	}

	// 导数mask
	qDebug() << "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_term
	cv::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_ext
	cv::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;
}

运行结果：