java 滤波 iir_IIR型高斯滤波的原理及实现

#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;

}