1. 边缘检测的概念
边缘检测是图像处理与计算机视觉中极为重要的一种分析图像的方法，至少在我做图像分析与识别时，边缘是我最喜欢的图像特征。边缘检测的目的就是找到图像中亮度变化剧烈的像素点构成的集合，表现出来往往是轮廓。在对现实世界的图像采集中，有下面4种情况会表现在图像中时形成一个边缘。
- 深度的不连续（物体处在不同的物平面上）；
- 表面方向的不连续（如正方体的不同的两个面）；
- 物体材料不同（这样会导致光的反射系数不同）；
- 场景中光照不同（如被树萌投向的地面）；
上面的图像是图像中水平方向7个像素点的灰度值显示效果，我们很容易地判断在第4和第5个像素之间有一个边缘，因为它俩之间发生了强烈的灰度跳变。在实际的边缘检测中，边缘远没有上图这样简单明显，我们需要取对应的阈值来区分出它们。
| 原图 | Roberts边缘 | Prewitt边缘 |
 |
 |
 |
| Sobel边缘 | Log边缘 | Canny边缘 |
 |
 |
 |
2. 边缘检测的基本方法
2.1 一阶微分边缘算子
一阶微分边缘算子也称为梯度边缘算子，它是利用图像在边缘处的阶跃性，即图像梯度在边缘取得极大值的特性进行边缘检测。梯度是
一个矢量，它具有方向$\theta$和模$|\Delta I|$：
$$\Delta I = \begin{pmatrix} \frac{\partial I}{\partial x} \\ \frac{\partial I}{\partial y} \end{pmatrix}$$
$$|\Delta I| = \sqrt{(\frac{\partial I}{\partial x} )^2+(\frac{\partial I}{\partial y} )^2} = \sqrt{I_x^2+I_y^2} $$
$$ \theta = arctan(I_y / I_x)$$
梯度的方向提供了边缘的趋势信息，因为梯度方向始终是垂直于边缘方向，梯度的模值大小提供了边缘的强度信息。
在实际使用中，通常利用有限差分进行梯度近似。对于上面的公式，我们有如下的近似：
$$\frac{\partial I}{\partial x} = \lim_{h\to 0}\frac{I(x+\Delta x,y)-I(x,y)}{\Delta x} \approx I(x+1,y)-I(x,y),(\Delta x = 1)$$
$$\frac{\partial I}{\partial y} = \lim_{h\to 0}\frac{I(x,y+\Delta xy)-I(x,y)}{\Delta y} \approx I(x,y+1)-I(x,y),(\Delta y = 1)$$
2.2 Roberts边缘检测算子
1963年，Roberts提出了这种寻找边缘的算子。Roberts边缘算子是一个2x2的模板，采用的是对角方向相邻的两个像素之差。从图像处理的实际效果来看，边缘定位较准，对噪声敏感。在Roberts检测算子中：
$$\frac{\partial I}{\partial x} = I(i,j)-I(i+1,j+1)$$
$$\frac{\partial I}{\partial y} = I(i+1,j)-I(i,j+1)$$
可以导出Roberts在点$(i+1/2,j+1/2)$处的水平与竖直边缘检测卷积核为：
$$m_x = \begin{bmatrix}1&0 \\ 0&-1\end{bmatrix}, m_y=\begin{bmatrix}0&-1 \\ 1&0\end{bmatrix}$$
2.3 Prewitt边缘检测算子
Prewitt利用周围邻域8个点的灰度值来估计中心的梯度，它的梯度计算公式如下：
$$\frac{\partial I}{\partial x} = I(i-1,j+1)+I(i,j+1)+I(i+1,j+1) – I(i-1, j-1)-I(i,j-1)-I(i+1,j-1)$$
$$\frac{\partial I}{\partial y} =-I(i+1,j-1)+I(i+1,j)+I(i+1,j+1) - I(i-1,j-1)+I(i-1,j)+I(i-1,j+1) $$
所以，Prewitt的卷积核为：
$$m_x = \begin{bmatrix}-1&0&+1 \\ -1&0&+1 \\ –1&0&+1\end{bmatrix}, m_y=\begin{bmatrix}-1&-1&-1 \\ 0&0&0 \\ +1&+1&+1\end{bmatrix}$$
2.4 Sobel边缘检测算子
比起Prewitt算子，Sobel也是用周围8个像素来估计中心像素的梯度，但是Sobel算子认为靠近中心像素的点应该给予更高的权重，所以Sobel算子把与中心像素4邻接的像素的权重设置为2或-2。
Sobel边缘检测算子的卷积核为：
$$m_x = \begin{bmatrix}-1&0&+1 \\ -2&0&+2\\ –1&0&+1\end{bmatrix}, m_y=\begin{bmatrix}-1&-2&-1 \\ 0&0&0 \\ +1&+2&+1\end{bmatrix}$$
Sobel进行边缘检测的实现可以参考我原来写的一篇博文：图像特征检测：sobel边缘检测，重要的是梯度图像计算后的阈值的确定与边缘的非极大值抑制算法，Roberts与Prewitt原理与sobel一致。
3. 二阶微分算子
学过微积分我们都知道，边缘即是图像的一阶导数局部最大值的地方，那么也意味着该点的二阶导数为零。二阶微分边缘检测算子就是利用图像在边缘处的阶跃性导致图像二阶微分在边缘处出现零值这一特性进行边缘检测的。
对于图像的二阶微分可以用拉普拉斯算子来表示：
$$\nabla^2I = \frac{\partial^2I}{\partial x^2}+\frac{\partial^2I}{\partial y^2}$$
我们在像素点$(i,j)$的$3\times 3$的邻域内，可以有如下的近似：
$$\frac{\partial^2I}{\partial x^2} = I(i,j+1)-2I(i,j)+I(i,j-1)$$
$$\frac{\partial^2I}{\partial y^2}=I(i+1,j)-2I(i,j)+I(i-1,j)$$
$$\nabla^2I = –4I(i,j)+I(i,j+1)+I(i,j-1)+I(i+1,j)+I(i-1,j)$$
对应的二阶微分卷积核为：
$$\boldsymbol{m} = \begin{bmatrix}0&1&0 \\ 1&4&1\\ 0&1&0\end{bmatrix}$$
所以二阶微分检测边缘的方法就分两步：1）用上面的Laplace核与图像进行卷积；2）对卷积后的图像，取得那些卷积结果为0的点。
虽然上述使用二阶微分检测边缘的方法简单，但它的缺点是对噪声十分敏感，同时也没有能够提供边缘的方向信息。为了实现对噪声的抑制，Marr等提出了LOG的方法。
为了减少噪声对边缘的影响，首先图像要进行低通滤波，LOG采用了高斯函数作为低通滤波器。高斯函数为：
$$G(x,y)=\frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}$$
上面的公式中$\sigma$决定了对图像的平滑程度。高斯函数生成的滤波模板尺寸一般设定为$6\sigma +1$（加1是会了使滤波器的尺寸为奇数）。使用高斯函数对图像进行滤波并对图像滤波结果进行二阶微分运算的过程，可以转换为先对高斯函数进行二阶微分，再利用高斯函数的二阶微分结果对图像进行卷积运算：
$$\nabla^2[G(x,y)*f(x,y)] = \nabla^2[G(x,y)]*f(x,y)$$
$$\nabla^2G(x,y) = -\frac{1}{2\pi\sigma^4}[1-\frac{x^2+y^2}{\sigma^2}]\cdot exp(-\frac{x^2+y^2}{2\sigma^2})$$
关于LOG算子的计算在上一篇文章：图像特征提取：斑点检测中有实现的代码。
4. Canny边缘检测
canny边缘检测实际上是一种一阶微分算子检测算法，但为什么这里拿出来说呢，因为它几乎是边缘检测算子中最为常用的一种，也是个人认为现在最优秀的边缘检测算子。Canny提出了边缘检测算子优劣评判的三条标准：
- 高的检测率。边缘检测算子应该只对边缘进行响应，检测算子不漏检任何边缘，也不应该将非边缘标记为边缘。
- 精确定位。检测到的边缘与实际边缘之间的距离要尽可能的小。
- 明确的响应。对每一条边缘只有一次响应，只得到一个点。
Canny边缘检测之所以优秀是因为它在一阶微分算子的基础上，增加了非最大值抑制和双阈值两项改进。利用非极大值抑制不仅可以有效地抑制多响应边缘，而且还可以提高边缘的定位精度；利用双阈值可以有效减少边缘的漏检率。
Canny边缘检测主要分四步进行：
- 去噪声；
- 计算梯度与方向角；
- 非最大值抑制；
- 滞后阈值化；
其中前两步很简单，先用一个高斯滤波器对图像进行滤波，然后用Sobel水平和竖直检测子与图像卷积，来计算梯度和方向角。
非极大值抑制
图像梯度幅值矩阵中的元素值越大，说明图像中该点的梯度值越大，但这不不能说明该点就是边缘（这仅仅是属于图像增强的过程）。在Canny算法中，非极大值抑制是进行边缘检测的重要步骤，通俗意义上是指寻找像素点局部最大值，将非极大值点所对应的灰度值置为0，这样可以剔除掉一大部分非边缘的点。
根据上图可知，要进行非极大值抑制，就首先要确定像素点C的灰度值在其8值邻域内是否为最大。图中蓝色的线条方向为C点的梯度方向，这样就可以确定其局部的最大值肯定分布在这条线上，也即出了C点外，梯度方向的交点dTmp1和dTmp2这两个点的值也可能会是局部最大值。因此，判断C点灰度与这两个点灰度大小即可判断C点是否为其邻域内的局部最大灰度点。如果经过判断，C点灰度值小于这两个点中的任一个，那就说明C点不是局部极大值，那么则可以排除C点为边缘。这就是非极大值抑制的工作原理。
在理解的过程中需要注意以下两点：
1）中非最大抑制是回答这样一个问题：“当前的梯度值在梯度方向上是一个局部最大值吗？” 所以,要把当前位置的梯度值与梯度方向上两侧的梯度值进行比较；
2）梯度方向垂直于边缘方向。
但实际上，我们只能得到C点邻域的8个点的值，而dTmp1和dTmp2并不在其中，要得到这两个值就需要对该两个点两端的已知灰度进行线性插值，也即根据图中的g1和g2对dTmp1进行插值，根据g3和g4对dTmp2进行插值，这要用到其梯度方向。
滞后阈值化
由于噪声的影响，经常会在本应该连续的边缘出现断裂的问题。滞后阈值化设定两个阈值：一个为高阈值$T_h$，一个为低阈值$T_l$。如果任何像素边缘算子的影响超过高阈值，将这些像素标记为边缘；响应超过低阈值（高低阈值之间）的像素，如果与已经标记为边缘的像素4-邻接或8-邻接，则将这些像素也标记为边缘。所以不整个过程描述如下：
- 如果该像素的梯度值小于$T_l$，则该像素为非边缘像素；
- 如果该像素的梯度值大于$T_h$，则该像素为边缘像素；
- 如果该像素的梯度值介于$T_l$与$T_h$之间，需要进一步检测该像素的$3\times 3$邻域内的8个点，如果这8个点内有一个或以上的点梯度超过了$T_h$，则该像素为边缘像素，否则不是边缘像素。
Canny边缘检测的实现
#include <iostream>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
using namespace std;
using namespace cv;
void getCannyEdge(const Mat& imgSrc, Mat& imgDst, double lowThresh = -1, double highThresh = -1, double sigma = 1);
double getGaussianThresh(const Mat& inputArray, double percentage);
int main(int argc, char** argv)
{
Mat image = imread(argv[1]);
Mat imgCanny;
getCannyEdge(image,imgCanny);
t = ((double)getTickCount() - t)/getTickFrequency();
imwrite("imgCanny.png", imgCanny);
return 0;
}
void getCannyEdge(const Mat& imgSrc, Mat& imgDst, double lowThresh, double highThresh, double sigma)
{
Mat gray;
if (imgSrc.channels() == 3)
{
cvtColor(imgSrc, gray, CV_BGR2GRAY);
}
else
{
gray = imgSrc.clone();
}
gray.convertTo(gray, CV_64F);
gray = gray / 255;
double gaussianDieOff = .0001;
double percentOfPixelsNotEdges = .7; // Used for selecting thresholds
double thresholdRatio = .4; // Low thresh is this fraction of the high
int possibleWidth = 30;
double ssq = sigma * sigma;
for (int i = 1; i <= possibleWidth; i++)
{
if (exp(-(i * i) / (2* ssq)) < gaussianDieOff)
{
possibleWidth = i - 1;
break;
}
}
if (possibleWidth == 30)
{
possibleWidth = 1; // the user entered a reallly small sigma
}
// get the 1D gaussian filter
int winSz = 2 * possibleWidth + 1;
Mat gaussKernel1D(1, winSz, CV_64F);
double* kernelPtr = gaussKernel1D.ptr<double>(0);
for (int i = 0; i < gaussKernel1D.cols; i++)
{
kernelPtr[i] = exp(-(i - possibleWidth) * (i - possibleWidth) / (2 * ssq)) / (2 * CV_PI * ssq);
}
// get the derectional derivatives of gaussian kernel
Mat dGaussKernel(winSz, winSz, CV_64F);
for (int i = 0; i < dGaussKernel.rows; i++)
{
double* linePtr = dGaussKernel.ptr<double>(i);
for (int j = 0; j< dGaussKernel.cols; j++)
{
linePtr[j] = - (j - possibleWidth) * exp(-((i - possibleWidth) * (i - possibleWidth) + (j - possibleWidth) * (j - possibleWidth)) / (2 * ssq)) / (CV_PI * ssq);
}
}
/* smooth the image out*/
Mat imgSmooth;
filter2D(gray, imgSmooth, -1, gaussKernel1D);
filter2D(imgSmooth, imgSmooth, -1, gaussKernel1D.t());
/*apply directional derivatives*/
Mat imgX, imgY;
filter2D(imgSmooth, imgX, -1, dGaussKernel);
filter2D(imgSmooth, imgY, -1, dGaussKernel.t());
Mat imgMag;
sqrt(imgX.mul(imgX) + imgY.mul(imgY), imgMag);
double magMax;
minMaxLoc(imgMag, NULL, &magMax, NULL, NULL);
if (magMax > 0 )
{
imgMag = imgMag / magMax;
}
if (lowThresh == -1 || highThresh == -1)
{
highThresh = getGaussianThresh(imgMag, percentOfPixelsNotEdges);
lowThresh = thresholdRatio * highThresh;
}
Mat imgStrong = Mat::zeros(imgMag.size(), CV_8U);
Mat imgWeak = Mat::zeros(imgMag.size(), CV_8U);
for (int dir = 1; dir <= 4; dir++)
{
Mat gradMag1(imgMag.size(), imgMag.type());
Mat gradMag2(imgMag.size(), imgMag.type());
Mat idx = Mat::zeros(imgMag.size(), CV_8U);
if (dir == 1)
{
Mat dCof = abs(imgY / imgX);
idx = (imgY <= 0 & imgX > -imgY) | (imgY >= 0 & imgX < -imgY);
idx.row(0).setTo(Scalar(0));
idx.row(idx.rows - 1).setTo(Scalar(0));
idx.col(0).setTo(Scalar(0));
idx.col(idx.cols - 1).setTo(Scalar(0));
for (int i = 1; i < imgMag.rows - 1; i++)
{
for (int j = 1; j < imgMag.cols - 1; j++)
{
gradMag1.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i,j + 1) + dCof.at<double>(i,j) * imgMag.at<double>(i - 1,j + 1);
gradMag2.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i,j - 1) + dCof.at<double>(i,j) * imgMag.at<double>(i + 1,j - 1);
}
}
}
else if(dir == 2)
{
Mat dCof = abs(imgX / imgY);
idx = (imgX > 0 & -imgY >= imgX) | (imgX < 0 & -imgY <= imgX);
for (int i = 1; i < imgMag.rows - 1; i++)
{
for (int j = 1; j < imgMag.cols - 1; j++)
{
gradMag1.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i - 1,j) + dCof.at<double>(i,j) * imgMag.at<double>(i - 1,j + 1);
gradMag2.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i + 1,j) + dCof.at<double>(i,j) * imgMag.at<double>(i + 1,j - 1);
}
}
}
else if(dir == 3)
{
Mat dCof = abs(imgX / imgY);
idx = (imgX <= 0 & imgX > imgY) | (imgX >= 0 & imgX < imgY);
for (int i = 1; i < imgMag.rows - 1; i++)
{
for (int j = 1; j < imgMag.cols - 1; j++)
{
gradMag1.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i - 1,j) + dCof.at<double>(i,j) * imgMag.at<double>(i - 1,j - 1);
gradMag2.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i + 1,j) + dCof.at<double>(i,j) * imgMag.at<double>(i + 1,j + 1);
}
}
}
else
{
Mat dCof = abs(imgY / imgX);
idx = (imgY <0 & imgX <= imgY) | (imgY > 0 & imgX >= imgY);
for (int i = 1; i < imgMag.rows - 1; i++)
{
for (int j = 1; j < imgMag.cols - 1; j++)
{
gradMag1.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i,j - 1) + dCof.at<double>(i,j) * imgMag.at<double>(i - 1,j - 1);
gradMag2.at<double>(i,j) = (1 - dCof.at<double>(i,j)) * imgMag.at<double>(i,j + 1) + dCof.at<double>(i,j) * imgMag.at<double>(i + 1,j + 1);
}
}
}
Mat idxLocalMax = idx & ((imgMag >= gradMag1) & (imgMag >= gradMag2));
imgWeak = imgWeak | ((imgMag > lowThresh) & idxLocalMax);
imgStrong= imgStrong| ((imgMag > highThresh) & imgWeak);
}
imgDst = Mat::zeros(imgWeak.size(),imgWeak.type());
for (int i = 1; i < imgWeak.rows - 1; i++)
{
uchar* pWeak = imgWeak.ptr<uchar>(i);
uchar* pDst = imgDst.ptr<uchar>(i);
uchar* pStrPre = imgStrong.ptr<uchar>(i - 1);
uchar* pStrMid = imgStrong.ptr<uchar>(i);
uchar* pStrAft = imgStrong.ptr<uchar>(i + 1);
for (int j = 1; j < imgWeak.cols - 1; j++)
{
if (!pWeak[j])
{
continue;
}
if (pStrMid[j])
{
pDst[j] = 255;
}
if (pStrMid[j-1] || pStrMid[j+1] || pStrPre[j-1] || pStrPre[j] || pStrPre[j+1] || pStrAft[j-1] || pStrAft[j] ||pStrAft[j+1])
{
pDst[j] = 255;
}
}
}
}
double getGaussianThresh(const Mat& inputArray, double percentage)
{
double thresh = -1.0;
// compute the 64-hist of inputArray
int nBins = 64;
double minValue, maxValue;
minMaxLoc(inputArray, &minValue, &maxValue, NULL, NULL);
double step = (maxValue - minValue) / nBins;
vector<unsigned> histBin(nBins,0);
for (int i = 0; i < inputArray.rows; i++)
{
const double* pData = inputArray.ptr<double>(i);
for(int j = 0; j < inputArray.cols; j++)
{
int index = (pData[j] - minValue) / step;
histBin[index]++;
}
}
unsigned cumSum = 0;
for (int i = 0; i < nBins; i++)
{
cumSum += histBin[i];
if (cumSum > percentage * inputArray.rows * inputArray.cols)
{
thresh = (i + 1) / 64.0;
break;
}
}
return thresh;
}
5. 参考资料
1. 《图像局部不变特征与描述》王永明，王贵锦。