Traditional Gaussian blur and optimization algorithm (with complete C + + code)

Source: Internet
Author: User

Gaussian Blur (English: Gaussian Blur), also known as Gaussian smoothing, is a widely used processing effect in image processing software such as Adobe Photoshop, GIMP, and paint.net, often used to reduce image noise and reduce levels of detail. The image generated by this blur technique is visually similar to a translucent screen that looks at the image, which is significantly different from the bokeh of the lens and the shadows of ordinary lighting. Gaussian smoothing is also used in the pre-processing phase of computer vision algorithms to enhance image performance at different proportions (see scale space representation and scale space implementations). From the point of view of mathematics, the Gaussian blur process of image is the convolution of image and normal distribution. Since the normal distribution is also called Gaussian distribution, this technique is called Gaussian blur. Blurring the image with the circular box will produce a more accurate out-of-focus imaging effect. Because the Fourier transform of the Gaussian function is another Gaussian function, the Gaussian blur is a low-pass filter for the image.

Gaussian Blur is an image fuzzy filter that calculates the transformation of each pixel in an image with a normal distribution. n-dimensional normal distribution equation is

In a two-dimensional space defined as

where R is the fuzzy radius (), Σ is the standard deviation of the normal distribution. In a two-dimensional space, the contours of the surface generated by this formula are concentric circles that begin with a normal distribution from the center. A convolution matrix that distributes nonzero pixels is transformed from the original image. The value of each pixel is a weighted average of the surrounding neighboring pixel values. The value of the original pixel has the largest Gaussian distribution value, so there is a maximum weight, and the neighboring pixels are getting farther away from the original pixels, and their weights are getting smaller. In this way, the blurring process retains the edge effect more than other equalization fuzzy filters, see scale space implementation.

Theoretically, the distribution of each point in the image is nonzero, which means that each pixel needs to contain an entire image. In practical applications, when calculating the discrete approximation of a Gaussian function, pixels outside the approximate 3σ distance can be considered as ineffective, and the calculation of these pixels can be ignored. In general, the image processing program only needs the computed matrix to guarantee the relevant pixel influence. For points on the boundary, it is common to copy the surrounding points to the other side and then perform the weighted averaging.

In addition to circular symmetry, Gaussian blur can also be used to calculate two independent one-dimensional spaces on a two-dimensional image, which is called linear . This means that the effect of using a two-dimensional matrix transformation can also be obtained by a Gaussian matrix transformation in the horizontal direction plus a Gaussian matrix transformation in the vertical direction. From a computational point of view, this is a useful feature, because it requires only a second calculation, and the irreducible matrix requires a secondary calculation, which is the dimension of the image that needs to be filtered, and the dimension of the filter.

The effect of multiple continuous Gaussian blur on an image can produce the same effect as a larger Gaussian blur, and the radius of a large Gaussian blur is the square root of the sum of squares of multiple Gaussian blur radii. For example, using a two-time Gaussian blur transform with a radius of 6 and 8 is equivalent to a Gaussian blur effect with a radius of 10. Depending on this relationship, using more than one continuous, smaller Gaussian blur will be less than a single Gaussian larger processing time.

Gaussian blur is often used in cases where the image size is reduced. In the cases of under-sampling, the image is usually treated with low-pass filtering prior to sampling. This ensures that false high-frequency information is not present in the sampled image. Gaussian Blur has good properties, such as no apparent boundary, so that it does not form a shock in the filtered image.

The above information is excerpted from Wikipedia (Gaussian blur entry):

Https://zh.wikipedia.org/wiki/%E9%AB%98%E6%96%AF%E6%A8%A1%E7%B3%8A

So how exactly?

The code Presents:

inline int* buildgausskern (int winsize, int sigma) {int wincenter, x;float sum = 0.0f;wincenter = Winsize/2;float *kern = (float*) malloc (winsize*sizeof (float)), int *ikern = (int*) malloc (winsize*sizeof (int)), float sqrt_2pi = 2.506628274631f;float sigmamul2pi = 1.0f/(Sigma * sqrt_2pi), float divSigmaPow2 = 1.0f/(2.0f * sigma * sigma); for (x = 0; X < Wincenter + 1; X + +) {Kern[wincenter-x] = kern[wincenter + x] = exp (-(x * x) * divSigmaPow2) * sigmamul2pi;sum + = Kern[wincenter-x] + (( X! = 0)? Kern[wincenter + x]: 0.0);} sum = 1.0f/sum;for (x = 0; x < winsize; + +) {Kern[x] *= sum;ikern[x] = kern[x] * 256.0f;} Free (kern); return Ikern;} void Gaussblur (unsigned char* pixels, unsigned int width, unsigned int height, unsigned int channels, int sigma) {WID th = 3 * width;if ((width% 4)! = 0) Width + = (4-(width% 4)); unsigned int winsize = (1 + (((int) ceil (3 * sigma) * 2))  ; int *gausskern = Buildgausskern (Winsize, sigma); winsize *= 3;unsigned int halfsize = WINSIZE/2; UnSigned char *tmpbuffer = (unsigned char*) malloc (Width * height* sizeof (unsigned char)); for (unsigned int h = 0; h < Hei Ght h++) {unsigned int rowwidth = h * width;for (unsigned int w = 0; w < width; w + = Channels) {unsigned int rowr = 0;unsign ed int rowg = 0;unsigned int ROWB = 0;int * Gausskernptr = gausskern;int whalfsize = w + width-halfsize;unsigned int cu RpoS = rowwidth + w;for (unsigned int k = 1; k < winsize; k + = Channels) {unsigned int pos = Rowwidth + ((k + whalfsize  )% width); int Fkern = *gausskernptr++;rowr + = (Pixels[pos] * fkern); ROWG + = (pixels[pos + 1] * Fkern); ROWB + = (pixels[pos + 2] * fkern);} Tmpbuffer[curpos] = ((unsigned char) (ROWR >> 8)), Tmpbuffer[curpos + 1] = ((unsigned char) (ROWG >> 8)); Tmpbuffe R[curpos + 2] = ((unsigned char) (ROWB >> 8);}} Winsize/= 3;halfsize = winsize/2;for (unsigned int w = 0; w < width; w++) {for (unsigned int h = 0; h < height; +) {unsigned int col_all = 0;int hhalfsize = h + height-halfsize;for (unsigned int k = 0; k < winsize; k++) {Col_all + = tmpbuffer[((k + hhalfsize)% height) * Width + W] * gausskern[k];} Pixels[h * Width + W] = (unsigned char) (Col_all >> 8);}} Free (Tmpbuffer), free (gausskern); }

Note:

To the original algorithm, I made a few minor changes, mainly to consider a little bit of performance problems.

Sometimes I write too many comments instead of being verbose, so I will look at ha.

This code, the measured speed is very bad, processing a 5000x3000 in the radius of 5 or so will take 10 seconds to a few 10 seconds, it is difficult to accept.

Because of the speed problem, there are many optimization algorithms on the net.

Before I also sent a "fast Gaussian fuzzy algorithm", under the same conditions, the algorithm has been like the above algorithm more than 10 times times faster.

Since this code is hard to read, I have made further adjustments and optimizations.

void Gaussianblur (Unsigned char* img, unsigned int width, unsigned int height, unsigned int channels, unsigned int radius {radius = min (max (1, RADIUS), 248), unsigned int kernelsize = 1 + radius * 2;unsigned int* kernel = (unsigned int*) malloc ( kernelsize* sizeof (unsigned int)); memset (kernel, 0, kernelsize* sizeof (unsigned int)); int (*mult) [n] = (int (*) [256]) malloc (kernelsize * *) sizeof (int)); memset (mult, 0, Kernelsize * () sizeof (int)); Intxstart = 0;intystart = 0;width = Xstart + Width-max (0, (xstart + width)-width), height = ystart + height-max (0, (y Start + height)-height), int imageSize = Width*height;int Widthstep = width*channels;if (Channels = = 3 | | channels = 4) { unsigned char * cacheimg = nullptr;    Cacheimg = (unsigned char *) malloc (sizeof (unsigned char) * imageSize * 6); if (cacheimg = = nullptr) return;unsigned char * Rcache = cacheimg;unsigned char * gcache = cacheimg + imagesize;unsigned char * bCache = cacheimg + imageSize * 2 ; unsigned char * R2CAChe = cacheimg + imageSize * 3;unsigned char * g2cache = cacheimg + imageSize * 4;unsigned char * B2cache = cacheimg + imageSize * 5;int sum = 0;for (int K = 1; K < radius; k++) {unsigned int szi = Radius-k;kernel[radius + K] = Kernel[szi] = szi*szi;sum + = Kernel[szi] + kernel[szi];for (int J = 0; J < 256; J + +) {Mult[radius + k][j] = mult[szi][j] = Kernel[szi] * j;}} Kernel[radius] = radius*radius;sum + = kernel[radius];for (int j = 0; J < N; j + +) {Mult[radius][j] = Kernel[radius] * j ;} for (int Y = 0; Y < height;  ++y) {unsigned char* lineps = img + y*widthstep;unsigned char* linepr = rcache + y*width;unsigned char* LINEPG = Gcache + y*width;unsigned char* LINEPB = BCache + y*width;for (int X = 0; X < width; ++X) {int p2 = x*channels; LINEPR[X] = lineps[p2]; LINEPG[X] = lineps[p2 + 1]; LINEPB[X] = lineps[p2 + 2];}} int kernelsum = 0;for (int K = 0; K < kernelsize; k++) {kernelsum + = Kernel[k];} float fkernelsum = 1.0f/kernelsum;for (int Y = Ystart;Y < height; y++) {int heightstep = Y * width;unsigned char* linepr = rcache + heightstep;unsigned char* linepg = Gcache + heigh tstep;unsigned char* LINEPB = BCache + heightstep;for (int X = Xstart; X < width; X + +) {int cb = 0;int CG = 0;int CR = 0;for (int K = 0; K < kernelsize; k++) {unsigned int readpos = ((X-radius + K + width)% width); int * Pmult = MULT[K];CR + = Pmult[linepr[readpos]];c G + = PMULT[LINEPG[READPOS]];CB + = Pmult[linepb[readpos]];} unsigned int p = heightstep + x;r2cache[p] = cr* Fkernelsum;g2cache[p] = cg* Fkernelsum;b2cache[p] = cb* fkernelsum;}} for (int X = Xstart; X < width;      X + +) {int Widthcomp = X*channels;int Widthstep = width*channels;unsigned char* lineps = img + x*channels;unsigned char* Linepr = R2cache + x;unsigned char* linepg = g2cache + x;unsigned char* LINEPB = B2cache + x;for (int Y = YSt Art Y < height; y++) {int cb = 0;int CG = 0;int CR = 0;for (int K = 0; K < kernelsize; k++) {unsigned int readpos = ((Y-radius + K + height)% height) * Width;int * Pmult = MULT[K];CR + = PMULT[LINEPR[READPOS]];CG + Pmult[linepg[readpos] ];CB + = Pmult[linepb[readpos]];} int p = y*widthstep; LINEPS[P] = (unsigned char) (CR * fkernelsum); Lineps[p + 1] = (unsigned char) (CG * fkernelsum); Lineps[p + 2] = (unsigned char) (cb* fkernelsum); }}free (cacheimg);} else if (channels = = 1) {unsigned char * cacheimg = nullptr;    Cacheimg = (unsigned char *) malloc (sizeof (unsigned char) * imageSize * 2); if (cacheimg = = nullptr) return;unsigned char * Rcache = cacheimg;unsigned char * r2cache = cacheimg + imagesize;int sum = 0;for (int K = 1; K < radius; k++) {unsigned int szi = Radius-k;kernel[radius + K] = Kernel[szi] = szi*szi;sum + = Kernel[szi] + kernel[szi];for (int J = 0; J < 256; J + +) {Mult[radius + k][j] = mult[szi][j] = Kernel[szi] * j;}} Kernel[radius] = radius*radius;sum + = kernel[radius];for (int j = 0; J < N; j + +) {Mult[radius][j] = Kernel[radius] * j ;} for (int Y = 0; Y < height; ++y) {Unsigned char* lineps = img + y*widthstep;unsigned char* linepr = Rcache + y*width;for (int X = 0; X < width; ++X) {Linepr[x] = lineps[x];}} int kernelsum = 0;for (int K = 0; K < kernelsize; k++) {kernelsum + = Kernel[k];} float fkernelsum = 1.0f/kernelsum;for (int Y = Ystart; Y < height; y++) {int heightstep = Y * width;unsigned char* linepr = Rcache + heightstep;for (int X = Xstart; X < width; X + +) {int cb = 0;int CG = 0;int CR = 0;for (int K = 0; K < kernelsize; k++) {unsigned int readpos = ((X-radius + k+width)%width); int * Pmult = MULT[K];CR + = Pmult[linepr[readpos]];} unsigned int p = heightstep + x;r2cache[p] = cr * FKERNELSUM;}} for (int X = Xstart; X < width;      X + +) {int Widthcomp = X*channels;int Widthstep = width*channels;unsigned char* lineps = img + x*channels;unsigned char* Linepr = R2cache + x;for (int Y = Ystart; Y < height; y++) {int cb = 0;int CG = 0;int CR = 0;for (int K = 0; K < kernelsize; k++) {unsigned int readpos = (y-rAdius + k+height)%height) * Width;int * Pmult = MULT[K];CR + = Pmult[linepr[readpos]];} int p = y*widthstep; LINEPS[P] = (unsigned char) (cr* fkernelsum);}} Free (cacheimg);} Free (kernel); free (mult);}

Some of the algorithm optimization techniques, want to also play a bit of a role.

Put a:

This article is only to give a try, if there are other related questions or needs can also be contacted by email to discuss.

e-mail address is:

[Email protected]

Traditional Gaussian blur and optimization algorithm (with complete C + + code)

Related Article

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.