c++ - 如何访问 cv::canny 阈值的梯度幅值图像

Question

我正在尝试通过计算梯度幅值图像的一些统计数据来设置两个 Canny 阈值(这似乎比像许多人所做的那样在灰度图像上计算阈值(如 Otsu)更好，因为这些阈值与实际应用阈值的梯度幅值图像相关，但其值大不相同)。然而，计算的阈值需要从与 Canny 最终在内部进行阈值化的完全相同的梯度幅值图像计算，否则结果将不会如预期。也就是说，cv::canny 在内部进行一些平滑处理(其参数未公开)，应用 Sobel 运算符，执行快速或完整的梯度幅值计算等，然后应用用户在进行细化/链接等之前指定阈值。在计算我的统计数据之前，我必须在外部执行完全相同的步骤，这样我传递给 cv::canny 的阈值才真正有意义。

有没有办法访问算法内部使用的这张图片？

Answer 1

您不能直接获取 OpenCV Canny 函数的内部状态，但您可以提取 OpenCV 代码并制作您自己的函数。

这是一个自动选择 Canny 阈值的函数(基于 egonSchiele implementation )。

注意，在这个函数中:

• 会输出Sobel梯度sobel_x和sobel_y的结果，避免用Sobel重新计算想稍后使用图像渐变。 (如果不需要，您可以轻松地重构它)

• 此代码始终使用 L1 梯度来计算统计数据。然后根据输入参数使用 L1 或 L2 进行实际幅度计算。

• 这里的魔数(Magic Number)是固定的。您可以轻松地重构代码以将它们作为输入参数传递。这些神奇的数字是:

  • NUM_BINS:用于计算统计数据的直方图的 bin 数
  • percent_of_pixels_not_edges:估计较高的 Canny 阈值
  • threshold_ratio:恢复较低的 Canny 阈值。

---

关于 Otsu 在灰度图像上的使用以恢复 Canny 阈值...好吧，这对我来说意义不大，因为“灰度”图像和“梯度幅值”图像具有不同的语义和值范围。

---

代码:

#include<opencv2/opencv.hpp>
using namespace cv;


// Based on https://gist.github.com/egonSchiele/756833
void cvCanny3(const void* srcarr, void* dstarr,
    void* dxarr, void* dyarr,
    int aperture_size)
{
    cv::AutoBuffer<char> buffer;
    std::vector<uchar*> stack;
    uchar **stack_top = 0, **stack_bottom = 0;

    CvMat srcstub, *src = cvGetMat(srcarr, &srcstub);
    CvMat dststub, *dst = cvGetMat(dstarr, &dststub);

    CvMat dxstub, *dx = cvGetMat(dxarr, &dxstub);
    CvMat dystub, *dy = cvGetMat(dyarr, &dystub);


    CvSize size;
    int flags = aperture_size;
    int low, high;
    int* mag_buf[3];
    uchar* map;
    ptrdiff_t mapstep;
    int maxsize;
    int i, j;
    CvMat mag_row;

    if (CV_MAT_TYPE(src->type) != CV_8UC1 ||
        CV_MAT_TYPE(dst->type) != CV_8UC1 ||
        CV_MAT_TYPE(dx->type) != CV_16SC1 ||
        CV_MAT_TYPE(dy->type) != CV_16SC1)
        CV_Error(CV_StsUnsupportedFormat, "");

    if (!CV_ARE_SIZES_EQ(src, dst))
        CV_Error(CV_StsUnmatchedSizes, "");

    aperture_size &= INT_MAX;
    if ((aperture_size & 1) == 0 || aperture_size < 3 || aperture_size > 7)
        CV_Error(CV_StsBadFlag, "");


    size.width = src->cols;
    size.height = src->rows;

    //aperture_size = -1; //SCHARR
    cvSobel(src, dx, 1, 0, aperture_size);
    cvSobel(src, dy, 0, 1, aperture_size);


    //% Calculate Magnitude of Gradient
    //magGrad = hypot(dx, dy);

    Mat1f magGrad(size.height, size.width, 0.f);
    float maxGrad(0);
    float val(0);
    for (i = 0; i<size.height; ++i)
    {
        float* _pmag = magGrad.ptr<float>(i);
        const short* _dx = (short*)(dx->data.ptr + dx->step*i);
        const short* _dy = (short*)(dy->data.ptr + dy->step*i);
        for (j = 0; j<size.width; ++j)
        {
            val = float(abs(_dx[j]) + abs(_dy[j]));
            _pmag[j] = val;
            maxGrad = (val > maxGrad) ? val : maxGrad;
        }
    }

    //% Normalize for threshold selection
    //normalize(magGrad, magGrad, 0.0, 1.0, NORM_MINMAX);

    //% Determine Hysteresis Thresholds

    // -------------------------------------------------
    //% Set magic numbers
    const int NUM_BINS = 64;
    const double percent_of_pixels_not_edges = 0.9;
    const double threshold_ratio = 0.25;
    // -------------------------------------------------

    //% Compute histogram
    int bin_size = cvFloor(maxGrad / float(NUM_BINS) + 0.5f) + 1;
    if (bin_size < 1) bin_size = 1;
    int bins[NUM_BINS] = { 0 };
    for (i = 0; i<size.height; ++i)
    {
        float *_pmag = magGrad.ptr<float>(i);
        for (j = 0; j<size.width; ++j)
        {
            int hgf = int(_pmag[j]);
            bins[int(_pmag[j]) / bin_size]++;
        }
    }




    //% Select the thresholds
    float total(0.f);
    float target = float(size.height * size.width * percent_of_pixels_not_edges);
    int low_thresh, high_thresh(0);

    while (total < target)
    {
        total += bins[high_thresh];
        high_thresh++;
    }
    high_thresh *= bin_size;
    low_thresh = cvFloor(threshold_ratio * float(high_thresh));

    if (flags & CV_CANNY_L2_GRADIENT)
    {
        Cv32suf ul, uh;
        ul.f = (float)low_thresh;
        uh.f = (float)high_thresh;

        low = ul.i;
        high = uh.i;
    }
    else
    {
        low = cvFloor(low_thresh);
        high = cvFloor(high_thresh);
    }


    buffer.allocate((size.width + 2)*(size.height + 2) + (size.width + 2) * 3 * sizeof(int));
    mag_buf[0] = (int*)(char*)buffer;
    mag_buf[1] = mag_buf[0] + size.width + 2;
    mag_buf[2] = mag_buf[1] + size.width + 2;
    map = (uchar*)(mag_buf[2] + size.width + 2);
    mapstep = size.width + 2;

    maxsize = MAX(1 << 10, size.width*size.height / 10);
    stack.resize(maxsize);
    stack_top = stack_bottom = &stack[0];

    memset(mag_buf[0], 0, (size.width + 2)*sizeof(int));
    memset(map, 1, mapstep);
    memset(map + mapstep*(size.height + 1), 1, mapstep);

    /* sector numbers
    (Top-Left Origin)

    1   2   3
    *  *  *
    * * *
    0*******0
    * * *
    *  *  *
    3   2   1
    */

#define CANNY_PUSH(d)    *(d) = (uchar)2, *stack_top++ = (d)
#define CANNY_POP(d)     (d) = *--stack_top

    mag_row = cvMat(1, size.width, CV_32F);

    // calculate magnitude and angle of gradient, perform non-maxima supression.
    // fill the map with one of the following values:
    //   0 - the pixel might belong to an edge
    //   1 - the pixel can not belong to an edge
    //   2 - the pixel does belong to an edge
    for (i = 0; i <= size.height; i++)
    {
        int* _mag = mag_buf[(i > 0) + 1] + 1;
        float* _magf = (float*)_mag;
        const short* _dx = (short*)(dx->data.ptr + dx->step*i);
        const short* _dy = (short*)(dy->data.ptr + dy->step*i);
        uchar* _map;
        int x, y;
        ptrdiff_t magstep1, magstep2;
        int prev_flag = 0;

        if (i < size.height)
        {
            _mag[-1] = _mag[size.width] = 0;

            if (!(flags & CV_CANNY_L2_GRADIENT))
                for (j = 0; j < size.width; j++)
                    _mag[j] = abs(_dx[j]) + abs(_dy[j]);

            else
            {
                for (j = 0; j < size.width; j++)
                {
                    x = _dx[j]; y = _dy[j];
                    _magf[j] = (float)std::sqrt((double)x*x + (double)y*y);
                }
            }
        }
        else
            memset(_mag - 1, 0, (size.width + 2)*sizeof(int));

        // at the very beginning we do not have a complete ring
        // buffer of 3 magnitude rows for non-maxima suppression
        if (i == 0)
            continue;

        _map = map + mapstep*i + 1;
        _map[-1] = _map[size.width] = 1;

        _mag = mag_buf[1] + 1; // take the central row
        _dx = (short*)(dx->data.ptr + dx->step*(i - 1));
        _dy = (short*)(dy->data.ptr + dy->step*(i - 1));

        magstep1 = mag_buf[2] - mag_buf[1];
        magstep2 = mag_buf[0] - mag_buf[1];

        if ((stack_top - stack_bottom) + size.width > maxsize)
        {
            int sz = (int)(stack_top - stack_bottom);
            maxsize = MAX(maxsize * 3 / 2, maxsize + 8);
            stack.resize(maxsize);
            stack_bottom = &stack[0];
            stack_top = stack_bottom + sz;
        }

        for (j = 0; j < size.width; j++)
        {
#define CANNY_SHIFT 15
#define TG22  (int)(0.4142135623730950488016887242097*(1<<CANNY_SHIFT) + 0.5)

            x = _dx[j];
            y = _dy[j];
            int s = x ^ y;
            int m = _mag[j];

            x = abs(x);
            y = abs(y);
            if (m > low)
            {
                int tg22x = x * TG22;
                int tg67x = tg22x + ((x + x) << CANNY_SHIFT);

                y <<= CANNY_SHIFT;

                if (y < tg22x)
                {
                    if (m > _mag[j - 1] && m >= _mag[j + 1])
                    {
                        if (m > high && !prev_flag && _map[j - mapstep] != 2)
                        {
                            CANNY_PUSH(_map + j);
                            prev_flag = 1;
                        }
                        else
                            _map[j] = (uchar)0;
                        continue;
                    }
                }
                else if (y > tg67x)
                {
                    if (m > _mag[j + magstep2] && m >= _mag[j + magstep1])
                    {
                        if (m > high && !prev_flag && _map[j - mapstep] != 2)
                        {
                            CANNY_PUSH(_map + j);
                            prev_flag = 1;
                        }
                        else
                            _map[j] = (uchar)0;
                        continue;
                    }
                }
                else
                {
                    s = s < 0 ? -1 : 1;
                    if (m > _mag[j + magstep2 - s] && m > _mag[j + magstep1 + s])
                    {
                        if (m > high && !prev_flag && _map[j - mapstep] != 2)
                        {
                            CANNY_PUSH(_map + j);
                            prev_flag = 1;
                        }
                        else
                            _map[j] = (uchar)0;
                        continue;
                    }
                }
            }
            prev_flag = 0;
            _map[j] = (uchar)1;
        }

        // scroll the ring buffer
        _mag = mag_buf[0];
        mag_buf[0] = mag_buf[1];
        mag_buf[1] = mag_buf[2];
        mag_buf[2] = _mag;
    }

    // now track the edges (hysteresis thresholding)
    while (stack_top > stack_bottom)
    {
        uchar* m;
        if ((stack_top - stack_bottom) + 8 > maxsize)
        {
            int sz = (int)(stack_top - stack_bottom);
            maxsize = MAX(maxsize * 3 / 2, maxsize + 8);
            stack.resize(maxsize);
            stack_bottom = &stack[0];
            stack_top = stack_bottom + sz;
        }

        CANNY_POP(m);

        if (!m[-1])
            CANNY_PUSH(m - 1);
        if (!m[1])
            CANNY_PUSH(m + 1);
        if (!m[-mapstep - 1])
            CANNY_PUSH(m - mapstep - 1);
        if (!m[-mapstep])
            CANNY_PUSH(m - mapstep);
        if (!m[-mapstep + 1])
            CANNY_PUSH(m - mapstep + 1);
        if (!m[mapstep - 1])
            CANNY_PUSH(m + mapstep - 1);
        if (!m[mapstep])
            CANNY_PUSH(m + mapstep);
        if (!m[mapstep + 1])
            CANNY_PUSH(m + mapstep + 1);
    }

    // the final pass, form the final image
    for (i = 0; i < size.height; i++)
    {
        const uchar* _map = map + mapstep*(i + 1) + 1;
        uchar* _dst = dst->data.ptr + dst->step*i;

        for (j = 0; j < size.width; j++)
        {
            _dst[j] = (uchar)-(_map[j] >> 1);
        }
    }
};

void Canny3(InputArray image, OutputArray _edges,
    OutputArray _sobel_x, OutputArray _sobel_y,
    int apertureSize = 3, bool L2gradient = false)
{
    Mat src = image.getMat();
    _edges.create(src.size(), CV_8U);
    _sobel_x.create(src.size(), CV_16S);
    _sobel_y.create(src.size(), CV_16S);


    CvMat c_src = src, c_dst = _edges.getMat();
    CvMat c_dx = _sobel_x.getMat();
    CvMat c_dy = _sobel_y.getMat();


    cvCanny3(&c_src, &c_dst,
        &c_dx, &c_dy,
        apertureSize + (L2gradient ? CV_CANNY_L2_GRADIENT : 0));
};

int main()
{
    Mat3b img = imread("path_to_image");
    Mat1b gray;
    cvtColor(img, gray, COLOR_BGR2GRAY);

    Mat1b edges;
    Mat1s sobel_x, sobel_y;
    Canny3(gray, edges, sobel_x, sobel_y);

    imshow("edges", edges);
    waitKey();

    return 0;
}