1 #include<opencv2/opencv.hpp>
  2 #include<iostream>
  3 using namespace  std ;
  4 using namespace  cv ;
  5 Mat gaussianlbrf( Mat scr, float sigma);//高斯低通滤波器函数
  6 Mat freqfilt( Mat scr,  Mat blur);//频率域滤波函数
  7 static void help(char* progName)
  8 {
  9 cout <<  endl
 10 << "This program demonstrated the use of the discrete Fourier
 11 transform (DFT). " <<  endl
 12 << "The dft of an image is taken and it's power spectrum is
 13 displayed." <<  endl
 14 << "Usage:" <<  endl
 15 << progName << " [image_name -- default ../data/lena.jpg] " <<  endl <<
 16 endl ;
 17 }
 18 int  main ()
 19 {
 20 const char* filename = "1.jpg";
 21 Mat input =  imread (filename,  IMREAD_GRAYSCALE );
 22 if (input. empty ())
 23 return -1;
 24 imshow ("input", input);//显示原图
 25 int w =  getOptimalDFTSize (input. cols ); //获取进行dtf的最优尺寸
 26 int h =  getOptimalDFTSize (input. rows ); //获取进行dtf的最优尺寸
 27 Mat padded;
 28 copyMakeBorder (input, padded, 0, h - input. rows , 0, w - input. cols ,
 29 BORDER_CONSTANT ,  Scalar :: all (0));//边界填充
 30 padded. convertTo (padded,  CV_32FC1 ); //将图像转换为flaot型
 31 Mat gaussianKernel = gaussianlbrf(padded, 45);//高斯低通滤波器
 32 Mat out = freqfilt(padded, gaussianKernel);//频率域滤波
 33 out = out( Rect (0, 0, input. cols , input. rows ));
 34 imshow ("结果图 sigma=40", out);
 35 waitKey (0);
 36 return 0;
 37 }
 38 //*****************高斯低通滤波器***********************
 39 Mat gaussianlbrf( Mat scr, float sigma)
 40 {
 41 Mat gaussianBlur(scr. size (),  CV_32FC1 ); //,CV_32FC1
 42 //高斯函数参数,越小,频率高斯滤波器越窄,滤除高频成分越多,图像就越平滑
 43 
 44 float d0 = 2 * sigma*sigma;
 45     for (int i = 0; i < scr. rows ; i++)
 46     {
 47     for (int j = 0; j < scr. cols ; j++)
 48     {
 49     float d =  pow (float(i - scr. rows / 2), 2) +    pow (float(j -
 50     scr. cols / 2), 2);//分子,计算pow必须为float型
 51     gaussianBlur. at <float>(i, j) =  expf (-d / d0);//expf为以e为底求幂
 52     (必须为float型)
 53     }
 54     }
 55     imshow ("高斯低通滤波器", gaussianBlur);
 56     return gaussianBlur;
 57     }
 58     //*****************频率域滤波*******************
 59     Mat freqfilt( Mat scr,    Mat blur)
 60     {
 61     //***********************DFT*******************
 62     Mat plane[] = { scr,  Mat :: zeros (scr. size () ,    CV_32FC1 ) }; //创建通道,存储    dft后的实部与虚部(CV_32F,必须为单通道数)
 63     Mat complexIm;
 64     merge (plane, 2, complexIm);//合并通道 (把两个矩阵合并为一个2通道的Mat类容
 65     器)
 66     dft (complexIm, complexIm);//进行傅立叶变换,结果保存在自身
 67     //***************中心化********************
 68     split (complexIm, plane);//分离通道(数组分离)
 69     int cx = plane[0]. cols / 2; int cy = plane[0]. rows / 2;//以下的操作是移动
 70     图像 (零频移到中心)
 71     Mat part1_r(plane[0],  Rect (0, 0, cx, cy));//元素坐标表示为(cx,cy)
 72     Mat part2_r(plane[0],  Rect (cx, 0, cx, cy));
 73     Mat part3_r(plane[0],  Rect (0, cy, cx, cy));
 74     Mat part4_r(plane[0],  Rect (cx, cy, cx, cy));
 75     Mat temp;
 76     part1_r. copyTo (temp);//左上与右下交换位置(实部)
 77     part4_r. copyTo (part1_r);
 78     temp. copyTo (part4_r);
 79     part2_r. copyTo (temp);//右上与左下交换位置(实部)
 80     part3_r. copyTo (part2_r);
 81     temp. copyTo (part3_r);
 82     Mat part1_i(plane[1],  Rect (0, 0, cx, cy));//元素坐标(cx,cy)
 83     Mat part2_i(plane[1],  Rect (cx, 0, cx, cy));
 84     Mat part3_i(plane[1],  Rect (0, cy, cx, cy));
 85     Mat part4_i(plane[1],  Rect (cx, cy, cx, cy));
 86     part1_i. copyTo (temp);//左上与右下交换位置(虚部)
 87     part4_i. copyTo (part1_i);
 88     temp. copyTo (part4_i);
 89     part2_i. copyTo (temp);//右上与左下交换位置(虚部)
 90     part3_i. copyTo (part2_i);
 91     temp. copyTo (part3_i);
 92     //*****************滤波器函数与DFT结果的乘积****************
 93     Mat blur_r, blur_i, BLUR;
 94     multiply (plane[0], blur, blur_r); //滤波(实部与滤波器模板对应元素相乘)
 95     multiply (plane[1], blur, blur_i);//滤波(虚部与滤波器模板对应元素相乘)
 96     Mat plane1[] = { blur_r, blur_i };
 97     merge (plane1, 2, BLUR);//实部与虚部合并
 98     //*********************得到原图频谱图***********************************
 99     magnitude (plane[0], plane[1], plane[0]);//获取幅度图像,0通道为实部通道,1    为虚部,因为二维傅立叶变换结果是复数
100     plane[0] +=  Scalar :: all (1);//傅立叶变o换后的图片不好分析,进行对数处理,结
101     果比较好看
102     log (plane[0], plane[0]);// float型的灰度空间为[0,1])
103     normalize (plane[0], plane[0], 1, 0,  CV_MINMAX );//归一化便于显示
104     imshow ("原图像频谱图", plane[0]);
105     //******************IDFT*******************************
106     /*
107     Mat part111(BLUR,Rect(0,0,cx,cy)); //元素坐标(cx,cy)
108     Mat part222(BLUR,Rect(cx,0,cx,cy));
109     Mat part333(BLUR,Rect(0,cy,cx,cy));
110     Mat part444(BLUR,Rect(cx,cy,cx,cy));
111     part111.copyTo(temp); //左上与右下交换位置(虚部)
112     part444.copyTo(part111);
113     temp.copyTo(part444);
114     part222.copyTo(temp); //右上与左下交换位置
115     part333.copyTo(part222);
116     temp.copyTo(part333);
117     */
118     idft (BLUR, BLUR);//idft结果也为复数
119     split (BLUR, plane);//分离通道,主要获取通道
120     magnitude (plane[0], plane[1], plane[0]);//求幅值(模)
121     normalize (plane[0], plane[0], 1, 0,  CV_MINMAX );//归一化便于显示
122     return plane[0];//返回参数
123     }