openCV入门----边缘检测（三）：canny算法

        对于canny算法，这个应用非常广泛算法，我着实想不到什么很好的开场词来介绍它。那么就套用《Learning openCV》中文版中的一段好了：

    “在图像边缘检测中，抑制噪声和边缘精确定位是无法同时满足的，一些边缘检测算法通过平滑滤波去除噪声的同事，也增加了边缘定位的不确定性；而提高边缘检测算子对边缘的敏感性的同事，也提高了对噪声的敏感性。canny算子力求在抗噪声干扰和精确定位之间寻求最佳折中方案。”

      回顾一下各种边缘检测的算法，正如上述所论。比如说，一阶算子检测边缘，如Robert，虽然对边缘定位精度较高，但是容易丢失边缘，对噪声无抑制能力，又如Sobel，对噪声仅有一点抑制能力，但无法完全排除虚假边缘；又如二阶算子-拉普拉斯算子，对噪声的处理更让人烦恼无比（对噪声响应很高）；又或是改进之后的高斯-拉普拉斯算法，虽然经过了高斯模糊，但是这些模糊不能完全去除噪声，没处理的好，反而会影响边缘检测的效果（也就是降低真实边缘的响应）···这种种都说明我们需要一个更牛逼的算法来完成边缘检测的艰巨任务，不用多说，那就是-------------canny算法！

     正如上文所说，canny旨在抗噪和精确定位之间寻求最佳折中方案，这就想数学建模中的多约束求最优解问题，那么首先要明确的是我们的约束是什么：

    （1）      信噪比准则： 信噪比越高越佳，错误率也就越小，噪声干扰越小；

    （2）      定位精度准则：检测出的边缘要尽可能接近真实边缘；

    （3）      单边缘响应准备：想拉普拉斯算子对于阶跃边缘的处理，往往会出现两条响应边缘，这是我们不希望出现的；

此处我略掉了一些数学公式（至少我现在看上去没那么有用），我认为明白这些概念就可以了，接下来我将想大家介绍canny处理的步骤，然后再回过头来对照准则，加深理解，也算是弥补理论上的不足吧！

Canny算法的步骤如下：

    （1）用高斯滤波器平滑图像；

    （2）用一阶偏导的有限差分计算梯度的幅值和方向；

                                            Hx = [ 1, -1; 1, -1 ]              Hy = [ -1, -1; 1, 1 ]

             设图像数据用 f 表示，则卷积后得到结果：

                                            phiX = f * Hx                        phiY = f * Hy

              以此来计算梯度的幅值和方向：

              梯度幅值：               sqrt（phiX^2 + phiY^2）

              梯度方向：               arctan（phiY / phiX）

    （3）对梯度幅值进行非极大值抑制。什么是非极大值抑制？从字面上理解，就是对非极大值的数据进行抑制，也可以理解成对非极大值数据排除其是边缘的可能性。那么，“极大值”这个概念如何理解呢？首先，我们要打破常规的思维----极大值往往是运用在数学中的连续函数中，这是一个相对的概念，实际上是指一定区域上的最大值或者最小值。而图像是由一个一个像素点构成，宏观上看就是一个离散的矩阵，那么我们要在矩阵中使用极大值的概念，那么我们首先要划分一个区域，有区域才有极值而言------8领域是最好的选择。

 如上图所示，以C点为中心的8领域（即C点为我们要讨论的对象）。首先明确，我们讨论的“极大值”是针对梯度幅值的，并且是在算出来的梯度方向上的梯度幅值极大值！所以，这个问题实际上被描述成：中心点C点的梯度幅值在其梯度方向上且在八领域范围内是否是最大的？否则抑制！

而很容易知道我们算出来的梯度方向应该是0°~360°的，并且是连续的角度。但是我们的八领域提供的是八个固定的离散角度（0°（360°），45°，90°，135°，180°，225°，270°分别对应上图的4,3,2,1,8,7,6,5区域）。所以，要算出梯度方向上的、且在中心点C八领域内的梯度幅值，那么就需要使用插值运算。举个例子来说，对于上图C点，若其梯度幅值算出为30°，要得到该梯度方向上、C点八领域内的梯度幅值，实际上是在区域3与区域4之间进行插值（相对的还有区域7和区域8之间进行插值，此图省略），具体如何插值呢？下文程序处将详细说明！

           （4）双阈值算法检测和连接边缘。

            首先也要明确，阈值是针对梯度幅值而言的！双阈值言下之意就是设置两个阈值TH,TL，经验上来说TH:TL = 2：1最佳。在高阈值TH的截取下，将得到一高阈值图像，配合步骤（3）中非极大值抑制的结果，于是我们可以得到较为一些边缘曲线。由于双重判断、并且阈值高，判断出来的边缘一般不会存在假边缘，但是正是因为阈值高，会出现边缘间断的现象。此时就要进行基于低阈值的边缘修补。总的来说，双阈值的作用是：

                         高阈值：判断真边缘的位置在哪里！

                         低阈值：修补判断出来的真边缘！

具体步骤：

 高阈值处理得高阈值图像  -> 搜寻高阈值图像中真边缘端点 -> 基于低阈值修补真边缘端点附近边缘             

 

接下来，懒得啰嗦了，直接上关键程序。程序中有很详细地注释。必要地方我会再做解释：

目标函数：IplImage* UseCanny(IplImage* img, int iSize, double sigma, double Th, double Tl)

输入参数： img    ------     待处理图像

                  iSize   ------     高斯滤波模板尺寸

                  sigma ------     高斯滤波标准差

                  Th      ------     双阈值中的高阈值

                   Tl       ------    双阈值中的低阈值

输出参数：处理完的图像

 

按照步骤，该函数各部分如下：

 /*
 step 1 ： 高斯滤波,得到处理好的图像guass
 */

 IplImage* MyImg = cvCreateImage(cvGetSize(img), 8, 1);
 cvConvertImage(img, MyImg, CV_BGR2GRAY);
 IplImage* guass = cvCreateImage(cvGetSize(MyImg), 8, 1);
 cvSmooth(MyImg, guass, CV_GAUSSIAN, iSize, iSize);

注：这里要注意，canny算法的输入图像必须是一个单通道图像，也就是灰度图像，考虑到大部分读进来的图像是BGR格式的，所以这里要先转成灰度图像。

 

/*
 step 2 ： 计算梯度幅值、方向
 这里使用到一阶水平、竖直模板  H1 = [-1， -1;  1， 1];            H2 = [1， -1 ;  1， -1];
 用公式表示：  H1----   f(i+1,j)-f(i,j)+f(i+1,j+1)-f(i,j+1)
               H2----   f(i,j)-f(i+1,j)+f(i,j+1)-f(i+1,j+1)
 */
 //初始化二维数组，准备存放结果
 double **difX;
 double **difY;
 double **value;
 double **dir;
 difX = new double*[guass->height];
 for (int i = 0; i < guass->height; ++i)
 {
  difX[i] = new double[guass->width];
 }
 difY = new double*[guass->height];
 for (int i = 0; i < guass->height; ++i)
 {
  difY[i] = new double[guass->width];
 }
 value = new double*[guass->height];
 for (int i = 0; i < guass->height; ++i)
 {
  value[i] = new double[guass->width];
 }
 dir = new double*[guass->height];
 for (int i = 0; i < guass->height; ++i)
 {
  dir[i] = new double[guass->width];
 }
 //正式计算过程
 for (int i = 0; i < guass->height - 1; ++i)
 {
  for (int j = 0; j < guass->width - 1; ++j)
  {
   difX[i][j] = ((char *)(guass->imageData + guass->widthStep * i))[j+1] - ((char *)(guass->imageData + guass->widthStep * i))[j]
    + ((uchar *)(guass->imageData + guass->widthStep * (i + 1)))[j + 1] - ((uchar *)(guass->imageData + guass->widthStep * (i + 1)))[j];//使用H1：f(i+1,j)-f(i,j)+f(i+1,j+1)-f(i,j+1)
   difY[i][j] = ((char *)(guass->imageData + guass->widthStep * (i + 1)))[j] - ((char *)(guass->imageData + guass->widthStep * i))[j]
    + ((uchar *)(guass->imageData + guass->widthStep * (i + 1)))[j + 1] - ((uchar *)(guass->imageData + guass->widthStep * i))[j + 1];//使用H1：f(i+1,j)-f(i,j)+f(i+1,j+1)-f(i,j+1)
   value[i][j] = sqrt(difX[i][j]*difX[i][j] + difY[i][j]*difY[i][j]);//计算梯度幅值
   dir[i][j] = atan2(difY[i][j], difX[i][j]);//计算梯度,返回-pi~pi
   if (dir[i][j] < 0)
   {
    dir[i][j] = dir[i][j] + 2 * P;//统一到0~2*pi
   }
  }
 }

 

注：这里计算梯度幅值、方向并没有什么技术含量，直接按照卷积模板进行运算即可。这里的P在程序前定义成了3.1415926常量。

 

/*
 step 3 ：非极大值抑制
 */
 double **MyResult;//定义一个灰度结果数组，放入非极大值抑制结果
 int g1 = 0, g2 = 0, g3 = 0, g4 = 0;//插值的边界值
 double weight = 0;//插值权重系数,这个系数是由角度的偏向确定
 double Aresult = 0, Bresult = 0;//将产生上下两个结果，A为上，B为下
 MyResult = new double*[guass->height];
 for (int i = 0; i < guass->height; ++i)
 {
  MyResult[i] = new double[guass->width];
 }
 //图像边缘不可能是边缘线，故认为清零
 for (int i = 1; i < guass->height - 1; ++i)
 {
  MyResult[i][0] = 0;
  MyResult[i][guass->width - 1] = 0;
 }
 for (int i = 1; i < guass->width - 1; ++i)
 {
  MyResult[0][i] = 0;
  MyResult[guass->height - 1][i] = 0;
 }
 
 for (int i = 1; i < guass->height - 1; ++i)
 {
  for (int j = 1; j < guass->width - 1; ++j)
  {
   /*
   情况1：  g1  g2    很明显，角度在（P/4 , p/2）||（5*P/4 , 3*P/2）之间
                    c
            g3  g4
   */
   if (((dir[i][j] >= P / 4) && (dir[i][j] <= P / 2)) || ((dir[i][j] >= 5 * P / 4) && (dir[i][j]) <= 3 * P / 2))
   {
    g1 = value[i - 1][j];
    g2 = value[i - 1][j + 1];
    g3 = value[i + 1][j + 1];
    g4 = value[i + 1][j];
    weight = fabs(difX[i][j]) / fabs(difY[i][j]);//1/正切值，保证权重小于1
    Aresult = (1 - weight) * g1 + weight * g2;
    Bresult = (1 - weight) * g4 + weight * g3;
   }
   /*
   情况2：g1  g2    很明显，角度在（P/2 , 3*P/4）||（3*P/2 , 7*P/4）之间
                        c
                        g3  g4
   */
   else if (((dir[i][j] >= P / 2) && (dir[i][j] <= 3 * P / 4)) || ((dir[i][j] >= 3 * P / 2) && (dir[i][j]) <= 7 * P / 4))
   {
    g1 = value[i - 1][j - 1];
    g2 = value[i - 1][j];
    g3 = value[i + 1][j];
    g4 = value[i + 1][j + 1];
    weight = fabs(difX[i][j]) / fabs(difY[i][j]);//正切值
    Aresult = (1 - weight) * g2 + weight * g1;
    Bresult = (1 - weight) * g3 + weight * g4;
   }
   /*
   情况3：     g1      很明显，角度在（0 , P/4）||（P , 5*P/4）之间
            g3 c   g2
            g4
   */
   else if (((dir[i][j] >= 0) && (dir[i][j] <= P / 4)) || ((dir[i][j] >= P) && (dir[i][j]) <= 5 * P / 4))
   {
    g1 = value[i - 1][j + 1];
    g2 = value[i][j + 1];
    g3 = value[i][j - 1];
    g4 = value[i + 1][j - 1];
    weight = fabs(difY[i][j]) / fabs(difX[i][j]);//正切值
    Aresult = (1 - weight) * g2 + weight * g1;
    Bresult = (1 - weight) * g3 + weight * g4;
   }
   /*
   情况4：g1         很明显，角度在（3*P/4 , P）||（7*P/4 , 2*P）之间
                 g2 c g3
                         g4
   */
   else if (((dir[i][j] >= 3 * P / 4) && (dir[i][j] <= P)) || ((dir[i][j] >= 7 * P / 4) && (dir[i][j]) <= 2 * P))
   {
    g1 = value[i - 1][j - 1];
    g2 = value[i][j - 1];
    g3 = value[i][j + 1];
    g4 = value[i + 1][j + 1];
    weight = fabs(difY[i][j]) / fabs(difX[i][j]);//正切值
    Aresult = (1 - weight) * g2 + weight * g1;
    Bresult = (1 - weight) * g3 + weight * g4;
   }
   if ((value[i][j] >= Aresult) && (value[i][j] >= Bresult))
   {
    MyResult[i][j] = 128;
   }
   else
   {
    MyResult[i][j] = 0;
   }
  }
 }

 

注：这里需要好好理解，四种情况如注释解释得很清楚了。这里唯一需要好好注意的就是weight权重系数的确定。首先，我们要保证weight在[0,1]区间内，读者一定要仔细注意weight的表达式，四种情况的表达式是不一样的：当梯度方向的tan值大于1，weight就要等于1/tan；当梯度方向的tan小于1，则weight直接等于tan即可（tan即为Y/X）。

接下来，又一个看似简单的问题了，就比如情况1中的下式：

                                        Aresult = (1 - weight) * g1 + weight * g2;
                                        Bresult = (1 - weight) * g4 + weight * g3;

读者自然而然就会想到，这里weight为什么是和g2、g3相乘，而（1-weight）为什么是和g1、g4相乘？

这是由g1、g2、g3、g4的位置以及weight的表达式所决定的

                                         weight = fabs(difX[i][j]) / fabs(difY[i][j]);

也就是说是X值比上Y值。再观察情况一的个g的位置：
                                        g1  g2    
                                         c
                                 g3  g4
可以这么假设，如果weight变大，则从表达式可以看出X变大，于是梯度方向就往X变大的方向偏，也就是往g2、g3的方向偏，实际上就是离g2、g3更近些，那么自然而然g2、g3的比重就要大些。这就说明了weight与g2、g3是呈正比的关系，故weight是和g2、g3相乘，同理说明（1-weight）。

 

不知你是否也和我一样，考虑过为什么要将weight定义成这个函数，大家网上找到的估计都是这个版本，我想答案应该很简单，对，就是很简单。插值只要能放映趋势就行，离谁近谁的比重就大，那取法实际上有很多。例如可以拿角度来做文章，梯度线离g1有30°，离g2有15°，那么g1的比重为30°/45°，g2的比重为15°/45°，这样也是可以的，只是做起来相对来说更加复杂。

 

 /*
  step 4 ：双阈值算法检测和连接边缘
  这里忽略双阈值的正确与否，TH>TL，并且理论上Th：Tl = 2：1，3：1为好.注：这个阈值是对梯度而言的
  */
  //首先使用高阈值得到初步图像，确定边缘，排除假边缘
  for (int i = 0; i < guass->height; ++i)
  {
   for (int j = 0; j < guass->width; ++j)
   {
    if ((MyResult[i][j] == 128) && (value[i][j] >= Th))
    {
     MyResult[i][j] = 255;
    }
   }
  }
  //再遍历结果，查看边缘端点，并利用低阈值修补边缘
  for (int i = 1; i < guass->height - 1; ++i)
  {
   for (int j = 1; j < guass->width - 1; ++j)
   {
    if (isSingle(MyResult, i, j))
    {
     TrackEdge(MyResult, value, Tl, i, j);
    }
   }
  }
  //将非边缘点置0
  for (int i = 0; i < guass->height; ++i)
  {
   for (int j = 0; j < guass->width; ++j)
   {
    if (MyResult[i][j] != 255)
    {
     MyResult[i][j] = 0;
    }
   }
  }

 这里使用到了两个函数isSingle和TrackEdge，分别定义如下：

/*
这是一个检测间断点的函数
*/
bool isSingle(double** temp, int i, int j)
{
 //对8邻域像素进行查询 
 int xNum[8] = { 1,1,0,-1,-1,-1,0,1 };
 int yNum[8] = { 0,1,1,1,0,-1,-1,-1 };
 int xx, yy;
 for (int k = 0; k < 8; k++)
 {
  yy = i + yNum[k];
  xx = j + xNum[k];
  if (temp[yy][xx] != 0)
  {
   return false;
  }
 }
 return true;
}

/*
这是一个通过低阈值判断边缘的函数
嵌套函数，当无边缘时停止
*/
void  TrackEdge(double** result, double** tidu, int TL, int i, int j)
{
 //对8邻域像素进行查询 
 int xNum[8] = { 1,1,0,-1,-1,-1,0,1 };
 int yNum[8] = { 0,1,1,1,0,-1,-1,-1 };
 int xx, yy;
 for (int k = 0; k < 8; k++)
 {
  yy = i + yNum[k];
  xx = j + xNum[k];
  if ((result[yy][xx] == 128) && (tidu[yy][xx] >= TL))
  {
   result[yy][xx] = 255;
   TrackEdge(result, tidu, TL, yy, xx);
  }
 }
 return;
}

 

注：一定要注意，步骤里说的很清楚，使用高阈值图像要遇到间断点或者说是端点的时候，再使用低阈值进行边缘修补，那就一定要严格按照此来执行。网上有些直接是走一步高阈值就马上走一步低阈值，当然，结果一定还是会出来，只是可以想想，某些地方是不是会出现噪声，canny算法的精髓是取最大梯度方向作为边缘判定！！！！

 

最后得到的MyResult数组，其中只含有0或者255，那就是我们想要的图像数据了！！！！！把他直接赋给图像对象并返回即可！看看效果吧！

 

实验结果：

 自己编写的程序：TH=200，TL=100，iSize=3，sigma=0.5
 

电脑的效果

 

感觉还是有差距的！这是为什么呢？参数选取几乎都是一样的。从图像上来看，应该是修补时出现的问题，导致边缘出现了一些断点、重复边缘之类的情况，后续有时间我会做一些修改和调整，争取让结果更加完美！

 

 

 

当然，故事还没有结束。很多人会想，canny算法中还有两个阈值如何确定呢？是的，双阈值的选取很大程度上决定了边缘检测的效果，高阈值取大了，可能会漏掉边缘，取低了，可能会带来噪声。。。这个问题，我会在后续进一步分析！