#include"stdafx.h"typedefstruct{floatB;float b[4];
} gauss_coefs;//参数计算
void compute_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型高斯滤波//**************************************************************//
//参考文献:Recursive implementation of the Gaussian filter//Src:输入图像//Dest:输出图像//sigma:高斯标准差//**************************************************************//IS_RET IIRGaussianFilter(TMatrix *Src,TMatrix *Dest,floatsigma)
{int X,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);//缓存
if(Buff==NULL){returnIS_RET_ERR_OUTOFMEMORY;}for(int i=0;i
{
LinePS=Src->Data+i;
LinePD=Dest->Data+i;//拷贝原图至缓存Buff
BPointer=Buff;for (Y=0;Y
{for (X=0;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){returnIS_RET_ERR_OUTOFMEMORY;}
w2=(float*)malloc(sizeof(float)*(Width+3));if(w2==NULL){returnIS_RET_ERR_OUTOFMEMORY;}for(Y=0;Y
{//前向滤波
w1[0]=w1[1]=w1[2]=BPointer[0];for(int n=3,i=0;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){returnIS_RET_ERR_OUTOFMEMORY;}
w2=(float*)realloc(w2,sizeof(float)*(Height+3));if(w2==NULL){returnIS_RET_ERR_OUTOFMEMORY;}for (X=0;X
{//前向滤波
w1[0]=w1[1]=w1[2]=BPointer[0];for(int n=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]);
}//后向滤波
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
{for(X=0;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;
}