一、前言
接触了图像的傅里叶变换，数学原理依旧不是很懂，因此不敢在这里妄言。下午用Opencv代码实现了这一变换，有一些经验心得
二、关键函数解析
2.1copyMakeBorder() 扩展图片尺寸
傅里叶变换的计算对图像的尺寸有一定要求，尺寸不满足要求的，可用copyMakeBorder() 函数进行扩展。函数定义如下：
void copyMakeBorder(InputArray src, //输入图像
OutputArray dst, //输出图像
int top, //上边界添加的像素行数
int bottom, //下边界添加的像素行数
int left, //左边界添加的像素列数
int right, //右边界添加的像素列数
int borderType, //表示边界的类型
const Scalar& value=Scalar()//表示如果边界的类型是BORDER_CONSTANT时边界的颜色值 )
borderType边界的类型有以下几种:
1)BORDER_REPLICATE:重复： aaaaaa|abcdefgh|hhhhhhh
2)BORDER_REFLECT:反射: fedcba|abcdefgh|hgfedcb
3)BORDER_REFLECT_101:反射101: gfedcb|abcdefgh|gfedcba
4)BORDER_WRAP:外包装： cdefgh|abcdefgh|abcdefg
5)BORDER_CONSTANT:常量复制： iiiiii|abcdefgh|iiiiiii(i的值由最后一个参数 const Scalar& value=Scalar()确定，如Scalar::all(0) )
2.2getOptimalDFTSize() 获取最佳计算尺寸
离散傅里叶变换的计算速度与图片的尺寸有关，当图片的尺寸是2，3，5的倍数时，计算速度最快。常与函数copyMakeBorder() 配合使用，保证较快的计算速度。
int getOptimalDFTSize(int vecsize)使用示例：
int m = getOptimalDFTSize(mImage.rows);//返回行的最佳尺寸
int n = getOptimalDFTSize(mImage.cols);//返回列的最佳尺寸
copyMakeBorder(mImage, mImage, 0, m - mImage.rows, 0, n - mImage.cols, BORDER_CONSTANT, Scalar(0));2.3dft()傅里叶变换计算
傅里叶变换计算函数。
void dft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0);参数解释：
InputArray src: 输入图像，可以是实数或虚数
OutputArray dst: 输出图像，其大小和类型取决于第三个参数flags
int flags = 0: 转换的标识符，有默认值0。其可取的值如下所示：
1）DFT_INVERSE: 用一维或二维逆变换取代默认的正向变换 ；
2）DFT_SCALE: 缩放比例标识符，根据数据元素个数平均求出其缩放结果，如有N个元素，则输出结果以1/N缩放输出，常与DFT_INVERSE搭配使用；
3）DFT_ROWS: 对输入矩阵的每行进行正向或反向的傅里叶变换；此标识符可在处理多种适量的的时候用于减小资源的开销，这些处理常常是三维或高维变换等复杂操作；
4）DFT_COMPLEX_OUTPUT: 对一维或二维的实数数组进行正向变换，这样的结果虽然是复数阵列，但拥有复数的共轭对称性（CCS），可以以一个和原数组尺寸大小相同的实数数组进行填充，这是最快的选择也是函数默认的方法。你可能想要得到一个全尺寸的复数数组（像简单光谱分析等等），通过设置标志位可以使函数生成一个全尺寸的复数输出数组；
5）DFT_REAL_OUTPUT: 对一维二维复数数组进行逆向变换，这样的结果通常是一个尺寸相同的复数矩阵，但是如果输入矩阵有复数的共轭对称性（比如是一个带有DFT_COMPLEX_OUTPUT标识符的正变换结果），便会输出实数矩阵。
int nonzeroRows = 0: 当这个参数不为0，函数会假设只有输入数组（没有设置DFT_INVERSE）的第一行或第一个输出数组（设置了DFT_INVERSE）包含非零值。这样的话函数就可以对其他的行进行更高效的处理节省一些时间，这项技术尤其是在采用DFT计算矩阵卷积时非常有效。
2.4magnitude()计算二维矢量幅值
void magnitude(InputArray x, InputArray y, OutputArray magnitude)；参数解释：
InputArray x: 浮点型数组的x坐标矢量，也就是实部
InputArray y: 浮点型数组的y坐标矢量，必须和x尺寸相同
OutputArray magnitude: 与x类型和尺寸相同的输出数组
其计算公式如下： 
2.5log()自然对数计算
log()函数的功能是计算每个数组元素绝对值的自然对数 。
void log(InputArray src,OutputArray dst)
参数解释：
InputArray src：为输入图像
OutputArray dst：为得到的对数值
其原理如下： 
2.6normalize()矩阵归一化
归一化就是把要处理的数据经过某种算法的处理限制在所需要的范围内。首先归一化是为了后面数据处理的方便，其次归一化能够保证程序运行时收敛加快。
void normalize(InputArray src, OutputArray dst,
double alpha=1, double beta=0, int norm_type=NORM_L2, int dtype=-1, InputArray mask=noArray() )参数解释：
InputArray src: 输入图像 ;
OutputArray dst: 输出图像，尺寸大小和src相同 ;
double alpha = 1: range normalization模式的最小值 ;
double beta = 0: range normalization模式的最大值，不用于norm normalization(范数归一化)模式 ;
int norm_type = NORM_L2: 归一化的类型，主要有 :
1)NORM_INF: 归一化数组的C-范数（绝对值的最大值）
2)NORM_L1: 归一化数组的L1-范数（绝对值的和）
3)NORM_L2: 归一化数组的L2-范数（欧几里得）
4)NORM_MINMAX: 数组的数值被平移或缩放到一个指定的范围，线性归一化，一般较常用。
int dtype = -1: 当该参数为负数时，输出数组的类型与输入数组的类型相同，否则输出数组与输入数组只是通道数相同，而depth = CV_MAT_DEPTH(dtype) ;
InputArray mask = noArray(): 操作掩膜版，用于指示函数是否仅仅对指定的元素进行操作。
归一化公式：
1)norm_type!=NORM_MINMAX(线性函数转换)：
if mask(i,j)!=0
dst(i,j)=(src(i,j)-min(src))*(b'-a')/(max(src)-min(src))+ a'
else
dst(i,j)=src(i,j)
其中b'=MAX(a,b), a'=MIN(a,b);
2)norm_type!=NORM_MINMAX:
if mask(i,j)!=0
dst(i,j)=src(i,j)*a/norm (src,norm_type,mask)
else
dst(i,j)=src(i,j)
其中，函数norm的功能是计算norm（范数）的绝对值
三、代码及结果分享
#include <opencv2/opencv.hpp>
#include <iostream>
using namespace std;
using namespace cv;
int main()
{
Mat mImage = imread("Lenna.jpg", 0);
if (mImage.data == 0)
{
cerr << "Image reading error" << endl;
system("pause");
return -1;
}
namedWindow("The original image", WINDOW_NORMAL);
imshow("The original image", mImage);
//Extending image
int m = getOptimalDFTSize(mImage.rows);
int n = getOptimalDFTSize(mImage.cols);
copyMakeBorder(mImage, mImage, 0, m - mImage.rows, 0, n - mImage.cols, BORDER_CONSTANT, Scalar(0));
//Fourier transform
Mat mFourier(mImage.rows + m, mImage.cols + n, CV_32FC2, Scalar(0, 0));
Mat mForFourier[] = { Mat_<float>(mImage), Mat::zeros(mImage.size(), CV_32F) };
Mat mSrc;
merge(mForFourier, 2, mSrc);
dft(mSrc, mFourier);
//channels[0] is the real part of Fourier transform,channels[1] is the imaginary part of Fourier transform
vector<Mat> channels;
split(mFourier, channels);
Mat mRe = channels[0];
Mat mIm = channels[1];
//Calculate the amplitude
Mat mAmplitude;
magnitude(mRe, mIm, mAmplitude);
//Logarithmic scale
mAmplitude += Scalar(1);
log(mAmplitude, mAmplitude);
//The normalized
normalize(mAmplitude, mAmplitude, 0, 255, NORM_MINMAX);
Mat mResult(mImage.rows, mImage.cols, CV_8UC1, Scalar(0));
for (int i = 0; i < mImage.rows; i++)
{
uchar* pResult = mResult.ptr<uchar>(i);
float* pAmplitude = mAmplitude.ptr<float>(i);
for (int j = 0; j < mImage.cols; j++)
{
pResult[j] = (uchar)pAmplitude[j];
}
}
Mat mQuadrant1 = mResult(Rect(mResult.cols / 2, 0, mResult.cols / 2, mResult.rows / 2));
Mat mQuadrant2 = mResult(Rect(0, 0, mResult.cols / 2, mResult.rows / 2));
Mat mQuadrant3 = mResult(Rect(0, mResult.rows / 2, mResult.cols / 2, mResult.rows / 2));
Mat mQuadrant4 = mResult(Rect(mResult.cols / 2, mResult.rows / 2, mResult.cols / 2, mResult.rows / 2));
Mat mChange1 = mQuadrant1.clone();
//mQuadrant1 = mQuadrant3.clone();
//mQuadrant3 = mChange1.clone();
mQuadrant3.copyTo(mQuadrant1);
mChange1.copyTo(mQuadrant3);
Mat mChange2 = mQuadrant2.clone();
//mQuadrant2 = mQuadrant4.clone();
//mQuadrant4 = mChange2.clone();
mQuadrant4.copyTo(mQuadrant2);
mChange2.copyTo(mQuadrant4);
namedWindow("The Fourier transform", WINDOW_NORMAL);
imshow("The Fourier transform", mResult);
waitKey();
destroyAllWindows();
return 0;
}
四、注意事项
4.1Mat对象类型问题
为方便运算，在计算过程中，Mat对象中数据都是float型，将其归一化至0-255范围后仍保留有小数部分，造成无法用imshow()正常显示。最后结果必须用uchar强制转换。源代码中以下片段实现了这一过程：
normalize(mAmplitude, mAmplitude, 0, 255, NORM_MINMAX);
Mat mResult(mImage.rows, mImage.cols, CV_8UC1, Scalar(0));
for (int i = 0; i < mImage.rows; i++)
{
uchar* pResult = mResult.ptr<uchar>(i);
float* pAmplitude = mAmplitude.ptr<float>(i);
for (int j = 0; j < mImage.cols; j++)
{
pResult[j] = (uchar)pAmplitude[j];
}
}4.2clone()与copyTo()的差异问题
clone 是完全的深拷贝，在内存中申请新的空间。copyTo 也是深拷贝，但是否申请新的内存空间，取决于dst矩阵头中的大小信息是否与src一至，若一致则只深拷贝并不申请新的空间，否则先申请空间后再进行拷贝。
Mat A = Mat::ones(4,5,CV_32F);
Mat B = A.clone() //clone 是完全的深拷贝，在内存中申请新的空间，与Ａ独立
Mat C;
A.copyTo(C) //此处的Ｃ矩阵大小与Ａ大小不一致，则申请新的内存空间，并完成拷贝，等同于clone()
Mat D = A.col(１);
A.col(0).copyTo(D) //此处D矩阵大小与Ａ.col(0)大小一致，因此不会申请空间，而是直接进行拷贝，相当于把Ａ的第１列赋值给第二列因此在进行不同象限数据对换时，用copyTo()能成功对换数据，而用clone()则不能起到作用，原因就在于mQuadrant.clone()时重新开辟了内存空间，数据对换时对换的并非是原矩阵mResult中的数据。所以源代码中，以下片段能够正常对换数据，而注释掉的部分不能正常工作。
Mat mChange1 = mQuadrant1.clone();
//mQuadrant1 = mQuadrant3.clone();
//mQuadrant3 = mChange1.clone();
mQuadrant3.copyTo(mQuadrant1);
mChange1.copyTo(mQuadrant3);
Mat mChange2 = mQuadrant2.clone();
//mQuadrant2 = mQuadrant4.clone();
//mQuadrant4 = mChange2.clone();
mQuadrant4.copyTo(mQuadrant2);
mChange2.copyTo(mQuadrant4);//*************
在学习信号与系统或通信原理等课程里面可能对傅里叶变换有了一定的了解。我们知道傅里叶变换是把一个信号从时域变换到其对应的频域进行分析。如果有小伙伴还对傅里叶变换处于很迷糊的状态，请戳这里，非常通俗易懂。而在图像处理中也有傅里叶分析的概念，我这里给出在其官方指导文件opencv_tutorials中给出的解释。
傅里叶变换可以将一幅图片分解为正弦和余弦两个分量，换而言之，他可以将一幅图像从其空间域（spatial domain）转换为频域（frequency domain）。这种变换的思想是任何函数可以很精确的接近无穷个sin()函数和cos()函数的和。傅里叶变换提供了这种方法来达到这种效果。对于二位图像其傅里叶变换公式如下：
式中f(i, j)是图像空间域的值而F是频域的值。傅里叶转换的结果是复数，这也显示出了傅里叶变换是一副实数图像（real image）和虚数图像（complex image）叠加或者是幅度图像（magitude image）和相位图像（phase image）叠加的结果。在实际的图像处理算法中仅有幅度图像（magnitude image）图像能够用到，因为幅度图像包含了我们所需要的所有图像几何结构的信息。但是，如果想通过修改幅度图像或者相位图像来间接修改原空间图像，需要保留幅度图像和相位图像来进行傅里叶逆变换，从而得到修改后图像。
1.dft()
首先看一下opencv提供的傅里叶变换函数dft()，其定义如下：
C++: void dft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0);
1
参数解释：
. InputArray src: 输入图像，可以是实数或虚数
. OutputArray dst: 输出图像，其大小和类型取决于第三个参数flags
. int flags = 0: 转换的标识符，有默认值0.其可取的值如下所示：
。DFT_INVERSE: 用一维或二维逆变换取代默认的正向变换
。DFT_SCALE: 缩放比例标识符，根据数据元素个数平均求出其缩放结果，如有N个元素，则输出结果以1/N缩放输出，常与DFT_INVERSE搭配使用。
。DFT_ROWS: 对输入矩阵的每行进行正向或反向的傅里叶变换；此标识符可在处理多种适量的的时候用于减小资源的开销，这些处理常常是三维或高维变换等复杂操作。
。DFT_COMPLEX_OUTPUT: 对一维或二维的实数数组进行正向变换，这样的结果虽然是复数阵列，但拥有复数的共轭对称性（CCS），可以以一个和原数组尺寸大小相同的实数数组进行填充，这是最快的选择也是函数默认的方法。你可能想要得到一个全尺寸的复数数组（像简单光谱分析等等），通过设置标志位可以使函数生成一个全尺寸的复数输出数组。
。DFT_REAL_OUTPUT: 对一维二维复数数组进行逆向变换，这样的结果通常是一个尺寸相同的复数矩阵，但是如果输入矩阵有复数的共轭对称性（比如是一个带有DFT_COMPLEX_OUTPUT标识符的正变换结果），便会输出实数矩阵。
. int nonzeroRows = 0: 当这个参数不为0，函数会假设只有输入数组（没有设置DFT_INVERSE）的第一行或第一个输出数组（设置了DFT_INVERSE）包含非零值。这样的话函数就可以对其他的行进行更高效的处理节省一些时间，这项技术尤其是在采用DFT计算矩阵卷积时非常有效。
2. getOptimalDFTSize()
返回给定向量尺寸经过DFT变换后结果的最优尺寸大小。其函数定义如下：
C++: int getOptimalDFTSize(int vecsize);
1
参数解释：
int vecsize: 输入向量尺寸大小(vector size)
DFT变换在一个向量尺寸上不是一个单调函数，当计算两个数组卷积或对一个数组进行光学分析，它常常会用0扩充一些数组来得到稍微大点的数组以达到比原来数组计算更快的目的。一个尺寸是2阶指数（2,4,8,16,32…）的数组计算速度最快，一个数组尺寸是2、3、5的倍数（例如：300 = 5*5*3*2*2）同样有很高的处理效率。
getOptimalDFTSize()函数返回大于或等于vecsize的最小数值N，这样尺寸为N的向量进行DFT变换能得到更高的处理效率。在当前N通过p,q,r等一些整数得出N = 2^p*3^q*5^r.
这个函数不能直接用于DCT（离散余弦变换）最优尺寸的估计，可以通过getOptimalDFTSize((vecsize+1)/2)*2得到。
3.magnitude()
计算二维矢量的幅值，其定义如下：
C++: void magnitude(InputArray x, InputArray y, OutputArray magnitude)；
1
参数解释：
. InputArray x: 浮点型数组的x坐标矢量，也就是实部
. InputArray y: 浮点型数组的y坐标矢量，必须和x尺寸相同
. OutputArray magnitude: 与x类型和尺寸相同的输出数组
其计算公式如下：
4. copyMakeBorder()
扩充图像边界，其函数定义如下：
C++: void copyMakeBorder(InputArray src, OutputArray dst, int top, int bottom, int left, int right, int borderType, const Scalar& value=Scalar() )；
1
参数解释：
. InputArray src: 输入图像
. OutputArray dst: 输出图像，与src图像有相同的类型，其尺寸应为Size(src.cols+left+right, src.rows+top+bottom)
. int类型的top、bottom、left、right: 在图像的四个方向上扩充像素的值
. int borderType: 边界类型，由borderInterpolate()来定义，常见的取值为BORDER_CONSTANT
. const Scalar& value = Scalar(): 如果边界类型为BORDER_CONSTANT则表示为边界值
5. normalize()
归一化就是把要处理的数据经过某种算法的处理限制在所需要的范围内。首先归一化是为了后面数据处理的方便，其次归一化能够保证程序运行时收敛加快。归一化的具体作用是归纳同意样本的统计分布性，归一化在0-1之间是统计的概率分布，归一化在某个区间上是统计的坐标分布，在机器学习算法的数据预处理阶段，归一化也是非常重要的步骤。其定义如下：
C++: void normalize(InputArray src, OutputArray dst, double alpha=1, double beta=0, int norm_type=NORM_L2, int dtype=-1, InputArray mask=noArray() )
1
参数解释：
. InputArray src: 输入图像
. OutputArray dst: 输出图像，尺寸大小和src相同
. double alpha = 1: range normalization模式的最小值
. double beta = 0: range normalization模式的最大值，不用于norm normalization(范数归一化)模式
. int norm_type = NORM_L2: 归一化的类型，主要有
。NORM_INF: 归一化数组的C-范数（绝对值的最大值）
。NORM_L1: 归一化数组的L1-范数（绝对值的和）
。NORM_L2: 归一化数组的L2-范数（欧几里得）
。NORM_MINMAX: 数组的数值被平移或缩放到一个指定的范围，线性归一化，一般较常用。
. int dtype = -1: 当该参数为负数时，输出数组的类型与输入数组的类型相同，否则输出数组与输入数组只是通道数相同，而depth = CV_MAT_DEPTH(dtype)
. InputArray mask = noArray(): 操作掩膜版，用于指示函数是否仅仅对指定的元素进行操作。
示例程序：
#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
using namespace std;
using namespace cv;
int main()
{
Mat I = imread("lena.jpg", IMREAD_GRAYSCALE); //读入图像灰度图
//判断图像是否加载成功
if (I.empty())
{
cout << "图像加载失败!" << endl;
return -1;
}
else
cout << "图像加载成功!" << endl << endl;
Mat padded; //以0填充输入图像矩阵
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//填充输入图像I，输入矩阵为padded，上方和左方不做填充处理
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
Mat complexI;
merge(planes, 2, complexI); //将planes融合合并成一个多通道数组complexI
dft(complexI, complexI); //进行傅里叶变换
//计算幅值，转换到对数尺度(logarithmic scale)
//=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(complexI, planes); //planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
//即planes[0]为实部,planes[1]为虚部
magnitude(planes[0], planes[1], planes[0]); //planes[0] = magnitude
Mat magI = planes[0];
magI += Scalar::all(1);
log(magI, magI); //转换到对数尺度(logarithmic scale)
//如果有奇数行或列，则对频谱进行裁剪
magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));
//重新排列傅里叶图像中的象限，使得原点位于图像中心
int cx = magI.cols / 2;
int cy = magI.rows / 2;
Mat q0(magI, Rect(0, 0, cx, cy)); //左上角图像划定ROI区域
Mat q1(magI, Rect(cx, 0, cx, cy)); //右上角图像
Mat q2(magI, Rect(0, cy, cx, cy)); //左下角图像
Mat q3(magI, 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-1之间的浮点数将矩阵变换为可视的图像格式
normalize(magI, magI, 0, 1, CV_MINMAX);
imshow("输入图像", I);
imshow("频谱图", magI);
waitKey(0);
return 0;
}
程序分析：
1.图像填充：
Mat padded;
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
根据前面的理论介绍可以知道当图像尺寸为2、3、5的倍数时可以得到最快的处理速度，所以通过getOptimalDFTSize()函数获取最佳DFT变换尺寸，之后再结合copyMakeBorder()函数对图像进行扩充。
2.为实部和虚部分配存储空间
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
傅里叶变换的结果是复数，这就意味着经过傅里叶变换每个图像值都会变成两个值，此外其频域(frequency domains)范围比空间域(spatial counterpart)范围大很多。我们通常以浮点型数据格式对结果进行存储。因此我们将输入图像转换为这种类型，通过另外的通道扩充图像。
3.傅里叶变换
dft(complexI, complexI);
1
傅里叶变换函数，对图像进行傅里叶变换。
4.将实数和复数的值转换为幅度值
split(complexI, planes);
magnitude(planes[0], planes[1], planes[0]);
Mat magI = planes[0];
1
2
3
复数包含实部和虚部两个部分，傅里叶变换的结果是一个复数，傅里叶变换的幅度计算公式是：
5.转换为对数尺度(Switch to a logarithmic scale)
magI += Scalar::all(1);
log(magI, magI);
1
2
之所以要进行对数转换是因为傅里叶变换后的结果对于在显示器显示来讲范围比较大，这样的话对于一些小的变化或者是高的变换值不能进行观察。因此高的变化值将会转变成白点，而较小的变化值则会变成黑点。为了能够获得可视化的效果，可以利用灰度值将我们的线性尺度(linear scale)转变为对数尺度(logarithmic scale)，其计算公式如下：
6.剪切和象限变换
magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));
int cx = magI.cols / 2;
int cy = magI.rows / 2;
Mat q0(magI, Rect(0, 0, cx, cy));
Mat q1(magI, Rect(cx, 0, cx, cy));
Mat q2(magI, Rect(0, cy, cx, cy));
Mat q3(magI, 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);
在进行傅里叶变换时，为了取得更快的计算效果，对图像进行了扩充，现在就需要对新增加的行列进行裁剪了。为了可视化的需要，我们同样需要对显示的结果图像像素进行调整，如果不进行调整，最后显示的结果是这样的：
可以看到四周的角上时高频分量，现在我们通过重新调整象限将高频分量调整到图像正中间。
7.归一化
normalize(magI, magI, 0, 1, CV_MINMAX);
1
对结果进行归一化处理同样是处于可视化的目的。现在我们得到了幅度值，但是这仍然超出了0-1的显示范围。这就需要利用normalize()函数对数据进行归一化处理。
程序运行结果如图：
————————————————
版权声明：本文为CSDN博主「梧桐栖鸦」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/keith_bb/article/details/53389819
//***********
傅里叶变换将图像分解成其正弦和余弦分量，它将图像由空域转换为时域。任何函数都可以近似的表示为无数正弦和余弦函数的和，傅里叶变换就是实现这一步的，数学上一个二维图像的傅里叶变换为：
公式中，f是图像在空域的值，F是频域的值。转换的结果是复数，但是不可能通过一个真实图像和一个复杂的图像或通过大小和相位图像去显示这样的一个图像。然而，在整个图像处理算法只对大小图像是感兴趣的，因为这包含了所有我们需要的图像几何结构的信息。
可通过以下几步显示一副傅里叶变换后的图像
1、将图像扩展到它的最佳尺寸，DFT(直接傅里叶变换)的性能依赖于图片的尺寸，当图像是2，3，5的倍数时往往是最快的。因此，为了达到最优性能通常采用垫边界值的方法，得到一个最佳的尺寸。
2、为傅立叶变换结果的实部和虚部分配存储空间。傅里叶变换的结果是一个复数，这意味着每幅图的结果都有一个实部和虚部，此外，频域范围远远大于它对应的空间范围。因此，我们这些通常至少以一个浮点数格式存储这些数值。因此，我们会将我们的输入图像转换为这种类型并且扩展它与另一通道存放复数值
3、进行傅里叶变换。
4、将复数转换为幅值，DFT的幅值由以下公式得出：
正在上传…重新上传取消转存失败重新上传取消
5、切换到对数刻度。对图像进行对数尺度的缩放，结果证明，傅立叶系数矩阵的动态范围太大，无法显示在屏幕上，我们无法通过这样去观察一些小的和高的变化值。因此那些高的数值将转化成白点而小的数值会变成黑点，使用灰度值进行可视化，我们可以将线性刻度转换为对数刻度，以便于观察。
6、剪切和重分布幅度图象，第一步我们扩展了图像，这里我们去掉扩展的那部分值，基于可视化的目的，我们还可以重新排列结果的象限，使原点(0，0)对应于与图像中心
7、归一化。目前得到的幅值图像仍然太大，超出了显示的范围，归一化这范围内的值，可以进一步达到可视化的目的
实现程序
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 |
|
傅里叶变换后的图片
Opencv 实现图像的离散傅里叶变换（DFT）、卷积运算（相关滤波）
我是做Tracking 的，对于速度要求非常高。发现傅里叶变换能够使用。
于是学习之。
核心： 最根本的一点就是将时域内的信号转移到频域里面。这样时域里的卷积能够转换为频域内的乘积！
在分析图像信号的频率特性时，对于一幅图像，直流分量表示预想的平均灰度。低频分量代表了大面积背景区域和缓慢变化部分，高频部分代表了它的边缘，细节，跳跃部分以及颗粒噪声. 因此，我们能够做对应的锐化和模糊的处理：提出当中的高频分量做傅里叶逆变换得到的就是锐化的结果。
提出当中的低频分量做傅里叶逆变换得到的就是模糊的结果。
最不能理解的应该是：截取频域图中的不论什么一个区域相应的都是原来的整张图的区域。而不是相应的局部。
由于频域内的各个点都反映的是整张图的一个状态。
我们能够用时间和频率来理解：当你走完一段单位路程的时候。如果你花了100秒，那么你的频率就是0.01HZ。
这个0.01HZ显然体现的是一个总体的结果。而不是局部。
我们再由公式来看：
能够非常明显的知道频域内的每个点的值都是由整个图像求出来的。当然以上得出的结果，我们一般仅仅关注幅值频谱图。
也就是说真正起作用的就是前面的那个cos x而已. 于是我们能够知道。在整个范围内（0<k <N, 0<l <N）,低频分量集中于四个角。
且其它地方的值仅仅可能比这个小。
在原点的傅里叶变换即等于图像的平均灰度级。由于 在原点处经常为零，F(0，0)有时称做 频率谱的直流成分。
使用：
当图像的尺寸是2，3，5的整数倍时，计算速度最快。因此opencv里面有一个函数：
int m = getOptimalDFTSize( I.rows ); int n = getOptimalDFTSize( I.cols ); // 在边缘加入0
它能够使得图片的尺寸能够满足这个要求。
可是这样就须要对原来的图像进行大小的处理，因此使用函数：CopyMakeBorder复制图像而且制作边界。
（处理边界卷积）
Mat padded;
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
将原始的图像I 扩充为理想的大小放在padded里面。
接下来我们须要给计算出来的结果分配空间：
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI); // 为延扩后的图像增添一个初始化为0的通道
然后便能够进行傅里叶变换了：
dft(complexI, complexI); // 变换结果非常好的保存在原始矩阵中
得到的结果有两部分。实数部分和虚数部分，你能够分别对这两部分进行操作：
split(complexI, planes); // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I)) magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude Mat magI = planes[0];
当然还能够进行：归一化：
normalize(magI, magI, 0, 1, CV_MINMAX); // 将float类型的矩阵转换到可显示图像范围
// (float [0。 1]).
另外重要的一个应用是： convolveDFT。
当中的 *代表的是 卷积。我认为这也是我们进行离散傅里叶变换的目的。
使得计算的速度大大的添加。
先来说一下卷积在图像中的意义：
如果图像f(x),模板是g(x),然后将模版g(x)在模版中移动,每到一个位置,就把f(x)与g(x)的定义域相交的元素进行乘积而且求和,得出新的图像一点,就是被卷积后的图像. 模版又称为卷积核.卷积核做一个矩阵的形状.（当然边缘点可能须要特殊的处理，同一时候这个操作和滤波也非常像，或许就是一回事）。
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
using namespace cv;
using namespace std;
//http://docs.opencv.org/modules/core/doc/operations_on_arrays.html#dft[2]
void convolveDFT(Mat A, Mat B, Mat& C)
{
// reallocate the output array if needed
C.create(abs(A.rows - B.rows)+1, abs(A.cols - B.cols)+1, A.type());
Size dftSize;
// calculate the size of DFT transform
dftSize.width = getOptimalDFTSize(A.cols + B.cols - 1);
dftSize.height = getOptimalDFTSize(A.rows + B.rows - 1);
// allocate temporary buffers and initialize them with 0's
Mat tempA(dftSize, A.type(), Scalar::all(0));//initial 0
Mat tempB(dftSize, B.type(), Scalar::all(0));
// copy A and B to the top-left corners of tempA and tempB, respectively
Mat roiA(tempA, Rect(0,0,A.cols,A.rows));
A.copyTo(roiA);
Mat roiB(tempB, Rect(0,0,B.cols,B.rows));
B.copyTo(roiB);
// now transform the padded A & B in-place;
// use "nonzeroRows" hint for faster processing
dft(tempA, tempA, 0, A.rows);
dft(tempB, tempB, 0, B.rows);
// multiply the spectrums;
// the function handles packed spectrum representations well
mulSpectrums(tempA, tempB, tempA, DFT_COMPLEX_OUTPUT);
//mulSpectrums(tempA, tempB, tempA, DFT_REAL_OUTPUT);
// transform the product back from the frequency domain.
// Even though all the result rows will be non-zero,
// you need only the first C.rows of them, and thus you
// pass nonzeroRows == C.rows
dft(tempA, tempA, DFT_INVERSE + DFT_SCALE, C.rows);
// now copy the result back to C.
tempA(Rect(0, 0, C.cols, C.rows)).copyTo(C);
// all the temporary buffers will be deallocated automatically
}
int main(int argc, char* argv[])
{
const char* filename = argc >=2 ? argv[1] : "Lenna.png";
Mat I = imread(filename, CV_LOAD_IMAGE_GRAYSCALE);
if( I.empty())
return -1;
Mat kernel = (Mat_<float>(3,3) << 1, 1, 1, 1, 1, 1, 1, 1, 1);
cout << kernel;
Mat floatI = Mat_<float>(I);// change image type into float
Mat filteredI;
convolveDFT(floatI, kernel, filteredI);
normalize(filteredI, filteredI, 0, 1, CV_MINMAX); // Transform the matrix with float values into a
// viewable image form (float between values 0 and 1).
imshow("image", I);
imshow("filtered", filteredI);
waitKey(0);
}
当中：
C.create(abs(A.rows - B.rows)+1, abs(A.cols - B.cols)+1, A.type());
C 为什么是这种勒？想想一个特殊的样例就知道了：当A，B尺寸相等的时候，这个时候的高斯滤波得到的也就是中心点的那一个值（卷积核滤波的区别在于须要绕中心180度旋转）。
MulSpectrums 是对于两张频谱图中每个元素的乘法。
void cvMulSpectrums( const CvArr* src1, const CvArr* src2, CvArr* dst, int flags ); src1 第一输入数组 src2 第二输入数组 dst 输出数组，和输入数组有同样的类型和大小。 flags 以下列举的值的组合: CV_DXT_ROWS - 把数组的每一行视为一个单独的频谱 (參见 cvDFT 的參数讨论). CV_DXT_MUL_CONJ - 在做乘法之前取第二个输入数组的共轭.
第四个參数flag值没有指定，应指定为DFT_COMPLEX_OUTPUT或是DFT_REAL_OUTPUT.
參考资料：
http://blog.sina.com.cn/s/blog_4bdb170b01019atv.html
http://www.cnblogs.com/xianglan/archive/2010/12/30/1922386.html
http://www.cnblogs.com/tornadomeet/archive/2012/07/26/2610414.html
http://blog.csdn.net/ubunfans/article/details/24787569
http://blog.csdn.net/lichengyu/article/details/18848281
714
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
