opencv的頻域濾波

下面是頻域濾波示例程序：

在本程序中，共有五個自定義函數，分別是：

1. myMagnitude（），在該函數中封裝了Opencv中的magnitude函數，實現對於復數圖像的幅值計算。

2. dftshift()，該函數實現對圖像四個象限的對角互換，相當於MatLab中 fftshift（），將頻譜的原點（0，0）移到圖像中心。示例1中采用了該函數實現了頻譜圖中心化。

3. srcCentralized()用於傅里葉變換前的預處理，以便得到傅里葉頻譜的原點（0，0）位於圖像的中心。

該函數與dftshift（）目的一致，實現方法不同，一個是變換前預處理，一個是變換后處理。

示例2中采用了該函數實現了頻譜圖中心化。

4. displayDftSpectrum()，該函數用於顯示復數圖像。

5. makeFilter（），用於制作頻域濾波器，該函數利用了ptr（）

   指針遍歷圖像的方法，實現了圓等濾波函數。

下面是兩個示例，分別采用了dftshift()、srcCentralized()實現頻譜圖的中心化。示例2與岡薩雷斯的《數字圖像處理（第3版）》的4.6.7小結中的流程一致。

在主程序內容，參見注釋。
示例1：在該程序中采用了第2個函數dftshift()，實現了頻譜圖中心化。

#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;
}

 

 示例2：在該示例中采用了第3個函數srcCentralized()，實現頻譜圖中心化。程序流程即

在本教材的6.2.3頻域濾波步驟小結中，或者在岡薩雷斯的《數字圖像處理》4.6.7節中，頻域濾波流程總結如下：

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)。

/*求取復數矩陣的幅值*/
void myMagnitude(Mat & complexImg,Mat & mI)
{
    Mat planes[2];
    split(complexImg,planes);
    magnitude(planes[0],planes[1],mI);
}
/*傅里葉變換后的頻譜圖后處理，將傅里葉普的原點(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)
{
    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)
                mv[0]=-mv[0];//如果i+j為奇數，該像素值取負值
        }
    }
}
/*在窗口中顯示復數圖像，如果是正向傅里葉矩陣，需要取log才能顯示更多頻譜信息
 *如果是逆傅里葉變換，通過normalize歸一化后，顯示頻譜圖*/
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\\interferometer2.jpg",0);
    int ny=src.rows,nx=src.cols;
    imshow("original image",src);
    src.convertTo(src,CV_32FC1);
    srcCentralized(src);
    cv::copyMakeBorder(src,src,0,ny,0,nx,BORDER_CONSTANT);

    /***************傅里葉變換*************/
    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,30,true);
    Mat fftTemp=dftDst.mul(filter);//復數頻譜圖與濾波器相乘，實現頻域濾波
    srcCentralized(fftTemp);

    Mat srcFiltered;
    dft(fftTemp,srcFiltered,DFT_INVERSE);//逆傅里葉變換
    //  顯示經過濾波后的圖像
    displayDftSpectrum(srcFiltered,"filtered image",true);//顯示逆傅里葉頻譜圖
    waitKey();
    return 0;
}