Canny边缘检测的c++&&OpenCV实现

Canny边缘检测算法是一种比较常见且效果较好的边缘检测算法，其优点在于得到的边缘为单像素边缘，但是速度较Sobel等边缘检测算法较慢。网上Canny算法的算法介绍很多，此处不仔细写了，主要分为五步，分别是：高斯去噪，梯度计算、梯度方向计算、非极大值抑制和双阈值边缘抑制。其中，非极大值抑制保证了边缘检测结果为单像素边缘，双阈值边缘抑制的作用为去掉孤立的低阈值边缘。以下代码为Canny算子的简单实现，其中包含了Sobel算子的实现和高斯去噪的简单版本。

高斯去噪：

cv::Mat GaussBlur(__in__ const cv::Mat& src) {
	int row = src.rows;
	int col = src.cols;
	cv::Mat Gauss_x(row, col, CV_8UC1);
	cv::Mat des(row, col, CV_8UC1);

	for (int i = 0; i < row; ++i) {
		for (int j = 1; j < col - 1; ++j) {
			int temp = ((int)src.data[i * col + j] << 1) + src.data[i * col + j - 1] + src.data[i * col + j + 1];
			Gauss_x.data[i * col + j] = (temp >> 2);
		}
	}
	for (int i = 1; i < row - 1; ++i) {
		for (int j = 1; j < col - 1; ++j) {
			int temp = ((int)Gauss_x.data[i * col + j] << 1) + Gauss_x.data[i * col + j - col] + Gauss_x.data[i * col + j + col];
			des.data[i * col + j] = (temp >> 2);
		}
	}

	return des;
}

 

Sobel算子：

cv::Mat SobelOperator(__in__ const cv::Mat& src, __out__ cv::Mat& Gx, __out__ cv::Mat& Gy) {

	int row = src.rows;
	int col = src.cols;

	cv::Mat des(row, col, CV_16UC1);
	Gx = cv::Mat(row, col, CV_16SC1);
	Gy = cv::Mat(row, col, CV_16SC1);

	for (int i = 1; i < row - 1; ++i) {
		for (int j = 1; j < col - 1; ++j) {
			int temp = 0;
			auto pdata = &src.data[i * col + j];

			temp = *(pdata - col - 1) + (((int)*(pdata - col)) << 1) + *(pdata - col + 1);
			temp -= (*(pdata + col - 1) + (((int)*(pdata + col)) << 1) + *(pdata + col + 1));
			Gy.at<short>(i, j) = temp;
			auto gy = temp;

			temp = *(pdata - col - 1) + (((int)*(pdata - 1)) << 1) + *(pdata + col - 1);
			temp -= *(pdata - col + 1) + (((int)*(pdata + 1)) << 1) + *(pdata + col + 1);
			Gx.at<short>(i, j) = temp;
			auto gx = temp;

			des.at<ushort>(i, j) = (ushort)sqrt(gy * gy + gx * gx);
		}
	}

	return des;
}

Canny算子：

cv::Mat CannyEdgeDetector(__in__ const cv::Mat& src, __in__ const float up, __in__ const float down) {
	
	const int row = src.rows;
	const int col = src.cols;
	const double A2R = 1.0 * 180.0 / 3.1415926535;
	//Step 1: Gaussian filtering to remove noise
	cv::Mat Gauss = GaussBlur(src);

	//Step 2: Sobel Operator
	cv::Mat Gx;//CV_16UC1
	cv::Mat Gy;//CV_16UC1
	cv::Mat G = SobelOperator(Gauss, Gx, Gy);//CV_16SC1

	//Step 3: Finding Gradient angle theta and tracing the edge in the image using theta
	cv::Mat theta = cv::Mat(row, col, CV_32FC1);
	for (int i = 1; i < row; ++i) {
		for (int j = 1; j < col; ++j) {
			 theta.at<float>(i, j) = (float)atan2(Gy.at<short>(i, j), Gx.at<short>(i, j) );
		}
	}

	uint max_pixel_value = 0;
	cv::Mat non_max = G.clone();
	//Step 4: Non-maximum suppression
	for (int i = 1; i < row - 1; ++i) {
		for (int j = 1; j < col - 1; ++j) {
			ushort* pdata = &G.at<ushort>(i, j);
			ushort* pnew = &non_max.at<ushort>(i, j);
			float th = theta.at<float>(i, j);
			if (th < 0)th += CV_PI;
			float angle = th * A2R;
			if (angle < 22.5 || angle >= 157.5) {
				if (*pdata < *(pdata + 1) || *pdata < *(pdata - 1))*pnew = 0;
			}
			else if (angle >= 22.5 && angle < 67.5){
				if (*pdata < *(pdata + col - 1) || *pdata < *(pdata - col + 1))*pnew = 0;
			}
			else if (angle >= 67.5 && angle < 112.5) {
				if (*pdata < *(pdata + col) || *pdata < *(pdata - col))*pnew = 0;
			}
			else if (angle >= 112.5 && angle < 157.5) {
				if (*pdata < *(pdata + 1 + col) || *pdata < *(pdata - 1 - col))*pnew = 0;
			}
			else *pnew = *pdata;
			
			if (max_pixel_value < *pdata)max_pixel_value = *pdata;
		}
	}

	cv::Mat result(row, col, CV_8UC1, cv::Scalar(0));
	//Step 5: Hysteresis Thresholding using two thresholds
	uint* account = new uint[max_pixel_value + 1]();
	for (int i = 1; i < row - 1; ++i) {
		for (int j = 1; j < col - 1; ++j) {
			account[non_max.at<ushort>(i, j)]++;
		}
	}
	//chack the thresholds' legality
	if (up < down) {
		cout << "up is bigger than down\n";
		throw param_exception("up is bigger than down");
	}
	if (up > 1) {
		cout << "up is bigger than 1\n";
		throw param_exception("up is bigger than 1");
	}
	if (down < 0) {
		cout << "down is lower than 0\n";
		throw param_exception("down is lower than 0");
	}
	
	int all = (row - 2) * (col - 2) - account[0];
	const int up_account = all * up;
	const int down_account = all * down;
	int up_num, down_num;
	while (all > up_account) {
		all -= account[max_pixel_value--];
	}
	up_num = ++max_pixel_value;
	while (all > down_account) {
		all -= account[max_pixel_value--];
	}
	down_num = ++max_pixel_value;

	for (int i = 1; i < row - 1; ++i) {
		for (int j = 1; j < col - 1; ++j) {
			if (non_max.at<ushort>(i, j) >= up_num)result.data[i * col + j] = 255;
			else if(non_max.at<ushort>(i, j) <= down_num)result.data[i * col + j] = 0;
			else result.data[i * col + j] = 127;
		}
	}
	for (int i = 1; i < row - 1; ++i) {
		for (int j = 1; j < col - 1; ++j) {
			if (127 == result.data[i * col + j]) {
				uchar* pdata = &result.data[i * col + j];
				if (*(pdata - 1) || *(pdata + 1) || *(pdata - 1 - col) || *(pdata - 1 + col) || *(pdata + 1 - col) || *(pdata + 1 + col) || *(pdata - col) || *(pdata + col))*pdata = 255;
			}
		}
	}

	return result;
}

使用这段代码需要#include <opencv.hpp>, 其中涉及了一段异常抛出，各位看官可以直接删除掉那段代码，或者直接用c++的exception。

需要注意几个坑：

第一：Sobel算子的计算结果并非CV_8UC1

第二：注意atan2的计算，不要把Gx和Gy搞反

第三，非极大值抑制的时候不要在G的原图上改，否则画出来的边缘和素描一样，并非单边缘

第四，openCV里面的双阈值为像素具体值，而本程序里面用的是百分比，测试中使用的值up为0.93， down为0.9，效果还可以。

此外，非极大值抑制部分有另一个版本的形式，下面为另一个版本，对比来看，前一个版本应该更好用一点

//Step 4: Non-maximum suppression
	for (int i = 1; i < row - 1; ++i) {
		for (int j = 1; j < col - 1; ++j) {
			ushort* pdata = &G.at<ushort>(i, j);
			ushort* pnew = &non_max.at<ushort>(i, j);
			float th = theta.at<float>(i, j);
			if (th < 0)th += CV_PI;
			float angle = th * A2R;
			double tanth = tan(th);
			double inv_tanth = 1. / tanth;

			if (th <= 45.0) {
				double data1 = *(pdata - col + 1) * tanth + (1.0 - tanth) * (*(pdata + 1));
				double data2 = *(pdata + col - 1) * tanth + (1.0 - tanth) * (*(pdata - 1));
				if (*pdata < data1 || *pdata < data2)*pnew = 0;
			}
			else if (th <= 90.0) {
				double data1 = *(pdata - col + 1) * inv_tanth + (1.0 - inv_tanth) * (*(pdata - col));
				double data2 = *(pdata + col - 1) * inv_tanth + (1.0 - inv_tanth) * (*(pdata + col));
				if (*pdata < data1 || *pdata < data2)*pnew = 0;
			}
			else if (th <= 135.0) {
				double data1 = *(pdata - col - 1) * (-inv_tanth) + (1.0 + inv_tanth) * (*(pdata - col));
				double data2 = *(pdata + col + 1) * (-inv_tanth) + (1.0 + inv_tanth) * (*(pdata + col));
				if (*pdata < data1 || *pdata < data2)*pnew = 0;
			}
			else {
				double data1 = *(pdata - col - 1) * (-tanth) + (1.0 + tanth) * (*(pdata - 1));
				double data2 = *(pdata + col + 1) * (-tanth) + (1.0 + tanth) * (*(pdata + 1));
				if (*pdata < data1 || *pdata < data2)*pnew = 0;
			}
			
			if (max_pixel_value < *pdata)max_pixel_value = *pdata;
		}
	}