在本程序中，共有六個自定義函數，分別是：
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時的濾波結果：
