The principle and realization of IIR Gaussian filter

Second, the realization

The implementation of IIR Gaussian filter in GIMP, the code is located in Contrast-retinex.c, the reader can see for themselves. Here is the core code I implemented:

#include"stdafx.h"typedefstruct{    floatB; floatb[4];} Gauss_coefs;//parameter CalculationvoidCOMPUTE_COEFS3 (Gauss_coefs *c,floatSigma) {    floatQ, Q2, Q3; if(Sigma >=2.5) {Q=0.98711* Sigma-0.96330; }    Else if(Sigma >=0.5) && (Sigma <2.5) ) {Q=3.97156-4.14554* (float) sqrt (float)1-0.26891*Sigma); }    Else{Q=0.1147705018520355224609375; } Q2= Q *Q; Q3= Q *Q2; C->b[0] =1/(1.57825+(2.44413*Q) + (1.4281*Q2) + (0.422205*Q3)); C->b[1] = (        (2.44413*Q) + (2.85619*Q2) + (1.26661*Q3)) *c->b[0]; C->b[2] = (                   -((1.4281*Q2) + (1.26661*Q3)) *c->b[0]; C->b[3] = (                                 (0.422205*Q3)) *c->b[0]; C->b =1.0-(c->b[1]+c->b[2]+c->b[3]);}//IIR Type Gaussian filter//**************************************************************////References: Recursive implementation of the Gaussian filter//SRC: Input Image//Dest: Output Image//Sigma: Gaussian standard deviation//**************************************************************//Is_ret iirgaussianfilter (Tmatrix *src,tmatrix *dest,floatSigma) {    intX,y,width=src->width,height=src->height,stride=src->widthstep,channel=src->Channel;    Gauss_coefs C; COMPUTE_COEFS3 (&C,sigma); unsignedChar*lineps,*LINEPD; float*buff,*bpointer,*w1,*W2; Buff=(float*)malloc(sizeof(float) *height*width);//Cache    if(Buff==null) {returnIs_ret_err_outofmemory;}  for(intI=0; i<channel;i++) {lineps=src->data+i; LINEPD=dest->data+i; //copy original to cache buffBpointer=Buff;  for(y=0; y)        {             for(x=0; x<width; X + +) {bpointer[0]=float(lineps[0]); Bpointer++; Lineps+=Channel; } lineps+=stride-width*Channel; }        //Transverse filteringBpointer=Buff; W1=(float*)malloc(sizeof(float) * (width+3)); if(W1==null) {returnIs_ret_err_outofmemory;} W2=(float*)malloc(sizeof(float) * (width+3)); if(W2==null) {returnIs_ret_err_outofmemory;}  for(y=0; y)        {            //Forward Filteringw1[0]=w1[1]=w1[2]=bpointer[0];  for(intn=3, i=0; i<width;n++,i++) {W1[n]=c.b*bpointer[i]+ (c.b[1]*w1[n-1]+c.b[2]*w1[n-2]+c.b[3]*w1[n-3]); }            //back-to-filterw2[width]=w2[width+1]=w2[width+2]=w1[width+2];  for(intn=width-1; n>=0; n--) {Bpointer[n]=w2[n]=c.b*w1[n+3]+ (c.b[1]*w2[n+1]+c.b[2]*w2[n+2]+c.b[3]*w2[n+3]); } bpointer+=Width; }        //longitudinal filteringBpointer=Buff; W1=(float*)realloc(W1,sizeof(float) * (height+3)); if(W1==null) {returnIs_ret_err_outofmemory;} W2=(float*)realloc(W2,sizeof(float) * (height+3)); if(W2==null) {returnIs_ret_err_outofmemory;}  for(x=0; x<width; X + +)        {            //Forward Filteringw1[0]=w1[1]=w1[2]=bpointer[0];  for(intn=3, i=0; i) {W1[n]=c.b*bpointer[i*width]+ (c.b[1]*w1[n-1]+c.b[2]*w1[n-2]+c.b[3]*w1[n-3]); }            //back-to-filterw2[height]=w2[height+1]=w2[height+2]=w1[height+2];  for(intn=height-1; n>=0; n--) {Bpointer[n*width]=w2[n]=c.b*w1[n+3]+ (c.b[1]*w2[n+1]+c.b[2]*w2[n+2]+c.b[3]*w2[n+3]); } bpointer++; }        //Copy Cache to resultsBpointer=Buff;  for(y=0; y)        {             for(x=0; x<width; X + +) {linepd[0]=bpointer[0]; if(bpointer[0]>255) linepd[0]=255; if(bpointer[0]<0) linepd[0]=0; Bpointer++; LINEPD+=Channel; } LINEPD+=stride-width*Channel; }         Free(W1);  Free(W2); }    returnIs_ret_ok;}

Experimental results: For a 1024*1024 color image, the algorithm is time consuming 175ms.

Reference Documents :

Young I T, Van Vliet L J. Recursive implementation of the Gaussian Filter[j]. Signal processing, 1995, 44 (2): 139-151.

The principle and realization of IIR Gaussian filter