Canny的C++实现

好的图像边缘检测应该满足

1. 尽可能的标记处所有的边缘；
2. 标记出的边缘就是实际图像内容的边缘；
3. 图像中的边缘只标记一次。

Canny边缘检测的步骤

1. 图像灰度化；
2. 高斯滤波（不限制）降噪；
3. 使用简单的算子（如Sobel、Prewitt等）检测图像水平$grad\_x$、垂直梯度$grad\_y$，并求出梯度$\sqrt{grad\_x^2+grad\_y^2}$和梯度方向（角度）$grad\_y/grad\_x$；
4. 对梯度进行NMS；
5. 使用双阈值法确定强边缘、非边缘、弱边缘，以及弱边缘的二次判定。

C++实现示例

auxiliary.hpp

#pragma once

#include <cmath>
#include <numeric>   // std::accumulate
#include <vector>
#include <opencv2/imgproc.hpp>

/*  f(x,y) = exp(-(x^2 + y^2) / (2 * sigma^2)) / (2 * pi * sigma ^ 2) */
int GaussianKernel(int kernel_size, std::vector<std::vector<float>> &kernel)
{
    kernel.clear();
    kernel.resize(kernel_size);
    for (auto &it : kernel)
    {
        it.resize(kernel_size);
    }

    std::vector<int> coord_val(kernel_size, -kernel_size / 2);  //  坐标取值范围,一般取0左右对称的相反数
    for (int i = 1; i < kernel_size; ++i)
    {
        coord_val[i] = coord_val[i - 1] + 1;
    }

    const float kSigma = 0.5f;     // sigma,σ,Standard Deviation
    float val1 = 1 / (2 * M_PI * kSigma * kSigma);
    float val2 = -1 / (2 * kSigma * kSigma);
    for (int i = 0; i < kernel_size; ++i)
    {
        for (int j = 0; j < kernel_size; ++j)
        {
            kernel[i][j] = val1 * exp(val2 * (pow(coord_val[i], 2) + pow(coord_val[j], 2)));
        }
    }

    float sum = 0.0f;
    for (int i = 0; i < kernel_size; ++i)
    {
        sum += std::accumulate(kernel[i].begin(), kernel[i].end(), 0.0);
    }

    for (int i = 0; i < kernel_size; ++i)
    {
        for (int j = 0; j < kernel_size; ++j)
        {
            kernel[i][j] /= sum;
        }
    }

    return 0;
}

int GaussianBlur(const std::vector<std::vector<float>> &kernel, const cv::Mat &src, cv::Mat &formula_gaussian)
{
    int kernel_size = kernel.size();

    /* 对原图灰度图补pad */
    int pad = kernel_size >> 1;
    cv::Mat formula_pad = cv::Mat(src.rows + 2 * pad, src.cols + 2 * pad, CV_32FC1);
    memset(formula_pad.data, 0, formula_pad.rows * formula_pad.cols * sizeof(float));
    cv::Rect rect(pad, pad, src.cols, src.rows);
    src.copyTo(formula_pad(rect));

    /* 执行高斯滤波 */
    formula_gaussian = cv::Mat(src.rows, src.cols, CV_32FC1);
    memset(formula_gaussian.data, 0, formula_gaussian.rows * formula_gaussian.cols * sizeof(float));
    for (int i = 0; i < src.rows; ++i)
    {
        for (int j = 0; j < src.cols; ++j)
        {
            float tmp = 0.0f;
            for (int m = 0; m < kernel_size; ++m)
            {
                for (int n = 0; n < kernel_size; ++n)
                {
                    tmp += kernel[m][n] * formula_pad.at<float>(i + m, j + n);
                }
            }

            formula_gaussian.at<float>(i, j) = tmp;
        }
    }

    return 0;
}

int DoSobel(const cv::Mat &formula_gaussian, cv::Mat &sobel_xy, cv::Mat &sobel_angle, float &mean_sobel_xy)
{
    /* Sobel滤波 */
    const int kSobelKernel = 3;
    int pad = kSobelKernel >> 1;   // Sobel的kernel大小是3
    std::vector<std::vector<int> > kernel_x{ {-1, 0, 1}, {-2, 0, 2}, {-1, 0, 1} };
    std::vector<std::vector<int> > kernel_y{ {1, 2, 1}, {0, 0, 0}, {-1, -2, -1} };
    cv::Mat gaussian_pad(formula_gaussian.rows + 2 * pad, formula_gaussian.cols + 2 * pad, CV_32FC1);
    memset(gaussian_pad.data, 0, gaussian_pad.rows * gaussian_pad.cols * sizeof(float));
    cv::Rect rect = cv::Rect(pad, pad, formula_gaussian.cols, formula_gaussian.rows);
    formula_gaussian.copyTo(gaussian_pad(rect));

    cv::Mat sobel_x(formula_gaussian.rows, formula_gaussian.cols, CV_32FC1);
    cv::Mat sobel_y(formula_gaussian.rows, formula_gaussian.cols, CV_32FC1);
    sobel_xy = cv::Mat(formula_gaussian.rows, formula_gaussian.cols, CV_32FC1);
    sobel_angle = cv::Mat(formula_gaussian.rows, formula_gaussian.cols, CV_32FC1);  // 这里角度用tan值表示
    const float kEps = 0.000001f;
    float sum_sobel_xy = 0.0f;
    mean_sobel_xy = 0.0f;
    for (int i = 0; i < formula_gaussian.rows; ++i)
    {
        for (int j = 0; j < formula_gaussian.cols; ++j)
        {
            float gradient_x = 0.0f;
            float gradient_y = 0.0f;
            float gradient_xy = 0.0f;
            for (int m = 0; m < kSobelKernel; ++m)
            {
                for (int n = 0; n < kSobelKernel; ++n)
                {
                    gradient_x += kernel_x[m][n] * gaussian_pad.at<float>(i + m, j + n);
                    gradient_y += kernel_y[m][n] * gaussian_pad.at<float>(i + m, j + n);
                }
            }
            gradient_xy = sqrt(pow(gradient_x, 2) + pow(gradient_y, 2));
            sobel_x.at<float>(i, j) = gradient_x;
            sobel_y.at<float>(i, j) = gradient_y;
            sobel_xy.at<float>(i, j) = gradient_xy;
            sobel_angle.at<float>(i, j) = gradient_y / (fabs(gradient_x) > kEps ? gradient_x : kEps);

            sum_sobel_xy += gradient_xy;
        }
    }
    mean_sobel_xy = sum_sobel_xy / (sobel_xy.rows * sobel_xy.cols);

    return 0;
}

/* Efficient Non-Maximum Suppression */
// 在梯度图每个像素的八邻域内，求该像素的梯度正反两个方向对应的虚拟像素点的梯度,
// 如果该像素的梯度比两个虚拟像素的梯度大则认为是真边缘,否则是假边缘,舍弃。
// 八邻域eight-neighbor标识说明,en4标识中心点像素,其它8个是其邻域: 
//                               en0 en1 en2
//                               en3 en4 en5
//                               en6 en7 en8 
int DoNMS(const cv::Mat &sobel_xy, const cv::Mat &sobel_angle, cv::Mat &nms_sobel_xy)
{
    nms_sobel_xy = cv::Mat(sobel_xy.rows, sobel_xy.cols, CV_32FC1);
    memset(nms_sobel_xy.data, 0, nms_sobel_xy.rows * nms_sobel_xy.cols * sizeof(float));

    for (int i = 1; i < nms_sobel_xy.rows - 1; ++i)        // 图像第一行、最后一行不认为是边缘
    {
        for (int j = 1; j < nms_sobel_xy.cols - 1; ++j)    // 图像第一列、最后一列不认为是边缘
        {
            // 8邻域各像素的梯度
            float en0 = sobel_xy.at<float>(i - 1, j - 1);
            float en1 = sobel_xy.at<float>(i - 1, j);
            float en2 = sobel_xy.at<float>(i - 1, j + 1);
            float en3 = sobel_xy.at<float>(i, j - 1);
            float en4 = sobel_xy.at<float>(i, j);          // 当前像素梯度
            float en5 = sobel_xy.at<float>(i, j + 1);
            float en6 = sobel_xy.at<float>(i + 1, j - 1);
            float en7 = sobel_xy.at<float>(i + 1, j);
            float en8 = sobel_xy.at<float>(i + 1, j + 1);

            // 插值interpolation
            float grad_inter1 = 0.0f;
            float grad_inter2 = 0.0f;
            float ratio = 0.0f;   // 插值时的权重比例,[0.0, 1.0f]
            float angle = sobel_angle.at<float>(i, j);  // 注意角度使用dy/dx表示的,因此ratio是|angle|或1/|angle|
            if (angle >= 0)  // 第一、三象限
            {
                if (angle >= 1)         // [45°, 90°]使用en1、en2得到插值一，[225°, 270°]使用en6、en7得到插值二
                {
                    ratio = 1.0f / angle;  // angle越大越靠近Y轴,en1、en7的占比就越大,ratio越小,1-ratio越大;
                                           // 当angele=1时是45°或225°方向此时正好en2或en6的占比是1;
                                           // 当angle无穷大时是90°或270°方向此时正好en1、en7的占比是1
                    grad_inter1 = (1 - ratio) * en1 + ratio * en2;
                    grad_inter2 = (1 - ratio) * en7 + ratio * en6;
                }
                else  // angle∈[0,1)   // [0°, 45°)使用en2、en5得到插值一，[180°, 225°)使用en3、en6得到插值二
                {
                    ratio = angle;         // angle越小越靠近X轴,en5、en3的占比就越大,ratio越小,1-ratio越大;
                                           // 当angle=0时是0°或180°方向此时正好en5、en3的占比是1;
                                           // 当angle=1时是45°或225°方向此时正好en2、en6的占比是1
                    grad_inter1 = (1 - ratio) * en5 + ratio * en2;
                    grad_inter2 = (1 - ratio) * en3 + ratio * en6;
                }
            }
            else             // 第二、四象限
            {
                if (fabs(angle) >= 1)  // (90°, 135°]使用en0、en1得到插值一，(270°, 315°]使用en7、en8得到插值二
                {
                    ratio = 1.0f / fabs(angle);  // fabs(angle)越大越靠近Y轴,en1、en7的占比就越大,ratio越小,1-ratio越大;
                                                 // 当fabs(angle)=1时是135°或315°方向此时正好en0或en8的占比是1;
                                                 // 当fabs(angle)无穷大时是90°或270°方向此时正好en1、en7的占比是1
                    grad_inter1 = (1 - ratio) * en1 + ratio * en0;
                    grad_inter2 = (1 - ratio) * en7 + ratio * en8;
                }
                else  // fabs(angle)∈[0,1)   // (135°, 180°]使用en0、en3得到插值一，(315°, 360°]使用en5、en8得到插值二
                {
                    ratio = fabs(angle);         // fabs(angle)越小越靠近X轴,en3、en5的占比就越大,ratio越小,1-ratio越大;
                                                 // 当fabs(angle)=0时是180°或360°方向此时正好en3或en5的占比是1;
                                                 // 当fabs(angle)=1时是135°或315°方向此时正好en0、en8的占比是1
                    grad_inter1 = (1 - ratio) * en3 + ratio * en0;
                    grad_inter2 = (1 - ratio) * en5 + ratio * en8;
                }
            }
            if (en4 > grad_inter1 && en4 > grad_inter2)
            {
                nms_sobel_xy.at<float>(i, j) = en4;
            }
        }
    }

    return 0;
}

/* 双阈值检测
   如果边缘像素梯度值大于高阈值则认为是强边缘;
   如果边缘像素梯度值小于低阈值则不认为是边缘;
   如果边缘像素梯度介于高低阈值之间,则视为弱边缘,
   需要进一步判断若边缘的8邻域是否有强边缘,8邻域内只要有一个强边缘则认为该弱边缘是真实的边缘,否则不认为是边缘 */
int DoBinaryThresh(const cv::Mat &nms_sobel_xy, float mean_sobel_xy, cv::Mat &img_canny)
{
    cv::Mat binary_thresh_canny = cv::Mat(nms_sobel_xy.rows, nms_sobel_xy.cols, CV_32FC1);
    memset(binary_thresh_canny.data, 0, binary_thresh_canny.rows * binary_thresh_canny.cols * sizeof(float));

    float low_thresh = mean_sobel_xy * 0.5f;
    float high_thresh = low_thresh * 3.0f;
    for (int i = 1; i < binary_thresh_canny.rows - 1; ++i)      // 图像第一行、最后一行不认为是边缘
    {
        for (int j = 1; j < binary_thresh_canny.cols - 1; ++j)  // 图像第一列、最后一列不认为是边缘
        {
            float en0 = nms_sobel_xy.at<float>(i - 1, j - 1);
            float en1 = nms_sobel_xy.at<float>(i - 1, j);
            float en2 = nms_sobel_xy.at<float>(i - 1, j + 1);
            float en3 = nms_sobel_xy.at<float>(i, j - 1);
            float en4 = nms_sobel_xy.at<float>(i, j);              // 当前像素梯度
            float en5 = nms_sobel_xy.at<float>(i, j + 1);
            float en6 = nms_sobel_xy.at<float>(i + 1, j - 1);
            float en7 = nms_sobel_xy.at<float>(i + 1, j);
            float en8 = nms_sobel_xy.at<float>(i + 1, j + 1);
            if (en4 >= high_thresh)       // 强边缘
            {
                binary_thresh_canny.at<float>(i, j) = 255.0f;
            }
            else if (en4 <= low_thresh)  // 不是边缘
            {
                binary_thresh_canny.at<float>(i, j) = 0.0f;
            }
            else                        // 弱边缘,继续8邻域的判断
            {
                if (en0 >= high_thresh || en1 >= high_thresh || en2 >= high_thresh
                    || en3 >= high_thresh || en5 >= high_thresh
                    || en6 >= high_thresh || en7 >= high_thresh || en8 >= high_thresh)
                {
                    binary_thresh_canny.at<float>(i, j) = 255.0f;
                }
            }
        }
    }

    binary_thresh_canny.convertTo(img_canny, CV_8UC1);

    return 0;
}

Canny.cpp

#include "auxiliary.hpp"
#include <iostream>
#include <opencv2/highgui.hpp>

int main()
{
    //cv::Mat img = cv::imread("../data/carton1.jpg", 0);
    cv::Mat img = cv::imread("../data/推导2.jpg", 0);
    cv::Mat imgf;
    img.convertTo(imgf, CV_32FC1, 1 / 255.0);

    /* 计算高斯滤波器 */
    int kKernelSize = 5;   // 滤波器大小,可自定义,但应当是奇数
    std::vector<std::vector<float>> kernel;
    int ret = GaussianKernel(kKernelSize, kernel);
    if (0 != ret)
    {
        system("pause");
        return ret;
    }

    /* 高斯滤波 */
    cv::Mat formula_gaussian;
    ret = GaussianBlur(kernel, imgf, formula_gaussian);
    if (0 != ret)
    {
        system("pause");
        return ret;
    }

    cv::Mat sobel_xy, sobel_angle;
    float mean_sobel_xy = 0.0f;
    ret = DoSobel(formula_gaussian, sobel_xy, sobel_angle, mean_sobel_xy);
    if (0 != ret)
    {
        system("pause");
        return ret;
    }

    cv::Mat nms_sobel_xy;
    ret = DoNMS(sobel_xy, sobel_angle, nms_sobel_xy);
    if (0 != ret)
    {
        system("pause");
        return ret;
    }

    cv::Mat img_canny;
    ret = DoBinaryThresh(nms_sobel_xy, mean_sobel_xy, img_canny);
    if (0 != ret)
    {
        system("pause");
        return ret;
    }

    cv::imshow("img_canny", img_canny);
    cv::waitKey(0);

    return 0;
}

非极大值抑制时双线性插值示意图

图1 非极大值抑制时双线性插值示意图

实验图像及结果

以原图推导2.jpg为示例

图2 推导2.jpg

gaussian.jpg

图3 gaussian.jpg

sobel_x.jpg

图4 sobel_x.jpg

sobel_y.jpg

图5 sobel_y.jpg

sobel_xy.jpg

图6 sobel_xy.jpg

sobel_angle.jpg在NMS时插值计算梯度正反两个方向虚拟像素点梯度时使用。

图7 sobel_angle.jpg

nms_sobel_xy.jpg

图8 nms_sobel_xy.jpg

最终结果img_canny.jpg

图9 最终结果 img_canny.jpg

Reference

Canny算子中的非极大值抑制（Non-Maximum Suppression）分析