當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

opencv学之图像傅里叶变换dft

發布時間：2023/12/18 编程问答 33 豆豆

生活随笔收集整理的這篇文章主要介紹了 opencv学之图像傅里叶变换dft 小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

一、前言

接觸了圖像的傅里葉變換，數學原理依舊不是很懂，因此不敢在這里妄言。下午用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 imageint 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 transformMat 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 amplitudeMat mAmplitude;magnitude(mRe, mIm, mAmplitude);//Logarithmic scalemAmplitude += Scalar(1);log(mAmplitude, mAmplitude);//The normalizednormalize(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、歸一化。目前得到的幅值圖像仍然太大，超出了顯示的范圍，歸一化這范圍內的值，可以進一步達到可視化的目的

實現程序

void _DFT(){

//1以灰度模式讀取原圖像并顯示

Mat srcImage = imread("miFan.jpg",0);

if (!srcImage.data){ cout << "Error\n"; }

imshow("原圖像", srcImage);

//2將輸入圖像擴展到最佳尺寸,邊界用0補充

int m = getOptimalDFTSize(srcImage.rows);

int n = getOptimalDFTSize(srcImage.cols);

//將添加的像素初始化為0

Mat padded;

copyMakeBorder(srcImage, padded, 0, m - srcImage.rows,

0, n - srcImage.cols, BORDER_CONSTANT, Scalar::all(0));

//3為傅里葉變換的結果(實部和虛部)分配存儲空間

//將數組組合合并為一個多通道數組

Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };

Mat complexI;

merge(planes, 2, complexI);

//4進行傅里葉變換

dft(complexI, complexI);

//5將復數轉換為幅值，即=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))

//將多通道數組分離為幾個單通道數組

split(complexI, planes);//planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))

magnitude(planes[0], planes[1], planes[0]);

Mat magImage = planes[0];

//6進行對數尺度縮放

magImage += Scalar::all(1);

log(magImage, magImage);//求自然對數

//7剪切和重分布幅度圖象限

//若有奇數行或奇數列，進行頻譜剪裁

magImage = magImage(Rect(0, 0, magImage.cols&-2, magImage.rows&-2));

//重新排列傅立葉圖像中的象限，使得原點位于圖像中心

int cx = magImage.cols / 2;

int cy = magImage.rows / 2;

Mat q0(magImage, Rect(0, 0, cx, cy));

Mat q1(magImage, Rect(cx, 0, cx, cy));

Mat q2(magImage, Rect(0,cy,cx,cy));

Mat q3(magImage, 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);

//8歸一化，用0到1的浮點值將矩陣變換為可視的圖像格式

normalize(magImage, magImage, 0, 1, CV_MINMAX);

//9顯示

imshow("頻譜增幅", magImage);

waitKey();

}

正在上傳…重新上傳取消轉存失敗重新上傳取消

傅里葉變換后的圖片

正在上傳…重新上傳取消轉存失敗重新上傳取消

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 neededC.create(abs(A.rows - B.rows)+1, abs(A.cols - B.cols)+1, A.type());Size dftSize;// calculate the size of DFT transformdftSize.width = getOptimalDFTSize(A.cols + B.cols - 1);dftSize.height = getOptimalDFTSize(A.rows + B.rows - 1);// allocate temporary buffers and initialize them with 0'sMat tempA(dftSize, A.type(), Scalar::all(0));//initial 0Mat tempB(dftSize, B.type(), Scalar::all(0));// copy A and B to the top-left corners of tempA and tempB, respectivelyMat 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 processingdft(tempA, tempA, 0, A.rows);dft(tempB, tempB, 0, B.rows);// multiply the spectrums;// the function handles packed spectrum representations wellmulSpectrums(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.rowsdft(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 floatMat 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

總結

以上是生活随笔為你收集整理的opencv学之图像傅里叶变换dft的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。

上一篇：解决:Unable to identif
下一篇：解决： Linux – git: com

编程问答

opencv学 之图像傅里叶变换dft

一、前言

二、關鍵函數解析

2.1copyMakeBorder() 擴展圖片尺寸

2.2getOptimalDFTSize() 獲取最佳計算尺寸

2.3dft()傅里葉變換計算

2.4magnitude()計算二維矢量幅值

2.5log()自然對數計算

2.6normalize()矩陣歸一化

三、代碼及結果分享

?

四、注意事項

4.1Mat對象類型問題

4.2clone()與copyTo()的差異問題

Opencv 實現圖像的離散傅里葉變換（DFT）、卷積運算（相關濾波）

總結

opencv学之图像傅里叶变换dft