opencv的频域滤波

频域滤波流程总结如下：

1. 给定一幅大小为M×N的输入图像f(x,y)，从式（6.1-25）和式（6.1-26）得到填充参数P和Q。典型地，我们选择P=2M和Q=2N；
2. 对f(x,y)添加必要数量的0，形成大小为P×Q填充后的图像
3. 用（-1）^(x+y)乘以f_p(x,y)，进行频谱中心化的预处理；
4. 计算中心化预处理过的f_p(x,y)的傅里叶变换，得到F_p(u,v);
5. 生成一个实的、对称的滤波函数H(u,v)，其大小为P×Q，频谱零点位于(P/2,Q/2)处。用阵列相乘形成乘积G(u,v)=F(u,v)H(u,v)；
6. 经过频域滤波后的图像：
         g_p(x,y)={real[fft^-1[G(u,v)]]}（-1）^(x+y)
   其中，为忽略由于计算不准确导致的寄生复分量，选择了实部，下标P指出我们处理的是填充后的阵列；
7. 通过从 g_p(x,y)的左上象限提取M×N区域，得到最终处理结果g(x,y)。

下面是频域滤波示例程序：

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

1. myMagnitude（），在该函数中封装了Opencv中的

magnitude函数，实现对于复数图像的幅值计算。

2. dftshift()，该函数实现对图像四个象限的对角互换，

相当于MatLab中 fftshift（）。

3. displayDftSpectrum()，该函数用于显示复数图像。

4. makeFilter（），用于制作频域滤波器，该函数利用了ptr（）

   指针遍历图像的方法，实现了圆等滤波函数。

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

#include <iostream>
#include<opencv2/opencv.hpp>
using namespace cv;
using namespace std;

void myMagnitude(Mat & complexImg,Mat & mI)
{
    Mat planes[2];
    split(complexImg,planes);
    magnitude(planes[0],planes[1],mI);
}
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);
}
void displayDftSpectrum(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);
}
void makeFilter(Mat&filter,int R,bool isLowPassFilter)
{
    int cx=filter.cols/2,cy=filter.rows/2;
    int R2=R*R;

    for(int i=0;i<filter.rows;i++)
    {
        Vec2d* pf=filter.ptr<Vec2d>(i);
        for(int j=0;j<filter.cols;j++)
        {
            int r2=(j-cx)*(j-cx)+(i-cy)*(i-cy);
            if(r2<R2)
            {
                pf[j][0]=1;
                pf[j][1]=pf[j][0];
            }
        }
    }
    if(!isLowPassFilter)
    {
     filter =Scalar::all(1)-filter;
    }
    Mat displayFilter;
    extractChannel(filter,displayFilter,0);
    imshow("filter image",displayFilter);
}

int main()
{
    ///读入灰度图像
    Mat src=imread("D:\\Qt\\MyImage\\baboon.jpg",0);
    int ny=src.rows,nx=src.cols;
    cv::copyMakeBorder(src,src,0,ny,0,nx,BORDER_CONSTANT);
    imshow("original image",src);
    /***************傅里叶变换*************/
    src.convertTo(src,CV_64FC1);//转成浮点数据类型
    Mat dftDst(src.size(),CV_64FC2);//预设dft的输出结果矩阵
    dft(src,dftDst,DFT_COMPLEX_OUTPUT,0);//离散傅立叶变换
    dftshift(dftDst);//将傅里叶变换的结果，四象限对角互换
    displayDftSpectrum(dftDst,"dftSpectrum display",false);//显示傅里叶频谱图

    /***************频域滤波*************/
    Mat filter(src.size(),CV_64FC2,Scalar::all(0));
    makeFilter(filter,80,false);
    Mat fftTemp=dftDst.mul(filter);//复数频谱图与滤波器相乘，实现频域滤波
    Mat srcFiltered;
    dft(fftTemp,srcFiltered,DFT_INVERSE);//逆傅里叶变换
    //显示经过滤波后的图像
    displayDftSpectrum(srcFiltered,"filtered image",true);//显示逆傅里叶频谱图
    waitKey();
    return 0;
}