DFT離散傅立葉變換
1.1傅立葉影象
對一張影象使用傅立葉變換就是將它分解成正弦和餘弦兩部分。也就是將影象從空間域(spatial domain)轉換到頻域(frequency domain)。 這一轉換的理論基礎來自於以下事實：任一函式都可以表示成無數個正弦和餘弦函式的和的形式。傅立葉變換就是一個用來將函式分解的工具。 2維影象的傅立葉變換可以用以下數學公式表達:
式中 f 是空間域(spatial domain)值， F 則是頻域(frequencydomain)值。 轉換之後的頻域值是複數，因此，顯示傅立葉變換之後的結果需要使用實數影象(real image) 加虛數影象(complex image), 或者幅度影象(magitude image)加相點陣圖像(phase image)。 在實際的影象處理過程中，僅僅使用了幅度影象，因為幅度影象包含了原影象的幾乎所有我們需要的幾何資訊。然而，如果你想通過修改幅度影象或者相點陣圖像的方法來間接修改原空間影象，你需要使用逆傅立葉變換得到修改後的空間影象，這樣你就必須同時保留幅度影象和相點陣圖像了。
在此示例中，我將展示如何計算以及顯示傅立葉變換後的幅度影象。由於數字影象的離散性，畫素值的取值範圍也是有限的。比如在一張灰度影象中，畫素灰度值一般在0到255之間。 因此，我們這裡討論的也僅僅是離散傅立葉變換(DFT)。 如果你需要得到影象中的幾何結構資訊，那你就要用到它了。請參考以下步驟(假設輸入影象為單通道的灰度影象 I)
1.2函式
(1)int getOptimalDFTSize(int vecsize)
該函式是為了獲得進行DFT計算的最佳尺寸。因為在進行DFT時，如果需要被計算的數字序列長度vecsize為2的n次冪的話，那麼其執行速度是非常快的。如果不是2的n次冪，但能夠分解成2,3,5的乘積，則運算速度也非常快。這裡的getOptimalDFTSize()函式就是為了獲得滿足分解成2,3,5的最小整數尺寸。很顯然，如果是多維矩陣需要進行DFT，則每一維單獨用這個函式獲得最佳DFT尺寸。
(2)void copyMakeBorder(InuptArray src, OutputArray dst, int top ,int bottom, int left, int right, int borderType, const Scalar&value=Scalar())
該函式是用來擴充套件一個影象的邊界的，第3～6個引數分別為原始影象的上下左右各擴充套件的畫素點的個數，第7個引數表示邊界的型別，如果其為BORDER_CONSTANT，則擴充的邊界畫素值則用第8個引數來初始化。將src影象擴充邊界後的結果儲存在dst影象中。
(3)merge()函式是把多個但通道陣列連線成1個多通道陣列，而split()函式則相反，把1個多通道函式分解成多個但通道函式。
(4)void magnitude(InputArray x, InputArray y, OutPutArray magnitude)
該函式是計算輸入矩陣x和y對應該的每個畫素平方求和後開根號儲存在輸出矩陣magnitude中。
(5)函式log(InputArraysrc, OutputArray dst)是對輸入矩陣src中每個畫素點求log，儲存在輸出矩陣dst的相對應的位置上。
1.3過程描述
將影象延擴到最佳尺寸。離散傅立葉變換的執行速度與圖片的尺寸息息相關。當影象的尺寸是2， 3，5的整數倍時，計算速度最快。因此，為了達到快速計算的目的，經常通過添湊新的邊緣畫素的方法獲取最佳影象尺寸。函式 getOptimalDFTSize() 返回最佳尺寸，而函式copyMakeBorder() 填充邊緣畫素:
Mat padded; //將輸入影象延擴到最佳的尺寸
int m = getOptimalDFTSize( I.rows );
int n = getOptimalDFTSize( I.cols ); // 在邊緣新增0
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));新增的畫素初始化為0.
為傅立葉變換的結果(實部和虛部)分配儲存空間。傅立葉變換的結果是複數，這就是說對於每個原影象值，結果是兩個影象值。此外，頻域值範圍遠遠超過空間值範圍，因此至少要將頻域儲存在 float 格式中。結果我們將輸入影象轉換成浮點型別，並多加一個額外通道來儲存複數部分：
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI); // 為延擴後的影象增添一個初始化為0的通道進行離散傅立葉變換。支援影象原地計算 (輸入輸出為同一影象):
dft(complexI, complexI); // 變換結果很好的儲存在原始矩陣中將複數轉換為幅度.複數包含實數部分(Re)和複數部分 (imaginary - Im)。 離散傅立葉變換的結果是複數，對應的幅度可以表示為:
轉化為OpenCV程式碼:
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];對數尺度(logarithmicscale)縮放. 傅立葉變換的幅度值範圍大到不適合在螢幕上顯示。高值在螢幕上顯示為白點，而低值為黑點，高低值的變化無法有效分辨。為了在螢幕上凸顯出高低變化的連續性，我們可以用對數尺度來替換線性尺度:
轉化為OpenCV程式碼:
magI += Scalar::all(1); // 轉換到對數尺度
log(magI, magI);剪下和重分佈幅度圖象限. 還記得我們在第一步時延擴了影象嗎? 那現在是時候將新新增的畫素剔除了。為了方便顯示，我們也可以重新分佈幅度圖象限位置(注：將第五步得到的幅度圖從中間劃開得到四張1/4子影象，將每張子影象看成幅度圖的一個象限，重新分佈即將四個角點重疊到圖片中心)。 這樣的話原點(0,0)就位移到影象中心。
int cx = magI.cols/2;
int cy = magI.rows/2;
Mat q0(magI, Rect(0, 0, cx, cy)); // Top-Left - 為每一個象限建立ROI
Mat q1(magI, Rect(cx, 0, cx, cy)); // Top-Right
Mat q2(magI, Rect(0, cy, cx, cy)); // Bottom-Left
Mat q3(magI, Rect(cx, cy, cx, cy)); //Bottom-Right
Mat tmp; // 交換象限 (Top-Left with Bottom-Right)
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp); // 交換象限 (Top-Right with Bottom-Left)
q2.copyTo(q1);
tmp.copyTo(q2)7.歸一化. 這一步的目的仍然是為了顯示。 現在我們有了重分佈後的幅度圖，但是幅度值仍然超過可顯示範圍[0,1] 。我們使用 normalize()函式將幅度歸一化到可顯示範圍。
normalize(magI, magI, 0, 1, CV_MINMAX); // 將float型別的矩陣轉換到可顯示影象範圍
// (float [0， 1]).1.4完整程式
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
#include <fstream>
using namespace cv;
using namespace std;
int main(int argc, char ** argv)
{
//以灰度影象模式讀入原始影象
const char* filename ="4.jpg";
Mat I = imread(filename, CV_LOAD_IMAGE_GRAYSCALE);
if( I.empty())
return -1;
//將影象擴充套件到最佳尺寸,邊界用0補充
Mat padded; //expand input image to optimal size
//將新增的畫素初始化為0
int m = getOptimalDFTSize( I.rows );//計算最佳擴充尺寸
int n = getOptimalDFTSize( I.cols ); // on the border add zero values
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0)); //擴充影象
//為傅立葉變換的結果(實部和虛部)分配儲存空間
//將planes陣列組合合併成一個多通道的陣列complexI
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)}; //新增維度，用於儲存傅立葉變換的結果
Mat complexI;
merge(planes, 2, complexI); // Add to the expanded another plane with zeros //合併通道
//進行傅立葉變換
dft(complexI, complexI); // this way the result may fit in the source matrix //離散傅立葉變換
//將複數轉換為幅值，即=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
// compute the magnitude and switch to 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)) //將儲存在complexI的結果分解到planes[0],planes[1]中 // 將多通道陣列complexI分離成幾個單通道陣列，planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude //計算複製
Mat magI = planes[0];
//進行對數尺度縮放
magI += Scalar::all(1); // switch to logarithmic scale //用對數表示
log(magI, magI); //求自然對數
//剪下和充分不幅度圖象限
//若有奇數行或者奇數列，進行頻譜裁剪
// crop the spectrum, if it has an odd number of rows or columns
//按位與，x & -2大概意思就是取不大於x 的最大偶數
magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));
//重新排列傅立葉影象中的象限，使得原點位於影象的中心
// rearrange the quadrants of Fourier image so that the origin is at the image center
int cx = magI.cols/2;
int cy = magI.rows/2;
Mat q0(magI, Rect(0, 0, cx, cy)); // Top-Left - Create a ROI per quadrant
Mat q1(magI, Rect(cx, 0, cx, cy)); // Top-Right
Mat q2(magI, Rect(0, cy, cx, cy)); // Bottom-Left
Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right
//交換象限(左上與右下)
Mat tmp; // swap quadrants (Top-Left with Bottom-Right)
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
//交換象限(右上與左下）
q1.copyTo(tmp); // swap quadrant (Top-Right with Bottom-Left)
q2.copyTo(q1);
tmp.copyTo(q2);
//歸一化
normalize(magI, magI, 0, 1, CV_MINMAX); // Transform the matrix with float values into a
// viewable image form (float between values 0 and 1).
Mat invDFT,invDFTcvt;
idft(complexI, invDFT, DFT_SCALE | DFT_REAL_OUTPUT );//離散傅立葉逆變換
invDFT.convertTo(invDFTcvt, CV_8U);
imshow("invDFTcvt", invDFTcvt);
imshow("Input Image" , I ); // Show the result
imshow("spectrum magnitude", magI);
waitKey();
return 0;
}