OpenCV 離散傅里葉變換

離散傅里葉變換（DFT）

定義

離散傅里葉變換（Discrete Fourier Transform，縮寫爲DFT），是傅里葉變換在時域和頻域上都呈離散的形式，將信號的時域採樣變換爲其DFT的頻域採樣。

對於N點序列{X[n]}(0 <= n <= N),它的離散傅里葉變換爲:

dft() 函數

dft()函數的做用是對一維或二維的浮點數數組進行正向或反向的離散傅里葉變換。

函數原型

void dft( InoutArray src, //輸入矩陣 OutputArray dst, ///輸出矩陣 int flags = 0, //轉換的標識符 int onozeroRows, ) 複製代碼

第三個參數，轉換的標識符分爲:

• DFT_INVERSE 用一維或二維逆變換代替默認的正向變換。
• DFT_SCALE 縮放比例標識符，輸出的結果都會以1 / N進行縮放。
• DFT_CMPLEX_OUTPUT、DFT_REAL_OUTPUT 進行一維或二維的數組反變換。

返回DFT最優尺寸大小：getOptimalDFTSize()函數

void int getOptimalDFTSize(int vecsize);
複製代碼

擴充圖像邊界：copyMakeBorder()函數

void copyMakeBorder( InputArray src, //輸入圖像 OutputArray dst, //輸出圖像 int top, //在圖像上方擴充的像素值 int bottom, //在圖像下方擴充的像素值 int left, //在圖像左方擴充的像素值 int right, //在圖像右方擴充的像素值 int borderType, //邊界類型· const Scalars, ) 複製代碼

計算二維矢量的幅值：magnitude()函數

用於計算二維矢量的幅值

void magnitude( InputArray x, //表示矢量的浮點型X座標值,即實部 InputArray y, //表示矢量的浮點型Y座標值，即虛部 OutputArray magnitude, //輸出的幅值 ) 複製代碼

計算天然對數：log() 函數

log()函數的功能是計算每一個數組元素絕對值的天然對數

void log( InputArray src, OutputArray dst, );

原理即爲：

if(src(I) != 0)
    log|src(I)|
else
    C
複製代碼

矩陣歸一化：normalize()函數

void normalize(
    InputArray src,
    OutputArray dst,
    double alpha = 1,  //歸一化以後的最大值，有默認值1
    double beta = 0,  //歸一化以後的最大值，有默認值0
    int norm_type = NORM_L2,  //歸一化類型
    int dtype = -1,   //爲負數時輸出矩陣和src有一樣的類型，不然，它和src有一樣的通道數，深度爲CV_MAT_DEPTH
    InputArray mask=noArray(),  //可選的操做掩模
)
複製代碼

getOptimalDFTSize函數

該函數返回給定向量尺寸的傅里葉最優尺寸大小。爲了提升離散傅里葉變換的運行速度，須要擴充圖像。

使用dft函數計算兩個二維實矩陣卷積的示例核心片斷

• 其中MulSpectrums的做用是計算兩個傅里葉頻譜的每一個元素的乘法

void convolveDft(InputArray A, InputArray B, OutputArray C)
{
    //初始化輸出矩陣
    C.create(abs(A.rows - B.rows) + 1, abs(A.cols - B.cols) + 1, A.type);

    //計算DFT變換的尺寸
    dftSize.width = getOptimalDFTSize(A.cols + B.cols - 1);
    dftSize.height = getOptimalDFTSize(A.rows + B.rows - 1);

    //分配臨時緩衝區並初始化置0
    Mat tempA(dftSize,A.type(),Scalar::all(0));
    Mat tempB(dftSize,B.type(),Scalar::all(0));

    //分別複製A和B到tempA和tempB的左上角
    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);

    //就地操做，進行快速傅里葉變換，並將nonzeroRows 參數置爲零，以進行更加快速的處理。
    dft(tempA,tempA,0,A.rows);
    dft(tempB,tempB,0,B.rows);

    //將獲得的頻譜相乘，結果存放於tempA當中
    mulSpectrums(tempA,tempB,tempA);  

    //將結果變換爲頻域，儘管結果行（result.rows）都爲非零，咱們只需其中的C.rows的第一行，因此採用nonzeroRows == C.rows
    dft(tempA,tempA,EFT_INVERSE + EFT_SCALE,C.rows);

    //將結果複製到C當中
    tempA(Rect(0,0,C.cols,C.rows)).copyTo(C);

    //全部的臨時緩存區將被自動釋放，因此無需收尾操做
}
複製代碼

例示程序：離散傅里葉變換

【1】載入原始圖像

//【1】ui灰度模式讀取原始圖像並顯示
Mat srcImage = imread("D:\\Desktop\\lena.jpg",0);
if (!srcImage.data)
{
    cout << "lena圖片讀取錯誤" << endl;
    return false;
}
imshow("原始圖像", srcImage);
複製代碼

【2】將圖像擴大到合適的尺寸

//【2】將輸入的圖像延擴到最佳的尺寸，邊界用0補充
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//將添加的像素初始化爲0
Mat padded;
copyMakeBorder(I,padded,0,m - I,rows,0,n - I.cols,BORDER_CONSTANT,Scalar::all(0));
複製代碼

【3】爲傅里葉變換的結果分配儲存空間

傅里葉變換的結果是複數，每一個原圖像值，結果會有兩個圖像值。

爲傅里葉變換的結果（實部和虛部）分配儲存空間
//將planes數組組合合併成一個多通道的數組complexI
Mat planes[] = (Mat_<float>(padded),Mat::zeros(padded.size(),CV_32F));
Mat complexI;
merge(planes,2,complexI);
複製代碼

【4】進行離散傅里葉變換

def(complexI,complexI);
複製代碼

【5】將複數轉換爲幅值

離散傅里葉變換的結果是複數，對應的幅值可表示爲：

//將複數轉換爲幅值，即=> log(1 + sqrt(Re(DFT(I)) ^ 2 + Im(DFT(I)) ^ 2))
split(complexI,planes);  //將多通道complexI分離成幾個單通道數組
planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
magnitude(planes[0],planes[1],planes[0]);
Mat magnitudeImage = planes[0];
複製代碼

【6】進行對數尺度縮放

傅里葉變換的幅度之範圍大到不合適在屏幕上顯示。高值在屏幕上顯示爲白點，而低值爲黑點，高低值的變化沒法有效判斷。爲了在屏幕上凸顯出高低變化的連續性，能夠用對數尺度來代替線性尺度，公式以下：

//進行對數尺度縮放
magnitudeImage += Scalar::all(1);
log(magnitudeImage,magnitudeImage);  //求天然對數
複製代碼

【7】剪切和重分佈幅度圖像象限

由於在第二步中延擴了圖像，那如今是時候將新添加的像素剔除了。爲了方便顯示，能夠從新分佈圖像象限位置。

//如有奇數行或技術列，進行譜寫裁剪
magnitudeImage = magnitudeImage(Rect(0,0,magnitudeImage.cols & -2,magnitudeImage.rows & -2));
//從新排列傅里葉圖像中的象限，使得原點位於圖像中心
int cx = magnitudeImage.cols / 2;
int cy = magnitudeImage.rows / 2;
Mat q0(magnitudeImage,Rect(0,0,cx,cy)); //ROI區域的左上
Mat q1(magnitudeImage,Rect(cx,0,cx,cy));  //ROI區域的右上
Mat q2(magnitudeImage,Rect(0,cy,cx,cy));  //ROI區域的左下
Mat q3(magnitudeImage,Rect(cx,cy,cx,cy));  //ROI區域的右下

//交換象限(左上與右下)
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);

//交換象限(右上與左下)
q1.copyTo(tmp);
q4.copyTo(q1);
tmp.copyTo(q4);
複製代碼

【8】歸一化

如今有了重分佈後的幅度圖，可是幅度值仍然超過可顯示範圍[0,1],咱們可使用 normalize()函數歸一化到可顯示範圍。

normalize(magnitudeImage,magnitudeImage,0,1,NORM_MINMAX);
複製代碼

【9】顯示效果

imshow("頻譜幅值"，magnitudeImage);
複製代碼

效果