Opencv 实现图像的离散傅里叶变换（DFT）、卷积运算（相关滤波）

 http://blog.csdn.net/ikerpeng/article/details/41845545

我是做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)的定义域相交的元素进行乘积并且求和,得出新的图像一点,就是被卷积后的图像. 模版又称为卷积核.卷积核做一个矩阵的形状.（当然边缘点可能需要特殊的处理，同时这个操作和滤波也很像，也许就是一回事）。

[cpp] view
plain copy

 print?





#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);  

  

}  

其中：

[cpp] view
plain copy

 print?





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.opencv.org.cn/opencvdoc/2.3.2/html/doc/tutorials/core/discrete_fourier_transform/discrete_fourier_transform.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