图像分割之cvSnakeImage

基本概念

        OpenCV之前有个cvSnakeImage函数实现了Snakes活动轮廓模型，但是cvSnakeImage函数在opencv2中已经被去掉。由于没有无法得知其依据的论文，这里根据源代码简单介绍它的原理。遍历每个轮廓的点，计算该点邻域内每个点的能量E，将能量最小的点代替当前点。E= alpha * Econt + beta * Ecurv + gamma * Eimg。Econ表示轮廓点的距离信息，Ecurv表示轮廓点的曲率信息，Eimg表示像素或者像素梯度。
        用OpenCV3.x重新实现后，发现貌似只对二值化图像分割有作用。

示例演示

        用OpenCV3.x重新实现cvSnakeImage函数。完整工程代码。

void SnakeImage::segment(std::vector<cv::Point> &points)
{
    drawCurve(points);
	if (points.size() < 3)
		return;

	int i = 0, j = 0, k = 0;
	//win_ is a domain, in which  search minimization of energy for each point
	//neighbors is the number of points in win_
	int neighbors = win_.height * win_.width;
	//the center of win_
	int centerx = win_.width >> 1;
	int centery = win_.height >> 1;

	//energy terms of neighbors of every points
	float *Econt = new float[neighbors];
	float *Ecurv = new float[neighbors];
	float *Eimg = new float[neighbors];
	float *E = new float[neighbors];

	//for two order difference
	int tx = 0, ty = 0;

	//for calculating distance
	float averagedistance = 0, diffx = 0, diffy = 0;

	//for calculating moving
	int offsetx = 0, offsety = 0;

	//temporary variables
	float tmp = 0.0, energy = 0.0;

	//if SNAKE_GRAD
	cv::Mat dx, dy, src;
	src_.convertTo(src, CV_16SC1);
	Sobel(src, dx, src.depth(), 1, 0, 1);
	Sobel(src, dy, src.depth(), 0, 1, 1);

	// the number of moved points in each iteration
	int movecount = 0;

	bool converged = false;
	int iteration = 0;
	while (!converged)
	{
		movecount = 0; // the number of moved points
		iteration++;

		//the average of distance between points
		averagedistance = 0;
		for (i = 0; i < points.size(); i++)
		{
			if (i != 0)
			{
				diffx = points[i - 1].x - points[i].x;
				diffy = points[i - 1].y - points[i].y;
			}
			else //n-1->0
			{
				diffx = points[points.size() - 1].x - points[i].x;
				diffy = points[points.size() - 1].y - points[i].y;
			}
			averagedistance += std::sqrt(std::pow(diffx, 2) + std::pow(diffy, 2));
		}
		averagedistance /= points.size();

		//search each point
		for (i = 0; i < points.size(); i++)
		{
			//calculate Econt
			float maxEcont = 0;
			float minEcont = SNAKE_BIG;
			//calculate reasonable search boundaries 
			int left = std::min(points[i].x, centerx);
			int right = std::min(src_.cols - 1 - points[i].x, centerx);
			int upper = std::min(points[i].y, centery);
			int bottom = std::min(src_.rows - 1 - points[i].y, centery);
			for (j = -upper; j <= bottom; j++)
			{
				for (k = -left; k <= right; k++)
				{
					if (i == 0)
					{
						diffx = points[points.size() - 1].x - (points[i].x + k);
						diffy = points[points.size() - 1].y - (points[i].y + j);
					}
					else
					{
						diffx = points[i - 1].x - (points[i].x + k);
						diffy = points[i - 1].y - (points[i].y + j);
					}
					energy = std::abs(averagedistance - std::sqrt(std::pow(diffx, 2) + std::pow(diffy, 2)));
					Econt[(j + centery) * win_.width + k + centerx] = energy;
					maxEcont = std::max(maxEcont, energy);
					minEcont = std::min(minEcont, energy);
				}
			}
			//if maxEcont == minEcont, set Econt all to zero
			//else normalize Econt
			tmp = maxEcont - minEcont;
			tmp = (tmp == 0.0f) ? 0.0f : (1.0 / tmp);
			for (k = 0; k < neighbors; k++)
			{
				Econt[k] = (Econt[k] - minEcont) * tmp;
			}

			//calculate Ecurv
			float maxEcurv = 0.0f;
			float minEcurv = SNAKE_BIG;
			for (j = -upper; j <= bottom; j++)
			{
				for (k = -left; k <= right; k++)
				{
					if (i == 0)
					{
						tx = points[points.size() - 1].x - 2 * (points[i].x + k) + points[i + 1].x;
						ty = points[points.size() - 1].y - 2 * (points[i].y + j) + points[i + 1].y;
 					}
					else if (i == points.size() - 1)
					{
						tx = points[i - 1].x - 2 * (points[i].x + k) + points[0].x;
						ty = points[i - 1].y - 2 * (points[i].y + j) + points[0].y;
					}
					else
					{
						tx = points[i - 1].x - 2 * (points[i].x + k) + points[i + 1].x;
						ty = points[i - 1].y - 2 * (points[i].y + j) + points[i + 1].y;
					}
					energy = std::pow(tx, 2) + std::pow(ty, 2);
					Ecurv[(j + centery) * win_.width + k + centerx] = energy;
					maxEcurv = std::max(maxEcurv, energy);
					minEcurv = std::min(minEcurv, energy);
				}
			}
			// normalize as above
			tmp = maxEcurv - minEcurv;
			tmp = (tmp == 0.0f) ? 0.0f : (1.0 / tmp);
			for (k = 0; k < neighbors; k++)
			{
				Ecurv[k] = (Ecurv[k] - minEcurv) * tmp;
			}

			//calculate Eimg
			float maxEimg = 0.0f;
			float minEimg = SNAKE_BIG;
			for (j = -upper; j <= bottom; j++)
			{
				for (k = -left; k <= right; k++)
				{
					if (scheme_ == SNAKE_GRAD)
					{
						energy = std::pow(dx.at<short>(points[i].y + j, points[i].x + k), 2);
						energy += std::pow(dy.at<short>(points[i].y + j, points[i].x + k), 2);
						Eimg[(j + centery) * win_.width + k + centerx] = energy;
					}
					else
					{
						energy = src_.at<uchar>(points[i].y + j, points[i].x + k)  ;
						Eimg[(j + centery) * win_.width + k + centerx] = energy;
					}
					maxEimg = std::max(maxEimg, energy);
					minEimg = std::min(minEimg, energy);
				}
			}
			// normalize as above 
			tmp = (maxEimg - minEimg);
 			tmp = (tmp == 0.0f) ? 0.0f : (1.0 / tmp);
			for (k = 0; k < neighbors; k++)
			{
				Eimg[k] = (minEimg - Eimg[k]) * tmp;
			}

			//calculate energy of neighbors
			for (k = 0; k < neighbors; k++)
			{
				E[k] = alpha_ * Econt[k] + beta_ * Ecurv[k] + gamma_ * Eimg[k];
			}
			//find minimum energy of neighbors
			float Emin = SNAKE_BIG;
			for (j = -upper; j <= bottom; j++)
			{
				for (k = -left; k <= right; k++)
				{
					if (E[(j + centery) * win_.width + k + centerx] < Emin)
					{
						Emin = E[(j + centery) * win_.width + k + centerx];
						offsetx = k;
						offsety = j;
					}
				}
			}
			//if moved
			if (offsetx || offsety)
			{
				points[i].x += offsetx;
				points[i].y += offsety;
				movecount++;
			}
		} // for each point

		//judge convergence
		converged = (movecount == 0);
		if ((termcriteria_.type & CV_TERMCRIT_ITER) && (iteration >= termcriteria_.maxCount))
			converged = 1;
		if ((termcriteria_.type & CV_TERMCRIT_EPS) && (movecount <= termcriteria_.epsilon))
			converged = 1;

		drawCurve(points);
	} //for each iteration

	//free resource
	delete[] Econt;
	delete[] Ecurv;
	delete[] Eimg;
	delete[] E;   
}

运行结果

相机

圆