图像边缘检测

一 一阶微分

   函数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);

    }