#include"stdafx.h"
typedef struct
{float B;float b[4];
} gauss_coefs;//参数计算
void compute_coefs3(gauss_coefs *c,float sigma)
{float q, 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型高斯滤波
//**************************************************************//
//参考文献：Recursive implementation of the Gaussian filter
//Src:输入图像
//Dest:输出图像
//sigma:高斯标准差
//**************************************************************//
IS_RET  IIRGaussianFilter(TMatrix *Src,TMatrix *Dest,float sigma)
{int X,Y,Width=Src->Width,Height=Src->Height,Stride=Src->WidthStep,Channel=Src->Channel;gauss_coefs  c;compute_coefs3(&c,sigma);unsigned char *LinePS,*LinePD;float *Buff,*BPointer,*w1,*w2;    Buff=(float*)malloc(sizeof(float)*Height*Width);//缓存if(Buff==NULL){return IS_RET_ERR_OUTOFMEMORY;}for(int i=0;i<Channel;i++){LinePS=Src->Data+i;LinePD=Dest->Data+i;//拷贝原图至缓存BuffBPointer=Buff;for (Y=0;Y<Height;Y++){for (X=0;X<Width;X++){BPointer[0]=float(LinePS[0]);BPointer++;LinePS+=Channel;}LinePS+=Stride-Width*Channel;}//横向滤波BPointer=Buff;w1=(float*)malloc(sizeof(float)*(Width+3));if(w1==NULL){return IS_RET_ERR_OUTOFMEMORY;}w2=(float*)malloc(sizeof(float)*(Width+3));if(w2==NULL){return IS_RET_ERR_OUTOFMEMORY;}for(Y=0;Y<Height;Y++){//前向滤波w1[0]=w1[1]=w1[2]=BPointer[0];for(int n=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]);}//后向滤波w2[Width]=w2[Width+1]=w2[Width+2]=w1[Width+2];for(int n=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;}//纵向滤波BPointer=Buff;w1=(float*)realloc(w1,sizeof(float)*(Height+3));if(w1==NULL){return IS_RET_ERR_OUTOFMEMORY;}w2=(float*)realloc(w2,sizeof(float)*(Height+3));if(w2==NULL){return IS_RET_ERR_OUTOFMEMORY;}for (X=0;X<Width;X++){//前向滤波w1[0]=w1[1]=w1[2]=BPointer[0];for(int n=3,i=0;i<Height;n++,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]);}//后向滤波w2[Height]=w2[Height+1]=w2[Height+2]=w1[Height+2];for(int n=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++;}//拷贝缓存至结果BPointer=Buff;for(Y=0;Y<Height;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);}return IS_RET_OK;
}