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.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