傅立叶变换与逆变换opencv实现，及频谱居中对idft后处理的影响

一、背景说明

在图像处理的滤波器操作中，经常需要使用傅立叶变换将空间域转换为频率域。同时需要对傅立叶变换后的频谱图进行居中处理，以方便与滤波器进行操作。

在滤波器的操作结束后，还需要进行傅立叶逆变换转换为空间域以得到空间域中的结果图。

本文核心探讨的就是逆变换为空间域所需要注意的事项，如频谱图做了象限换位的居中操作之后是否需要将象限转换回去再进行逆变换。

单纯从理论分析怕有错漏，故以下直接实践出真知。

二、傅立叶变换与逆变换调用的函数

void dft(InputArray src, OutputArray dst, int flags = 0, int nonzeroRows = 0);
void idft(InputArray src, OutputArray dst, int flags = 0, int nonzeroRows = 0);

实际在idft中也是调用了dft。

参数说明

InputArray类型的src，输入图像，如Mat类型。 OutputArray类型的dst，输出图像。
int类型的flags，转换的标识符，默认为0，也就是使用标准的正向变换。
int类型的nonzeroRows，默认值为0。如果该值设为非零值，dft函数会将该值作为非零行的有效区间长度，只对非零行进行处理，提高计算效率。

第三个参数flags，有几种常见的标识符：

DFT_INVERSE：一维或二维逆变换。 DFT_SCALE：缩放比例标识符，输出的结果会以1/N进行缩放。
DFT_ROWS：对输入矩阵的每行进行正向或反向的变换，适用于处理多种矢量的场景，可减小资源开销。
DFT_COMPLEX_OUTPUT：一维或二维实数数组正变换。 DFT_REAL_OUTPUT：一维或二维复数数组反变换。

DFT_REAL_OUTPUT - 返回IDFT后的Mat数组序列为0的Mat。

三、实现方式

1.dft与idft最简单的处理方式：dft后不改变频谱直接idft

//first,border preprocess
Mat img = imread("/home/seeking/Pictures/222.png",0);
int M = getOptimalDFTSize( img.rows );
int N = getOptimalDFTSize( img.cols ); 
copyMakeBorder(img, padded, 0, M - img.rows, 0, N - img.cols, BORDER_CONSTANT, Scalar::all(0));

//second,dft
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexImg;
merge(planes, 2, complexImg);
dft(complexImg, complexImg);   

//third,idft
Mat ifft;  
idft(complexImg,ifft,DFT_REAL_OUTPUT);  

//last,get result and input image,check.
normalize(ifft,ifft,0,1,CV_MINMAX);
normalize(ifft,result,0,255,CV_MINMAX);
result.convertTo(result,CV_8UC1);
normalize(padded,input,0,255,CV_MINMAX);
input.convertTo(input,CV_8UC1);

2.如果使用模的结果作为输出，居中的频谱可以不用转换回去

//first,border preprocess
Mat input = imread("/home/seeking/Pictures/222.png",0);
int M = getOptimalDFTSize( input.rows );
int N = getOptimalDFTSize( input.cols );
Mat padded;  
copyMakeBorder(input, padded, 0, M - input.rows, 0, N - input.cols, BORDER_CONSTANT, Scalar::all(0));

//second,dft
input = padded;
input.convertTo(input,CV_32FC1);
Mat plane[] = {input , Mat::zeros(input.size() , CV_32FC1)};
Mat dftMat;
merge(plane,2,dftMat);
dft(dftMat,dftMat);

//third, shift frequency domain
split(dftMat,plane);
int cx = plane[0].cols/2;
int cy = plane[0].rows/2;
Mat part1_r(plane[0],Rect(0,0,cx,cy));
Mat part2_r(plane[0],Rect(cx,0,cx,cy));
Mat part3_r(plane[0],Rect(0,cy,cx,cy));
Mat part4_r(plane[0],Rect(cx,cy,cx,cy));

Mat temp;
part1_r.copyTo(temp);
part4_r.copyTo(part1_r);
temp.copyTo(part4_r);

part2_r.copyTo(temp);
part3_r.copyTo(part2_r);
temp.copyTo(part3_r);

Mat part1_i(plane[1],Rect(0,0,cx,cy));
Mat part2_i(plane[1],Rect(cx,0,cx,cy));
Mat part3_i(plane[1],Rect(0,cy,cx,cy));
Mat part4_i(plane[1],Rect(cx,cy,cx,cy));

part1_i.copyTo(temp);
part4_i.copyTo(part1_i);
temp.copyTo(part4_i);

part2_i.copyTo(temp);
part3_i.copyTo(part2_i);
temp.copyTo(part3_i);

Mat BLUR;
merge(plane,2,BLUR);

//last, from frequency domain to spatial domain
//use magnitude for result
idft( BLUR, BLUR);
split(BLUR,plane);
magnitude(plane[0],plane[1],plane[0]);
normalize(plane[0],plane[0],1,0,CV_MINMAX);

normalize(plane[0], result, 0, 1, CV_MINMAX);
result.convertTo(result,CV_8UC1);
normalize(input, input, 0, 1, CV_MINMAX);
input.convertTo(input,CV_8UC1);

3.如果将居中的频率图转换回去则可以使用plane[0]，或者idft(complexImg,ifft,DFT_REAL_OUTPUT);

注意，需要用size为复数的图像作为输入。如果size为单数的话，在IDFT阶段就不能使用模去求解，其结果会不同，亲测。

但是，如果转换回去再进行IDFT结果就还是相同。

size转为复数可以使用：

input = input(Rect(0, 0, input.cols & -2, input.rows & -2));

不管怎样，最好对边界做前处理，能够获得最适合傅立叶变换的size，实现方法如下：

int M = getOptimalDFTSize( input.rows );
int N = getOptimalDFTSize( input.cols );
Mat padded;  
copyMakeBorder(input, padded, 0, M - input.rows, 0, N - input.cols, BORDER_CONSTANT, Scalar::all(0));