Canny原理

Canny的原理就不细说了，冈萨雷斯的《数字图像处理》（第三版）P463~465讲解的比较清楚，主要就四个步骤：
1. 对图像进行高斯滤波
2. 计算梯度大小和梯度方向
3. 对梯度幅值图像进行非极大抑制
4. 双阈值处理和连接性分析（通常这一步与非极大抑制并行，详见下面的代码）

下面重点说一下非极大抑制。

非极大抑制

对一幅图像计算梯度大小和梯度方向后，需要进行非极大抑制，一般都是通过计算梯度方向，沿着梯度方向，判断该像素点的梯度大小是否是极大值。这里主要说一下方向的判断。

图像坐标系

《数字图像处理》（第三版）这本书中图像坐标顺序与上诉坐标系是相反的，不知道作者为什么这么写，有知道的朋友告诉我一下哈，看书的时候，注意一下坐标的顺序。

图像坐标系中，从x轴正方向往y轴正方向旋转为正（我没有查到权威文献资料，不知道这句话对不对）。

边缘方向区间

非极大抑制中，通常将边缘量化为4个方向，水平，垂直，45°和-45°，实际中，通过定义一个方向角方位，在该方位内认为是某一方向的边缘，实现中，我们通过计算梯度方向的范围从而判断边缘的方向（边缘的方向与梯度方向垂直）。

由于边缘方向没有正负，所以梯度方向$\theta$只考虑180°区间内。如上图所示：
1. 当$-22.5°<\theta<22.5°$的时候，梯度方向为0°，也就是垂直边缘
2. 当$22.5°<\theta<67.5°$的时候，梯度方向为45°，也就是135°边缘
3. 其他两个方向依次类推

实现中，通过三角函数的性质计算，比如第一种情况，使用下面条件来判断
$-22.5°<\theta<22.5°$ =>
$tan(-22.5°)<tan\theta<tan(22.5°)$ =>
$tan(-22.5°)<{fy \over fx } <tan(22.5°)$=>
$fx*tan(-22.5°)<{fy } <fx* tan(22.5°)$
其中$fx和fy$表示x方向和y方向的偏导数。

其他三种情况读者可以自行推导。

判断出梯度方向后，就可以进行非极大值抑制了。还是以第一种情况为例，比如我计算出了P5这个像素点处的梯度方向为0°，则这个时候，我们要判断的就是M(P5)>M(P4)&&M(P5)>M(P6),也就是P5是否是极大值，其中M表示该像素点的梯度大小。

下面是我根据上面的分析写出来的Canny的实现

Canny算法的实现（1.0版）

#define CANNY_SHIFT 16
#define TAN_225  (int)(0.4142135623730950488016887242097*(1 << CANNY_SHIFT));
#define TAN_675  (int)(2.4142135623730950488016887242097*(1 << CANNY_SHIFT));
void Canny1(const Mat &srcImage, Mat &dstImage, double lowThreshold, double highThreshold, int sizeOfAperture, bool L2)
{
    // 只支持灰度图
    CV_Assert(srcImage.type() == CV_8UC1);
    dstImage.create(srcImage.size(), srcImage.type());

    // L2范数计算边缘强度的时候，距离采用平方的方式，所以阈值也需要采用平方
    if (L2)
    {
        lowThreshold = std::min(32767.0, lowThreshold);
        highThreshold = std::min(32767.0, highThreshold);

        if (lowThreshold > 0) lowThreshold *= lowThreshold;
        if (highThreshold > 0) highThreshold *= highThreshold;
    }

    // 计算fx,fy,强度图
    Mat fx(srcImage.size(), CV_32SC1);
    Mat fy(srcImage.size(), CV_32SC1);
    Mat enlargedImage;
    Mat magnitudeImage(srcImage.rows + 2, srcImage.cols + 2, CV_32SC1);
    magnitudeImage.setTo(Scalar(0));
    copyMakeBorder(srcImage, enlargedImage, 1, 1, 1, 1, cv::BORDER_REPLICATE);
    int stepOfEnlargedImage = enlargedImage.cols;
    int stepOffx = fx.cols;
    int height = srcImage.rows;
    int width = srcImage.cols;
    uchar *rowOfEnlargedImage = enlargedImage.data + stepOfEnlargedImage + 1;
    int *rowOffx = (int *)fx.data;
    int *rowOffy = (int *)fy.data;
    int *rowOfMagnitudeImage = (int *)magnitudeImage.data + stepOfEnlargedImage + 1;
    for (int y = 0; y <= height - 1; ++y, rowOfEnlargedImage += stepOfEnlargedImage, rowOfMagnitudeImage += stepOfEnlargedImage, rowOffx += stepOffx, rowOffy += stepOffx)
    {
        uchar *colOfEnlargedImage = rowOfEnlargedImage;
        int *colOffx = rowOffx;
        int *colOffy = rowOffy;
        int *colOfMagnitudeImage = rowOfMagnitudeImage;
        for (int x = 0; x <= width - 1; ++x, ++colOfEnlargedImage, ++colOffx, ++colOffy, ++colOfMagnitudeImage)
        {
            // fx
            colOffx[0] = colOfEnlargedImage[1 - stepOfEnlargedImage] + 2 * colOfEnlargedImage[1] + colOfEnlargedImage[1 + stepOfEnlargedImage] -
                colOfEnlargedImage[-1 - stepOfEnlargedImage] - 2 * colOfEnlargedImage[-1] - colOfEnlargedImage[-1 + stepOfEnlargedImage];

            // fy
            colOffy[0] = colOfEnlargedImage[stepOfEnlargedImage - 1] + 2 * colOfEnlargedImage[stepOfEnlargedImage] + colOfEnlargedImage[stepOfEnlargedImage + 1] -
                colOfEnlargedImage[-stepOfEnlargedImage - 1] - 2 * colOfEnlargedImage[-stepOfEnlargedImage] - colOfEnlargedImage[-stepOfEnlargedImage + 1];

            // 计算边缘强度，由于只是用于比较，为了加快速度，只计算平方和
            if (L2)
            {
                colOfMagnitudeImage[0] = colOffx[0] * colOffx[0] + colOffy[0] * colOffy[0];
            }
            else
            {
                colOfMagnitudeImage[0] = std::abs(colOffx[0]) + std::abs(colOffy[0]);
            }

        }
    }


    // 非极大抑制，同时标记图做标记,双阈值处理
    //   0 - 可能是边缘
    //   1 - 不是边缘
    //   2 - 一定是边缘
    Mat labelImage(srcImage.rows + 2, srcImage.cols + 2, CV_8UC1);
    memset(labelImage.data, 1, labelImage.rows*labelImage.cols);
    rowOffx = (int *)fx.data;
    rowOffy = (int *)fy.data;
    rowOfMagnitudeImage = (int *)magnitudeImage.data + stepOfEnlargedImage + 1;
    uchar *rowOfLabelImage = labelImage.data + stepOfEnlargedImage + 1;
    queue<uchar*> queueOfEdgePixel;
    for (int y = 0; y <= height - 1; ++y, rowOfMagnitudeImage += stepOfEnlargedImage, rowOffx += stepOffx, rowOffy += stepOffx, rowOfLabelImage += stepOfEnlargedImage)
    {
        int *colOffx = rowOffx;
        int *colOffy = rowOffy;
        int *colOfMagnitudeImage = rowOfMagnitudeImage;
        uchar *colOfLabelImage = rowOfLabelImage;
        for (int x = 0; x <= width - 1; ++x, ++colOffx, ++colOffy, ++colOfMagnitudeImage, ++colOfLabelImage)
        {
            int fx = colOffx[0];
            int fy = colOffy[0];

            // 该像素点才有可能是边缘点
            if (colOfMagnitudeImage[0] > lowThreshold)
            {
                // 非极大抑制
                fy = fy*(1 << CANNY_SHIFT);

                int tan225 = fx * TAN_225;
                int tan675 = fx * TAN_675;

                // 梯度方向0
                if (fy>-1 * tan225 && fy < tan225)
                {
                    // 极大值，有可能是边缘
                    if (colOfMagnitudeImage[0] >= colOfMagnitudeImage[-1] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[1])
                    {
                        // 大于高阈值，是边缘，标记为2
                        if (colOfMagnitudeImage[0] > highThreshold)
                        {
                            // 进入队列，并设置标记
                            colOfLabelImage[0] = 2;
                            queueOfEdgePixel.push(colOfLabelImage);
                        }
                        else
                        {
                            // 有可能是边缘，标记为0
                            colOfLabelImage[0] = 0;

                        }
                    }
                }

                // 梯度方向45
                if (fy >= tan225&&fy <= tan675)
                {
                    // 极大值，有可能是边缘
                    if (colOfMagnitudeImage[0] > colOfMagnitudeImage[stepOfEnlargedImage + 1] && colOfMagnitudeImage[0]  > colOfMagnitudeImage[-stepOfEnlargedImage - 1])
                    {
                        // 大于高阈值，是边缘，标记为2
                        if (colOfMagnitudeImage[0] > highThreshold)
                        {
                            // 进入队列，并设置标记
                            colOfLabelImage[0] = 2;
                            queueOfEdgePixel.push(colOfLabelImage);
                        }
                        else
                        {
                            // 有可能是边缘，标记为0
                            colOfLabelImage[0] = 0;

                        }

                    }

                }

                // 梯度方向90
                if (fy>tan675 || fy<-tan675)
                {
                    // 极大值，有可能是边缘
                    if (colOfMagnitudeImage[0] >= colOfMagnitudeImage[stepOfEnlargedImage] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[-stepOfEnlargedImage])
                    {
                        // 大于高阈值，是边缘，标记为2
                        if (colOfMagnitudeImage[0] > highThreshold)
                        {
                            // 进入队列，并设置标记
                            colOfLabelImage[0] = 2;
                            queueOfEdgePixel.push(colOfLabelImage);
                        }
                        else
                        {
                            // 有可能是边缘，标记为0
                            colOfLabelImage[0] = 0;

                        }

                    }

                }

                // 梯度方向135
                if (fy >= -1 * tan675&&fy <= -1 * tan225)
                {
                    // 极大值，有可能是边缘
                    if (colOfMagnitudeImage[0] > colOfMagnitudeImage[stepOfEnlargedImage - 1] && colOfMagnitudeImage[0]  > colOfMagnitudeImage[-stepOfEnlargedImage + 1])
                    {
                        // 大于高阈值，是边缘，标记为2
                        if (colOfMagnitudeImage[0] > highThreshold)
                        {
                            // 进入队列，并设置标记
                            colOfLabelImage[0] = 2;
                            queueOfEdgePixel.push(colOfLabelImage);
                        }
                        else
                        {
                            // 有可能是边缘，标记为0
                            colOfLabelImage[0] = 0;

                        }

                    }

                }

            }
        }

    }

    // 连接性分析,这里采用队列实现（广度优先遍历）
    // 连接性分析也可以采用栈实现(深度优先遍历，OpenCV的做法)
    while (!queueOfEdgePixel.empty())
    {
        uchar *m = queueOfEdgePixel.front();
        queueOfEdgePixel.pop();

        // 在8领域搜索
        if (!m[-1])
        {
            m[-1] = 2;
            queueOfEdgePixel.push(m - 1);
        }

        if (!m[1])
        {
            m[1] = 2;
            queueOfEdgePixel.push(m + 1);
        }
        if (!m[-stepOfEnlargedImage - 1])
        {
            m[-stepOfEnlargedImage - 1] = 2;
            queueOfEdgePixel.push(m - stepOfEnlargedImage - 1);
        }
        if (!m[-stepOfEnlargedImage])
        {
            m[-stepOfEnlargedImage] = 2;
            queueOfEdgePixel.push(m - stepOfEnlargedImage);
        }
        if (!m[-stepOfEnlargedImage + 1])
        {
            m[-stepOfEnlargedImage + 1] = 2;
            queueOfEdgePixel.push(m - stepOfEnlargedImage + 1);
        }
        if (!m[stepOfEnlargedImage - 1])
        {
            m[stepOfEnlargedImage - 1] = 2;
            queueOfEdgePixel.push(m + stepOfEnlargedImage - 1);
        }
        if (!m[stepOfEnlargedImage])
        {
            m[stepOfEnlargedImage] = 2;
            queueOfEdgePixel.push(m + stepOfEnlargedImage);
        }
        if (!m[stepOfEnlargedImage + 1])
        {
            m[stepOfEnlargedImage + 1] = 2;
            queueOfEdgePixel.push(m + stepOfEnlargedImage + 1);
        }
    }

    // 最后生成边缘图
    rowOfLabelImage = labelImage.data + stepOfEnlargedImage + 1;
    uchar *rowOfDst = dstImage.data;
    for (int y = 0; y <= height - 1; ++y, rowOfLabelImage += stepOfEnlargedImage, rowOfDst += stepOffx)
    {
        uchar *colOfLabelImage = rowOfLabelImage;
        uchar *colOfDst = rowOfDst;
        for (int x = 0; x <= width - 1; ++x, ++colOfDst, ++colOfLabelImage)
        {
            if (colOfLabelImage[0] == 2)
                 colOfDst[0] = 255;
            else
            {
                colOfDst[0] = 0;
            }
        }
    }

}

Canny算法的实现(2.0版)

我看了OpenCV源码之后，将角度的判断修改为OpenCV的方式.

void Canny2(const Mat &srcImage, Mat &dstImage, double lowThreshold, double highThreshold, int sizeOfAperture, bool L2)
{
    // 只支持灰度图
    CV_Assert(srcImage.type() == CV_8UC1);
    dstImage.create(srcImage.size(), srcImage.type());

    // L2范数计算边缘强度的时候，距离采用平方的方式，所以阈值也需要采用平方
    if (L2)
    {
        lowThreshold = std::min(32767.0, lowThreshold);
        highThreshold = std::min(32767.0, highThreshold);

        if (lowThreshold > 0) lowThreshold *= lowThreshold;
        if (highThreshold > 0) highThreshold *= highThreshold;
    }

    // 计算fx,fy,强度图
    Mat fx(srcImage.size(), CV_32SC1);
    Mat fy(srcImage.size(), CV_32SC1);
    Mat enlargedImage;
    Mat magnitudeImage(srcImage.rows + 2, srcImage.cols + 2, CV_32SC1);
    magnitudeImage.setTo(Scalar(0));
    copyMakeBorder(srcImage, enlargedImage, 1, 1, 1, 1, cv::BORDER_REPLICATE);
    int stepOfEnlargedImage = enlargedImage.cols;
    int stepOffx = fx.cols;
    int height = srcImage.rows;
    int width = srcImage.cols;
    uchar *rowOfEnlargedImage = enlargedImage.data + stepOfEnlargedImage + 1;
    int *rowOffx = (int *)fx.data;
    int *rowOffy = (int *)fy.data;
    int *rowOfMagnitudeImage = (int *)magnitudeImage.data + stepOfEnlargedImage + 1;
    for (int y = 0; y <= height - 1; ++y, rowOfEnlargedImage += stepOfEnlargedImage, rowOfMagnitudeImage += stepOfEnlargedImage, rowOffx += stepOffx, rowOffy += stepOffx)
    {
        uchar *colOfEnlargedImage = rowOfEnlargedImage;
        int *colOffx = rowOffx;
        int *colOffy = rowOffy;
        int *colOfMagnitudeImage = rowOfMagnitudeImage;
        for (int x = 0; x <= width - 1; ++x, ++colOfEnlargedImage, ++colOffx, ++colOffy, ++colOfMagnitudeImage)
        {
            // fx
            colOffx[0] = colOfEnlargedImage[1 - stepOfEnlargedImage] + 2 * colOfEnlargedImage[1] + colOfEnlargedImage[1 + stepOfEnlargedImage] -
                colOfEnlargedImage[-1 - stepOfEnlargedImage] - 2 * colOfEnlargedImage[-1] - colOfEnlargedImage[-1 + stepOfEnlargedImage];

            // fy
            colOffy[0] = colOfEnlargedImage[stepOfEnlargedImage - 1] + 2 * colOfEnlargedImage[stepOfEnlargedImage] + colOfEnlargedImage[stepOfEnlargedImage + 1] -
                colOfEnlargedImage[-stepOfEnlargedImage - 1] - 2 * colOfEnlargedImage[-stepOfEnlargedImage] - colOfEnlargedImage[-stepOfEnlargedImage + 1];

            // 计算边缘强度，由于只是用于比较，为了加快速度，只计算平方和
            if (L2)
            {
                colOfMagnitudeImage[0] = colOffx[0] * colOffx[0] + colOffy[0] * colOffy[0];
            }
            else
            {
                colOfMagnitudeImage[0] = std::abs(colOffx[0]) + std::abs(colOffy[0]);
            }

        }
    }
#define CANNY_SHIFT 15
#define TG22  (int)(0.4142135623730950488016887242097*(1 << CANNY_SHIFT) + 0.5);

    // 遍历强度图，计算角度，并使用非极大抑制，同时标记图做标记
    //   0 - 可能是边缘
    //   1 - 不是边缘
    //   2 - 一定是边缘
    Mat labelImage(srcImage.rows + 2, srcImage.cols + 2, CV_8UC1);
    memset(labelImage.data, 1, labelImage.rows*labelImage.cols);
    rowOffx = (int *)fx.data;
    rowOffy = (int *)fy.data;
    rowOfMagnitudeImage = (int *)magnitudeImage.data + stepOfEnlargedImage + 1;
    uchar *rowOfLabelImage = labelImage.data + stepOfEnlargedImage + 1;
    queue<uchar*> queueOfEdgePixel;
    for (int y = 0; y <= height - 1; ++y, rowOfMagnitudeImage += stepOfEnlargedImage, rowOffx += stepOffx, rowOffy += stepOffx, rowOfLabelImage += stepOfEnlargedImage)
    {
        int *colOffx = rowOffx;
        int *colOffy = rowOffy;
        int *colOfMagnitudeImage = rowOfMagnitudeImage;
        uchar *colOfLabelImage = rowOfLabelImage;
        for (int x = 0; x <= width - 1; ++x, ++colOffx, ++colOffy, ++colOfMagnitudeImage, ++colOfLabelImage)
        {
            int xs = colOffx[0];
            int ys = colOffy[0];

            // 该像素点才有可能是边缘点
            if (colOfMagnitudeImage[0] > lowThreshold)
            {
                // 非极大抑制
                int x = std::abs(xs);
                int y = std::abs(ys) << CANNY_SHIFT;

                int tg22x = x * TG22;

                // 梯度方向0
                // |dy|/|dx|<0.414,计算出来-22.5<theta<22.5
                if (y < tg22x)
                {
                    // 极大值，有可能是边缘
                    if (colOfMagnitudeImage[0] > colOfMagnitudeImage[-1] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[1])
                    {
                         // 大于高阈值，是边缘，标记为2
                        if (colOfMagnitudeImage[0] > highThreshold)
                        {
                            // 进入队列，并设置标记
                            colOfLabelImage[0] = uchar(2);
                            queueOfEdgePixel.push(colOfLabelImage);
                        }
                        else
                        {
                             // 有可能是边缘，标记为0
                            colOfLabelImage[0] = 0;

                        }
                    }
                }
                else
                {
                    // 梯度方向90
                    int tg67x = tg22x + (x << (CANNY_SHIFT + 1));

                    // 水平边缘|dy|/|dx|>tan67.5=2.414,注意tan函数曲线计算出来67.5<theta<112.5
                    if (y > tg67x)
                    {
                        // 极大值，有可能是边缘
                        if (colOfMagnitudeImage[0] > colOfMagnitudeImage[stepOfEnlargedImage] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[-stepOfEnlargedImage])
                        {
                            // 大于高阈值，是边缘，标记为2
                            if (colOfMagnitudeImage[0] > highThreshold)
                            {
                                // 进入队列，并设置标记
                                colOfLabelImage[0] = 2;
                                queueOfEdgePixel.push(colOfLabelImage);
                            }
                            else
                            {
                                // 有可能是边缘，标记为0
                                colOfLabelImage[0] = 0;

                            }

                        }
                    }
                    else
                    {
                        // 梯度方向 +45°/-45°
                        int s = (xs ^ ys) < 0 ? -1 : 1;// ^:异或

                        // 极大值，有可能是边缘
                        if (colOfMagnitudeImage[0] > colOfMagnitudeImage[-stepOfEnlargedImage - s] && colOfMagnitudeImage[0]  > colOfMagnitudeImage[stepOfEnlargedImage + s])
                        {
                            // 大于高阈值，是边缘，标记为2
                            if (colOfMagnitudeImage[0] > highThreshold)
                            {
                                // 进入队列，并设置标记
                                colOfLabelImage[0] = 2;
                                queueOfEdgePixel.push(colOfLabelImage);
                            }
                            else
                            {
                                // 有可能是边缘，标记为0
                                colOfLabelImage[0] = 0;

                            }

                        }
                    }
                }
            }
        }

    }

    // 连接性分析,这里采用队列实现（广度优先遍历）
    // 连接性分析也可以采用栈实现(深度优先遍历，OpenCV的做法)
    while (!queueOfEdgePixel.empty())
    {
        uchar *m = queueOfEdgePixel.front();
        queueOfEdgePixel.pop();

        // 在8领域搜索
        if (!m[-1])
        {
            m[-1] = 2;
            queueOfEdgePixel.push(m - 1);
        }

        if (!m[1])
        {
            m[1] = 2;
            queueOfEdgePixel.push(m + 1);
        }
        if (!m[-stepOfEnlargedImage - 1])
        {
            m[-stepOfEnlargedImage - 1] = 2;
            queueOfEdgePixel.push(m - stepOfEnlargedImage - 1);
        }
        if (!m[-stepOfEnlargedImage])
        {
            m[-stepOfEnlargedImage] = 2;
            queueOfEdgePixel.push(m - stepOfEnlargedImage);
        }
        if (!m[-stepOfEnlargedImage + 1])
        {
            m[-stepOfEnlargedImage + 1] = 2;
            queueOfEdgePixel.push(m - stepOfEnlargedImage + 1);
        }
        if (!m[stepOfEnlargedImage - 1])
        {
            m[stepOfEnlargedImage - 1] = 2;
            queueOfEdgePixel.push(m + stepOfEnlargedImage - 1);
        }
        if (!m[stepOfEnlargedImage])
        {
            m[stepOfEnlargedImage] = 2;
            queueOfEdgePixel.push(m + stepOfEnlargedImage);
        }
        if (!m[stepOfEnlargedImage + 1])
        {
            m[stepOfEnlargedImage + 1] = 2;
            queueOfEdgePixel.push(m + stepOfEnlargedImage + 1);
        }
    }

    // 最后生成边缘图
    rowOfLabelImage = labelImage.data + stepOfEnlargedImage + 1;
    uchar *rowOfDst = dstImage.data;
    for (int y = 0; y <= height - 1; ++y, rowOfLabelImage += stepOfEnlargedImage, rowOfDst += stepOffx)
    {
        uchar *colOfLabelImage = rowOfLabelImage;
        uchar *colOfDst = rowOfDst;
        for (int x = 0; x <= width - 1; ++x, ++colOfDst, ++colOfLabelImage)
        {
            if (colOfLabelImage[0] == 2)
                colOfDst[0] = 255;
            else
            {
                colOfDst[0] = 0;
            }
        }
    }
}

实验结果

实验代码

int main(int argc, char *argv[])
{
    Mat srcImage = imread("D:/Image/Gray/Lena512.bmp", -1);
    Mat canny1,canny,canny2,canny3;
    Canny1(srcImage, canny1, 50, 150, 3, false);
    Canny2(srcImage, canny2, 50, 150, 3, false);
    Canny(srcImage, canny, 50, 150);
    imwrite("D:/Canny.bmp", canny);
    imwrite("D:/Canny1.bmp", canny1);
    imwrite("D:/Canny2.bmp", canny2);
    return 0;
}

使用的是标准的Lena图

OpenCV 的处理结果为：

Canny2 方法的结果:

Canny1 方法的结果:

结果表明，OpenCV的方式效果比自己实现的要好，具体原因还不太清楚，希望知道的朋友能留言，不胜感激。

2016-6-19 01:54:58