图像处理之其他杂项（四）之cvSnakeImage()函数代码升级，从C接口到C++接口：snakeImage()

cvSnakeImage()函数代码升级，从C接口到C++接口：snakeImage()

1.snake算法存在于老版本的opencv中，采用的是C接口，函数为cvSnakeImage()，新版本opencv2中取消了这个函数，现对于网上存在的老代码进行升级为C++接口，适合opencv2中使用，对于算法的具体实现原理未具体明白，只是在原有代码的基础上按葫芦画瓢，对于相应代码、函数进行了替换，提供给大家参考。

2.因为对于SNAKE未真正深入领悟，更改代码肯定有错误之处，且现阶段运行结果存在问题，未实现最终理想效果，具体原因未知。

3.其中有些注释为老源代码带有。

4.老源代码：蚂蚁搬家　http://blog.csdn.net/hongxingabc/article/details/51606520　　《主动轮廓线模型Snake模型简介&openCV中cvSnakeImage()函数代码分析》

5.在此之前有个版本   对于老源代码做了些许兼容性改动，使代码可以在opencv2中运行，但总体还是C接口    http://blog.csdn.net/coming_is_winter/article/details/73010825      《cvSnakeImage改进升级兼容》。

/*
1.snake算法存在于老版本的opencv中，采用的是C接口，函数为cvSnakeImage()，新版本opencv2中取消了这个函数，现对于网上存在的老代码进行升级为C++接口，适合opencv2中使用，对于算法的具体实现原理未具体明白，只是在原有代码的基础上按葫芦画瓢，对于相应代码、函数进行了替换，提供给大家学习参考。
2.因为对于SNAKE未真正深入领悟，更改代码肯定有错误之处，且现阶段运行结果存在问题，未实现最终理想效果，具体原因未知。
3.其中有些注释为老源代码带有。
4.老源代码：蚂蚁搬家　http://blog.csdn.net/hongxingabc/article/details/51606520　　《主动轮廓线模型Snake模型简介&openCV中cvSnakeImage()函数代码分析》
5.在此之前有个版本   对于老源代码做了些许兼容性改动，使代码可以在opencv2中运行，但总体还是C接口    http://blog.csdn.net/coming_is_winter/article/details/73010825      《cvSnakeImage改进升级兼容》。
*/

#include<iostream> 
#include<opencv.hpp>
using namespace std;
using namespace cv;
#define SNAKE_BIG 2.e+38f  
#define SNAKE_IMAGE 1  
#define SNAKE_GRAD  2  
#define VALUE 30

int icvSnake8uC1R(unsigned char *src,   //原始图像数据  
	int srcStep,         //每行的字节数  
	Size roi,         //图像尺寸  
	Point * pt,       //轮廓点(变形对象)  
	int n,            //轮廓点的个数        //γ的值，同α    
	Size win,       //每个点用于搜索的最小的领域大小,宽度为奇数  
	TermCriteria criteria,   //递归迭代终止的条件准则  
	int scheme)         //确定图像能量场的数据选择，1为灰度，2为灰度梯度  
{
	
	int i, j, k;
	int neighbors = win.height*win.width;    //当前点领域中点的个数  											  
	int centerx = win.width >> 1;
	int centery = win.height >> 1;
	float invn;        //n 的倒数  
	int iteration = 0;     //迭代次数  
	int converged = 0;      //收敛标志，0为非收敛  

							  
	float *Econt=new float;    //能量 
	float *Ecurv= new float;   //轮廓曲线能量  
	float *Eimg= new float;    //图像能量  
	float *E= new float;      //  


	float *gradient = NULL;
	uchar *map = NULL;
	int map_width = ((roi.width - 1) >> 3) + 1;
	int map_height = ((roi.height - 1) >> 3) + 1;
#define WTILE_SIZE 8  
#define TILE_SIZE (WTILE_SIZE + 2)         
	short dx[TILE_SIZE*TILE_SIZE], dy[TILE_SIZE*TILE_SIZE];
	Mat _dx = Mat(TILE_SIZE, TILE_SIZE, CV_16SC1, dx);
	Mat _dy = Mat(TILE_SIZE, TILE_SIZE, CV_16SC1, dy);
	Mat _src = Mat(roi.height, roi.width, CV_8UC1, src);

	invn = 1 / ((float)n);        //轮廓点数n的倒数，用于求平均  

	if (scheme ==SNAKE_GRAD)
	{
	
		memset((void *)map, 0, map_width * map_height);
	}

	while (!converged)     
	{
		float ave_d = 0;  
		int moved = 0;     

		converged = 0;       
		iteration++;       

							
							
							 
		for (i = 1; i < n; i++)
		{
			int diffx = pt[i - 1].x - pt[i].x;
			int diffy = pt[i - 1].y - pt[i].y;
			ave_d += sqrt((float)(diffx * diffx + diffy * diffy));
					
		}
		
		//再加上从点n-1到点0的距离,形成回路轮廓  
		ave_d += sqrt((float)((pt[0].x - pt[n - 1].x) *
			(pt[0].x - pt[n - 1].x) +
			(pt[0].y - pt[n - 1].y) * (pt[0].y - pt[n - 1].y)));
		//求平均，得出平均距离  
		ave_d *= invn;
	
		
		//对于每个轮廓点进行特定循环迭代求解  
		for (i = 0; i < n; i++)
		{
			
			//初始化各个能量  
			float maxEcont = 0;
			float maxEcurv = 0;
			float maxEimg = 0;
			float minEcont =SNAKE_BIG;
			float minEcurv =SNAKE_BIG;
			float minEimg =SNAKE_BIG;
			float Emin =SNAKE_BIG;
			//初始化变形后轮廓点的偏移量  
			int offsetx = 0;
			int offsety = 0;
			float tmp;

			//计算边界  
			/* compute bounds */
			//计算合理的搜索边界，以防领域搜索超过ROI图像的范围  
			int left = MIN(pt[i].x, win.width >> 1);
			int right = MIN(roi.width - 1 - pt[i].x, win.width >> 1);
			int upper = MIN(pt[i].y, win.height >> 1);
			int bottom = MIN(roi.height - 1 - pt[i].y, win.height >> 1);
			//初始化Econt  
			maxEcont = 0;
			minEcont =SNAKE_BIG;
			//在合理的搜索范围内进行Econt的计算  
			for (j = -upper; j <= bottom; j++)
			{
				for (k = -left; k <= right; k++)
				{
					int diffx, diffy;
					float energy;
					//在轮廓点集的首尾相接处作相应处理,求轮廓点差分  
					if (i == 0)
					{
						diffx = pt[n - 1].x - (pt[i].x + k);
						diffy = pt[n - 1].y - (pt[i].y + j);
					}
					else
						//在其他地方作一般处理  

					{
						diffx = pt[i - 1].x - (pt[i].x + k);
						diffy = pt[i - 1].y - (pt[i].y + j);
					}
					//将邻域陈列坐标转成Econt数组的下标序号，计算邻域中每点的Econt  
					//Econt的值等于平均距离和此点和上一点的距离的差的绝对值（这是怎么来的？）  
					Econt[(j + centery) * win.width + k + centerx] = energy =
						(float)fabs(ave_d -
							sqrt((float)(diffx * diffx + diffy * diffy)));
					//求出所有邻域点中的Econt的最大值和最小值  
					maxEcont = MAX(maxEcont, energy);
					minEcont = MIN(minEcont, energy);
				}
			}
			
			//求出邻域点中最大值和最小值之差,并对所有的邻域点的Econt进行标准归一化,若最大值最小  
			//相等，则邻域中的点Econt全相等，Econt归一化束缚为0  
			tmp = maxEcont - minEcont;
			tmp = (tmp == 0) ? 0 : (1 / tmp);
			for (k = 0; k < neighbors; k++)
			{
				Econt[k] = (Econt[k] - minEcont) * tmp;
			}


			//计算每点的Ecurv  
			/*  Calculate Ecurv */
			maxEcurv = 0;
			minEcurv =SNAKE_BIG;
			for (j = -upper; j <= bottom; j++)
			{
				for (k = -left; k <= right; k++)
				{
					int tx, ty;
					float energy;
					//第一个点的二阶差分  
					if (i == 0)
					{
						tx = pt[n - 1].x - 2 * (pt[i].x + k) + pt[i + 1].x;
						ty = pt[n - 1].y - 2 * (pt[i].y + j) + pt[i + 1].y;
					}
					//最后一个点的二阶差分  
					else if (i == n - 1)
					{
						tx = pt[i - 1].x - 2 * (pt[i].x + k) + pt[0].x;
						ty = pt[i - 1].y - 2 * (pt[i].y + j) + pt[0].y;
					}
					//其余点的二阶差分  
					else
					{
						tx = pt[i - 1].x - 2 * (pt[i].x + k) + pt[i + 1].x;
						ty = pt[i - 1].y - 2 * (pt[i].y + j) + pt[i + 1].y;
					}
					//转换坐标为数组序号，并求各点的Ecurv的值,二阶差分后取平方  
					Ecurv[(j + centery) * win.width + k + centerx] = energy =
						(float)(tx * tx + ty * ty);
					//取最小的Ecurv和最大的Ecurv  
					maxEcurv = MAX(maxEcurv, energy);
					minEcurv = MIN(minEcurv, energy);
				}
			}
			
			//对Ecurv进行标准归一化  
			tmp = maxEcurv - minEcurv;
			tmp = (tmp == 0) ? 0 : (1 / tmp);
			for (k = 0; k < neighbors; k++)
			{
				Ecurv[k] = (Ecurv[k] - minEcurv) * tmp;
			}

				//求Eimg  
			/* Calculate Eimg */
			for (j = -upper; j <= bottom; j++)
			{
				for (k = -left; k <= right; k++)
				{
					float energy;
					//若采用灰度梯度数据  
					if (scheme ==SNAKE_GRAD)
					{
						/* look at map and check status */
						int x = (pt[i].x + k) / WTILE_SIZE;
						int y = (pt[i].y + j) / WTILE_SIZE;
						//若此处的图像能量还没有获取，则对此处对应的图像分块进行图像能量的求解  
						if (map[y * map_width + x] == 0)
						{
							int l, m;

							/* evaluate block location */
							//计算要进行梯度算子处理的图像块的位置  
							int upshift = y ? 1 : 0;
							int leftshift = x ? 1 : 0;
							int bottomshift = MIN(1, roi.height - (y + 1)*WTILE_SIZE);
							int rightshift = MIN(1, roi.width - (x + 1)*WTILE_SIZE);
							//图像块的位置大小（由于原ROI不一定是8的倍数，所以图像块会大小不一）  
							Rect g_roi = { x*WTILE_SIZE - leftshift, y*WTILE_SIZE - upshift,
								leftshift + WTILE_SIZE + rightshift, upshift + WTILE_SIZE + bottomshift };
							Mat _src1;
							//cvGetSubArr(&_src, &_src1, g_roi);  //得到图像块的数据  
																//分别对图像的X方向和Y方向进行梯度算子  
							Sobel(_src1, _dx, 1, 0, 1);
							Sobel(_src1, _dy, 1, 0, 1);
							//pX.process(&_src1, &_dx);
							//pY.process(&_src1, &_dy);
							//求分块区域中的每个点的梯度  
							for (l = 0; l < WTILE_SIZE + bottomshift; l++)
							{
								for (m = 0; m < WTILE_SIZE + rightshift; m++)
								{
									gradient[(y*WTILE_SIZE + l) * roi.width + x*WTILE_SIZE + m] =
										(float)(dx[(l + upshift) * TILE_SIZE + m + leftshift] *
											dx[(l + upshift) * TILE_SIZE + m + leftshift] +
											dy[(l + upshift) * TILE_SIZE + m + leftshift] *
											dy[(l + upshift) * TILE_SIZE + m + leftshift]);
								}
							}
							//map相应位置置1表示此处图像能量已经获取  
							map[y * map_width + x] = 1;
						}
						//以梯度数据作为图像能量  
						Eimg[(j + centery) * win.width + k + centerx] = energy =
							gradient[(pt[i].y + j) * roi.width + pt[i].x + k];
					}
					else
					{
						//以灰度作为图像能量  
						Eimg[(j + centery) * win.width + k + centerx] = energy =
							src[(pt[i].y + j) * srcStep + pt[i].x + k];
					}
					//获得邻域中最大和最小的图像能量  
					maxEimg = MAX(maxEimg, energy);
					minEimg = MIN(minEimg, energy);
				}
			}
			
			//Eimg的标准归一化  
			tmp = (maxEimg - minEimg);
			tmp = (tmp == 0) ? 0 : (1 / tmp);

			for (k = 0; k < neighbors; k++)
			{
				Eimg[k] = (minEimg - Eimg[k]) * tmp;
			}
			//加入系数  
			/* locate coefficients */
			
			
			/* Find Minimize point in the neighbors */
			//求得每个邻域点的Snake能量  
			for (k = 0; k < neighbors; k++)
			{
				E[k] = (float)( 40 * Econt[k] + 50 * Ecurv[k] +60 * Eimg[k]);
			
			}
			

			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 (offsetx || offsety)
			{
				pt[i].x += offsetx;
				pt[i].y += offsety;
				moved++;
			}
		}
		
		//各个轮廓点迭代计算完成后，如果没有移动的点了，则收敛标志位有效，停止迭代  
		converged = (moved == 0);
		//达到最大迭代次数时，收敛标志位有效，停止迭代  
		if ((criteria.type & 1) && (iteration >= criteria.maxCount))
			converged = 1;
		//到大相应精度时，停止迭代（与第一个条件有相同效果）  
		if ((criteria.type & 2) && (moved <= criteria.epsilon))
			converged = 1;
	}
	
	return 0;
}


void snakeImage( Mat& src,Point* points,
	int length,	 Size win,
	TermCriteria criteria, int calcGradient)
{

	Size size;
	int step;

	uchar* data = src.data;
	step = src.cols;
	size.height = src.rows;
	size.width = src.cols;

	icvSnake8uC1R(data, step, size, points, length, win, criteria,
		1);
}

/*
以上代码为核心SNAKE算法代码，涉及到能量场问题，感觉头大，没有研究透，主要是根据设定，将预先找出的轮廓点集根据算法进行特定移动，能量场稳定后轮廓固定为收缩后效果，运算结束
*/

Mat image;  
Mat image2;  
Mat image3;
using namespace std;
int main()

{
	image3 = imread("E://素材//圆.jpg", 1);
	image2 = imread("E://素材//圆.jpg", 1); 
	image = imread("E://素材//圆.jpg", 0);

	threshold(image, image, 147, 205, CV_THRESH_BINARY_INV); //请选取适当的阈值，进行操作     
	imshow("二值图", image);
	vector<vector<Point>> contours;

	findContours(image, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_SIMPLE);//查找轮廓
	int a = 0;
	int j = 0;
	for (int i = 0; i < contours.size(); i++)//这里的主要目的是选取出图像中的最大轮廓
	{
		
		if (a < contours[i].size())
		{
			j = i;
			a = contours[i].size();
		}
	}

	int length = contours[j].size();//最大轮廓的点集数
	for (int i = 0; i < length; i++)//画出原始轮廓
	{
		int l = (i + 1) % length;
		line(image2, contours[j][i], contours[j][l], Scalar(255, 0, 0), 2, 8, 0);

	}
	
	
	imshow("蓝色为原始轮廓", image2);

	
	cout << "轮廓点的数量是 "<<length<<"个。程序已继续....\n\r";
    waitKey(5000);
	Point* points = new Point[length];
	for (int i = 0; i < length; i++)
	{
		points[i] = contours[j][i];
		cout <<"这是snake前第"<<i<<"个点:"<<"  "<< points[i]<<".\n\r";
		waitKey(20);
	}

	Size size;
	size.width =5;//这里我理解是后续能量场操作的窗口大小
	size.height =5;
	TermCriteria criteria;
	criteria.type = CV_TERMCRIT_ITER;
	criteria.maxCount = 1000;
	criteria.epsilon = 0.1;

	snakeImage(image, points,length, size, criteria, 0);//主要操作
	//int length1 = points.
	for (int i = 0; i < length; i++)
	{
		cout << "这是snake后第" << i << "个点:" << "  " << points[i] << ".\n\r";
		waitKey(20);
     }
		
	for (int i = 0; i < length; i++)//画出处理后的轮廓
	{
		int j = (i + 1) % length;
		line(image3, points[i], points[j], Scalar( 0,255, 0), 2, 8, 0);
	
	}
	imshow("绿色为SNAKE轮廓", image3);

	waitKey(0);
	//return 0;
}

最终实验结果如下图所示，有点怪，最理想的情况应该是比原轮廓更小的圆才对。