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[0]+x*src.step[1])

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

在主程序内容，参见注释。
**********************************************/
/*1.求取复数矩阵的幅值*/

void myMagnitude(Mat & complexImg,Mat & magnitudeImage)
{
    Mat planes[2];
    split(complexImg,planes);
    magnitude(planes[0],planes[1],magnitudeImage);
}

/*傅里叶变换后的频谱图后处理，将傅里叶普的原点(0,0)平移到图像的中心*/
void dftshift(Mat& ds)
{
    int cx=ds.cols/2;//图像的中心点x 坐标
    int cy=ds.rows/2;//图像的中心点y 坐标
    Mat q0=ds(Rect(0,0,cx,cy));//左上
    Mat q1=ds(Rect(cx,0,cx,cy));//右上
    Mat q2=ds(Rect(0,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=0;i<src.rows;i++)
        {
            for(int j=0;j<src.cols;j++)
            {
                float* mv=CV_MAT_ELEM2(src,float,i,j);
                if((i+j)%2!=0)
                {
                    for(int c=0;c<src.channels();c++)//对所有通道同样操作
                        mv[c]=-mv[c];//如果i+j为奇数，该像素值取负值
                }
            }
        }
        break;
    case CV_64F:
        for(int i=0;i<src.rows;i++)
        {
            for(int j=0;j<src.cols;j++)
            {
                double* mv=CV_MAT_ELEM2(src,double,i,j);
                if((i+j)%2!=0)
                {
                    for(int c=0;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(1);
        log(magI,magI);
    }
    normalize(magI,magI,0,1,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/2;
    int cy=filter.rows/2;
    switch(filterform)
    {
    case LOW_PASS_FILTER:
        for(int i=0;i<filter.rows;i++)
        {
            Vec2d* pf=filter.ptr<Vec2d>(i);
            for(int j=0;j<filter.cols;j++)
            {
                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);//rs表示r的平方
                pf[j][0]=1./(1.+pow(rs/Rs,n));//Rs是R的平方,
                pf[j][1]=pf[j][0];
            }
        }
        break;
    case HIGH_PASS_FILTER:
        for(int i=0;i<filter.rows;i++)
        {
            Vec2d* pf=filter.ptr<Vec2d>(i);
            for(int j=0;j<filter.cols;j++)
            {
                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);
                double Lp=1./(1.+pow(rs/Rs,n));//巴特沃斯公式
                pf[j][0]=1.0-Lp;
                pf[j][1]=pf[j][0];
            }
        }
        break;
    case BAND_STOP_FILTER:
        for(int i=0;i<filter.rows;i++)
        {
            Vec2d* pf=filter.ptr<Vec2d>(i);
            for(int j=0;j<filter.cols;j++)
            {
                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);
                double r=std::sqrt(rs);//r相当于书上的D
                pf[j][0]=1./(1.+pow(r*W/(rs-Rs),2*n));
                pf[j][1]=pf[j][0];
            }
        }
        break;
    case BAND_PASS_FILTER:
        for(int i=0;i<filter.rows;i++)
        {
            Vec2d* pf=filter.ptr<Vec2d>(i);
            for(int j=0;j<filter.cols;j++)
            {
                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);
                double r=std::sqrt(rs);//r相当于书上的D
                pf[j][0]=1-1./(1.+pow(r*W/(rs-Rs),2*n));
                pf[j][1]=pf[j][0];
            }
        }
        break;
    }
///***************【显示滤波器，如果不需要显示，将代码注销】*************/
    Mat displayFilter;
    extractChannel(filter,displayFilter,1);
    imshow("btw filter image",displayFilter);

}


void myDft(Mat&src,Mat&dst,bool isProCentralized=false,bool doubleSizeOrNot=false)
{
    CV_Assert(src.channels()==1);//验证src是否是单通道
    int ny=src.rows,nx=src.cols;
    Mat srcPadded;
    if(doubleSizeOrNot)//如果doubleOrNot为真，对src尺度扩展一倍
    {
        cv::copyMakeBorder(src,srcPadded,0,ny,0,nx,BORDER_CONSTANT);
    }
    else//否则，将src填补为最优傅里叶变换尺寸
    {
        int padx=getOptimalDFTSize(nx);//获取最优傅里叶变换列数
        int pady=getOptimalDFTSize(ny);//获取最优傅里叶变换行数数
        cv::copyMakeBorder(src,srcPadded,0,pady-ny,0,padx-nx,BORDER_CONSTANT);
    }
    if(srcPadded.type()!=CV_64FC1)
        srcPadded.convertTo(srcPadded,CV_64FC1);//转成浮点数据类型
    //如果isProCentralized为真，则在傅里叶变换前，对src做中心化预处理,
    //这样在傅里叶变换后的频谱图的（0，0）点就会位于频谱图的中心
    if(isProCentralized)
        srcCentralized(srcPadded);

    Mat planes[2]={srcPadded,Mat::zeros(srcPadded.rows,srcPadded.cols,srcPadded.type())};
    Mat srcComplex;
    merge(planes,2,srcComplex);
    dft(srcComplex,dst,DFT_COMPLEX_OUTPUT,0);//离散傅立叶变换
    cout<<"srcPadded.rows="<<srcPadded.rows<<"  "<<"src.cols="<<srcPadded.cols<<endl;
}

int main()
{
    ///读入灰度图像
    Mat src=imread("D:\Qt\MyImage\Fig0333(a).tif",0);

    ///***************【调用自定义的傅里叶变换myDft()】*************/
    Mat dftDst;//预声明dft的输出结果矩阵
    myDft(src,dftDst,0,0);//调用自定义的dft函数，对src执行傅里叶变换，dftDst是傅里叶变换结果
    dftshift(dftDst);//将傅里叶变换的结果，四象限对角互换
    imshowComplexMat(dftDst,"dftSpectrum display",false);//显示傅里叶频谱图

    ///***************【创建滤波器】*************/
    Mat filter(dftDst.rows,dftDst.cols,CV_64FC2,Scalar(0,0));
    //巴特沃斯滤波器的阶数n=2，半径分别取R=10，30，60，160，460
    createFilterButterworth(filter,2,100,60,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",1);
    imshow("srcimage",src);

    waitKey();
    return 0;
}

       

      

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

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