【openCV】Canny边缘检测

手动实现了Canny边缘检测算法的部分步骤。

#include "opencv2/opencv.hpp"
#include <math.h>
using namespace cv;

const int dir[8][2] = {{-1,-1},{-1,0},{-1,1},{0,-1},{0,1},{1,-1},{1,0},{1,1}};

inline int sqr(int x){
	return x*x;
}

void grad_low_thresholding(Mat &nms,const Mat &image,float threshold_l,int x,int y){
//在每个点周围探索，满足条件的点设为真边界，并以那个点为中心继续探索
	for (int k = 0;k < 8;k++){
		int newx = x + dir[k][0];
		int newy = y + dir[k][1];
		if (x < 0 || y < 0 || x >= nms.cols || y>=nms.rows) continue;
		if (nms.at<uchar>(newx,newy) == 128 && image.at<float>(newx,newy) >= threshold_l){
			nms.at<uchar>(newx,newy) = 255;
			grad_low_thresholding(nms,image,threshold_l,newx,newy);
		}
	}
}

int main(){
	
	// 从文件中加载原图
	Mat srcImage = imread("image.jpg");
	if(srcImage.empty())
	{//如果读入图像失败
         fprintf(stderr, "Can not load image\n");
         return 0;
    }
	
	//步骤1：高斯滤波
	Mat gauImage;
	GaussianBlur(srcImage,gauImage,Size(5,5),0,0);
	srcImage.release();
	
	//步骤2：转为灰度图 
	Mat grayImage;
	cvtColor(gauImage, grayImage, CV_BGR2GRAY);
	//cvShowImage("Original_Gray_Level_Image",grayImage);
	//cvWaitKey(0);
	//cvDestroyWindow("Original_Gray_Level_Image");
	gauImage.release();

	
	//步骤3：利用Sobel算子求x,y方向上偏导（grad_x,gray_y），再求梯度(grad)、梯度方向(theta)
	Mat grad_x,grad_y;//,abs_grad_x,abs_grad_y;
	Sobel(grayImage,grad_x,CV_32FC1,1,0);
	Sobel(grayImage,grad_y,CV_32FC1,0,1);
	//convertScaleAbs(grad_x,abs_grad_x);
	//convertScaleAbs(grad_y,abs_grad_y);

	//求梯度，角度
	//Mat grad(gauImage.rows,gauImage.cols,CV_16S),theta(gauImage.rows,gauImage.cols,CV_16S);
	Mat grad,theta;
	magnitude(grad_x,grad_y,grad);
	phase(grad_x,grad_y,theta,true);

	
	//步骤4：非极大值抑制
	//对每个点，跟其梯度方向上相邻的两个点进行比较，如果不是极大值则设为0，是则设为128。
	//处理结果存放在nms中（Non-maximum suppression）
	//Deg为梯度方向的角度值，根据角度值判断四种情况（为注释中的情况），
	//分别用g1和g2、g3和g4加上角度估算出该点在梯度方向上相邻点的梯度tmp1、tmp2，并进行比较。

	Mat nms(grad.rows,grad.cols,CV_8UC1);
	for (int j = 0;j < grad.cols;j++){
		nms.at<uchar>(0,j) = 0;
	}

	for (int i = 1; i < grad.rows-1;i++){
		nms.at<uchar>(i,0) = 0;
		for (int j = 1;j < grad.cols-1;j++){
			float deg = theta.at<float>(i,j);
			float g1,g2,g3,g4,pt;
			float weight,tmp1,tmp2;
			pt = grad.at<float>(i,j);

			//while (deg>360){deg-=360;}
			//while (deg<0){deg+=360;}
			if ( (deg >= 0 && deg < 45) || (deg >=180 && deg < 225) ){
				//***g1**tmp1**g2*************
				//****************************
				//*************pt*************
				//****************************
				//*************g3**tmp2**g4***
				g1 = grad.at<float>(i-1,j-1);
				g2 = grad.at<float>(i-1,j);
				g3 = grad.at<float>(i+1,j);
				g4 = grad.at<float>(i+1,j+1);
				weight = abs(grad_y.at<float>(i,j) / grad_x.at<float>(i,j));
				tmp1 = g1*weight + g2*(1-weight);
				tmp2 = g4*weight + g3*(1-weight);
			}
			else if ( (deg >= 45 && deg < 90) || (deg >=225 && deg < 270) ){
				//***g1***********************
				//**tmp1**********************
				//***g2********pt********g3***
				//**********************tmp2**
				//***********************g4***
				g1 = grad.at<float>(i-1,j-1);
				g2 = grad.at<float>(i,j-1);
				g3 = grad.at<float>(i,j+1);
				g4 = grad.at<float>(i+1,j+1);
				weight = abs(grad_x.at<float>(i,j) / grad_y.at<float>(i,j));
				tmp1 = g1*weight + g2*(1-weight);
				tmp2 = g4*weight + g3*(1-weight);
			}
			else if ( (deg >= 90 && deg < 135) || (deg >=270 && deg < 315) ){
				//***********************g3***
				//**********************tmp2**
				//****g1*******pt********g4***
				//****tmp1********************
				//****g2**********************
				g1 = grad.at<float>(i,j-1);
				g2 = grad.at<float>(i+1,j-1);
				g3 = grad.at<float>(i-1,j+1);
				g4 = grad.at<float>(i,j+1);
				weight = abs(grad_x.at<float>(i,j) / grad_y.at<float>(i,j));
				tmp1 = g2*weight + g1*(1-weight);
				tmp2 = g3*weight + g4*(1-weight);
			}
			else if ( (deg >= 135 && deg < 180) || (deg >=315 && deg < 360 ) ){
				//*************g1**tmp1**g2***
				//****************************
				//*************pt*************
				//****************************
				//***g3**tmp2**g4*************
				g1 = grad.at<float>(i-1,j);
				g2 = grad.at<float>(i-1,j+1);
				g3 = grad.at<float>(i+1,j-1);
				g4 = grad.at<float>(i+1,j);
				weight = abs(grad_y.at<float>(i,j) / grad_x.at<float>(i,j));
				tmp1 = g2*weight + g1*(1-weight);
				tmp2 = g3*weight + g4*(1-weight);
			}
			
			nms.at<uchar>(i,j) = ((pt >= tmp1 && pt >= tmp2)?128:0);
		}
		nms.at<uchar>(i,grad.cols-1) = 0;
	}

	for (int j = 0;j < grad.cols;j++){
		nms.at<uchar>(grad.rows-1,j) = 0;
	}
	
	imwrite("nms.jpg",nms);

	//步骤5：双阀值检测
	float threshold_h = 65;
	float threshold_l = 15;

	//首先判断高阀值，如果灰度大于高阀值，即为真边界，再在其周围根据低阀值搜索
	for (int i = 0; i < grad.rows;i++){
		for (int j = 0;j < grad.cols;j++){
			if (nms.at<uchar>(i,j) == 128 && grad.at<float>(i,j) >= threshold_h){
				nms.at<uchar>(i,j) = 255;
				grad_low_thresholding(nms,grad,threshold_l,i,j);
				//在每个点周围探索，满足条件的点设为真边界，并以那个点为中心继续探索
			}
		}
	}

	for (int i = 0; i < grad.rows;i++){
		for (int j = 0;j < grad.cols;j++){
			if (nms.at<uchar>(i,j) != 255){
				nms.at<uchar>(i,j) = 0;
			}
		}
	}

	//cvCanny(&grayImage,&nms,150,30);
	imwrite("gray.jpg",grayImage);
	imwrite("output.jpg",nms);

	grayImage.release();
	grad_x.release();
	grad_y.release();
	//abs_grad_x.release();
	//abs_grad_y.release();
	grad.release();
	theta.release();
	nms.release();
	return 0;
}

另外，openCV中有库函数cvCanny()实现了Canny边缘检测算法。