Canny 边缘检测

算法代码如下：

#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <iostream>
#include <cmath>

using namespace cv;
using namespace std;

void TraceEdge(int y, int x, int nThrLow, Mat &edgeDefualt, Mat &gradient);

int main()
{
    double duration;
    Mat src = imread("stuff.bmp", 0);
    Mat srcBorder( src.rows+2, src.cols+2, CV_8UC1, Scalar(128));
    Mat dst( src.rows, src.cols, CV_8UC1, Scalar(0));
    Mat gradientX( src.rows, src.cols, CV_32FC1, Scalar(0));
    Mat gradientY( src.rows, src.cols, CV_32FC1, Scalar(0));
    Mat gradient( src.rows, src.cols, CV_8UC1, Scalar(0));
    Mat gradientDegree( src.rows, src.cols, CV_32FC1, Scalar(0));
    Mat edgeDefualt( src.rows, src.cols, CV_8UC1, Scalar(0));
    double temp = 0;
    double temp1 = 0;
    double temp2 = 0;
duration = static_cast<double>(getTickCount());
    for ( int i = 0; i < src.rows; i ++ ) //边界扩充
    {
        uchar *srcPtr = src.ptr<uchar>(i);
        uchar *srcBorderPtr = srcBorder.ptr<uchar>(i+1);
        for ( int j = 0; j < src.cols; j ++ )
        {
            srcBorderPtr[j+1] = srcPtr[j];
        }
    }

    for ( int i = 1; i < srcBorder.rows-1; i ++ ) //高斯滤波
    {
        uchar *srcBorderPtrPrevious = srcBorder.ptr<uchar>(i-1);
        uchar *srcBorderPtrCurrent  = srcBorder.ptr<uchar>(i);
        uchar *srcBorderPtrNext     = srcBorder.ptr<uchar>(i+1);
        uchar *dstPtr               = dst.ptr<uchar>(i-1);
        for ( int j = 1; j < srcBorder.cols-1; j ++ )
        {
            temp = 0;
            temp = srcBorderPtrPrevious[j-1] + 2*srcBorderPtrPrevious[j] +  srcBorderPtrPrevious[j+1]
                + 2*srcBorderPtrCurrent[j-1] + 4* srcBorderPtrCurrent[j] + 2*srcBorderPtrCurrent[j+1]
                +      srcBorderPtrNext[j-1] + 2*    srcBorderPtrNext[j] +      srcBorderPtrNext[j+1];
            temp /= 16;
            dstPtr[j-1] = temp;
        }
    }

    for ( int i = 1; i < dst.rows-1; i ++ ) //计算梯度幅值和梯度方向角
    {
        uchar *dstPtrPrevious = dst.ptr<uchar>(i-1);
        uchar *dstPtrCurrent  = dst.ptr<uchar>(i);
        uchar *dstPtrNext     = dst.ptr<uchar>(i+1);
        float *gradientXPtr         = gradientX.ptr<float>(i-1);
        float *gradientYPtr         = gradientY.ptr<float>(i-1);
        uchar *gradientPtr          = gradient.ptr<uchar>(i-1);
        float *gradientDegreePtr    = gradientDegree.ptr<float>(i-1);
        for ( int j = 1; j < dst.cols-1; j ++ )
        {
            //X方向梯度
            temp1 = - dstPtrPrevious[j-1] +  dstPtrPrevious[j+1]
                   - 2*dstPtrCurrent[j-1] + 2*dstPtrCurrent[j+1]
                   -      dstPtrNext[j-1] +      dstPtrNext[j+1];
            gradientXPtr[j-1] = temp1;
            //Y方向梯度
            temp2 = dstPtrPrevious[j-1] + 2*dstPtrPrevious[j] + dstPtrPrevious[j+1]
                      - dstPtrNext[j-1] -     2*dstPtrNext[j] -     dstPtrNext[j+1];
            gradientYPtr[j-1] = temp2;
            //计算梯度幅值
            temp = sqrt( temp1*temp1 + temp2*temp2) + 0.5;
            gradientPtr[j-1] = temp;
            //计算梯度方向角度
            temp = atan2(temp2, temp1) * 57.3;
            if(temp < 0) //将角度转化成 0~360° 范围内
                temp += 360;
            gradientDegreePtr[j-1] = temp;
        }
    }
//    Sobel( dst, gradientX, CV_32FC1, 1, 0, 3 );
//    Sobel( dst, gradientY, CV_32FC1, 0, 1, 3 );
//    addWeighted( gradientX, 0.5, gradientY, 0.5, 0, gradient );//求缩小后目标图片边缘梯度

    int g1=0, g2=0, g3=0, g4=0;                            //用于进行插值，得到亚像素点坐标值
    double dTmp1=0.0, dTmp2=0.0;                           //保存两个亚像素点插值得到的灰度数据
    double dWeight=0.0;                                    //插值的权重
    for ( int i = 1; i < edgeDefualt.rows-1; i ++ )
    {
        float *gradientXPtr        = gradientX.ptr<float>(i);
        float *gradientYPtr        = gradientY.ptr<float>(i);
        uchar *gradientPtrPrevious = gradient.ptr<uchar>(i-1);
        uchar *gradientPtrCurrent  = gradient.ptr<uchar>(i);
        uchar *gradientPtrNext     = gradient.ptr<uchar>(i+1);
        float *gradientDegreePtr   = gradientDegree.ptr<float>(i);
        uchar *edgeDefualtPtr      = edgeDefualt.ptr<uchar>(i);
        for ( int j = 1; j < edgeDefualt.cols-1; j ++ )
        {
            if(gradientPtrCurrent[j] == 0)
                edgeDefualtPtr[j] = 0; //如果当前梯度幅值为0，则不是局部最大对该点赋为0
            else
            {
                首先判断属于那种情况，然后根据情况插值///
                第一种情况///
                /       g1  g2                  /
                /           C                   /
                /           g3  g4              /
                /
                if( ((gradientDegreePtr[j]>=90)&&(gradientDegreePtr[j]<135)) ||
                    ((gradientDegreePtr[j]>=270)&&(gradientDegreePtr[j]<315)))
                {
                //根据斜率和四个中间值进行插值求解
                    g1 = gradientPtrPrevious[j-1];
                    g2 = gradientPtrPrevious[j];
                    g3 = gradientPtrNext[j];
                    g4 = gradientPtrNext[j+1];
                    dWeight = fabs(gradientXPtr[j])/fabs(gradientYPtr[j]);   //反正切
                    dTmp1 = g1*dWeight+g2*(1-dWeight);
                    dTmp2 = g4*dWeight+g3*(1-dWeight);
                }
                第二种情况///
                /       g1                      /
                /       g2  C   g3              /
                /               g4              /
                /
                else if( ((gradientDegreePtr[j]>=135)&&(gradientDegreePtr[j]<180)) ||
                    ((gradientDegreePtr[j]>=315)&&(gradientDegreePtr[j]<360)))
                {
                    g1 = gradientPtrPrevious[j-1];
                    g2 = gradientPtrCurrent[j-1];
                    g3 = gradientPtrCurrent[j+1];
                    g4 = gradientPtrNext[j+1];
                    dWeight = fabs(gradientYPtr[j])/fabs(gradientXPtr[j]);   //正切
                    dTmp1 = g2*dWeight+g1*(1-dWeight);
                    dTmp2 = g4*dWeight+g3*(1-dWeight);
                }
                第三种情况///
                /           g1  g2              /
                /           C                   /
                /       g4  g3                  /
                /
                else if( ((gradientDegreePtr[j]>=45)&&(gradientDegreePtr[j]<90)) ||
                    ((gradientDegreePtr[j]>=225)&&(gradientDegreePtr[j]<270)))
                {
                    g1 = gradientPtrPrevious[j];
                    g2 = gradientPtrPrevious[j+1];
                    g3 = gradientPtrNext[j];
                    g4 = gradientPtrNext[j-1];
                    dWeight = fabs(gradientXPtr[j])/fabs(gradientYPtr[j]);   //反正切
                    dTmp1 = g2*dWeight+g1*(1-dWeight);
                    dTmp2 = g3*dWeight+g4*(1-dWeight);
                }
                第四种情况///
                /               g1              /
                /       g4  C   g2              /
                /       g3                      /
                /
                else if( ((gradientDegreePtr[j]>=0)&&(gradientDegreePtr[j]<45)) ||
                    ((gradientDegreePtr[j]>=180)&&(gradientDegreePtr[j]<225)))
                {
                    g1 = gradientPtrPrevious[j+1];
                    g2 = gradientPtrCurrent[j+1];
                    g3 = gradientPtrCurrent[j-1];
                    g4 = gradientPtrNext[j-1];
                    dWeight = fabs(gradientYPtr[j])/fabs(gradientXPtr[j]);   //正切
                    dTmp1 = g1*dWeight+g2*(1-dWeight);
                    dTmp2 = g3*dWeight+g4*(1-dWeight);
                }
            }
            //进行局部最大值判断，并写入检测结果
            if((gradientPtrCurrent[j]>=dTmp1) && (gradientPtrCurrent[j]>=dTmp2))
                edgeDefualtPtr[j]= 128;
            else
                edgeDefualtPtr[j] = 0;
        }
    }

    int nHist[1024];
    int nEdgeNum;             //可能边界数
    int nMaxMag = 0;          //最大梯度数
    int nHighCount;
    for ( int i = 0; i < 1024; i ++ )
        nHist[i] = 0;
    for( int i = 0; i < edgeDefualt.rows; i ++ )
    {
        uchar *edgeDefualtPtr = edgeDefualt.ptr<uchar>(i);
        uchar *gradientPtrCurrent  = gradient.ptr<uchar>(i);
        for( int j = 0; j < edgeDefualt.cols; j ++ )
        {
            if( edgeDefualtPtr[j] == 128 )
                nHist[gradientPtrCurrent[j]]++;
        }
    }
    nEdgeNum = nHist[0];
    nMaxMag = 0;                    //获取最大的梯度值
    for( int i = 1; i < 1024; i ++ )           //统计经过“非最大值抑制”后有多少像素
    {
        if(nHist[i] != 0)       //梯度为0的点是不可能为边界点的
        {
            nMaxMag = i;
        }
        nEdgeNum += nHist[i];   //经过non-maximum suppression后有多少像素
    }
    double  dRatHigh = 0.9;
    double  dThrHigh;
    double  dThrLow;
    double  dRatLow = 0.5;
    nHighCount = (int)(dRatHigh * nEdgeNum + 0.5);
    int n = 1;
    nEdgeNum = nHist[1];
    while( (n < (nMaxMag-1)) && (nEdgeNum < nHighCount) )
    {
           n ++;
           nEdgeNum += nHist[n];
    }
    dThrHigh = n;                                   //高阈值
    dThrLow = (int)((dThrHigh) * dRatLow + 0.5);    //低阈值


    for( int i = 0; i < edgeDefualt.rows; i ++ ) //边界追踪
    {
        uchar *edgeDefualtPtr = edgeDefualt.ptr<uchar>(i);
        uchar *gradientPtrCurrent  = gradient.ptr<uchar>(i);
        for( int j = 0; j < edgeDefualt.cols; j ++ )
        {
            if((edgeDefualtPtr[j] == 128) && (gradientPtrCurrent[j] >= dThrHigh))
            {
                edgeDefualtPtr[j] = 255;
                TraceEdge(i, j, dThrLow, edgeDefualt, gradient);
            }
        }
    }
    for( int i = 0; i < edgeDefualt.rows; i ++ )
    {
        uchar *edgeDefualtPtr = edgeDefualt.ptr<uchar>(i);
        for(int j = 0; j < edgeDefualt.cols; j ++ )
        {
            if(edgeDefualtPtr[j] != 255)
            {
                edgeDefualtPtr[j]  = 0 ;   // 设置为非边界点
            }
        }
    }

duration = static_cast<double>(getTickCount()) - duration;
duration /= getTickFrequency();
cout << duration << endl;
    namedWindow("src", 0);
    imshow("src", src);
    namedWindow("dst", 0);//高斯滤波后的图像
    imshow("dst", dst);
    namedWindow("gradient", 0);//梯度图像
    imshow("gradient", gradient);
//    namedWindow("gradientDegree", 0);
//    imshow("gradientDegree", gradientDegree);
    namedWindow("edgeDefualt", 0);//边缘图像
    imshow("edgeDefualt", edgeDefualt);
    Mat dstt;
    Canny(src, dstt, 100, 200, 3);
    namedWindow("dstt", 0);
    imshow("dstt", dstt);
    waitKey(0);
    return 0;
}

void TraceEdge(int y, int x, int nThrLow, Mat &edgeDefualt, Mat &gradient)
{
    //对8邻域像素进行查询
    int xNum[8] = {1,1,0,-1,-1,-1,0,1};
    int yNum[8] = {0,1,1,1,0,-1,-1,-1};
    int yy,xx,k;
    for(k=0;k<8;k++)
    {
        yy = y+yNum[k];
        xx = x+xNum[k];
        uchar *edgeDefualtPtr = edgeDefualt.ptr<uchar>(yy);
        uchar *gradientPtrCurrent  = gradient.ptr<uchar>(yy);
        if(edgeDefualtPtr[xx] == 128 && gradientPtrCurrent[xx] >= nThrLow )
        {
            //该点设为边界点
            edgeDefualtPtr[xx] = 255;
            //以该点为中心再进行跟踪
            TraceEdge( yy, xx, nThrLow, edgeDefualt, gradient );
        }
    }
}