IIR型高斯滤波的原理及实现

#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;
        //拷贝原图至缓存Buff
        BPointer=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;
}