6.3.2巴特沃斯（butterworth）低通滤波器

在本程序中，共有六个自定义函数，分别是：

1. myMagnitude（Mat & complexImg,Mat & magnitudeImage），在该函数中封装了Opencv中的

   magnitude函数，实现对于复数图像的幅值计算。该函数共有两个参数：

   complexImg--输入的复数阵列，或复数图像

   magnitudeImage--输出的幅值阵列，或幅值图像

2. dftshift(Mat& ds)，该函数实现对图像四个象限的对角互换，相当于MatLab中 fftshift（），将频谱的原点（0，0）移到图像中心。

3. srcCentralized(Mat& src)用于傅里叶变换前的预处理，以便得到傅里叶频谱的原点（0，0）位于图像的中心。

   该函数与dftshift（）目的一致，实现方法不同，一个是变换前预处理，一个是变换后处理。

4. imshowComplexMat(Mat&dftDst,String winName,bool inverseSpectrum)，该函数用于显示复数图像或双通道矩阵，共有三个参数：

   dftDst--待显示的复数矩阵

   winName--显示复数矩阵的窗口名字

   inverseSpectrum-输入的dftDst是正向傅里叶变换的结果，还是逆傅里叶变换的结果

5. createFilterButterworth(Mat&filter,int n,int R,int W,FilterForm filterform)，用于制作Butterworth频域滤波器，该函数利用了ptr（）

   指针遍历图像的方法。该函数可以实现低通、高通、带通、带阻滤波器。目前该函数共有五个参数：

   filter--输入的矩阵，要求数据类型为CV_64FC2；

   n--巴特沃斯阶数

   R--截止频率半径，如果小于0，则返回一个全口径滤波器，否则返回一个口径受限的滤波器

   W--带宽

   filterform--滤波器形式，它是个枚举类型数据，enum FilterForm{LOW_PASS_FILTER,HIGH_PASS_FILTER,BAND_PASS_FILTER,BAND_STOP_FILTER};

6.void myDft(Mat&src,Mat&dst,bool isProCentralized,bool doubleSizeOrNot)，该函数更有四个参数

  src--是输入的原图像

  dst--是傅里叶变换的输出图像

  isProCentralized--表示是否调用SRCCentralized函数，对src进行中心化预处理

  doubleSizeOrNot--表示是否需要将原图像尺寸扩展为两倍，以便解决卷积缠绕问题

#include <iostream>

#include<opencv2/opencv.hpp>

using namespace cv;

using namespace std;

#define CV_MAT_ELEM2(src,dtype,y,x) \

    (dtype*)(src.data+y*src.step[]+x*src.step[])

/**************程序说明****************

在主程序内容，参见注释。

**********************************************/

/*1.求取复数矩阵的幅值*/

void myMagnitude(Mat & complexImg,Mat & magnitudeImage)

{

    Mat planes[];

    split(complexImg,planes);

    magnitude(planes[],planes[],magnitudeImage);

}

/*傅里叶变换后的频谱图后处理，将傅里叶普的原点(0,0)平移到图像的中心*/

void dftshift(Mat& ds)

{

    int cx=ds.cols/;//图像的中心点x 坐标

    int cy=ds.rows/;//图像的中心点y 坐标

    Mat q0=ds(Rect(,,cx,cy));//左上

    Mat q1=ds(Rect(cx,,cx,cy));//右上

    Mat q2=ds(Rect(,cy,cx,cy));//左下

    Mat q3=ds(Rect(cx,cy,cx,cy));//右下

    Mat tmp;

    q0.copyTo(tmp);

    q3.copyTo(q0);

    tmp.copyTo(q3);

    q1.copyTo(tmp);

    q2.copyTo(q1);

    tmp.copyTo(q2);

}

/*傅里叶变换前的预处理，以便频谱图的原点（0，0）移动到图像的中心*/

void srcCentralized(Mat& src)

{

    switch (src.depth())

    {

    case CV_32F:

        for(int i=;i<src.rows;i++)

        {

            for(int j=;j<src.cols;j++)

            {

                float* mv=CV_MAT_ELEM2(src,float,i,j);

                if((i+j)%!=)

                {

                    for(int c=;c<src.channels();c++)//对所有通道同样操作

                        mv[c]=-mv[c];//如果i+j为奇数，该像素值取负值

                }

            }

        }

        break;

    case CV_64F:

        for(int i=;i<src.rows;i++)

        {

            for(int j=;j<src.cols;j++)

            {

                double* mv=CV_MAT_ELEM2(src,double,i,j);

                if((i+j)%!=)

                {

                    for(int c=;c<src.channels();c++)//遍历各个通道

                        mv[c]=-mv[c];//如果i+j为奇数，该像素值取负值

                }

            }

        }

        break;

    default:

        break;

    }

}

/*在窗口中显示复数图像，如果是正向傅里叶矩阵，需要取log才能显示更多频谱信息

 *如果是逆傅里叶变换，通过normalize归一化后，显示频谱图*/

void imshowComplexMat(Mat&dftDst,String winName,bool inverseSpectrum)

{

    Mat magI;

    myMagnitude(dftDst,magI);

    if(!inverseSpectrum)//如果是正向傅里叶变换谱

    {

        magI+=Scalar::all();

        log(magI,magI);

    }

    normalize(magI,magI,,,NORM_MINMAX);

    imshow(winName,magI);

}

///

enum FilterForm{LOW_PASS_FILTER,HIGH_PASS_FILTER,BAND_PASS_FILTER,BAND_STOP_FILTER};

void createFilterButterworth(Mat&filter,int n,int R,int W,FilterForm filterform)

{

    double Rs=R*R;//R1_square

    int cx=filter.cols/;

    int cy=filter.rows/;

    switch(filterform)

    {

    case LOW_PASS_FILTER:

        for(int i=;i<filter.rows;i++)

        {

            Vec2d* pf=filter.ptr<Vec2d>(i);

            for(int j=;j<filter.cols;j++)

            {

                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);//rs表示r的平方

                pf[j][]=./(.+pow(rs/Rs,n));//Rs是R的平方,

                pf[j][]=pf[j][];

            }

        }

        break;

    case HIGH_PASS_FILTER:

        for(int i=;i<filter.rows;i++)

        {

            Vec2d* pf=filter.ptr<Vec2d>(i);

            for(int j=;j<filter.cols;j++)

            {

                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);

                double Lp=./(.+pow(rs/Rs,n));//巴特沃斯公式

                pf[j][]=1.0-Lp;

                pf[j][]=pf[j][];

            }

        }

        break;

    case BAND_STOP_FILTER:

        for(int i=;i<filter.rows;i++)

        {

            Vec2d* pf=filter.ptr<Vec2d>(i);

            for(int j=;j<filter.cols;j++)

            {

                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);

                double r=std::sqrt(rs);//r相当于书上的D

                pf[j][]=./(.+pow(r*W/(rs-Rs),*n));

                pf[j][]=pf[j][];

            }

        }

        break;

    case BAND_PASS_FILTER:

        for(int i=;i<filter.rows;i++)

        {

            Vec2d* pf=filter.ptr<Vec2d>(i);

            for(int j=;j<filter.cols;j++)

            {

                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);

                double r=std::sqrt(rs);//r相当于书上的D

                pf[j][]=-./(.+pow(r*W/(rs-Rs),*n));

                pf[j][]=pf[j][];

            }

        }

        break;

    }

///***************【显示滤波器，如果不需要显示，将代码注销】*************/

    Mat displayFilter;

    extractChannel(filter,displayFilter,);

    imshow("btw filter image",displayFilter);

}

void myDft(Mat&src,Mat&dst,bool isProCentralized=false,bool doubleSizeOrNot=false)

{

    CV_Assert(src.channels()==);//验证src是否是单通道

    int ny=src.rows,nx=src.cols;

    Mat srcPadded;

    if(doubleSizeOrNot)//如果doubleOrNot为真，对src尺度扩展一倍

    {

        cv::copyMakeBorder(src,srcPadded,,ny,,nx,BORDER_CONSTANT);

    }

    else//否则，将src填补为最优傅里叶变换尺寸

    {

        int padx=getOptimalDFTSize(nx);//获取最优傅里叶变换列数

        int pady=getOptimalDFTSize(ny);//获取最优傅里叶变换行数数

        cv::copyMakeBorder(src,srcPadded,,pady-ny,,padx-nx,BORDER_CONSTANT);

    }

    if(srcPadded.type()!=CV_64FC1)

        srcPadded.convertTo(srcPadded,CV_64FC1);//转成浮点数据类型

    //如果isProCentralized为真，则在傅里叶变换前，对src做中心化预处理,

    //这样在傅里叶变换后的频谱图的（0，0）点就会位于频谱图的中心

    if(isProCentralized)

        srcCentralized(srcPadded);

    Mat planes[]={srcPadded,Mat::zeros(srcPadded.rows,srcPadded.cols,srcPadded.type())};

    Mat srcComplex;

    merge(planes,,srcComplex);

    dft(srcComplex,dst,DFT_COMPLEX_OUTPUT,);//离散傅立叶变换

    cout<<"srcPadded.rows="<<srcPadded.rows<<"  "<<"src.cols="<<srcPadded.cols<<endl;

}

int main()

{

    ///读入灰度图像

    Mat src=imread("D:\\Qt\\MyImage\\Fig0333(a).tif",);

    ///***************【调用自定义的傅里叶变换myDft()】*************/

    Mat dftDst;//预声明dft的输出结果矩阵

    myDft(src,dftDst,,);//调用自定义的dft函数，对src执行傅里叶变换，dftDst是傅里叶变换结果

    dftshift(dftDst);//将傅里叶变换的结果，四象限对角互换

    imshowComplexMat(dftDst,"dftSpectrum display",false);//显示傅里叶频谱图

    ///***************【创建滤波器】*************/

    Mat filter(dftDst.rows,dftDst.cols,CV_64FC2,Scalar(,));
    //巴特沃斯滤波器的阶数n=2，半径分别取R=10，30，60，160，460

    createFilterButterworth(filter,,,,LOW_PASS_FILTER);

    ///***************【频域滤波】*************/

    Mat temp;

    temp=dftDst.mul(filter);

    ///***************【调用opencv中的dft，执行逆傅里叶变换】*************/

    Mat filteredResult;

    dft(temp,filteredResult,DFT_INVERSE);

    imshowComplexMat(filteredResult,"Frequency domain butterworth filtered image",);
    imshow("srcimage",src);

    waitKey();

    return ;

}

上面四幅图，分别是巴特沃斯滤波器，n=2，截止半径=10；第二副图像是src的频谱图；第三幅图像是src原图像；第四幅是滤波后的结果。

下面是截止半径R分别为30,60,160和460时的滤波结果：

巴特西

6.3.2巴特沃斯（butterworth）低通滤波器

最新文章

热门文章