1. 边缘检测的概念
边缘检测是图像处理与计算机视觉中极为重要的一种分析图像的方法,至少在我做图像分析与识别时,边缘是我最喜欢的图像特征。边缘检测的目的就是找到图像中亮度变化剧烈的像素点构成的集合,表现出来往往是轮廓。在对现实世界的图像采集中,有下面4种情况会表现在图像中时形成一个边缘。
- 深度的不连续(物体处在不同的物平面上);
- 表面方向的不连续(如正方体的不同的两个面);
- 物体材料不同(这样会导致光的反射系数不同);
- 场景中光照不同(如被树萌投向的地面);
原图 Roberts边缘 Prewitt边缘 Sobel边缘 Log边缘 Canny边缘
2. 边缘检测的基本方法
2.1 一阶微分边缘算子
一阶微分边缘算子也称为梯度边缘算子,它是利用图像在边缘处的阶跃性,即图像梯度在边缘取得极大值的特性进行边缘检测。梯度是
梯度的方向提供了边缘的趋势信息,因为梯度方向始终是垂直于边缘方向,梯度的模值大小提供了边缘的强度信息。
在实际使用中,通常利用有限差分进行梯度近似。对于上面的公式,我们有如下的近似:
2.2 Roberts边缘检测算子
可以导出Roberts在点 (i+1/2,j+1/2) 处的水平与竖直边缘检测卷积核为:
2.3 Prewitt边缘检测算子
Prewitt利用周围邻域8个点的灰度值来估计中心的梯度,它的梯度计算公式如下:
2.4 Sobel边缘检测算子
比起Prewitt算子,Sobel也是用周围8个像素来估计中心像素的梯度,但是Sobel算子认为靠近中心像素的点应该给予更高的权重,所以Sobel算子把与中心像素4邻接的像素的权重设置为2或-2。
3. 二阶微分算子
学过微积分我们都知道,边缘即是图像的一阶导数局部最大值的地方,那么也意味着该点的二阶导数为零。二阶微分边缘检测算子就是利用图像在边缘处的阶跃性导致图像二阶微分在边缘处出现零值这一特性进行边缘检测的。
我们在像素点 (i,j) 的 3×3 的邻域内,可以有如下的近似:
所以二阶微分检测边缘的方法就分两步:1)用上面的Laplace核与图像进行卷积;2)对卷积后的图像,取得那些卷积结果为0的点。
虽然上述使用二阶微分检测边缘的方法简单,但它的缺点是对噪声十分敏感,同时也没有能够提供边缘的方向信息。为了实现对噪声的抑制,Marr等提出了LOG的方法。
为了减少噪声对边缘的影响,首先图像要进行低通滤波,LOG采用了高斯函数作为低通滤波器。高斯函数为:
关于LOG算子的计算在上一篇文章:图像特征提取:斑点检测中有实现的代码。
4. Canny边缘检测
- 高的检测率。边缘检测算子应该只对边缘进行响应,检测算子不漏检任何边缘,也不应该将非边缘标记为边缘。
- 精确定位。检测到的边缘与实际边缘之间的距离要尽可能的小。
- 明确的响应。对每一条边缘只有一次响应,只得到一个点。
其中前两步很简单,先用一个高斯滤波器对图像进行滤波,然后用Sobel水平和竖直检测子与图像卷积,来计算梯度和方向角。
非极大值抑制
滞后阈值化
- 如果该像素的梯度值小于 Tl ,则该像素为非边缘像素;
- 如果该像素的梯度值大于 Th ,则该像素为边缘像素;
- 如果该像素的梯度值介于 Tl 与 Th 之间,需要进一步检测该像素的 3×3 邻域内的8个点,如果这8个点内有一个或以上的点梯度超过了 Th ,则该像素为边缘像素,否则不是边缘像素。
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;
} -
文章引子:
http://www.cnblogs.com/ronny/p/4001910.html