opencv-12-gaussian filter-bilateral filter (with C + + code implementation)

Before you start

In recent days, I haven't written for my own reasons. One is that I'm lazy. The other is that I don't want to write down some problems I have encountered here. Let's try to finish this chapter first. Before that, I introduced the more commonly used and well understood mean and median filtering. However, in the routine Smoothing Images There are other filtering methods, mainly Gaussian filtering and bilateral filtering,

We've finished filtering and smoothing this time. We've written a little more, but we don't want to write any more. Come on

Catalog

Objectives of this paper

This paper mainly introduces

  • Gauss filter
  • Bilateral filtering

As before, we still introduce the principle, give the specific form, then use opencv to carry out the implementation process, and finally use our previous images to test the algorithm analysis and summary

text

Gaussian filter

We introduced that the median filter is the result of statistical sorting and belongs to non-linear result. The mean filter is operated by using template kernel. We also mentioned that the weight of mean filter must be considered in the calculation process of mean filter, and then proposed the operation of weighted mean filter, such as the most common operation kernel of weighted mean filter

\[M = \frac{1}{16} \left [ \begin{array}{c} 1 & 2 & 1 \\ 2& 4 & 2 \\ 1 & 2 & 1 \end{array} \right ] \]

In fact, this kernel is the result value of discrete rounding of Gaussian filter in 3 * 3 window. The most obvious feature is that the coefficient of template changes with the distance from the center of template, which can effectively suppress noise and smooth image. Compared with mean filter, it can smooth image better and retain image edge

Gauss filter principle

Since our image is two-dimensional, but the Gaussian distribution is one-dimensional, let's first consider the one-dimensional Gaussian distribution, which is the normal distribution curve commonly used by us,

\[G(x) = \frac{1}{\sqrt{2\pi \sigma}} e^{-\frac{x^2}{2\sigma^2}} \]

For two-dimensional Gaussian distribution, we can consider the result of superposition of two direction operations

\[G(x,y) = \frac{1}{2\pi \sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}} = G(x)*G(y) \]

Considering that the image calculation is actually a discrete coordinate, for the template whose window size is \ ((2k + 1) \times (2k + 1)), we can express it as

\[G{i,j} = \frac{1}{2\pi \sigma ^ 2}e ^{-\frac{(i - k - 1)^2 + (j - k - 1)^2}{2 \sigma ^ 2}} \]

You can refer to Fundamentals of image processing (4): details of Gaussian filter
The method given in it, using

void generateGaussianTemplate(double window[][11], int ksize, double sigma)
{
    static const double pi = 3.1415926;
    int center = ksize / 2; // The center position of the template, that is, the origin of the coordinates
    double x2, y2;
    for (int i = 0; i < ksize; i++)
    {
        x2 = pow(i - center, 2);
        for (int j = 0; j < ksize; j++)
        {
            y2 = pow(j - center, 2);
            double g = exp(-(x2 + y2) / (2 * sigma * sigma));
            g /= 2 * pi * sigma*sigma;	// 
            window[i][j] = g;
        }
    }
    double k = 1 / window[0][0]; // Normalize the coefficient of the upper left corner to 1
    for (int i = 0; i < ksize; i++)
    {
        for (int j = 0; j < ksize; j++)
        {
            window[i][j] *= k;
        }
    }
}

The Gaussian template of \ (3 \ times 3, sigma = 0.8 \) is generated, and the corresponding rounding is obtained

\[M = \frac{1}{16} \left [ \begin{array}{c} 1 & 2 & 1 \\ 2& 4 & 2 \\ 1 & 2 & 1 \end{array} \right ] \]

The above article also introduces the significance of \ (\ sigma \) in statistics in detail, which can be referred to for learning
But from high school, we can see the curve of normal distribution

C + + implementation

As we mentioned earlier Fundamentals of image processing (4): details of Gaussian filter Here is the code implementation based on opencv. Here is the algorithm implementation of \ (O(m*n*k^2) \)

// Source: https://www.cnblogs.com/wangguchangqing/p/6407717.html
void GaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)
{
    CV_Assert(src.channels() || src.channels() == 3); // Only process single channel or three channel images
    const static double pi = 3.1415926;
    // Generating Gaussian filter template based on window size and sigma
    // Apply for a two-dimensional array to store the generated Gaussian template matrix
    double **templateMatrix = new double*[ksize];
    for (int i = 0; i < ksize; i++)
        templateMatrix[i] = new double[ksize];
    int origin = ksize / 2; // Take the center of the template as the origin
    double x2, y2;
    double sum = 0;
    for (int i = 0; i < ksize; i++)
    {
        x2 = pow(i - origin, 2);
        for (int j = 0; j < ksize; j++)
        {
            y2 = pow(j - origin, 2);
            // The constants before Gaussian function can be eliminated in the process of normalization without calculation
            double g = exp(-(x2 + y2) / (2 * sigma * sigma));
            sum += g;
            templateMatrix[i][j] = g;
        }
    }
    for (int i = 0; i < ksize; i++)
    {
        for (int j = 0; j < ksize; j++)
        {
            templateMatrix[i][j] /= sum;
            cout << templateMatrix[i][j] << " ";
        }
        cout << endl;
    }
    // Apply template to image
    int border = ksize / 2;
    copyMakeBorder(src, dst, border, border, border, border, BorderTypes::BORDER_REFLECT);
    int channels = dst.channels();
    int rows = dst.rows - border;
    int cols = dst.cols - border;
    for (int i = border; i < rows; i++)
    {
        for (int j = border; j < cols; j++)
        {
            double sum[3] = { 0 };
            for (int a = -border; a <= border; a++)
            {
                for (int b = -border; b <= border; b++)
                {
                    if (channels == 1)
                    {
                        sum[0] += templateMatrix[border + a][border + b] * dst.at<uchar>(i + a, j + b);
                    }
                    else if (channels == 3)
                    {
                        Vec3b rgb = dst.at<Vec3b>(i + a, j + b);
                        auto k = templateMatrix[border + a][border + b];
                        sum[0] += k * rgb[0];
                        sum[1] += k * rgb[1];
                        sum[2] += k * rgb[2];
                    }
                }
            }
            for (int k = 0; k < channels; k++)
            {
                if (sum[k] < 0)
                    sum[k] = 0;
                else if (sum[k] > 255)
                    sum[k] = 255;
            }
            if (channels == 1)
                dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
            else if (channels == 3)
            {
                Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
                dst.at<Vec3b>(i, j) = rgb;
            }
        }
    }
    // Release template array
    for (int i = 0; i < ksize; i++)
        delete[] templateMatrix[i];
    delete[] templateMatrix;
}

Then, the implementation of separation is also given. After the horizontal operation, the vertical operation is carried out, and the calculation time will be improved

// Source: https://www.cnblogs.com/wangguchangqing/p/6407717.html
// Calculation of separation
void separateGaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)
{
    CV_Assert(src.channels()==1 || src.channels() == 3); // Only process single channel or three channel images
    // Generating one-dimensional Gaussian filter template
    double *matrix = new double[ksize];
    double sum = 0;
    int origin = ksize / 2;
    for (int i = 0; i < ksize; i++)
    {
        // The constants before Gaussian function can be eliminated in the process of normalization without calculation
        double g = exp(-(i - origin) * (i - origin) / (2 * sigma * sigma));
        sum += g;
        matrix[i] = g;
    }
    // normalization
    for (int i = 0; i < ksize; i++)
        matrix[i] /= sum;
    // Apply template to image
    int border = ksize / 2;
    copyMakeBorder(src, dst, border, border, border, border, BorderTypes::BORDER_REFLECT);
    int channels = dst.channels();
    int rows = dst.rows - border;
    int cols = dst.cols - border;
    // horizontal direction
    for (int i = border; i < rows; i++)
    {
        for (int j = border; j < cols; j++)
        {
            double sum[3] = { 0 };
            for (int k = -border; k <= border; k++)
            {
                if (channels == 1)
                {
                    sum[0] += matrix[border + k] * dst.at<uchar>(i, j + k); // Row invariable, column changing; do horizontal convolution first
                }
                else if (channels == 3)
                {
                    Vec3b rgb = dst.at<Vec3b>(i, j + k);
                    sum[0] += matrix[border + k] * rgb[0];
                    sum[1] += matrix[border + k] * rgb[1];
                    sum[2] += matrix[border + k] * rgb[2];
                }
            }
            for (int k = 0; k < channels; k++)
            {
                if (sum[k] < 0)
                    sum[k] = 0;
                else if (sum[k] > 255)
                    sum[k] = 255;
            }
            if (channels == 1)
                dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
            else if (channels == 3)
            {
                Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
                dst.at<Vec3b>(i, j) = rgb;
            }
        }
    }
    // Vertical direction
    for (int i = border; i < rows; i++)
    {
        for (int j = border; j < cols; j++)
        {
            double sum[3] = { 0 };
            for (int k = -border; k <= border; k++)
            {
                if (channels == 1)
                {
                    sum[0] += matrix[border + k] * dst.at<uchar>(i + k, j); // Column invariance, row variation; vertical convolution
                }
                else if (channels == 3)
                {
                    Vec3b rgb = dst.at<Vec3b>(i + k, j);
                    sum[0] += matrix[border + k] * rgb[0];
                    sum[1] += matrix[border + k] * rgb[1];
                    sum[2] += matrix[border + k] * rgb[2];
                }
            }
            for (int k = 0; k < channels; k++)
            {
                if (sum[k] < 0)
                    sum[k] = 0;
                else if (sum[k] > 255)
                    sum[k] = 255;
            }
            if (channels == 1)
                dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
            else if (channels == 3)
            {
                Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
                dst.at<Vec3b>(i, j) = rgb;
            }
        }
    }
    delete[] matrix;
}

The algorithms here are all https://www.cnblogs.com/wangguchangqing/p/6407717.html. For details, please refer to the content

opencv Gaussian filtering

In fact, this article Image processing -- Gaussian filtering Well written
In fact, the main structure is the process given by him

In fact, the whole process of Gaussian filtering is to create a Gaussian kernel, and then use the filter2D method for filtering operation. If you want to go deep, you can see the function call graph, which is the same idea for implementation, very simple operation. We will test the effect later

// /modules\imgproc\src\smooth.dispatch.cpp:600
void GaussianBlur(InputArray _src, OutputArray _dst, Size ksize,
                  double sigma1, double sigma2,
                  int borderType)
  • src ‪ input image
  • dst output image
  • Odd size of ksize core
  • Sigma value in the direction of sigma x
  • sigma value in the direction of sign ‪ y
  • borderType boundary processing method

Gauss filter effect comparison

We still use the high salt and pepper noise image before, and then directly filter the algorithm, the calculation result is good, very similar to the previous test image, here

The four pictures here correspond to the high noise image, the result of direct Gaussian filtering, the result of separating xy direction for filtering, and the effect picture of opencv's own Gaussian filtering. Here is the preview image. The actual detection result is the parameter value given above. In fact, the effect can only be said to be general. We will compare the algorithm level later

image-noise: psnr:19.4727, mssim: B:0.353134 G:0.383638 R:0.629353
image-noise: psnr:26.3155, mssim: B:0.584585 G:0.617172 R:0.812303
image-noise: psnr:26.1721, mssim: B:0.574719 G:0.607494 R:0.809844
image-noise: psnr:26.4206, mssim: B:0.598176 G:0.630657 R:0.819658

Bilateral filter

Principle of bilateral filtering

The principle of Gaussian filtering is to give different weights for different distance from the center of the template, while the bilateral filtering not only considers the spatial distance of the image, but also the gray distance. The higher the weight of the point close to the middle gray value, the smaller the weight of the point with large gray value difference

The principle of bilateral filtering can be referred to Detailed explanation of Bilateral Filter,
You can refer to Bilateral Filtering for Gray and Color Images

In the article Fundamentals of image processing (5): bilateral filter Bilateral filtering is introduced in detail
In fact, it is consistent with the above filtering demonstration, which is to filter the noise under the condition of ensuring the image edge information

You can refer to Common understanding of bilateral filter The mathematical expression of bilateral filtering given

\[g(x,y) = \frac{\sum_{kl}f(k,l)w(i,j,k,l)}{\sum_{kl}w(i,j,k,l)} \]

For different template coefficients, there are two parts, mainly spatial domain template weight \ (w_d \) and gray domain template weight \ (w_r \),

\[\begin{array}{rl} w_d(i,j,k,l) &= e^{-\frac{(i-k)^2 +(j-l)^2}{2\sigma_d^2}} \\ w_r(i,j,k,l) &= e^{-\frac{\left \| f(i,j) - f(k,l) \right \|} {2\sigma_r^2}} \\ w &= w_d * w_r \end{array} \]

Where \ (q(i,j) \) is the coordinate of other coefficients of the template window, \ (f(i,j) \) is the pixel value of the image at point \ (q(i,j) \); \ (p(k,l) \) is the central coordinate point of the template window, and the corresponding pixel value is \ (f(k,l) \); \ (\ sigma_r \) is the standard deviation of the Gaussian function.

C + + implementation of bilateral filtering

I think it's very well written here Fundamentals of image processing (5): bilateral filter , manual implementation of bilateral filtering, we can refer to in detail, here

// Reference source: https://www.cnblogs.com/wangguchangqing/p/6416401.html
void myBilateralFilter(const Mat &src, Mat &dst, int ksize, double space_sigma, double color_sigma)
{
    int channels = src.channels();
    CV_Assert(channels == 1 || channels == 3);
    double space_coeff = -0.5 / (space_sigma * space_sigma);
    double color_coeff = -0.5 / (color_sigma * color_sigma);
    int radius = ksize / 2;
    Mat temp;
    copyMakeBorder(src, temp, radius, radius, radius, radius, BorderTypes::BORDER_REFLECT);
    vector<double> _color_weight(channels * 256); // Square of storage difference
    vector<double> _space_weight(ksize * ksize); // Space template coefficient
    vector<int> _space_ofs(ksize * ksize); // Coordinates of the template window
    double *color_weight = &_color_weight[0];
    double *space_weight = &_space_weight[0];
    int    *space_ofs = &_space_ofs[0];
    for (int i = 0; i < channels * 256; i++)
        color_weight[i] = exp(i * i * color_coeff);
    // Generate space template
    int maxk = 0;
    for (int i = -radius; i <= radius; i++)
    {
        for (int j = -radius; j <= radius; j++)
        {
            double r = sqrt(i*i + j * j);
            if (r > radius)
                continue;
            space_weight[maxk] = exp(r * r * space_coeff); // Coefficient of template storage
            space_ofs[maxk++] = i * temp.step + j * channels; // The position where the formwork is stored corresponds to the formwork coefficient
        }
    }
    // Filtering process
    for (int i = 0; i < src.rows; i++)
    {
        const uchar *sptr = temp.data + (i + radius) * temp.step + radius * channels;
        uchar *dptr = dst.data + i * dst.step;
        if (channels == 1)
        {
            for (int j = 0; j < src.cols; j++)
            {
                double sum = 0, wsum = 0;
                int val0 = sptr[j]; // Pixels in the center of the template
                for (int k = 0; k < maxk; k++)
                {
                    int val = sptr[j + space_ofs[k]];
                    double w = space_weight[k] * color_weight[abs(val - val0)]; // Template coefficient = space coefficient * gray value coefficient
                    sum += val * w;
                    wsum += w;
                }
                dptr[j] = (uchar)cvRound(sum / wsum);
            }
        }
        else if (channels == 3)
        {
            for (int j = 0; j < src.cols * 3; j+=3)
            {
                double sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
                int b0 = sptr[j];
                int g0 = sptr[j + 1];
                int r0 = sptr[j + 2];
                for (int k = 0; k < maxk; k++)
                {
                    const uchar *sptr_k = sptr + j + space_ofs[k];
                    int b = sptr_k[0];
                    int g = sptr_k[1];
                    int r = sptr_k[2];
                    double w = space_weight[k] * color_weight[abs(b - b0) + abs(g - g0) + abs(r - r0)];
                    sum_b += b * w;
                    sum_g += g * w;
                    sum_r += r * w;
                    wsum += w;
                }
                wsum = 1.0f / wsum;
                b0 = cvRound(sum_b * wsum);
                g0 = cvRound(sum_g * wsum);
                r0 = cvRound(sum_r * wsum);
                dptr[j] = (uchar)b0;
                dptr[j + 1] = (uchar)g0;
                dptr[j + 2] = (uchar)r0;
            }
        }
    }
}

Implementation of bilateral filtering with opencv

void bilateralFilter( InputArray _src, OutputArray _dst, int d,
                      double sigmaColor, double sigmaSpace,
                      int borderType )
  • InputArray src: input image. It can be of Mat type. The image must be an 8-bit or floating-point single channel or three channel image.
  • OutputArray dst: the output image has the same size and type as the original image.
  • int d: represents the diameter range of each pixel neighborhood during filtering. If the value is not a positive number, the function evaluates the value from the fifth parameter, sigma space.
  • double sigmaColor: the sigma value of the color space filter. The larger the value of this parameter is, the wider the color in the neighborhood of the pixel will be mixed together, resulting in a larger semi equal color region. (this parameter can be understood as value domain kernel)
  • Double sigma space: the Sigma value of the filter in the coordinate space. If the value is large, it means that the farther the pixels will affect each other, so that the more similar colors in the larger area can get the same color. When d > 0, D specifies the neighborhood size and has nothing to do with sigma space, otherwise D is proportional to sigma space
  • int borderType=BORDER_DEFAULT: a border mode used to infer pixels outside an image, with a default value of BORDER_DEFAULT

Comparison of bilateral filtering algorithms

At the beginning, I couldn't understand the bilateral filtering, and I didn't know what the purpose of doing this was, and what the final result represented. We tested our image according to the previous method, and the result is really the worst of several algorithms, but it just doesn't apply to our image results. In the actual use process, we still need to test before we can draw a conclusion

The test results are as follows: 3 windows are used for the results of original graph, manual implementation and opencv, and the value of sigma is 255
, this article https://blog.csdn.net/Jfuck/article/details/8932978 is very good. It introduces the influence of parameters on filtering, which can be learned

image-noise: psnr:19.4727, mssim: B:0.353134 G:0.383638 R:0.629353
image-noise: psnr:24.4502, mssim: B:0.538774 G:0.570666 R:0.776195
image-noise: psnr:24.4691, mssim: B:0.539177 G:0.571087 R:0.776461

summary

In fact, the use of bilateral filtering by individuals is not so much. I only learned a lot when I studied the guided filtering before. The writing here is poor. I can only barely read it. It is highly recommended https://www.cnblogs.com/wangguchangqing/category/740760.html This series, will be very detailed, many are the contents of the blog, you can refer to learning, Gaussian filtering is relatively simple, in fact, the complex filtering process is mainly to understand the algorithm, and then carry out the code implementation process according to the idea of the algorithm, and finally do some program optimization, understanding first, implementation second... Hope to give readers a little inspiration

I didn't prepare to write so much at the beginning, but the more I wrote, the more I couldn't stop. Many places I couldn't say I knew very well, this time I also had a deep understanding, but it was still very rigid and could only be used. Here I would like to recommend looking at the links I gave or looking up the relevant content myself. I just give a brief introduction here, If there is any mistake, please name it. Thank you very much

Reference link

  1. Fast Gauss filtering, Gauss blur, Gauss smoothing (two-dimensional convolution is divided into one-dimensional convolution step by step) AI Qingchengshan Xiaohe CSDN blog. See may 10, 2020 https://blog.csdn.net/qq_36359022/article/details/80188873.
  2. Bilateral filtering - Qitingshe Blog, May 10, 2020 Https://qitinghe.github.io/2018/06/14/bilateral filtering/.
  3. Detailed explanation of Bilateral Filter (CSDN blog), a column of artificial intelligence (Jfuck). See may 10, 2020 https://blog.csdn.net/Jfuck/article/details/8932978.
  4. Bilateral filter. Income Wikipedia, free encyclopedia, November 16, 2019 https://zh.wikipedia.org/w/index.php?title = bilateral filter & oldid = 56898678.
  5. CSDN blog, image processing - Gauss filtering, network, L-inYi, May 10, 2020 https://blog.csdn.net/L_inYi/article/details/8915116.
  6. Fundamentals of image processing (4): Gauss filter details - brook'icv - blog park. See may 10, 2020 https://www.cnblogs.com/wangguchangqing/p/6407717.html.
  7. Fundamentals of image processing (5): bilateral filters - brook'icv - blog park. 10 may 2020 https://www.cnblogs.com/wangguchangqing/p/6416401.html.
  8. Image processing - linear filtering - 3 Gauss filter - Tony Ma - blog park. 10 may 2020 https://www.cnblogs.com/pegasus/archive/2011/05/20/2052031.html.
  9. [transform] Gauss image filtering principle and its programming discretization implementation method (see Sina blog, May 10, 2020) http://blog.sina.com.cn/s/blog_640577ed0100yz8v.html.
  10. Common understanding of bilateral filter: CSDN blog of Pan Jinquan. See may 10, 2020 https://blog.csdn.net/guyuealian/article/details/82660826.
  11. Bilateral Filtering, May 10, 2020 http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/MANDUCHI1/Bilateral_Filtering.html.
  12. Image Processing, Computer Vision, image processing, OpenCV China, May 10, 2020 http://wiki.opencv.org.cn/index.php/Cv image processing.
  13. The principle, process, implementation and effect of o(1) complexity bilateral filtering algorithm. -Cloud + community - Tencent cloud. See may 10, 2020 https://cloud.tencent.com/developer/article/1011738.

This article is based on the platform of blog one article multiple sending OpenWrite release!

Tags: C++ OpenCV PHP Windows github

Posted on Sun, 10 May 2020 05:52:58 -0400 by stunna671