图像边缘检测
一 一阶微分
函数f(x, y)的一阶微分构成梯度grad(f)：
，梯度幅度mag： 
，梯度方向
为：
，梯度方向垂直于边缘方向。
在离散情况下，需要将一阶微分转换为一阶差分，具体如下：
考虑一维函数g(x)，其泰勒展开式为：
，
求解一阶导数为：
，其误差为：
；
使用
与联合求解得：
， 其误差为：
。
基于以上分析，函数f(x, y)的梯度可离散化为：
。
在实际图像边缘检测中，会对离散梯度函数进行适当扩展，这样就形成了Prewitt算子（前两组）及Sobel算子（后两组），具体如下：
，
。
Sobel算子在Prewitt算子基础上改变了中心像素的权重，在编码中，一般会优先使用Sobel算子，这里有必要给出更大尺寸的Sobel算子，如5*5， 7*7等；
如果对Prewitt算子进行更大尺寸扩展，直接扩展尺寸即可，但对于Sobel算子的扩展，由于各点权重不一致，还需要确定扩展后的权重值；仅考虑一个维度上的数据(1, 2, 1), 当扩展为更大尺寸时应该如何选择权重，如(x1, x2, x3, x4, x5)？
这里应该采用高斯函数来确定
系数；假设图像噪声为正态分布，高斯平滑被认为是最优的平滑方式；由于差分运算为线性运算，对各点应用差分后，其噪声也应该满足正态分布，使用高斯函数确定
系数可以突出信号并且最大程度抑制噪声。
高斯函数：
，其能量主要集中在距离原点
区间内，针对模板尺寸为5时，有
，带人
到高斯函数中，可求解系数。
则5*5Sobel算子为：
。
使用一阶微分算子，在图像边缘上会产生较大的响应，可以通过设定合适的阈值来检测图像边缘。
使用一阶微分检测边缘前需要使用高斯函数对图像进行滤波处理。高斯函数剔除图像高频成分保留低频部分，从而有效避免图像随机噪声干扰，如下图：
二 二阶微分
函数f(x)的一阶微分：
，
函数f(x)的二阶微分： 
，
函数f(x, y)的二阶偏微分构成laplace算子：
函数f(x,y)的二阶微分构成laplace算子：
，其矩阵形式如下：
。
二阶微分强调灰度突变，并不 强调灰度缓慢变化区域，利用这种特性，可以将灰度突变叠加到原始图像中以实现图像边缘锐化，其方法为：
。
二阶微分会产生过零点，可以通过检测二阶微分过零点来检测图像边缘。
三 Marr-Hildreth边缘检测
假设图像中噪声为信号无关噪声，满足平均值为0，符合正态分布，可使用如下步骤检测边缘：
1）使用高斯滤波可以减小加性噪声的影响；
2）在平滑后图像上应用laplace算子；
3）检测图像零交叉点，作为图像边缘。
将以上操作用数学表达式描述为：
，其中f(x, y)为图像函数，G(x, y)为高斯函数；
微分与卷积均为线性运算，上式可改写为：
，这样就可以提前计算出
（laplace of Gaussian）；
，
，
，
；
欲形成5*5拉普拉斯模板，根据公式有，带人（x, y)坐标可得出离散高斯拉普拉斯模板：
。
使用高斯拉普拉斯模板与图像卷积运算，完成1）2）步操作，接下来需要检测过零点以确定图像边缘，可能的方法有：
1）在3*3区域内检测图像左/右， 上/下或对角线符号变化，若符号发生改变，则确定为过零点；
2）将3*3区域划分为4个象限，每个象限4个像素（同一像素可能被划分到不同的象限），求每个象限的平均值，若存在两个平均值符号不同的象限，则确定为过零点。
四 Canny边缘检测
canny边缘检测思路如下：
1）使用高斯滤波剔除白噪声，使用一阶微分算子（如Sobel）与图像卷积（以上操作可先对高斯函数做一次卷积，再与图像卷积）；
2）计算卷积后图像梯度幅值与梯度方向；
3）在梯度方向上使用非极大值抑制，细化边缘；
4）使用滞后阈值方案（高低阈值）连接边缘，高阈值确定边缘，低阈值连接边缘。
void cannyTrace(int y,int x,int width,int widthstep,int thresh,uchar* pImageData,int* pMag) { int y_trace[8] = {-1,-1,-1, 0,0, 1,1,1}; int x_trace[8] = {-1, 0, 1,-1,1,-1,0,1}; for(int k = 0;k < 8;k++) { int yy = y + y_trace[k]; int xx = x + x_trace[k]; if(*(pImageData + yy * widthstep + xx) == 128 && pMag[yy * width + xx] > thresh) { *(pImageData + yy * widthstep + xx) = 255; cannyTrace(yy,xx,width,widthstep,thresh,pImageData,pMag); } } } void cannyGradHV(unsigned char* img, int* dx, int* dy, int* mag, float* theta, int width, int height) { int widthstep = (width + 3) / 4 * 4; for(short y = 1; y < height - 1; ++y) { for(short x = 1; x < width - 1; ++x) { *(dx + y * width + x) = *(img + (y - 1) * widthstep + x - 1) + *(img + (y ) * widthstep + x - 1) * 2 + *(img + (y + 1) * widthstep + x - 1) - *(img + (y - 1) * widthstep + x + 1) - *(img + (y ) * widthstep + x + 1) * 2 - *(img + (y + 1) * widthstep + x + 1); *(dy + y * width + x) = *(img + (y - 1) * widthstep + x - 1) + *(img + (y - 1) * widthstep + x ) * 2 + *(img + (y - 1) * widthstep + x + 1) - *(img + (y + 1) * widthstep + x - 1) - *(img + (y + 1) * widthstep + x ) * 2 - *(img + (y + 1) * widthstep + x + 1); *(mag + y * width + x) = abs(*(dx + y * width + x)) + abs(*(dy + y * width + x)); *(theta + y * width + x) = atan2((float)(*(dy + y * width + x)),(float)(*(dx + y * width + x))); } } }
void CannyCustomized(IplImage* pImage, IplImage* pCanny, float upper_threshold, float lower_threshold, bool hysteresis, CANNY_TPYE type) { IplImage* pFilterImage = cvCreateImage(cvGetSize(pImage),IPL_DEPTH_8U,1); uchar* pFilterImageData = (uchar*)pFilterImage->imageData; //cvSmooth(pImage,pFilterImage,CV_GAUSSIAN,5,1); cvCopy(pImage, pFilterImage); int width = pImage->width; int height = pImage->height; int widthstep = pImage->widthStep; // 梯度值与梯度方向 int* dx = new int[pImage->width * pImage->height]; int* dy = new int[pImage->width * pImage->height]; int* mag = new int[pImage->width * pImage->height]; float* theta = new float[pImage->width * pImage->height]; switch(type) { case CANNY_HORIZONTAL_VERTICAL: cannyGradHV(pFilterImageData, dx, dy, mag, theta, pImage->width, pImage->height); break; case CANNY_HORIZONTAL: cannyGradH(pFilterImageData, dx, dy, mag, theta, pImage->width, pImage->height); break; case CANNY_VERTICAL: cannyGradV(pFilterImageData, dx, dy, mag, theta, pImage->width, pImage->height); break; case CANNY_CUSTOM: cannyGradC(pFilterImageData, dx, dy, mag, theta, pImage->width, pImage->height); break; default: cannyGradHV(pFilterImageData, dx, dy, mag, theta, pImage->width, pImage->height); } // 非极大值抑制 IplImage* pRestrainImage = cvCreateImage(cvGetSize(pImage),IPL_DEPTH_8U,1); uchar* pRestrainImageData = (uchar*)pRestrainImage->imageData; cvZero(pRestrainImage); for(short y = 1; y < pFilterImage->height - 1; y++) { for(short x = 1; x < pFilterImage->width - 1; x++) { int index = y * width + x; if(mag[index] > 0) { int temp1 = 0; int temp2 = 0; /* 以下为图像梯度示意 * g1 g2 * c * g3 g4 * c为当前素 梯度角度从g1 g2 之间 到 g3 g4 之间 */ if((theta[index] >= CV_PI / 2 && theta[index] < CV_PI * 3 / 4) || (theta[index] >= -CV_PI / 2 && theta[index] < -CV_PI / 4 )){ int g1 = mag[index + width - 1]; int g2 = mag[index + width ]; int g3 = mag[index - width ]; int g4 = mag[index - width + 1]; double weight = fabs((double)dx[index] / (double)dy[index]); temp1 = g1 * weight + g2 * (1.0 - weight); temp2 = g4 * weight + g3 * (1.0 - weight); /* 以下为图像梯度示意 * g1 * g2 c g3 * g4 * c为当前素 梯度角度从g1 g2 之间 到 g3 g4 之间 */ }else if((theta[index] >= CV_PI * 3 / 4 && theta[index] < CV_PI) || (theta[index] >= -CV_PI / 4 && theta[index] < 0 )){ int g1 = mag[index + width - 1]; int g2 = mag[index - 1 ]; int g3 = mag[index + 1 ]; int g4 = mag[index - width + 1]; double weight = fabs((double)dy[index] / (double)dx[index]); temp1 = g1 * weight + g2 * (1.0 - weight); temp2 = g4 * weight + g3 * (1.0 - weight); /* 以下为图像梯度示意 * g1 * g3 c g2 * g4 * c为当前素 梯度角度从g1 g2 之间 到 g3 g4 之间 */ }else if((theta[index] >= 0 && theta[index] < CV_PI / 4) || (theta[index] >= -CV_PI && theta[index] < -CV_PI * 3 / 4)){ int g1 = mag[index + width + 1]; int g2 = mag[index + 1 ]; int g3 = mag[index - 1 ]; int g4 = mag[index - width - 1]; double weight = fabs((double)dy[index] / (double)dx[index]); temp1 = g1 * weight + g2 * (1.0 - weight); temp2 = g4 * weight + g3 * (1.0 - weight); /* 以下为图像梯度示意 * g2 g1 * c * g4 g3 * c为当前素 梯度角度从g1 g2 之间 到 g3 g4 之间 */ }else if((theta[index] >= CV_PI / 4 && theta[index] < CV_PI / 2) || (theta[index] >= -CV_PI * 3 / 4 && theta[index] < -CV_PI / 2)){ int g1 = mag[index + width + 1]; int g2 = mag[index + width ]; int g3 = mag[index - width ]; int g4 = mag[index - width - 1]; double weight = fabs((double)dx[index] / (double)dy[index]); temp1 = g1 * weight + g2 * (1.0 - weight); temp2 = g4 * weight + g3 * (1.0 - weight); } // 标记128为候选边缘 if(mag[index] >= temp1 && mag[index] >= temp2) *(pRestrainImageData + y * widthstep + x) = 128; } } } // 确定高低阈值 short upper = upper_threshold; short lower = lower_threshold; if(upper_threshold < 1.0 || lower_threshold < 1.0) { int hist[1024]; memset(hist, 0, sizeof(int) * 1024); for(short y = 1; y < pFilterImage->height - 1; y++) { for(short x = 1; x < pFilterImage->width - 1; x++) { ++hist[mag[y * width + x]]; } } int threshold = (pFilterImage->width - 2) * (pFilterImage->height - 2) * 0.85; if(upper_threshold > 0.01) threshold = (pFilterImage->width - 2) * (pFilterImage->height - 2) * upper_threshold; int sum = 0; for(short i = 0; i < 1024; ++i) { sum += hist[i]; if(sum > threshold) { upper = i; lower = upper * 0.5; if(lower_threshold > 0.01) lower = upper * lower_threshold; break; } } } // 根据高低阈值标记边缘 if(hysteresis) { for(short y = 0; y < pFilterImage->height; ++y) { for(short x = 0; x < pFilterImage->width; ++x) { if(*(pRestrainImageData + y * widthstep + x) == 128 && mag[y * width + x] > upper) { *(pRestrainImageData + y * widthstep + x) = 255; cannyTrace(y,x,width,widthstep,lower,pRestrainImageData,mag); } } } } else { for(short y = 0; y < pFilterImage->height; ++y) { for(short x = 0; x < pFilterImage->width; ++x) { if(*(pRestrainImageData + y * widthstep + x) == 128 && mag[y * width + x] > upper) *(pRestrainImageData + y * widthstep + x) = 255; } } } // 剔除无效边缘 for(short y = 0; y < pFilterImage->height; ++y) { for(short x = 0; x < pFilterImage->width; ++x) { if(*(pRestrainImageData + y * widthstep + x) != 255) *(pRestrainImageData + y * widthstep + x) = 0; } } cvCopy(pRestrainImage, pCanny); delete []dx; delete []dy; delete []mag; delete []theta; cvReleaseImage(&pFilterImage); cvReleaseImage(&pRestrainImage); }
