opencv代码分析--canny算子

一、canny算子优点
相比较于常见的robert算子、sobel算子、梯度求边缘等方法，使用canny算子可以提取单像素二值化的细边缘，这是其他方法所不具备的

二、opencv canny算子提取边缘基本原理
opencv中canny算子的基本原理主要有以下几个步骤：

1. 梯度或sobel算子求边缘图
下图是用梯度求的边缘图，梯度值取水平梯度和垂直梯度之和，可以看到求出的边缘具有一定的宽度，在某一象素沿切线的法线方向，梯度值先增大，在边缘中心位置达到最大后，依次减小

2. 对1求的边缘进行法线方向的非极大值抑制，并借助双阀值来确定最终的边缘点
(1)极大值点判定
由梯度或者sobel算子求出的边缘，在法线方向，边缘点的梯度值变化都是前半段大致单调递增，边缘中线位置梯度值最大，后半段大致单调递减，所以只需要比较法线方向上该点与相邻前后两点的值，就可以确定该点是不是极大值点
(2)法线方向判定
opencv中，将边缘点的法线方向粗略的分为4个方向。边缘点的法线方向是这样确定的，以该边缘点的正切值来判断法线的大致方向：图中A区域，法线方向为近似水平方向，这时求法线方向的极大值，只需比较4邻域内的左右两个象素。C区域，法线方向为近似垂直方向，这时求法线方向的极大值，只需比较4邻域内的上下两个象素，B和D区域，分别取8邻域中的另外4个象素比较

(3)最终边缘点判定
边缘法线方向极大值点且该点梯度值大于高阀值，则这个点是最终边缘点
边缘法线方向极大值点且该点梯度值大于低阀值，可能是边缘点，如果这个可能的边缘点和最终边缘点8邻域内是相邻的，则可以判断它也是最终边缘点
边缘法线方向非极值点或者梯度值小于低阀值的点肯定不是最终边缘点

三、opencv代码分析

以下以opencv2.2为例，代码在modules\imgproc\src\canny.cpp文件中

CV_IMPL void cvCanny( const void* srcarr, void* dstarr,
                      double low_thresh, double high_thresh,
                      int aperture_size )
//参数说明：
//srcarr  			输入图像CvArr结构指针
//dstarr  			输出图像CvArr结构指针
//low_thresh  	双阀值抑制中的低阀值
//high_thresh  	双阀值抑制中的高阀值
//aperture_size	sobel算子模板大小
{
	cv::Ptr<CvMat> dx, dy;
	cv::AutoBuffer<char> buffer;
	std::vector<uchar*> stack;
	uchar **stack_top = 0, **stack_bottom = 0;
	
	CvMat srcstub, *src = cvGetMat( srcarr, &srcstub );
	CvMat dststub, *dst = cvGetMat( dstarr, &dststub );
	CvSize size;
	int flags = aperture_size;
	int low, high;
	int* mag_buf[3];
	uchar* map;
	ptrdiff_t mapstep;
	int maxsize;
	int i, j;
	CvMat mag_row;
	
	if( CV_MAT_TYPE( src->type ) != CV_8UC1 ||
		CV_MAT_TYPE( dst->type ) != CV_8UC1 )
		CV_Error( CV_StsUnsupportedFormat, "" );
	
	if( !CV_ARE_SIZES_EQ( src, dst ))
		CV_Error( CV_StsUnmatchedSizes, "" );
	
	if( low_thresh > high_thresh )
	{
		double t;
		CV_SWAP( low_thresh, high_thresh, t );
	}
	
	aperture_size &= INT_MAX;
	if( (aperture_size & 1) == 0 || aperture_size < 3 || aperture_size > 7 )
		CV_Error( CV_StsBadFlag, "" );
	
	size = cvGetMatSize( src );
	
	dx = cvCreateMat( size.height, size.width, CV_16SC1 );
	dy = cvCreateMat( size.height, size.width, CV_16SC1 );
	
	//sobel算子求水平边缘和垂直边缘，aperture_size表示模板大小，aperture_size=3表示3x3模板
	cvSobel( src, dx, 1, 0, aperture_size );
	cvSobel( src, dy, 0, 1, aperture_size );

	if( flags & CV_CANNY_L2_GRADIENT )
	{
		Cv32suf ul, uh;
		ul.f = (float)low_thresh;
		uh.f = (float)high_thresh;
		
		low = ul.i;
		high = uh.i;
	}
	else
	{
		low = cvFloor( low_thresh );
		high = cvFloor( high_thresh );
	}
		
	//分配一个连续的buffer, buffer宽为图像宽度+2，高度为图像高度+2+3。其中前三行用于
	//保存象素点的梯度值，因为之后的边缘梯度非极大值抑制过程中，为确定极大值而进行梯度
	//值比较，都是在8邻域范围内进行的，所以循环使用3行的内存能满足使用，这种做法节省了
	//内存。map指向第3行开始的地址，用于保存边缘梯度非极大值抑制过程中确定的图像中各点
	//的状态
	buffer.allocate( (size.width+2)*(size.height+2) + (size.width+2)*3*sizeof(int) );
	
	mag_buf[0] = (int*)(char*)buffer;
	mag_buf[1] = mag_buf[0] + size.width + 2;
	mag_buf[2] = mag_buf[1] + size.width + 2;
	
	//用于标记图像中所有点的状态：0-像素点可能是最终边缘， 1-象素点不是最终边缘，2-象素点是最终边缘		
	map = (uchar*)(mag_buf[2] + size.width + 2);  
	mapstep = size.width + 2;
	
	//初始化栈大小，该栈用于保存所有边缘梯度非极大值抑制过程中确定的状态为2的点		
	maxsize = MAX( 1 << 10, size.width*size.height/10 ); 
	stack.resize( maxsize );
	stack_top = stack_bottom = &stack[0];
	
	memset( mag_buf[0], 0, (size.width+2)*sizeof(int) );
	memset( map, 1, mapstep );
	memset( map + mapstep*(size.height + 1), 1, mapstep );
	
	
	#define CANNY_PUSH(d)    *(d) = (uchar)2, *stack_top++ = (d)
	#define CANNY_POP(d)     (d) = *--stack_top

	mag_row = cvMat( 1, size.width, CV_32F );
	
	// calculate magnitude and angle of gradient, perform non-maxima supression.
	// fill the map with one of the following values:
	//   0 - the pixel might belong to an edge
	//   1 - the pixel can not belong to an edge
	//   2 - the pixel does belong to an edge
	for( i = 0; i <= size.height; i++ )
	{
		int* _mag = mag_buf[(i > 0) + 1] + 1;
		float* _magf = (float*)_mag;
		const short* _dx = (short*)(dx->data.ptr + dx->step*i);
		const short* _dy = (short*)(dy->data.ptr + dy->step*i);
		uchar* _map;
		int x, y;
		ptrdiff_t magstep1, magstep2;
		int prev_flag = 0;

		if( i < size.height )
		{
			//计算梯度值
			_mag[-1] = _mag[size.width] = 0;
		
			if( !(flags & CV_CANNY_L2_GRADIENT) )
				for( j = 0; j < size.width; j++ )
					_mag[j] = abs(_dx[j]) + abs(_dy[j]);
			
			else
			{
				for( j = 0; j < size.width; j++ )
				{
					x = _dx[j]; y = _dy[j];
					_magf[j] = (float)std::sqrt((double)x*x + (double)y*y);
				}
			}
		}
		else
			memset( _mag-1, 0, (size.width + 2)*sizeof(int) );

		// at the very beginning we do not have a complete ring
		// buffer of 3 magnitude rows for non-maxima suppression
		if( i == 0 )
			continue;
		
		_map = map + mapstep*i + 1;
		_map[-1] = _map[size.width] = 1;
		
		_mag = mag_buf[1] + 1; // 无论mag_buf中的3行内存怎样进行循环滚动，mag_buf[1]永远指向中间的那行
		_dx = (short*)(dx->data.ptr + dx->step*(i-1));
		_dy = (short*)(dy->data.ptr + dy->step*(i-1));
		
		magstep1 = mag_buf[2] - mag_buf[1];
		magstep2 = mag_buf[0] - mag_buf[1];

		if( (stack_top - stack_bottom) + size.width > maxsize ) //栈不够大时扩展
		{
			int sz = (int)(stack_top - stack_bottom);
			maxsize = MAX( maxsize * 3/2, maxsize + 8 );
			stack.resize(maxsize);
			stack_bottom = &stack[0];
			stack_top = stack_bottom + sz;
		}

		for( j = 0; j < size.width; j++ )
		{
			#define CANNY_SHIFT 15
			//角度22.5度的tan值,左移12位避免浮点运算
			#define TG22  (int)(0.4142135623730950488016887242097*(1<<CANNY_SHIFT) + 0.5) 
			
			x = _dx[j];
			y = _dy[j];
			int s = x ^ y;
			int m = _mag[j];
			
			x = abs(x);
			y = abs(y);
			if( m > low ) //大于低阀值的，有可能是最终边缘点
			{
				int tg22x = x * TG22;
				
				//tan67°简化计算
				//tan67x = xtan(22.5°+45°)
				//  = x(tan45+tan22.5°)/(1-tan45°*tan22.5°)
				//  = x(1+tan22.5°)/(1-tan22.5°)
				//  = x((1+tan22.5°)(1+tan22.5°)/1-tan22.5°^2)
				//  = x((tan22.5°^2 + 2*tan22.5° + 1)/2tan22.5°)	
				//  = x(2 + (tan22.5°^2 + 1 - 2tan22.5°)/2tan22.5°)
				//  = x(2 + (tan22.5°^2 + 1 - 1 + tan22.5°^2)/2tan22.5°  
				//  = x(2 + tan22.5°)
				//  = (x + x) + xtan22.5°
				int tg67x = tg22x + ((x + x) << CANNY_SHIFT);
				
				y <<= CANNY_SHIFT;
                
				//y/x < tan22说明法线方向在0-22.5度之间，这时近似取法线为水平方向，
				//因此当前点的梯度值只需要和邻域的左右两个像素比较
				if( y < tg22x )  
				{
					if( m > _mag[j-1] && m >= _mag[j+1] )  //当前点是法线方向的极大值点
					{
						//因为扫描顺序是从上到下，从左到右，如果该点左边和正上方的点已经确认
						//为最终边缘点, 那么即使该点满足最终边缘点的条件,该点也不需要入栈,这
						//样做可以节省栈空间。将该点暂时设置成可能的最终边缘点，在后续的出栈
						//操作中，可以将它更正成最终边缘点
						if( m > high && !prev_flag && _map[j-mapstep] != 2 )    
						{
							//所有确定的最终边缘点入栈保存,用于之后判断可能的最终边缘点
							CANNY_PUSH( _map + j );  
							prev_flag = 1;
						}
						else
							_map[j] = (uchar)0;
						continue;
					}
				}
				//y/x > tan67说明法线方向在66.5-90度之间，这时近似取法线为垂直方向，
				//因此当前点的梯度值只需要和邻域的上下两个像素比较
				else if( y > tg67x ) 
				{
					if( m > _mag[j+magstep2] && m >= _mag[j+magstep1] )  //当前点是法线方向的极大值
					{
						if( m > high && !prev_flag && _map[j-mapstep] != 2 )
						{
							//所有确定的最终边缘点入栈保存,用于之后判断可能的最终边缘点
							CANNY_PUSH( _map + j );  
							prev_flag = 1;
						}
						else
							_map[j] = (uchar)0;
						continue;
					}
				}
				else  //法线方向在B区和D区
				{
					s = s < 0 ? -1 : 1;  //根据x, y的值判断法线方向是B区还是D区
					if( m > _mag[j+magstep2-s] && m > _mag[j+magstep1+s] )
					{
						if( m > high && !prev_flag && _map[j-mapstep] != 2 )
						{
							//所有确定的最终边缘点入栈保存,用于之后判断可能的最终边缘点
							CANNY_PUSH( _map + j );  
							prev_flag = 1;
						}
						else
							_map[j] = (uchar)0;
						continue;
					}
				}
			}
			prev_flag = 0;
			_map[j] = (uchar)1;
		}

		//更新指示状态，循环使用mag_buf所指向的3行内存
		_mag = mag_buf[0];
		mag_buf[0] = mag_buf[1];
		mag_buf[1] = mag_buf[2];
		mag_buf[2] = _mag;
	}

	// 栈中保存了经过上述运算，已经确定的最终边缘点；接下来只需要判断哪些可能
	//的最终边缘点,到底是不是最终边缘点就可以了, 根据边缘点判定规则, 与最终边
	//缘点8邻域相邻的可能边缘点，肯定是最终边缘点
	while( stack_top > stack_bottom )
	{
		uchar* m;
		if( (stack_top - stack_bottom) + 8 > maxsize )
		{
			int sz = (int)(stack_top - stack_bottom);
			maxsize = MAX( maxsize * 3/2, maxsize + 8 );
			stack.resize(maxsize);
			stack_bottom = &stack[0];
			stack_top = stack_bottom + sz;
		}
		
		CANNY_POP(m);
		
		//8邻域内点判断
		if( !m[-1] )
			CANNY_PUSH( m - 1 );
		if( !m[1] )
			CANNY_PUSH( m + 1 );
		if( !m[-mapstep-1] )
			CANNY_PUSH( m - mapstep - 1 );
		if( !m[-mapstep] )
			CANNY_PUSH( m - mapstep );
		if( !m[-mapstep+1] )
			CANNY_PUSH( m - mapstep + 1 );
		if( !m[mapstep-1] )
			CANNY_PUSH( m + mapstep - 1 );
		if( !m[mapstep] )
			CANNY_PUSH( m + mapstep );
		if( !m[mapstep+1] )
			CANNY_PUSH( m + mapstep + 1 );
	}

	//二值化，最终边缘点设置成255
	for( i = 0; i < size.height; i++ )
	{
		const uchar* _map = map + mapstep*(i+1) + 1;
		uchar* _dst = dst->data.ptr + dst->step*i;
		
		for( j = 0; j < size.width; j++ )
		  _dst[j] = (uchar)-(_map[j] >> 1);
	}
}

四、问题

canny算法中，高低阀值确定是比较关键的问题，算法中只是建议取低阀值为0.4倍的高阀值，实际应用中，可根据对直方图的分析，或者借鉴类似OTSU算法的思想，自适应获取阀值