Canny边缘检测

1、概述

Canny边缘检测是多阶段优化的方法，该方法严格遵循以下三个准则：
1.1、高信噪比。
这一准则说明好的检测效果要具有高的信噪比，对于图像的边缘要能够实实在在地检测到，对于图像中非边缘的部分，不能够被检测成边缘。
1.2、控制对边缘检测的精度。
指算法所检测的边缘要尽可能的逼近实际的边缘，这是边缘检测算法的基本要求，检测精度直接反映了一个算法的好坏。
1.3、降低对同一边缘的相应次数。
对于检测的图像边缘，应该准确检测出边缘点，避免出现伪边缘。

2、步骤：

Canny边缘检测算法的具体的实现步骤为：
2.1、对图像进行高斯平滑
由于图像中的边缘信息对噪声比较敏感，因此使用滤波器来对噪声进行处理，这里使用的高斯平滑用来削弱图像中的噪声对边缘信息提取的影响，对于图像 ，有

其中，⊗表示卷积运算， ， ，σ表示高斯卷积的平滑因子。σ越大图像越平滑，对于噪声越不敏感，σ越小图像平缓程度越小，对噪声越敏感。
2.2、求梯度的方向和幅值
对经过高斯平滑后的图像再求其偏导，得到图像的偏导数和，对其求方向和幅值

2.3、非极大值抑制
图像中像素点导数的幅值大，但是不能够认定该点是边缘点。非极大值抑制就是指在图像中的局部最大值里，将非极大值点进行剔除，不对这些点进行进一步的检测，这样可以将一些由噪声引起的、出现局部梯度较大而可能会成为伪边缘的点剔除掉。实现起来就是将点与邻域的8个点按导数减小的方向进行幅值比较，若点小于该方向的两个邻域的点的值，那么该点不是极大值，在此处将该点梯度的幅值置0。
2.4、图像二值化
双阈值检测指的是首先使用高阈值进行对图像边缘的检测，检测到高于阈值的点就将置设为255，低于阈值的点就置为0。这时得到的图像边缘虽然比较真实，噪声也很少，但是图像边缘信息丢失比较严重，这时使用低阈值在高阈值的邻域进行二次检测，由于图像边缘一般都是连续的，所以检测到的新边缘可以填充高阈值检测到的图像边缘的缝隙中，完成边缘的提取。
阈值的取值关系到检测结果，阈值的取值对于图像边沿的提取至关重要，为了使选取的阈值能够适应不同场景的图像，选择使用图像直方图确定高阈值，统计梯度幅值大小为 的点的数量，得到函数

表示幅值大小为i的点的数量，则每个幅值大小相同的点在图像中的概率为

其中N表示图像中像素点的总数量，根据统计的梯度幅值大小的信息设置双阈值的高阈值，然后取高阈值的一半作为低阈值。
实验验证：
σ取0.3，高阈值取0.6
原图：

检测结果：

3、C语言实现：

/*/////////////////////////////////////////////////////////////////////函数名：int CannyEdge()结果：double Canny[IMGROWS][IMGCOLS];*//////////////////////////////////////////////////////////////////////int CannyEdge(){    int i,j,k,l;    int tempSize=0;    double m,n;    char name[20];    double sigma=0.3;//0.4;       double number;//中间变量    double gaussConall;    double conBuff=0.0;    double conBuffSum=0.0;    double g1=0, g2=0, g3=0, g4=0;     //中间变量                double dTmp1=0.0, dTmp2=0.0;    //中间变量              double dWeight=0.0;             //插值权重    int grayHist[50000]={0};    //只能用于灰度图像    double maxMag=0;          //最大梯度值    int notEdgeNum=0;    double  RatHigh = 0.6;  //直方图的比例    double  ThrHigh;      double  ThrLow;      double  RatLow = 0.6;      int HighCount=0;    FILE *fpTarget;    ////////////////////////高斯滤波  使用一维x2///////////////////////    tempSize=(int)(sigma*6+1.9);//保证模板尺寸    printf("\nsigma: %f\n",sigma);    if(!(tempSize&1))    {        tempSize++;    }    gaussConall=0;    number=sqrt(2.0*Pi)*sigma;    for(i=0;i<tempSize;i++)    {           m=i-tempSize/2;        gaussCon1[i]=(exp(-(m*m)/(2.0*sigma*sigma)))/number;        gaussConall+=gaussCon1[i];        //printf("gaussCon1   %d : %f\n",i,gaussCon1[i]);    }    for(i=0;i<tempSize;i++)//归一化    {        gaussCon1[i]/=gaussConall;        //printf("gaussCon1   %d : %f\n",i,gaussCon1[i]);    }    //printf("tempSize: %d\n",tempSize);    for(i=0;i<IMGROWS;++i)//横向    {        for(j=0;j<IMGCOLS;++j)        {                   conBuff=0;            conBuffSum=0;            for(k=0;k<tempSize;++k)            {                if(j-(tempSize-1)/2+k>=0&&j-(tempSize-1)/2+k<IMGCOLS)                {                    conBuff+=(gaussCon1[k]*(double)BMPHEAD.grayMatrix[i][j-(tempSize-1)/2+k]);                    conBuffSum+=gaussCon1[k];                }            }            imageMid[i][j]=conBuff/conBuffSum;              }    }/*for(i=0;i<IMGROWS;++i)*/    for(j=0;j<IMGCOLS;++j)//纵向    {        for(i=0;i<IMGROWS;++i)        {                   conBuff=0;            conBuffSum=0;            for(k=0;k<tempSize;++k)            {                if(i-(tempSize-1)/2+k>=0&&i-(tempSize-1)/2+k<IMGCOLS)                {                    conBuff+=(gaussCon1[k]*imageMid[i-(tempSize-1)/2+k][j]);                    conBuffSum+=gaussCon1[k];                }            }            image1[i][j]=conBuff/conBuffSum;                }    }/*for(i=0;i<IMGROWS;++i)*/    for(i=0;i<BMPHEAD.rows;i++)    {        for(j=0;j<BMPHEAD.cols;j++)        {            imageShow[i][j]=(unsigned char)image1[i][j];//(gradientMag[i][j]*100);    //      image1[i][j]=BMPHEAD.grayMatrix[i][j];        }    }    matrix2bmp8("gauss.bmp");    for(i=0; i<(IMGROWS-1); i++)//求偏导数    {          for(j=0; j<(IMGCOLS-1); j++)          {                dx[i][j] = (image1[i][j+1]-image1[i][j]+image1[i+1][j+1]-image1[i][j]+image1[i-1][j+1]-image1[i][j]+image1[i][j]-image1[i][j-1]+image1[i][j]-image1[i-1][j-1]+image1[i][j]-image1[i+1][j-1])/6.0;//+(image1[i+1][j+1]-image1[i][j]+image1[i+1][j-1]-image1[i][j]+image1[i-1][j+1]-image1[i][j]+image1[i-1][j-1]-image1[i][j])/8;                dy[i][j] = (image1[i+1][j]-image1[i][j]+image1[i+1][j+1]-image1[i][j]+image1[i+1][j-1]-image1[i][j]+image1[i][j]-image1[i-1][j]+image1[i][j]-image1[i-1][j-1]+image1[i][j]-image1[i-1][j+1])/6.0;//+(image1[i+1][j+1]-image1[i][j]+image1[i+1][j-1]-image1[i][j]+image1[i-1][j+1]-image1[i][j]+image1[i-1][j-1]-image1[i][j])/8;          }      }      for(i=0; i<IMGROWS; i++)//梯度幅值和方向    {          for(j=0; j<IMGCOLS; j++)          {                gradientMag[i][j] = sqrt(dx[i][j]*dx[i][j] + dy[i][j]*dy[i][j]);               gradientDir[i][j] = atan2(dy[i][j], dx[i][j]) * 180/Pi;                if(gradientDir[i][j] < 0)                      gradientDir[i][j] += 360;              //0--360        }      //  printf("%f\n",gradientDir[i][50]);    }      ////////////////极大值抑制//////////////////////    for(i=1; i<(IMGROWS-1); i++)      {          for(j=1; j<(IMGCOLS-1); j++)          {              if(gradientMag[i][j] == 0)              {                notGreat[i][j] = 0;         //如果当前梯度幅值为0，则不是局部最大，对该点赋为0              }            else              {              ////////首先判断属于那种情况，然后根据情况插值///////              ////////////////////第一种情况///////////////////////              /////////       g1  g2                  /////////////              /////////           C                   /////////////              /////////           g3  g4              /////////////              /////////////////////////////////////////////////////                  if( ((gradientDir[i][j]>=90)&&(gradientDir[i][j]<135)) ||                       ((gradientDir[i][j]>=270)&&(gradientDir[i][j]<315)))                  {                      //////根据斜率和四个中间值进行插值求解                      g1 = gradientMag[i-1][j-1];                      g2 = gradientMag[i-1][j];                      g3 = gradientMag[i+1][j];                      g4 = gradientMag[i+1][j+1];                      dWeight = fabs(dx[i][j])/fabs(dy[i][j]);   //反正切                      dTmp1 = g1*dWeight+g2*(1.0-dWeight);                      dTmp2 = g4*dWeight+g3*(1.0-dWeight);                  //  printf("90");                }              ////////////////////第二种情况///////////////////////              /////////       g1                      /////////////              /////////       g2  C   g3              /////////////              /////////               g4              /////////////              /////////////////////////////////////////////////////                  else if( ((gradientDir[i][j]>=135)&&(gradientDir[i][j]<180)) ||                       ((gradientDir[i][j]>=315)&&(gradientDir[i][j]<360)))                  {                      g1 = gradientMag[i-1][j-1];                      g2 = gradientMag[i][j-1];                      g3 = gradientMag[i][j+1];                      g4 = gradientMag[i+1][j+1];                      dWeight = fabs(dy[i][j])/fabs(dx[i][j]);   //正切                      dTmp1 = g2*dWeight+g1*(1.0-dWeight);                      dTmp2 = g4*dWeight+g3*(1.0-dWeight);                  //  printf("180");                }              ////////////////////第三种情况///////////////////////              /////////           g1  g2              /////////////              /////////           C                   /////////////              /////////       g4  g3                  /////////////              /////////////////////////////////////////////////////                  else if( ((gradientDir[i][j]>=45)&&(gradientDir[i][j]<90)) ||                       ((gradientDir[i][j]>=225)&&(gradientDir[i][j]<270)))                  {                      g1 = gradientMag[i-1][j];                      g2 = gradientMag[i-1][j+1];                      g3 = gradientMag[i+1][j];                      g4 = gradientMag[i+1][j-1];                      dWeight = fabs(dx[i][j])/fabs(dy[i][j]);   //反正切                      dTmp1 = g2*dWeight+g1*(1.0-dWeight);                      dTmp2 = g3*dWeight+g4*(1.0-dWeight);                  //  printf("270");                }                  ////////////////////第四种情况///////////////////////                  /////////               g1              /////////////                  /////////       g4  C   g2              /////////////                  /////////       g3                      /////////////                  /////////////////////////////////////////////////////                  else if( ((gradientDir[i][j]>=0)&&(gradientDir[i][j]<45)) ||                       ((gradientDir[i][j]>=180)&&(gradientDir[i][j]<225)))                  {                      g1 = gradientMag[i-1][j+1];                      g2 = gradientMag[i][j+1];                      g3 = gradientMag[i+1][j-1];                      g4 = gradientMag[i][j-1];                      dWeight = fabs(dy[i][j])/fabs(dx[i][j]);   //正切                      dTmp1 = g1*dWeight+g2*(1.0-dWeight);                      dTmp2 = g3*dWeight+g4*(1.0-dWeight);                     //printf("225");                }                }                 //////////进行局部最大值判断，并写入检测结果////////////////              if((gradientMag[i][j]>=dTmp1) && (gradientMag[i][j]>=dTmp2))              {                notGreat[i][j] = 128;                  //printf("1");            }            else              {                notGreat[i][j] = 0;                  //printf("\n");            }        }      }     test=0;    for(i=1; i<(IMGROWS-1); i++)      {        for(j=1; j<(IMGCOLS-1); j++)          {            if(gradientMag[i][j]>test)            {                test=gradientMag[i][j];            }        }       }    //////////////////////双阈值检测////////////////////////    for(i=1; i<(IMGROWS-1); i++)      {        for(j=1; j<(IMGCOLS-1); j++)          {            if(notGreat[i][j] == 128)            grayHist[(int)(gradientMag[i][j]*255.0/test)]++;        }    }/*  test=10000;    for(i=1; i<(IMGROWS-1); i++)      {        for(j=1; j<(IMGCOLS-1); j++)          {            if(gradientMag[i][j]<test)            {                test=gradientMag[i][j];            }        }       }    printf("////////////////////%f\n",test);*/    notEdgeNum=0;    for(i=1; i<50000; i++)       //统计有多少像素      {        if(grayHist[i] != 0)        {            maxMag = (double)(i)*test/255.0;          }        notEdgeNum += grayHist[i];     }    HighCount = (int)(RatHigh * notEdgeNum + 0.5);    printf("HighCount  :  %d\n",HighCount);    for(i=0;i<50000;++i)    {        HighCount-=grayHist[i];        printf("HighCount  :  %d\n",HighCount);        if(HighCount<0)        {            ThrHigh=(double)(i)*test/255.0;        //  printf("i  :  %d\n",i);        //              break;        }    }                                ThrLow = ThrHigh*RatLow;       printf("ThrHigh  :  %f\n",ThrHigh);    printf("ThrLow  :  %f\n",ThrLow);    for(i=1; i<(IMGROWS-1); i++)      {        for(j=1; j<(IMGCOLS-1); j++)          {            if((notGreat[i][j]==128) && (gradientMag[i][j] >= ThrHigh)) //强边界             {                  Canny[i][j] = 255;                //notGreat[i][j] = 255;              }          }    }    for(i=2; i<(IMGROWS-2); i++)      {        for(j=2; j<(IMGCOLS-2); j++)          {            if((Canny[i][j] != 255)&& (notGreat[i][j]==128) && (gradientMag[i][j]>= ThrLow) && (255==Canny[i+1][j]||255==Canny[i-1][j]||255==Canny[i][j+1]||255==Canny[i][j-1]||255==Canny[i+1][j+1]||255==Canny[i+1][j-1]||255==Canny[i-1][j+1]||255==Canny[i-1][j-1])) //弱边界                 {                      Canny[i][j] = 255;        //          printf("////////////////////%d\n",j);                }        }    //  printf("弱边界\n");    }    return 1;}