原理：

实现：

//阶乘
int factorial(int n)
{
	int fac = 1;
	if (n == 0)	return fac;
	for (int i = 1; i <= n; ++i)	fac *= i;
	return fac;
}

//获得Sobel平滑算子
Mat getSobelSmooth(int size)
{
	int n = size - 1;
	Mat SobelSmoothoper = Mat::zeros(size, 1, CV_32F);
	for (int k = 0; k <= n; k++)
	{
		float *pt = SobelSmoothoper.ptr<float>(0);
		pt[k] = factorial(n) / (factorial(k)*factorial(n - k));
	}
	return SobelSmoothoper;
}

//获得Sobel差分算子
Mat getSobeldiff(int size)
{
	Mat Sobeldiffoper = Mat::zeros(Size(size, 1), CV_32F);
	Mat SobelSmooth = getSobelSmooth(size - 1);
	for (int k = 0; k < size; k++)
	{
		if (k == 0)
			Sobeldiffoper.at<float>(0, k) = 1;
		else if (k == size - 1)
			Sobeldiffoper.at<float>(0, k) = -1;
		else
			Sobeldiffoper.at<float>(0, k) = SobelSmooth.at<float>(0, k) - SobelSmooth.at<float>(0, k - 1);
	}
	return Sobeldiffoper;
}

//卷积实现
void conv2D(Mat& src, Mat& dst, Mat kernel)
{
	flip(kernel, kernel, -1);
	filter2D(src, dst, CV_32F, kernel);
}

//可分离卷积———先水平方向卷积，后垂直方向卷积
void sepConv2D_X_Y(Mat& src, Mat& dst, Mat kernel_X, Mat kernel_Y)
{
	Mat dst_kernel_X;
	conv2D(src, dst_kernel_X, kernel_X); //水平方向卷积
	conv2D(dst_kernel_X, dst, kernel_Y); //垂直方向卷积
}

//可分离卷积———先垂直方向卷积，后水平方向卷积
void sepConv2D_Y_X(Mat& src, Mat& dst, Mat kernel_Y, Mat kernel_X)
{
	Mat dst_kernel_Y;
	conv2D(src, dst_kernel_Y, kernel_Y); //垂直方向卷积
	conv2D(dst_kernel_Y, dst, kernel_X); //水平方向卷积
}

//Sobel算子边缘检测
void sobel(Mat& src, Mat& dst, Mat& dst_X, Mat& dst_Y, int size)
{
	Mat SobelSmoothoper = getSobelSmooth(size); //平滑系数
	Mat Sobeldiffoper = getSobeldiff(size);		//差分系数

	sepConv2D_X_Y(src, dst_Y, SobelSmoothoper, Sobeldiffoper.t()); //得到水平边缘
	sepConv2D_Y_X(src, dst_X, SobelSmoothoper.t(), Sobeldiffoper); //得到垂直边缘	

	//边缘强度（近似）
	dst = abs(dst_X) + abs(dst_Y);
	convertScaleAbs(dst, dst);
	//convertScaleAbs(dst_X, dst_X);
	//convertScaleAbs(dst_Y, dst_Y);
}

// 确定一个点的坐标是否在图像内
bool checkInRange(int r, int c, int rows, int cols) 
{
	if (r >= 0 && r < rows && c >= 0 && c < cols)
		return true;
	else
		return false;
}

//从确定边缘点出发，延长边缘
void trace(Mat &edgeMag_noMaxsup, Mat &edge, float Th, int r, int c, int rows, int cols)
{
	if (edge.at<uchar>(r, c) == 0)
	{
		for (int i = -1; i <= 1; ++i)
		{
			for (int j = -1; j <= 1; ++j)
			{
				if (checkInRange(r + i, c + j, rows, cols) && edgeMag_noMaxsup.at<float>(r + i, c + j) > Th)
					edge.at<uchar>(r, c) = 255;
			}
		}
	}
}

//Canny边缘检测
void canny(Mat &src, Mat &dst, float Tl, float Th, int ksize = 3, bool L2graydient = false) 
{
	//高斯滤波
	GaussianBlur(src, src, Size(3, 3), 0);
	//sobel算子
	Mat dx, dy, sobel_dst;
	sobel(src, sobel_dst, dx, dy, ksize);

	//计算梯度幅值
	Mat edgeMag;
	if (L2graydient)
		magnitude(dx, dy, edgeMag);  //开平方
	else
		edgeMag = abs(dx) + abs(dy); //绝对值之和近似

	//计算梯度方向以及非极大值抑制
	Mat edgeMag_noMaxsup = Mat::zeros(src.size(), CV_32F);
	for (int i = 1; i < src.rows - 1; ++i)
	{
		for (int j = 1; j < src.cols - 1; ++j) 
		{
			float angle =  atan2f(dy.at<float>(i, j), dx.at<float>(i, j)) / CV_PI * 180; //当前位置梯度方向
			float cur = edgeMag.at<float>(i, j);  //当前位置梯度幅值
			
			//非极大值抑制 垂直边缘--梯度方向为水平方向--3*3邻域内左右方向比较
			if (abs(angle) < 22.5 || abs(angle) > 157.5)
			{
				float left = edgeMag.at<float>(i, j - 1);
				float right = edgeMag.at<float>(i, j + 1);
				if (cur >= left && cur >= right)
					edgeMag_noMaxsup.at<float>(i, j) = cur;
			}

			//水平边缘--梯度方向为垂直方向--3*3邻域内上下方向比较
			if ((angle >= 67.5 && angle <= 112.5) || (angle >= -112.5 && angle <= -67.5)) 
			{
				float top = edgeMag.at<float>(i - 1, j);
				float down = edgeMag.at<float>(i + 1, j);
				if (cur >= top && cur >= down)
					edgeMag_noMaxsup.at<float>(i, j) = cur;
			}

			//+45°边缘--梯度方向为其正交方向--3*3邻域内右上左下方向比较
			if ((angle>112.5 && angle <= 157.5) || (angle>-67.5 && angle <= -22.5)) 
			{
				float right_top = edgeMag.at<float>(i - 1, j + 1);
				float left_down = edgeMag.at<float>(i + 1, j - 1);
				if (cur >= right_top && cur >= left_down)
					edgeMag_noMaxsup.at<float>(i, j) = cur;
			}

			//+135°边缘--梯度方向为其正交方向--3*3邻域内右下左上方向比较
			if ((angle >= 22.5 && angle < 67.5) || (angle >= -157.5 && angle < -112.5)) 
			{
				float left_top = edgeMag.at<float>(i - 1, j - 1);
				float right_down = edgeMag.at<float>(i + 1, j + 1);
				if (cur >= left_top && cur >= right_down)
					edgeMag_noMaxsup.at<float>(i, j) = cur;
			}

		}
	}

	//双阈值处理及边缘连接
	dst = Mat::zeros(src.size(), CV_8U);
	for (int i = 1; i < src.rows - 1; ++i) 
	{
		for (int j = 1; j < src.cols - 1; ++j) 
		{
			float mag = edgeMag_noMaxsup.at<float>(i, j);
			//大于高阈值，为确定边缘点
			if (mag > Th)
				dst.at<uchar>(i, j) = 255;
			else if (mag < Tl)
				dst.at<uchar>(i, j) = 0;
			else
				trace(edgeMag_noMaxsup, dst, Th, i, j, src.rows, src.cols);
		}
	}

}