Opencv实现的陷波滤波器

Posted 2020-11-19 phoenixdsg

 在本示例中，共设计了三个函数，分别是巴特沃斯滤波器BLPF（）、巴特沃斯陷波滤波器notchFilter_BTW（）、高斯陷波滤波器notchFilter_GAUSS（）

巴特沃斯陷波滤波器参见书上6.4.5选择性滤波器，高斯陷波滤波器参照6.3.3节中的高斯低通滤波器。

参数说明：

rows——滤波器的行数；
cols——滤波器的列数；
D0——频域截止半径；
n——巴特沃斯的阶数；
cvtype——滤波器的数据类型和通道数，默认为双通道浮点数

 

Mat BLPF(int rows,int cols,float D0,int n=1,int cvtype=CV_32FC2)
{
    Mat filt(rows,cols,cvtype,Scalar::all(0));
    int cx=cols/2,cy=rows/2;
    float D02=D0*D0;
    for(int i=0;i<rows;i++)
    {
        for(int j=0;j<cols;j++)
        {
            int u=cx-j,v=cy-i;//中心坐标
            float Duv2=u*u+v*v;//距离中心半径的平方
            float H0=1-1/(1+pow(Duv2/D02,n));
            int u1=u-130,v1=v-130;
            float D2=u1*u1+v1*v1;
            float H1=1-1/(1+pow(D2/D02,n));
            //.data返回的是uchar*型指针，所以要强制转换成浮点数型
            float* p=(float*)(filt.data+i*filt.step[0]+j*filt.step[1]);

            for(int c=0;c<filt.channels();c++)
            {
                p[c]=H1*H0;
            }

        }
    }
    return filt;
}


Mat notchFilter_BTW(int rows,int cols,std::vector<cv::Point> np,
                float* D,int n=1,int cvtype=CV_32FC2)
{
    Mat filt(rows,cols,cvtype,Scalar::all(0));
    int cx=cols/2,cy=rows/2;
    int numNotch=np.size();
    float* D2=D;
    for(int i=0;i<numNotch;i++)
    {
        D2[i]=D[i]*D[i];
    }
    int uk[numNotch],vk[numNotch];//在画面上的实际陷波坐标点
    int u_k[numNotch],v_k[numNotch];//陷波共轭点
    float Dk[numNotch],D_k[numNotch];//陷波半径r
    float Hk[numNotch],H_k[numNotch];

    for(int i=0;i<rows;i++)
    {
        for(int j=0;j<cols;j++)
        {
            int u=cx-j,v=cy-i;//中心坐标
            for(int s=0;s<numNotch;s++)
            {
                uk[s]=u+np[s].x,vk[s]=v+np[s].y;
                Dk[s]=uk[s]*uk[s]+vk[s]*vk[s];//距离中心半径的平方
                Hk[s]=1-1/(1+pow(Dk[s]/D2[s],n));

                u_k[s]=u-np[s].x,v_k[s]=v-np[s].y;
                D_k[s]=u_k[s]*u_k[s]+v_k[s]*v_k[s];
                H_k[s]=1-1/(1+pow(D_k[s]/D2[s],n));
            }
            //.data返回的是uchar*型指针，所以要强制转换成浮点数型
            float* p=(float*)(filt.data+i*filt.step[0]+j*filt.step[1]);

            for(int c=0;c<filt.channels();c++)
            {
                p[c]=Hk[0]*H_k[0];
                for(int k=1;k<numNotch;k++)
                {
                    p[c]*=Hk[k]*H_k[k];
                }
            }

        }
    }
    return filt;
}
Mat notchFilter_GAUSS(int rows,int cols,std::vector<cv::Point> np,
                float* D,int cvtype=CV_32FC2)
{
    Mat filt(rows,cols,cvtype,Scalar::all(0));
    int cx=cols/2,cy=rows/2;

//    float D02=D0*D0;
    int numNotch=np.size();
    float* D2=D;
    for(int i=0;i<numNotch;i++)
    {
        D2[i]=D[i]*D[i];
    }
    int uk[numNotch],vk[numNotch];//在画面上的实际陷波坐标点
    int u_k[numNotch],v_k[numNotch];//陷波共轭点
    float Dk[numNotch],D_k[numNotch];//陷波半径r
    float Hk[numNotch],H_k[numNotch];

    for(int i=0;i<rows;i++)
    {
        for(int j=0;j<cols;j++)
        {
            int u=cx-j,v=cy-i;//中心坐标
            for(int s=0;s<numNotch;s++)
            {
                uk[s]=u+np[s].x,vk[s]=v+np[s].y;
                Dk[s]=uk[s]*uk[s]+vk[s]*vk[s];//距离中心半径的平方
                Hk[s]=1-exp(-Dk[s]/(D2[s]*2));

                u_k[s]=u-np[s].x,v_k[s]=v-np[s].y;
                D_k[s]=u_k[s]*u_k[s]+v_k[s]*v_k[s];
                H_k[s]=1-exp(-D_k[s]/(D2[s]*2));
            }
            //.data返回的是uchar*型指针，所以要强制转换成浮点数型
            float* p=(float*)(filt.data+i*filt.step[0]+j*filt.step[1]);

            for(int c=0;c<filt.channels();c++)
            {
                p[c]=Hk[0]*H_k[0];
                for(int k=1;k<numNotch;k++)
                {
                    p[c]*=Hk[k]*H_k[k];
                }
            }

        }
    }
    return filt;
}

int main()
{
    Point np[]={Point(130,130),Point(90,130),Point(130,100),Point(90,100)};//输入陷波坐标数组
    vector<Point> vnp(np,np+4);
    float D[4]={5,10,15,20};
//    Mat filt1=notchFilter_BTW(500,600,vnp,D,2);
    Mat filt1=notchFilter_GAUSS(500,600,vnp,D);
    Mat fc1;
    extractChannel(filt1,fc1,0);
    imshow("filter ",fc1);

    waitKey();
    return 0;
}

高斯陷波滤波器的演示实例结果如下：

下面是巴特沃斯陷波滤波器的演示结果：