图像的边缘是指灰度值发生急剧变换的位置。在某种程度上，边缘不随光照和视角的变化而变化。
边缘检测的目的是制作一个线图，在不会损害图像内容的情况下，同时又大大减少图像的数据量,提供了对图像数据的合理概述。
边缘通过检查每个像素的领域并对其灰度变化进行量化的，相当于微积分里连续函数中方向导数或者离散数列的查分。
8.1、Roberts算子
注意，为了方便，这里把锚点的位置标注在了卷积核上，这是不合适的，应该标注在卷积核逆时针翻转180°的结果上。这里标注的地方只是一个位置说明，即Roberts135翻转180°后，锚点的位置在第0行第0列，Roberts45翻转180°后，锚点的位置 在第0行第1列。
与Roberts核卷积，本质上是两个对角方向上的差分。
反映的是在垂直方向和水平方向上的边缘。
通常有四种方式来衡量最后多次卷积后输出的边缘强度：
其中取绝对值的最大值的方式，对边缘的走向有些敏感，而其他几种方式可以获得性能更一致的全方位响应。当然，取平方和的开方的方式效果一般是最好的，但是同时会更加耗时。
void roberts(InputArray src, OutputArray dst, int ddepth, int x, int y, int borderType)
{
if (x == 0 && y == 0)
{
return;
}
//roberts135=|1 0| roberts45=|0 1|
// |0 -1| |-1 0|
Mat roberts_1 = ( Mat_<float>(2, 2) << 1, 0, 0, -1 );
Mat roberts_2 = ( Mat_<float>(2, 2) << 0, 1, -1, 0 );
//当x不等于0时，src与robert_1卷积,135°卷积核,检测的是45°的边缘
if (x != 0 && y == 0)
{
conv2D(src, roberts_1, dst, ddepth, Point(0, 0), borderType);
}
//当y不等于0时，src与robert_2卷积，45°卷积核，检测的是135°的边缘
if (x == 0 && y != 0)
{
conv2D(src, roberts_2, dst, ddepth, Point(0, 0), borderType);
}
}
Mat img = imread("Koala.jpg", IMREAD_GRAYSCALE);
//图像矩阵和robert算子卷积
Mat roberts_1;
roberts(img, roberts_1, CV_32FC1, 1, 0);
Mat roberts_2;
roberts(img, roberts_2, CV_32FC1, 0, 1);
//两个卷积结果的灰度级显示
Mat abs_roberts_1, abs_roberts_2;
convertScaleAbs(roberts_1, abs_roberts_1, 1, 0);
convertScaleAbs(roberts_2, abs_roberts_2, 1, 0);
//边缘强度,这里采用平方根的方式
Mat roberts_1_2, roberts_2_2;
pow(roberts_1, 2.0, roberts_1_2);
pow(roberts_2, 2.0, roberts_2_2);
Mat edge;
sqrt(roberts_1_2 + roberts_2_2, edge);
edge.convertTo(edge, CV_8UC1);
Roberts边缘检测因为使用了很少的邻域像素来近似边缘强度， 因此对图像中的噪声具有高度敏感性。可以先对图像进行平滑处理， 然后再进行Roberts边缘检测。 下面介绍几种具有平滑作用 的边缘提取卷积核。
8.2、Prewitt算子
图像与prewittx卷积后可以反映图像垂直方向上的边缘，与prewitty卷积后可以反映图像水平方向上的边缘。 而且，这两个卷积核均是可分离的，其中
从分离的结果可以看出，prewittx算子实际上先对图像进行垂直方向上的非归一化的均值平滑，然后进行水平方向上的差分；而prewitty算子实际上先对图像进行水平方向上的非归一化的均值平滑，然后进行垂直方向上的差分。
反映的是在45°和135°方向上的边缘。遗憾的是，这两个卷积核是不可分离的。
void prewitt(InputArray src, OutputArray dst, int ddepth, int x = 1, int y = 0, int borderType )
{
if (x == 0 && y == 0)
{
return;
}
//prewitt_x=|1 0 -1| 分离 prewitt_x= | 1 | *| 1 0 -1|
// |1 0 -1| | 1 |
// |1 0 -1| | 1 |
//prewitt_y=| 1 1 1| 分离 prewitt_y=| 1 1 1|*| 1 |
// | 0 0 0| | 0 |
// |-1 -1 -1| | -1 |
Mat prewitt_x_y = (Mat_<float>(3,1) << 1, 1, 1);
Mat prewitt_x_x = (Mat_<float>(1,3) << 1, 0, -1);
Mat prewitt_y_y = prewitt_x_y.t();
Mat prewitt_y_x = prewitt_x_x.t();
//x 不为零的话，执行x方向的prewitt卷积
if (x != 0 && y == 0)
{
sepConv2D_Y_X(src, dst, ddepth, prewitt_x_y, prewitt_x_x, Point(-1, -1), borderType);
}
// y 不为零的话，执行y方向的prewitt卷积
if (x == 0 && y != 0)
{
sepConv2D_Y_X(src, dst, ddepth, prewitt_y_y, prewitt_y_x, Point(-1, -1), borderType);
}
}
Mat img = imread("Koala.jpg", IMREAD_GRAYSCALE);
//图像矩阵和prewitt算子卷积
Mat prewitt_1;
prewitt(img, prewitt_1, CV_32FC1, 1, 0);
Mat prewitt_2;
prewitt(img, prewitt_2, CV_32FC1, 0, 1);
//两个卷积结果的灰度级显示
Mat abs_prewitt_1, abs_prewitt_2;
convertScaleAbs(prewitt_1, abs_prewitt_1, 1, 0);
convertScaleAbs(prewitt_2, abs_prewitt_2, 1, 0);
//边缘强度,这里采用平方根的方式
Mat prewitt_1_2, prewitt_2_2;
pow(prewitt_1, 2.0, prewitt_1_2);
pow(prewitt_2, 2.0, prewitt_2_2);
Mat edge;
sqrt(prewitt_1_2 + prewitt_2_2, edge);
edge.convertTo(edge, CV_8UC1);
从Roberts和Prewitt边缘检测的效果图可以清晰地理解差分方向（或称梯度方向）与得到的边缘方向是垂直的，如水平差分方向上的卷积反映的是垂直方向上的边缘。
在图像的平滑处理中，高斯平滑的效果往往比均值平滑要好，因此把Prewitt算子的非归一化的均值卷积核替换成非归一化的高斯卷积核，就可以构建3阶的Sobel边缘检测算子。
8.3、Sobel算子
sobelx=| 1 |*|1 0 -1| = | 1 0 -1|
| 2 | | 2 0 -2|
| 1 | | 1 0 -1|
sobely=|1 2 1|*| 1 | = | 1 2 1|
| 0 | | 0 0 0|
|-1 | | -1 -2 -1|
Sobel算子是在一个坐标轴方向上进行非归一化的高斯平滑，在另一个坐标轴方向上 进行差分处理。
n×n的Sobel算子是由平滑算子和差分算子full卷积而得到的，其中n为奇 数。对于窗口大小为n的非归一化的高斯平滑算子等于n-1阶的二项式展开式的系数，
窗口大小为n的差分算子是在n-2阶的二 项式展开式的系数两侧补零， 然后后向差分得到的。
//阶乘
int factorial(int n)
{
int fac = 1;
if (n == 0)
{
return fac;
}
for (int i = 1; i <= n; i++)
{
fac *= i;
}
return fac;
}
//获取平滑算子
Mat getPascalSmooth(int n)
{
Mat pascalSmooth = Mat::zeros(Size(n, 1), CV_32FC1);
for (int i = 0; i < n; i++)
{
pascalSmooth.at<float>(0, i) = factorial(n - 1) / (factorial(i)*factorial(n - 1 - i));
}
return pascalSmooth;
}
//获取差分算子
Mat getPascalDiff(int n)
{
Mat pascalDiff = Mat::zeros(Size(n, 1), CV_32FC1);
Mat pascalSmooth_privious = getPascalSmooth(n - 1);
for (int i = 0; i < n; i++)
{
if (i == 0)
{
pascalDiff.at<float>(0, i) = 1;
}
else if (i == n - 1)
{
pascalDiff.at<float>(0, i) = -1;
}
else
{
pascalDiff.at<float>(0, i) = pascalSmooth_privious.at<float>(0, i) - pascalSmooth_privious.at<float>(0, i - 1);
}
}
return pascalDiff;
}
//sobel算子卷积运算
Mat sobel(Mat image, int x_flag, int y_flag, int winSize, int borderType)
{
//卷积核窗口大小应该为大于3的奇数
if (winSize < 3 || winSize % 2 == 1)
{
printf("error winSize!");
}
//平滑系数
Mat pascalSmooth = getPascalSmooth(winSize);
//差分系数
Mat pascalDiff = getPascalDiff(winSize);
//sobel卷积
Mat image_con_sobel;
if (x_flag != 0)
{
sepConv2D_Y_X(image, image_con_sobel, CV_32FC1, pascalSmooth.t(), pascalDiff, Point(-1, -1), borderType);
}
if(x_flag == 0 && y_flag != 0)
{
sepConv2D_Y_X(image, image_con_sobel, CV_32FC1, pascalSmooth, pascalDiff.t(), Point(-1, -1), borderType);
}
return image_con_sobel;
}
//main.cpp
Mat img = imread("Koala.jpg", IMREAD_GRAYSCALE);
//图像矩阵和sobel算子卷积
Mat sobel_x;
sobel_x = sobel(img, 1, 0, 5, BORDER_DEFAULT);
Mat sobel_y;
sobel_y = sobel(img, 0, 1, 5, BORDER_DEFAULT);
//两个卷积结果的灰度级显示
Mat abs_sobel_x, abs_sobel_y;
convertScaleAbs(sobel_x, abs_sobel_x, 1, 0);
convertScaleAbs(sobel_y, abs_sobel_y, 1, 0);
//边缘强度,这里采用平方根的方式
Mat sobel_x_1, sobel_y_2;
pow(sobel_x, 2.0, sobel_x_1);
pow(sobel_y, 2.0, sobel_y_2);
Mat edge;
sqrt(sobel_x_1 + sobel_y_2, edge);
edge.convertTo(edge, CV_8UC1);
一个有趣的应用， 就是对边缘检测的结果进行反色处理会呈现出铅笔素描的效果。
OpenCV提供的接口：
void Sobel( InputArray src, OutputArray dst, int ddepth,
int dx, int dy, int ksize = 3,
double scale = 1, double delta = 0,
int borderType = BORDER_DEFAULT );
8.4、Scharr算子
这两个卷积核均是不可分离的。
scharrx反映垂直方向的边缘强度，scharry反映水平方向的边缘强度。
scharr45反映135°方向的边缘强度，scharr135反映45°方向的边缘强度。
void scharr(InputArray src, OutputArray dst, int ddepth, int x, int y, int borderType)
{
if (x == 0 && y == 0)
{
printf("error x or y input");
}
Mat scharr_x = (Mat_<float>(3, 3) << 3, 0, -3, 10, 0, -10, 3, 0, -3);
Mat scharr_y = (Mat_<float>(3, 3) << 3, 10, 3, 0, 0, 0, -3, -10, -3);
if (x != 0 && y == 0)
{
conv2D(src, scharr_x, dst, ddepth, Point(-1, -1), borderType);
}
if (x == 0 && y != 0)
{
conv2D(src, scharr_y, dst, ddepth, Point(-1, -1), borderType);
}
}
OpenCV中提供的接口
void Scharr( InputArray src, OutputArray dst, int ddepth,
int dx, int dy, double scale = 1, double delta = 0,
int borderType = BORDER_DEFAULT );
8.5、Kirsch算子和Bobinson算子
8.5.1、Krisch算子
Kirsch算子由以下8个卷积核组成：
图像与每个卷积核进行卷积，然后取绝对值作为对应方向上的边缘强度的量化。对8个卷积结果取绝对值，然后在对应值位置取最大值作为最后输出的边缘强度。
8.5.2、Robinson算子
与Kirsch类似
因为Kirsch算子和Robinson算子使用了8个方向上的卷积核， 所以其检测到的边缘比标准的Prewitt算子和Sobel算子检测到的边缘会显得更加丰富。
8.6、Canny算子
基于卷积运算的边缘检测算法， 比如Sobel、 Prewitt等， 有如下两个缺点：
（1） 没有充分利用边缘的梯度方向。
（2） 最后输出的边缘二值图，只是简单地利用阈值进行处理，显然如果阈值过大，则会损失很多边缘信息； 如果阈值过小，则会有很多噪声。
而Canny边缘检测基于这两点做了改进，提出了：
（1） 基于边缘梯度方向的非极大值抑制。
（2） 双阈值的滞后阈值处理。
算法步骤如下：
第一步：图像矩阵I分别与水平方向上的卷积核sobelx和垂直方向上的卷积核sobely卷积得到dx和dy，然后利用平方和的开方得到边缘强度。这一步的过程和Sobel边缘检测一样，这里也可以将卷积核换为Prewitt核。
第二步：利用第一步计算出的dx和dy，计算出梯度方向angle=arctan2（dy，dx），即对每一个位置（r， c） ，angle（r， c）=arctan 2（dy（r， c），dx（r， c））代表该 位置的梯度方向，一般用角度表示，即angle（r， c）∈[0，180]∪[-180，0]。
第三步： 对每一个位置进行非极大值抑制的处理。magnitude（1， 1）分别与其右上方和左下方的值做比较， 这里912>292且912>276，如果它的值均大于梯度方向上邻域的值，则可看作极大值，令 nonMaxSup（1，1） =magni tude（1，1）；如果它的值不全大于梯度方向上邻域的值，则可看作非极大值，就需要抑制，令nonMaxSup（1，1） =0。
总结上述非极大值抑制的过程：
如果magnitude（r， c）在沿着梯度方向angle（r，c）上的邻域内是最大的则为极大值；否则，设置为0。 对于非极大值抑制的实现，将梯 度方向一般离散化为以下四种情况：
· angle(r,c)∈[0,22.5)∪(-22.5,0]∪(157.5,180]∪(-180,157.5)
· angle(r,c)∈[22.5,67.5)∪[-157.5,-112.5)
· angle(r,c)∈[67.5,112.5]∪[-112.5,-67.5]
· angle(r,c)∈(112.5,157.5]∪[-67.5.-22.5]
接下来介绍常用的非极大值抑制的第二种方式，它可以弥补这一点，没有舍弃任何信息， 而是用插值法拟合梯度方向上的边缘强度，这样会更加准确地衡量梯度方向上的边缘强度。
一般将梯度方向离散化为以下四种情况：
· angle(r,c)∈(45,90]∪(-135,-90]
· angle(r,c)∈(90,135]∪(-90,-45]
· angle(r,c)∈[0,45]∪[-180,-135]
· angle(r,c)∈(135,180]∪(-45,0)
如果angle（r， c）∈（45， 90]∪（-135， -90]，那么两个插值分别为：
其余不赘述
第四步： 双阈值的滞后阈值处理。
（1） 边缘强度大于高阈值的那些点作为确定边缘点。
（2） 边缘强度比低阈值小的那些点立即被剔除。
（3） 边缘强度在低阈值和高阈值之间的那些点， 按照以下原则进行处理——只有这 些点能按某一路径与确定边缘点相连时，才可以作为边缘点被接受。组成这一路径的所 有点的边缘强度都比低阈值要大。对这一过程可以理解为，首先选定边缘强度大于高阈 值的所有确定边缘点，然后在边缘强度大于低阈值的情况下尽可能延长边缘。
OpenCV提供的接口：
void Canny( InputArray image, OutputArray edges,
double threshold1, double threshold2,
int apertureSize = 3, bool L2gradient = false );
8.7、Laplacian算子
二维函数f （x， y）的Laplacian（拉普拉斯）变换， 由以下计算公式定义：
图像矩阵与拉普拉斯核的卷积本质上是计算任意位置的值与其在水平方向和垂直方向上四个相邻点平均值之间的差值（只是相差一个4的倍数）。
注意：没有平滑处理，对噪声敏感，误将噪声作为边缘，并且得不到有方向的边缘。
显然， 它无法像Sobel和Prewitt算子那样单独得到水平方向、垂直方向或者其他固定方向上的边缘。拉普拉斯算子的优点是它只有一个卷积核，所以其计算成本比其他算子要低。
其他形式如下
拉普拉斯核内所有值的和必须等于0，这样就使得在恒等灰度值区域不会产生错误的边缘，而且上述几种形式的拉普拉斯算子均是不可分离的。
void Laplacian( InputArray src, OutputArray dst, int ddepth,
int ksize = 1, double scale = 1, double delta = 0,
int borderType = BORDER_DEFAULT );
8.8、高斯拉普拉斯（LoG）算子
步骤：
第一步： 构建窗口大小为H×W、 标准差为σ的LoG卷积核。
其中H、 W 均为奇数且一般H=W， 卷积核锚点的位置在
第二步： 图像矩阵与LoGH×W 核卷积， 结果记为I_Cov_LoG。
第三步： 边缘二值化显示。
//构建分离的高斯卷积核,kernelX:水平方向；kernelY：垂直方向
void getSepLoGKernel(float sigma, int length, Mat &kernelX, Mat &kernelY)
{
//分配内存
kernelX.create(Size(length, 1), CV_32FC1);
kernelY.create(Size(1, length), CV_32FC1);
int center = (length - 1) / 2;
double sigma2 = pow(sigma, 2);
//构建可分离的高斯拉普拉斯核
for (int c = 0; c < length; c++)
{
float norm2 = pow(c - center, 2.0);
kernelY.at<float>(c, 0) = exp(-norm2 / (2 * sigma));
kernelX.at<float>(0, c) = (norm2 / sigma2 - 1.0)*kernelY.at<float>(c, 0);
}
}
//LoG算子
Mat LoG(InputArray image, float sigma, int win)
{
Mat kernelX, kernelY;
//得到两个分离核
getSepLoGKernel(sigma, win, kernelX, kernelY);
//先进行水平卷积，再进行垂直卷积
Mat convXY;
sepConv2D_Y_X(image, convXY, CV_32FC1, kernelX, kernelY, Point(-1, -1), BORDER_DEFAULT);
//卷积核转置
Mat kernelX_T = kernelX.t();
Mat kernelY_T = kernelY.t();
//先进行垂直卷积，再进行水平卷积
Mat convYX;
sepConv2D_Y_X(image, convYX, CV_32FC1, kernelX_T, kernelY_T, Point(-1, -1), BORDER_DEFAULT);
//计算两个卷积的结果和，得到高斯拉普拉斯卷积
Mat LoG;
add(convXY, convYX, LoG);
return LoG;
}
虽然高斯拉普拉斯核可分离，但是当核的尺寸较大时，计算量仍然很大，下面通过 高斯差分近似高斯拉普拉斯，从而进一步减少计算量。
8.9、高斯差分（DoG）算子
当 k=0.95时， 高斯拉普拉斯和高斯差分的值是近似相等的 。
高斯差分边缘检测的步骤如下：
第一步： 构建窗口大小为H×W 的高斯差分卷积核。
第二步： 图像矩阵与DoGH×W 核卷积， 结果记为I_Cov_DoG。
第三步： 与拉普拉斯边缘检测相同的二值化或者水墨效果的显示。
注意：为了减少计算量，可以不用创建高斯差分核，而是根据卷积的加法分配律和结合律的性质，图像矩阵分别与两个高斯核卷积，然后做差即可。
//高斯卷积
Mat gaussConv(Mat I, float sigma, int s)
{
//构建水平方向上的非归一化的高斯核
Mat xkernel = Mat::zeros(1, s, CV_32FC1);
//中心位置
int cs = (s - 1) / 2;
//方差
float sigma2 = pow(sigma, 2);
for (int c = 0; c < s; c++)
{
float norm2 = pow(float(c - cs), 2);
xkernel.at<float>(0, c) = exp(-norm2 / (2 * sigma2));
}
//将xkernel转置，得到垂直方向上的卷积核
Mat ykernel = xkernel.t();
//分离卷积核的卷积运算
Mat gauConv;
sepConv2D_Y_X(I, gauConv, CV_32F, xkernel, ykernel, Point(-1, -1), BORDER_DEFAULT);
gauConv.convertTo(gauConv, CV_32F, 1.0 / sigma2);
return gauConv;
}
//DoG 算子
Mat DoG(Mat I, float sigma, int s, float k)
{
//与标准差为sigma的非归一化的高斯核卷积
Mat Ig = gaussConv(I, sigma, s);
//与标准差为k*sigma的非归一化的高斯核卷积
Mat Igk = gaussConv(I, k*sigma, s);
//两个高斯卷积结果做差
Mat doG = Igk - Ig;
return doG;
}
8.10、Marr-Hildreth算子
对于LoG和DoG边缘检测，最后一步只是简单地进行阈值化处理，显然得到的边缘很粗略。那么Marr-Hildreth边缘检测可以简单地理解为对高斯差分和高斯拉普拉斯检 测到的边缘的细化，就像Canny对Sobel、Prewitt检测到的边缘的细化一样。
算法步骤如下：
第一步： 构建窗口大小为H×W 的高斯拉普拉斯或者高斯差分卷积核。
第二步： 图像矩阵与LoGH×W 核或者DoGH×W 核卷积。
第三步： 通过第二步得到的卷积结果寻找过零点位置，过零点位置即为边缘位置。
寻找过零点的方法
方法1：针对图像矩阵与高斯差分核（或者高斯拉普拉斯核）的卷积结果，对每一个位置判断以该位置为中心的3×3邻域内的上/下方向、左/右方向、左上/右下方向、右上/ 左下方向的值是否有异号出现。对于这四种情况，只要有一种情况出现异号，该位置（r， c）就是过零点，即为边 缘点。
方法2：与第一种方式类似，只是首先计算左上、右上、左下、右下的4个2×2 邻域内的均值对于这四个邻域内的均值，只要任意两个均值是异号的，该位置就是过零点，即为边缘点。
//第一种方法计算过零点
void zero_cross_default(InputArray _src, OutputArray _dst)
{
Mat src = _src.getMat();
_dst.create(src.size(), CV_8UC1);
Mat dst = _dst.getMat();
int rows = src.rows;
int cols = src.cols;
//零交叉点
for (int r = 1; r < rows - 2; r++)
{
for (int c = 1; c < cols - 2; c++)
{
//左上、右下
if (src.at<float>(r - 1, c - 1)*src.at<float>(r + 1, c + 1) < 0)
{
dst.at<uchar>(r, c) = 255;
continue;
}
//上、下
if (src.at<float>(r - 1, c)*src.at<float>(r + 1, c) < 0)
{
dst.at<uchar>(r, c) = 255;
continue;
}
//左、右
if (src.at<float>(r, c - 1)*src.at<float>(r, c + 1) < 0)
{
dst.at<uchar>(r, c) = 255;
continue;
}
//右上、左下
if (src.at<float>(r - 1, c + 1)*src.at<float>(r + 1, c - 1) < 0)
{
dst.at<uchar>(r, c) = 255;
continue;
}
}
}
}
//第二种方法计算过零点：均值
void zero_cross_mean(InputArray _src, OutputArray _dst)
{
Mat src = _src.getMat();
_dst.create(src.size(), CV_8UC1);
Mat dst = _dst.getMat();
int rows = src.rows;
int cols = src.cols;
double minValue;
double maxValue;
Mat temp(1, 4, CV_32FC1);
//零交叉点
for (int r = 1; r < rows - 1; r++)
{
for (int c = 1; c < cols - 1; c++)
{
//左上领域
Mat left_top(src, Rect(c - 1, r - 1, 2, 2));
temp.at<float>(0, 0) = mean(left_top)[0];
//右上领域
Mat right_top(src, Rect(c, r - 1, 2, 2));
temp.at<float>(0, 1) = mean(right_top)[0];
//左下领域
Mat left_bottom(src, Rect(c - 1, r, 2, 2));
temp.at<float>(0, 2) = mean(left_bottom)[0];
//右下领域
Mat right_bottom(src, Rect(c, r, 2, 2));
temp.at<float>(0, 3) = mean(right_bottom)[0];
//找出四个方向领域的最值
minMaxLoc(temp, &minValue, &maxValue);
if (minValue*minValue < 0)
{
dst.at<uchar>(r, c) = 255;
}
}
}
}
//Marr_Hildreth算子
Mat Marr_Hildreth(InputArray image, int win, float sigma, ZERO_CROSS_TYPE type)
{
//高斯拉普拉斯
Mat loG = LoG(image, sigma, win);
//过零点
Mat zeroCrossImage;
switch (type)
{
case ZERO_CROSS_DEFALUT:
zero_cross_default(loG, zeroCrossImage);
break;
case ZERO_CROSS_MEAN:
zero_cross_mean(loG, zeroCrossImage);
break;
default:
break;
}
return zeroCrossImage;
}
