边缘检测：canny

canny实现边缘检测
主要流程见主函数
主要参考：https://blog.csdn.net/dcrmg/article/details/52344902
https://www.cnblogs.com/love6tao/p/5152020.html
不过有些细节有问题，我改过来了

#include <iostream>
#include<opencv2/opencv.hpp>
#include <string>

using namespace std;
using namespace cv;

#define PI 3.141592653589793


void conv(char *G, Mat pic, Mat res)
{
	int i, j, i1, j1;
	int tmp;

	for (i = 1; i < pic.rows-1; i++)
	{
		for (j = 1; j < pic.cols - 1; j++)
		{
			tmp = 0;
			for (i1 = -1; i1 <= 1; i1++)
			{
				for (j1 = -1; j1 <= 1; j1++)
				{
					tmp += pic.at<uchar>(i + i1, j + j1)*G[(j1 + 1) * 3 + i1 + 1];
				}
			}
			res.at<short>(i, j) = tmp;
		}
	}
}

void merge1(Mat resX, Mat resY, Mat &res, Mat &theta)
{
	int tmp;
	double theTmp;
	for (int i = 0; i < res.rows; i++)
	{
		for (int j = 0; j < res.cols; j++)
		{
			tmp = abs(resX.at<short>(i, j)) + abs(resY.at<short>(i, j));
			//tmp = sqrt(resX.at<short>(i, j)*resX.at<short>(i, j) + resY.at<short>(i, j)*resY.at<short>(i, j));
			if (tmp > 255)
				tmp = 255;
			res.at<uchar>(i,j) = tmp;
			theTmp = atan2(resY.at<short>(i, j), resX.at<short>(i, j))*180/PI;
			if (theTmp < 0)
				theTmp += 360;
			theta.at<float>(i, j) = theTmp;
		}
	}
}

void mySobel(Mat pic, Mat &src, Mat &theta)
{
	char Gx[3 * 3] = {-1, 0, 1,
	                  -2, 0, 2,
	                  -1, 0, 1};
	char Gy[3 * 3] = { 1,  2,  1,
		               0,  0,  0,
		              -1, -2, -1};
	Mat resX = Mat::zeros(pic.rows, pic.cols, CV_16SC1);
	Mat resY = Mat::zeros(pic.rows, pic.cols, CV_16SC1);
	conv(Gx, pic, resX);
	conv(Gy, pic, resY);
	merge1(resX, resY, src, theta);


}

void BGR2Gray(Mat src, Mat &grayPic)
{

	for (int i = 0; i < src.rows; i++)
	{
		for (int j = 0; j < src.cols; j++)
		{
			//Gray = R*0.299 + G*0.587 + B*0.114
		
			grayPic.at<uchar>(i, j) = src.at<Vec3b>(i, j)[0] * 0.114 + src.at<Vec3b>(i, j)[1] * 0.587 + src.at<Vec3b>(i, j)[2] * 0.299;
		}

	}
}

void getGaussKer(double **K, int size, double sigma)
{
	int center = size / 2;
	double sum = 0;

	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			K[i][j] = exp(((i - center)*(i - center) + (j - center)*(j - center)) / (-2 * sigma*sigma));///最后还是需要归一化加不加系数无所谓
			sum += (K[i][j]);
		}
	}

	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			K[i][j] /= sum;
		}
	}

}

void GaussFilter(Mat pic, Mat &src, double **K, int size)
{
	int center = size / 2;
	for (int i = 0; i < pic.rows; i++)
	{
		for (int j = 0; j < pic.cols; j++)
		{
			double tmp = 0;
			int nowRow;
			int nowCol;

			for (int i1 = -size / 2; i1 <= size / 2; i1++)
			{
				for (int j1 = -size / 2; j1 <= size / 2; j1++)
				{
					nowRow = i + i1;
					if (nowRow < 0|| nowRow>=pic.rows)
						continue;
					nowCol = j + j1;
					if (nowCol < 0 || nowCol >= pic.cols)
						continue;
					tmp += K[center + i1][center + j1] * pic.at<uchar>(nowRow, nowCol);
				}
			}
			if (tmp < 0)
				tmp = 0;
			if (tmp > 255)
				tmp = 255;
			src.at<uchar>(i, j) = tmp;

		}
	}
}

void locMaxValue(Mat amp, Mat theta, Mat &src)
{
	Mat output = amp.clone();深拷贝

	for (int i = 1; i < amp.rows - 1; i++)
	{
		for (int j = 1; j < amp.cols - 1; j++)///判断每个像素点，最外一圈8领域不完整就不管了
		{
			float thetaNow = theta.at<float>(i, j);
			float k=abs(tan(thetaNow));
			float g1, g2, g3, g4;

			float dTmp1, dTmp2;
			float nowGray = amp.at<uchar>(i, j);

			if ((thetaNow >= 90 && thetaNow < 135) || (thetaNow >= 270 && thetaNow < 315))
			{
				g1 = output.at<uchar>(i - 1, j-1);
				g2 = output.at<uchar>(i - 1, j);
				g3 = output.at<uchar>(i + 1, j);
				g4 = output.at<uchar>(i + 1, j + 1);


				dTmp1 = g1 / k + (1 - 1 / k)*g2;
				dTmp2 = g4 / k + (1 - 1 / k)*g3;

			}
			else if ((thetaNow >= 135 && thetaNow < 180) || (thetaNow >= 315 && thetaNow < 360))
			{
				g1 = output.at<uchar>(i - 1, j - 1);
				g2 = output.at<uchar>(i, j - 1);
				g3 = output.at<uchar>(i, j + 1);
				g4 = output.at<uchar>(i + 1, j + 1);

				dTmp1 = g1*k + (1 - k)*g2;
				dTmp2 = g4*k + (1 - k)*g3;
				
			}
			else if ((thetaNow >= 45 && thetaNow < 90) || (thetaNow >= 225 && thetaNow < 270))
			{
				g1 = output.at<uchar>(i - 1, j );
				g2 = output.at<uchar>(i-1, j + 1);
				g3 = output.at<uchar>(i-1, j );
				g4 = output.at<uchar>(i + 1, j + 1);

				dTmp1 = g2 / k + (1 - 1 / k)*g1;
				dTmp2 = g3 / k + (1 - 1 / k)*g4;
			}
			else if (thetaNow >= 0 && thetaNow < 45 || (thetaNow >= 180 && thetaNow < 225))
			{
				g1 = output.at<uchar>(i - 1, j+1);
				g2 = output.at<uchar>(i , j + 1);
				g3 = output.at<uchar>(i , j-1);
				g4 = output.at<uchar>(i + 1, j - 1);

				dTmp1 = g1*k + (1 - k)*g2;
				dTmp2 = g4*k + (1 - k)*g3;
			}

			if (nowGray < dTmp1 || nowGray < dTmp2)
			{
				output.at<uchar>(i, j) = 0;
			}

		}
	}
	src = output.clone();

}


void traceLink(Mat &src, int minValue, int i, int j)
{
	int xn[8] = { 1,1,0,-1,-1,-1,0,1 };
	int yn[8] = { 0, 1, 1, 1, 0, -1, -1, -1 };

	int xx = 0;
	int yy = 0;
	bool change = true;
	while (change)
	{
		change = false;

		for (int k = 0; k < 8; k++)
		{
			xx = i + xn[k];
			yy = j + yn[k];

			if (src.at<uchar>(xx,yy)!=255&&src.at<uchar>(xx, yy) > minValue)
			{
				change = true;
				src.at<uchar>(xx, yy) = 255;
				i = xx;
				j = yy;
				break;
			}
		}
	}

}

void doubleThreshold(Mat &src, int maxValue, int minValue)
{
	for (int i = 1; i < src.rows-1; i++)
	{
		for (int j = 1; j < src.cols-1; j++)
		{
			int value = src.at<uchar>(i, j);
			//int t=value;
			if (src.at<uchar>(i, j) >= maxValue)
			{
				src.at<uchar>(i, j) = 255;

				traceLink(src, minValue, i, j);
			}
			
		}
	}

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



int main()
{
	int size = 3;
	double sigma = 1;

	double **k = new double*[size];
	for (int i = 0; i < size; i++)
	{
		k[i] = new double[size];
	}
	Mat pic = imread("C:\\Users\\Thinkpad\\Desktop\\c++work\\pic\\lena.jpg", 1);
	Mat grayPic=Mat::zeros(pic.rows, pic.cols, CV_8UC1); 
	Mat gaussPic=Mat::zeros(pic.rows, pic.cols, CV_8UC1);
	Mat amplitude = Mat::zeros(pic.rows, pic.cols, CV_8UC1);
	Mat theta = Mat::zeros(pic.rows, pic.cols, CV_32FC1);
	Mat src;

	BGR2Gray(pic, grayPic);转成灰度图

	getGaussKer(k, size, sigma);///计算高斯核

	GaussFilter(grayPic, gaussPic, k, 3);///计算高斯模糊后的图

	mySobel(gaussPic, amplitude, theta);///用sobel算子计算幅值和角度

	locMaxValue(amplitude, theta, src);非极大抑制，梯度方向不是极大值的进行抑制

	doubleThreshold(src, 100, 60);双阈值链接
	

	imshow("1", src);
	imshow("0", amplitude);
	waitKey();

	/*
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			cout << k[i][j] << ",";
		}
		cout << endl;
	}*/



	for (int i = 0; i < size; i++)
	{
		delete[]k[i];
	}
	delete[]k;
	
	return 0;
}

这里简单介绍下梯度算子的方向：通常OpenCV打开的的图像坐标原点在左上角，x轴方向向下，y轴方向向右，但是在选择算子进行梯度计算时比如sobel算子，其“X”方向的计算是右边减去左边，“Y”方向的计算是上边减去下边，这种梯度的计算方式和我们正常左手系坐标系的计算方式相同，也就是x轴方向向右，y轴方向向上，而dx,dy的计算与坐标系选取无关，但是由dx,dy计算得到的角度是与“X”轴正方向的夹角，所以无论dx,dy是基于哪个坐标系算的，sobel算子计算的梯度方向都是以右为x轴正方向的左手系坐标系的夹角，下图是基于此，我进行非极大抑制的计算过程草图