Gaussian filter is a linear filter, which can effectively suppress noise and smooth image. Its working principle is similar to the mean filter, which takes the mean value of pixels in the filter window as the output. The coefficients of the window template are different from those of the mean filter, and the template coefficients of the mean filter are the same, which is 1; The template coefficient of Gaussian filter decreases with the increase of distance from the template center. Therefore, Gaussian filter is less blurred than mean filter.

What is Gaussian filter

Since the name is Gaussian filter, it has a certain relationship with Gaussian distribution (normal distribution). A two-dimensional Gaussian function is as follows:

Where (x, y) (x, y) is the point coordinate, which can be considered as an integer in image processing; σσ Is the standard deviation. To get a Gaussian filter template, the Gaussian function can be discretized, and the obtained Gaussian function value is used as the coefficient of the template. For example: to generate a 3 × thirty-three × The Gaussian filter template of 3 takes the central position of the template as the coordinate origin for sampling. The coordinates of the template at each position are as follows (the x-axis is horizontal to the right and the y-axis is vertical to the down)

In this way, the coordinates of each position are brought into the Gaussian function, and the value obtained is the coefficient of the template.

The size of the window template is (2k + 1) × (2k + 1), the calculation formula of each element value in the template is as follows:

The template calculated in this way has two forms: decimal and integer.

The decimal template is the value calculated directly without any processing;

In integer form, normalization is required to normalize the value in the upper left corner of the template to 1, which will be described in detail below. When using integer templates, you need to add a coefficient in front of the template, which is the reciprocal of the sum of template coefficients.

Generation of Gaussian template

Knowing the principle of template generation is not difficult to implement

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;=””>

{

x2 = pow(i – center, 2);

for (int j = 0; j < ksize;=””>

{

y2 = pow(j – center, 2);

double g = exp(-(x2 + y2) / (2 * sigma * sigma));

g /= 2 * pi * sigma;

window[i][j] = g;

}

}

double k = 1 / window[0][0]; // Normalize the coefficient in the upper left corner to 1

for (int i = 0; i < ksize;=””>

{

for (int j = 0; j < ksize;=””>

{

window[i][j] *= k;

}

}

}

A two-dimensional array is required to store the generated coefficients (here, it is assumed that the maximum size of the template will not exceed 11); The second parameter is the size of the template (no more than 11); The third parameter is the standard deviation of Gaussian distribution.

In the generation process, first find the central position of the template ksize / 2 according to the size of the template. Then it is traversal, and the value of each coefficient in the template is calculated according to the function of Gaussian distribution.

It should be noted that in the final normalization process, the reciprocal of the coefficient in the upper left corner of the template is used as the normalized coefficient (the coefficient value in the upper left corner is normalized to 1), each coefficient in the template is multiplied by the value (the reciprocal of the coefficient in the upper left corner), and then the obtained value is rounded to obtain the integer Gaussian filter template.

The following screenshot is generated with a size of 3 × 3, σ= zero point eight three × 3, σ= 0.8 template

After rounding the above solution results, the following template is obtained:

This template is familiar. It is based on σ= 0.8 Gaussian function generated template.

The generation of decimal form is also relatively simple. The normalization process is removed, and after the solution process, each coefficient of the template is divided by the sum of all coefficients. The specific codes are as follows:

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;

double sum = 0;

for (int i = 0; i < ksize;=””>

{

x2 = pow(i – center, 2);

for (int j = 0; j < ksize;=””>

{

y2 = pow(j – center, 2);

double g = exp(-(x2 + y2) / (2 * sigma * sigma));

g /= 2 * pi * sigma;

sum += g;

window[i][j] = g;

}

}

//double k = 1 / window[0][0]; // Normalize the coefficient in the upper left corner to 1

for (int i = 0; i < ksize;=””>

{

for (int j = 0; j < ksize;=””>

{

window[i][j] /= sum;

}

}

}

three × 3, σ= 0.8 decimal template.

σσ Significance and selection of value

Through the above implementation process, it is not difficult to find that the most important parameter for the generation of Gaussian filter template is the standard deviation of Gaussian distribution σσ。 The standard deviation represents the degree of dispersion of the data, if σσ Smaller, the center coefficient of the generated template is larger, while the surrounding coefficient is smaller, so the smoothing effect on the image is not very obvious; conversely, σσ If it is large, the difference between the coefficients of the generated template is not very large. Compared with the mean template, the smoothing effect of the image is obvious.

Let’s take a look at the probability distribution density diagram of one-dimensional Gaussian distribution:

The horizontal axis represents the possible value x, and the vertical axis represents the probability distribution density f (x). It is not difficult to understand that the graphic area surrounded by such a curve and the X axis is 1. σσ (standard deviation) determines the width of this figure. It can be concluded that: σσ The larger the, the wider the figure, the smaller the peak, and the figure is more gentle; σσ The smaller, the narrower and more concentrated the figure, the sharper the middle part, and the figure changes violently. In fact, it is easy to understand that if the sigma, that is, the larger the standard deviation, it means that the density distribution must be relatively dispersed. Since the area is 1, the peak part decreases and the width is wider (the more dispersed the distribution is); Similarly, when σσ The smaller the, the more concentrated the density distribution, so the sharper the peak, the narrower the width!

The following conclusions can be drawn:

σσ The larger the, the more dispersed the distribution, and the proportion of each part has little difference, so the value of each element of the generated template has little difference, which is similar to the average template;

σσ The smaller the, the more concentrated the distribution, and the proportion of the middle part is much higher than that of other parts. Reflected in the Gaussian template, the value of the central element is much greater than that of other elements, so it is naturally equivalent to the operation of the middle value point.

Implementation based on OpenCV

When generating Gaussian template, its simple implementation is no different from other spatial filters. The specific code is as follows:

void GaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)

{

CV_ Assert(src.channels() || src.channels() == 3); // Only single channel or three channel images are processed

const staTIc double pi = 3.1415926;

//Gaussian filter template is generated according to 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;=””>

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;=””>

{

x2 = pow(i – origin, 2);

for (int j = 0; j < ksize;=””>

{

y2 = pow(j – origin, 2);

//The constant in front of the Gaussian function does not need to be calculated and will be eliminated in the process of normalization

double g = exp(-(x2 + y2) / (2 * sigma * sigma));

sum += g;

templateMatrix[i][j] = g;

}

}

for (int i = 0; i < ksize;=””>

{

for (int j = 0; j < ksize;=””>

{

templateMatrix[i][j] /= sum;

cout < templatematrix[i][j]=””>< “=””>

}

cout <>

}

//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;=””>

{

for (int j = border; j < cols;=””>

{

double sum[3] = { 0 };

for (int a = -border; a <= border;=””>=>

{

for (int b = -border; b <= border;=””>=>

{

if (channels == 1)

{

sum[0] += templateMatrix[border + a][border + b] * dst.at(i + a, j + b);

}

else if (channels == 3)

{

Vec3b rgb = dst.at(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;=””>

{

if (sum[k] <>

sum[k] = 0;

else if (sum[k] > 255)

sum[k] = 255;

}

if (channels == 1)

dst.at(i, j) = staTIc_ cast(sum[0]);

else if (channels == 3)

{

Vec3b rgb = { staTIc_cast(sum[0]), static_cast(sum[1]), static_cast(sum[2]) };

dst.at(i, j) = rgb;

}

}

}

//Release template array

for (int i = 0; i < ksize;=””>

delete[] templateMatrix[i];

delete[] templateMatrix;

}

Only single channel or three channel images are processed. After the template is generated, the filtering (convolution process) is relatively simple. However, for such a Gaussian filtering process, the number of cyclic operations is m × n × Ksize2, where m and N are the size of the image; Ksize is the size of the Gaussian filter. In this way, its time complexity is O (ksize2), which increases squarely with the size of the filter template. When the size of the Gaussian filter is large, its operation efficiency is very low. In order to improve the operation speed of filtering, the two-dimensional Gaussian filtering process can be decomposed.

Separation and Gaussian filtering

Due to the separability of Gaussian function, the larger Gaussian filter can be divided into two steps: firstly, the image is convoluted with one-dimensional Gaussian function in the horizontal (vertical) direction; Then, the convolution result is convoluted in the vertical (horizontal) direction using the template obtained by the same one-dimensional Gaussian function. The specific implementation code is as follows:

//Calculation of separation

void separateGaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)

{

CV_ Assert(src.channels()==1 || src.channels() == 3); // Only single channel or three channel images are processed

//Generate one-dimensional Gaussian filter template

double *matrix = new double[ksize];

double sum = 0;

int origin = ksize / 2;

for (int i = 0; i < ksize;=””>

{

//The constant in front of the Gaussian function does not need to be calculated and will be eliminated in the process of normalization

double g = exp(-(i – origin) * (i – origin) / (2 * sigma * sigma));

sum += g;

matrix[i] = g;

}

//Normalization

for (int i = 0; i < ksize;=””>

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;=””>

{

for (int j = border; j < cols;=””>

{

double sum[3] = { 0 };

for (int k = -border; k <= border;=””>=>

{

if (channels == 1)

{

sum[0] += matrix[border + k] * dst.at(i, j + k); // The row remains unchanged and the column changes; First do the horizontal convolution

}

else if (channels == 3)

{

Vec3b rgb = dst.at(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;=””>

{

if (sum[k] <>

sum[k] = 0;

else if (sum[k] > 255)

sum[k] = 255;

}

if (channels == 1)

dst.at(i, j) = static_ cast(sum[0]);

else if (channels == 3)

{

Vec3b rgb = { static_cast(sum[0]), static_cast(sum[1]), static_cast(sum[2]) };

dst.at(i, j) = rgb;

}

}

}

//Vertical direction

for (int i = border; i < rows;=””>

{

for (int j = border; j < cols;=””>

{

double sum[3] = { 0 };

for (int k = -border; k <= border;=””>=>

{

if (channels == 1)

{

sum[0] += matrix[border + k] * dst.at(i + k, j); // The column remains unchanged and the row changes; Convolution in vertical direction

}

else if (channels == 3)

{

Vec3b rgb = dst.at(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;=””>

{

if (sum[k] <>

sum[k] = 0;

else if (sum[k] > 255)

sum[k] = 255;

}

if (channels == 1)

dst.at(i, j) = static_ cast(sum[0]);

else if (channels == 3)

{

Vec3b rgb = { static_cast(sum[0]), static_cast(sum[1]), static_cast(sum[2]) };

dst.at(i, j) = rgb;

}

}

}

delete[] matrix;

}

The code is not refactored for a long time, but its implementation principle is relatively simple. Firstly, the template of one-dimensional Gaussian function is obtained. In the process of convolution (filtering), the row remains unchanged and the column changes, and the convolution operation is performed in the horizontal direction; Then, on the results obtained above, keep the columns not edge and the rows change, and perform convolution operation in the vertical direction. In this way, the time complexity of the algorithm is O (ksize) O (ksize), and the amount of computation and the template size of the filter increase linearly.

Gaussian blur is also encapsulated in opencv, and its declaration is as follows:

CV_ EXPORTS_ W void GaussianBlur( InputArray   src, OutputArray dst, Size ksize,

double sigmaX, double sigmaY = 0,

int borderType = BORDER_ DEFAULT );

The standard deviation of two-dimensional Gaussian function should have a standard deviation in the X and Y directions respectively. In the above code, it is always set that the standards in the X and Y directions are equal. In the Gaussian filter in opencv, different standard deviations can be set in the X and Y directions.

The following figure shows the comparison between the Gaussian filter implemented by ourselves and the Gaussian blur in OpenCV

The picture above is 5 × 5, σ= 0.8 Gaussian filter, it can be seen that the results obtained by the two implementations are not very different.

summary

Gaussian filter is a linear smoothing filter, and its filter template is obtained by discretizing two-dimensional Gaussian function. Because the center value of Gaussian template is the largest and decreases gradually around, the filtered result is better than the mean filter.

The most important parameter of Gaussian filter is the standard deviation of Gaussian distribution σσ， The standard deviation and Gaussian filter have great smoothing ability, σσ The larger the, the wider the frequency band of the Gaussian filter and the better the smoothness of the image. By adjusting σσ Parameter can balance the suppression of image noise and the blurring of image.