通过傅里叶变换方法求图像卷积-OpenCV实现

https://blog.csdn.net/lichengyu/article/details/18848281

在图像处理中经常会遇到各种滤波(平滑、锐化)的情况，基本方法都是将图像与一个核进行卷积实现。而卷积定理指出，两个函数的卷积的傅里叶变换等于各自的傅里叶变换的乘积，即：

[1]

那么，两个函数的卷积可以通过如下方式得到，对两个函数傅里叶变换的乘积做傅里叶反变换，即：

[1]

 

在进行卷积运算时，一般是将核沿着图像从左到右从上到下计算每一个像素处与核卷积后的值，这样的计算量较大，采用傅里叶变换的方法可以提高运算效率。

 

1.	#include "opencv2/core/core.hpp"  
2.	#include "opencv2/imgproc/imgproc.hpp"  
3.	#include "opencv2/highgui/highgui.hpp"  
4.	#include <iostream>  
5.	  
6.	using namespace cv;  
7.	using namespace std;  
8.	  
9.	//http://docs.opencv.org/modules/core/doc/operations_on_arrays.html#dft[2]  
10.	void convolveDFT(Mat A, Mat B, Mat& C)  
11.	{  
12.	    // reallocate the output array if needed  
13.	    C.create(abs(A.rows - B.rows)+1, abs(A.cols - B.cols)+1, A.type());  
14.	    Size dftSize;  
15.	    // calculate the size of DFT transform  
16.	    dftSize.width = getOptimalDFTSize(A.cols + B.cols - 1);  
17.	    dftSize.height = getOptimalDFTSize(A.rows + B.rows - 1);  
18.	  
19.	    // allocate temporary buffers and initialize them with 0's  
20.	    Mat tempA(dftSize, A.type(), Scalar::all(0));  
21.	    Mat tempB(dftSize, B.type(), Scalar::all(0));  
22.	  
23.	    // copy A and B to the top-left corners of tempA and tempB, respectively  
24.	    Mat roiA(tempA, Rect(0,0,A.cols,A.rows));  
25.	    A.copyTo(roiA);  
26.	    Mat roiB(tempB, Rect(0,0,B.cols,B.rows));  
27.	    B.copyTo(roiB);  
28.	  
29.	    // now transform the padded A & B in-place;  
30.	    // use "nonzeroRows" hint for faster processing  
31.	    dft(tempA, tempA, 0, A.rows);  
32.	    dft(tempB, tempB, 0, B.rows);  
33.	  
34.	    // multiply the spectrums;  
35.	    // the function handles packed spectrum representations well  
36.	    mulSpectrums(tempA, tempB, tempA, DFT_COMPLEX_OUTPUT);  
37.	    //mulSpectrums(tempA, tempB, tempA, DFT_REAL_OUTPUT);  
38.	  
39.	    // transform the product back from the frequency domain.  
40.	    // Even though all the result rows will be non-zero,  
41.	    // you need only the first C.rows of them, and thus you  
42.	    // pass nonzeroRows == C.rows  
43.	    dft(tempA, tempA, DFT_INVERSE + DFT_SCALE, C.rows);  
44.	  
45.	    // now copy the result back to C.  
46.	    tempA(Rect(0, 0, C.cols, C.rows)).copyTo(C);  
47.	  
48.	    // all the temporary buffers will be deallocated automatically  
49.	}  
50.	  
51.	  
52.	int main(int argc, char* argv[])  
53.	{  
54.	    const char* filename = argc >=2 ? argv[1] : "Lenna.png";  
55.	  
56.	    Mat I = imread(filename, CV_LOAD_IMAGE_GRAYSCALE);  
57.	    if( I.empty())  
58.	        return -1;  
59.	  
60.	    Mat kernel = (Mat_<float>(3,3) << 1, 1, 1, 1, 1, 1, 1, 1, 1);  
61.	    cout << kernel;  
62.	  
63.	    Mat floatI = Mat_<float>(I);// change image type into float  
64.	    Mat filteredI;  
65.	    convolveDFT(floatI, kernel, filteredI);  
66.	      
67.	    normalize(filteredI, filteredI, 0, 1, CV_MINMAX); // Transform the matrix with float values into a  
68.	                                            // viewable image form (float between values 0 and 1).  
69.	    imshow("image", I);  
70.	    imshow("filtered", filteredI);  
71.	    waitKey(0);  
72.	  
73.	}  

 

convolveDFT函数是从官方文档中抄录并做了修改，因为原来的程序有问题。一是输出Mat C应声明为引用；二是其中的mulSpectrums函数的第四个参数flag值没有指定，应指定为DFT_COMPLEX_OUTPUT或是DFT_REAL_OUTPUT.


main函数中首先按灰度图读入图像，然后创造一个平滑核kernel，将输入图像转换成float类型(注意这步是必须的，因为dft只能处理浮点数)，在调用convolveDFT求出卷积结果后，将卷积结果归一化方便显示观看。


需要注意的是，一般求法中，利用核游走整个图像进行卷积运算，实际上进行的是相关运算，真正意义上的卷积，应该首先把核翻转180度，再在整个图像上进行游走。OpenCV中的filter2D实际上做的也只是相关，而非卷积。"The function does actually compute correlation, not the convolution: ... That is, the kernel is not mirrored around the anchor point. If you need a real convolution, flip the kernel using flip() and set the new anchor to (kernel.cols - anchor.x - 1, kernel.rows - anchor.y - 1)"[3]


参考文献:
[1] http://zh.wikipedia.org/wiki/%E5%8D%B7%E7%A7%AF%E5%AE%9A%E7%90%86
[2] http://docs.opencv.org/modules/core/doc/operations_on_arrays.html#dft
[3] http://docs.opencv.org/modules/imgproc/doc/filtering.html#filter2d