離散傅立葉變換

目錄

• 目標
• 原始碼
• 解釋
  • 將影像擴充套件到最佳尺寸
  • 為複數值和實數值騰出空間
  • 進行離散傅立葉變換
  • 將實數值和複數值轉換為幅度
  • 切換到對數尺度
  • 裁剪和重排
  • Normalize
• 結果

上一個教程： 改變影像的對比度和亮度！
下一個教程： 使用XML / YAML / JSON檔案進行檔案輸入和輸出

原始作者	Bernát Gábor
相容性	OpenCV >= 3.0

目標

我們將尋找以下問題的答案

• 什麼是傅立葉變換，為什麼要使用它？
• 如何在OpenCV中實現？
• 使用以下函式： copyMakeBorder()、merge()、dft()、getOptimalDFTSize()、log() 和 normalize()。

原始碼

C++

你可以從此處下載或在OpenCV原始碼庫的 samples/cpp/tutorial_code/core/discrete_fourier_transform/discrete_fourier_transform.cpp 中找到它。

Java

你可以從此處下載或在OpenCV原始碼庫的 samples/java/tutorial_code/core/discrete_fourier_transform/DiscreteFourierTransform.java 中找到它。

Python

你可以從此處下載或在OpenCV原始碼庫的 samples/python/tutorial_code/core/discrete_fourier_transform/discrete_fourier_transform.py 中找到它。

以下是 dft() 的示例用法

C++

#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"
 
#include <iostream>
 
using namespace cv;
using namespace std;
 
static void help(char ** argv)
{
cout << endl
<< "本程式演示了離散傅立葉變換（DFT）的使用。 " << endl
<< "對影像進行dft變換並顯示其功率譜。" << endl << endl
<< "用法：" << endl
<< argv[0] << " [image_name -- 預設 lena.jpg]" << endl << endl;
}
 
int main(int argc, char ** argv)
{
help(argv);
 
    const char* filename = argc >=2 ? argv[1] : "lena.jpg";
 
    Mat I = imread( samples::findFile( filename ), IMREAD_GRAYSCALE);
    if( I.empty()){
cout << "開啟影像錯誤" << endl;
        return EXIT_FAILURE;
    }
 
    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));
 
    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexI;
merge(planes, 2, complexI); // 在擴充套件後的影像中新增一個全零的通道
 
dft(complexI, complexI); // 這樣結果可以儲存在源矩陣中
 
    // 計算幅度並切換到對數尺度
    // => 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]);// planes[0] = 幅度
    Mat magI = planes[0];
 
magI += Scalar::all(1); // 切換到對數尺度
log(magI, magI);
 
    // 如果頻譜的行數或列數為奇數，則進行裁剪
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);
 
normalize(magI, magI, 0, 1, NORM_MINMAX); // 將浮點值矩陣轉換為
                                            // 視覺化影像形式（浮點值介於0和1之間）。
 
imshow("Input Image" , I ); // 顯示結果
imshow("spectrum magnitude", magI);
cv::Mat::empty
 
    return EXIT_SUCCESS;
}

Java

import org.opencv.core.*;
import org.opencv.highgui.HighGui;
import org.opencv.imgcodecs.Imgcodecs;
 
import java.util.List;
import java.util.*;
 
class DiscreteFourierTransformRun{
    private void help() {
System.out.println("" +
                "本程式演示了離散傅立葉變換（DFT）的使用。 \n" +
                "對影像進行dft變換並顯示其功率譜。\n" +
                "用法:\n" +
                "./DiscreteFourierTransform [image_name -- 預設 ../data/lena.jpg]");
    }
 
    public void run(String[] args){
 
help();
 
String filename = ((args.length > 0) ? args[0] : "../data/lena.jpg");
 
Mat I = Imgcodecs.imread(filename, Imgcodecs.IMREAD_GRAYSCALE);
        if( I.empty() ) {
System.out.println("開啟影像錯誤");
System.exit(-1);
        }
 
Mat padded = new Mat(); //將輸入影像擴充套件到最佳尺寸
        int m = Core.getOptimalDFTSize( I.rows() );
        int n = Core.getOptimalDFTSize( I.cols() ); // 在邊界新增零值
Core.copyMakeBorder(I, padded, 0, m - I.rows(), 0, n - I.cols(), Core.BORDER_CONSTANT, Scalar.all(0));
 
List<Mat> planes = new ArrayList<Mat>();
padded.convertTo(padded, CvType.CV_32F);
planes.add(padded);
planes.add(Mat.zeros(padded.size(), CvType.CV_32F));
Mat complexI = new Mat();
Core.merge(planes, complexI); // 在擴充套件後的影像中新增一個全零的通道
 
Core.dft(complexI, complexI); // 這樣結果可以儲存在源矩陣中
 
        // 計算幅度並切換到對數尺度
        // => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
Core.split(complexI, planes); // planes.get(0) = Re(DFT(I)
                                                                    // planes.get(1) = Im(DFT(I))
Core.magnitude(planes.get(0), planes.get(1), planes.get(0));// planes.get(0) = 幅度
Mat magI = planes.get(0);
 
Mat matOfOnes = Mat.ones(magI.size(), magI.type());
Core.add(matOfOnes, magI, magI); // 切換到對數尺度
Core.log(magI, magI);
 
        // 如果頻譜的行數或列數為奇數，則進行裁剪
magI = magI.submat(new Rect(0, 0, magI.cols() & -2, magI.rows() & -2));
 
        // 重新排列傅立葉影像的象限，使原點位於影像中心
        int cx = magI.cols()/2;
        int cy = magI.rows()/2;
 
Mat q0 = new Mat(magI, new Rect(0, 0, cx, cy)); // 左上角 - 為每個象限建立一個感興趣區域
Mat q1 = new Mat(magI, new Rect(cx, 0, cx, cy)); // 右上角
Mat q2 = new Mat(magI, new Rect(0, cy, cx, cy)); // 左下角
Mat q3 = new Mat(magI, new Rect(cx, cy, cx, cy)); // 右下角
 
Mat tmp = new Mat(); // 交換象限（左上角與右下角）
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
 
q1.copyTo(tmp); // 交換象限（右上角與左下角）
q2.copyTo(q1);
tmp.copyTo(q2);
 
magI.convertTo(magI, CvType.CV_8UC1);
Core.normalize(magI, magI, 0, 255, Core.NORM_MINMAX, CvType.CV_8UC1); // 將浮點值矩陣
                                                                            // 轉換為視覺化影像形式（浮點值介於
                                                                            // 0和255之間）。
 
HighGui.imshow("Input Image" , I ); // 顯示結果
HighGui.imshow("Spectrum Magnitude", magI);
HighGui.waitKey();
 
System.exit(0);
    }
}
 
 
public class DiscreteFourierTransform {
    public static void main(String[] args) {
        // 載入本地庫。
System.loadLibrary(Core.NATIVE_LIBRARY_NAME);
        new DiscreteFourierTransformRun().run(args);
    }
}

Python

from __future__ import print_function
import sys
 
import cv2 as cv
import numpy as np
 
 
def print_help()
print('''
本程式演示了離散傅立葉變換（DFT）的使用。
對影像進行dft變換並顯示其功率譜。
用法
discrete_fourier_transform.py [image_name -- 預設 lena.jpg]''')
 
 
def main(argv)
 
print_help()
 
filename = argv[0] if len(argv) > 0 else 'lena.jpg'
 
I = cv.imread(cv.samples.findFile(filename), cv.IMREAD_GRAYSCALE)
    if I is None
print('開啟影像錯誤')
        return -1
    
rows, cols = I.shape
m = cv.getOptimalDFTSize( rows )
n = cv.getOptimalDFTSize( cols )
padded = cv.copyMakeBorder(I, 0, m - rows, 0, n - cols, cv.BORDER_CONSTANT, value=[0, 0, 0])
    
planes = [np.float32(padded), np.zeros(padded.shape, np.float32)]
complexI = cv.merge(planes) # 在擴充套件後的影像中新增一個全零的通道
    
    cv.dft(complexI, complexI) # 這樣結果可以儲存在源矩陣中
    
    cv.split(complexI, planes) # planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
    cv.magnitude(planes[0], planes[1], planes[0])# planes[0] = 幅度
magI = planes[0]
    
matOfOnes = np.ones(magI.shape, dtype=magI.dtype)
    cv.add(matOfOnes, magI, magI) # 切換到對數尺度
    cv.log(magI, magI)
    
magI_rows, magI_cols = magI.shape
    # 如果頻譜的行數或列數為奇數，則進行裁剪
magI = magI[0:(magI_rows & -2), 0:(magI_cols & -2)]
cx = int(magI_rows/2)
cy = int(magI_cols/2)
 
q0 = magI[0:cx, 0:cy] # 左上角 - 為每個象限建立一個感興趣區域
q1 = magI[cx:cx+cx, 0:cy] # 右上角
q2 = magI[0:cx, cy:cy+cy] # 左下角
q3 = magI[cx:cx+cx, cy:cy+cy] # 右下角
 
tmp = np.copy(q0) # 交換象限（左上角與右下角）
magI[0:cx, 0:cy] = q3
magI[cx:cx + cx, cy:cy + cy] = tmp
 
tmp = np.copy(q1) # 交換象限（右上角與左下角）
magI[cx:cx + cx, 0:cy] = q2
magI[0:cx, cy:cy + cy] = tmp
    
    cv.normalize(magI, magI, 0, 1, cv.NORM_MINMAX) # 將浮點值矩陣轉換為
    
    cv.imshow("Input Image" , I ) # 顯示結果
    cv.imshow("spectrum magnitude", magI)
    cv.waitKey()
 
if __name__ == "__main__"
    main(sys.argv[1:])

解釋

傅立葉變換將影像分解為正弦和餘弦分量。換句話說，它將影像從空間域轉換到頻率域。其思想是任何函式都可以透過無限多個正弦和餘弦函式的和來精確近似。傅立葉變換就是實現這一目的的一種方法。在數學上，二維影像的傅立葉變換是

\[F(k,l) = \displaystyle\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1} f(i,j)e^{-i2\pi(\frac{ki}{N}+\frac{lj}{N})}\]

\[e^{ix} = \cos{x} + i\sin {x}\]

這裡 f 是影像在空間域的值，F 是其在頻率域的值。變換結果是複數。顯示覆數可以透過*實部*影像和*虛部*影像，也可以透過*幅度*影像和*相位*影像。然而，在整個影像處理演算法中，只有*幅度*影像是令人感興趣的，因為它包含了我們所需的所有關於影像幾何結構的資訊。儘管如此，如果你打算以這些形式對影像進行一些修改，然後需要對其進行逆變換，則需要同時保留實部和虛部（或幅度和相位）。

在此示例中，我將展示如何計算並顯示傅立葉變換的*幅度*影像。數字影像是離散的。這意味著它們可以取給定域值中的一個值。例如，在基本的灰度影像中，值通常介於0到255之間。因此，傅立葉變換也需要是離散型別，從而產生離散傅立葉變換（DFT）。當你需要從幾何角度確定影像結構時，你都會希望使用它。以下是要遵循的步驟（針對灰度輸入影像I）

將影像擴充套件到最佳尺寸

DFT的效能取決於影像大小。對於尺寸是2、3和5倍數的影像，它的速度往往最快。因此，為了獲得最佳效能，通常最好在影像上填充邊界值，以獲得具有這些特性的尺寸。 getOptimalDFTSize() 函式返回此最佳尺寸，我們可以使用 copyMakeBorder() 函式來擴充套件影像的邊界（附加的畫素用零初始化）

C++

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

Java

Mat padded = new Mat(); //將輸入影像擴充套件到最佳尺寸
        int m = Core.getOptimalDFTSize( I.rows() );
        int n = Core.getOptimalDFTSize( I.cols() ); // 在邊界新增零值
Core.copyMakeBorder(I, padded, 0, m - I.rows(), 0, n - I.cols(), Core.BORDER_CONSTANT, Scalar.all(0));

Python

rows, cols = I.shape
m = cv.getOptimalDFTSize( rows )
n = cv.getOptimalDFTSize( cols )
padded = cv.copyMakeBorder(I, 0, m - rows, 0, n - cols, cv.BORDER_CONSTANT, value=[0, 0, 0])

為複數值和實數值騰出空間

傅立葉變換的結果是複數。這意味著對於每個影像值，結果是兩個影像值（每個分量一個）。此外，頻率域的範圍比其空間域對應物大得多。因此，我們通常至少以*浮點*格式儲存這些值。因此，我們會將輸入影像轉換為此型別，並用另一個通道來儲存複數值。

C++

    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexI;
    merge(planes, 2, complexI); // 在擴充套件後的影像中新增一個全零的通道

Java

List<Mat> planes = new ArrayList<Mat>();
padded.convertTo(padded, CvType.CV_32F);
planes.add(padded);
planes.add(Mat.zeros(padded.size(), CvType.CV_32F));
Mat complexI = new Mat();
Core.merge(planes, complexI); // 在擴充套件後的影像中新增一個全零的通道

Python

planes = [np.float32(padded), np.zeros(padded.shape, np.float32)]
complexI = cv.merge(planes) # 在擴充套件後的影像中新增一個全零的通道

進行離散傅立葉變換

可以進行原地計算（輸入和輸出相同）

C++

    dft(complexI, complexI); // 這樣結果可以儲存在源矩陣中

Java

Core.dft(complexI, complexI); // 這樣結果可以儲存在源矩陣中

Python

    cv.dft(complexI, complexI) # 這樣結果可以儲存在源矩陣中

將實數值和複數值轉換為幅度

複數有一個實部（Re）和一個復（虛部 - Im）部分。DFT的結果是複數。DFT的幅度是

\[M = \sqrt[2]{ {Re(DFT(I))}^2 + {Im(DFT(I))}^2}\]

轉換為OpenCV程式碼

C++

    split(complexI, planes); // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
    magnitude(planes[0], planes[1], planes[0]);// planes[0] = 幅度
    Mat magI = planes[0];

Java

Core.split(complexI, planes); // planes.get(0) = Re(DFT(I)
                                                                    // planes.get(1) = Im(DFT(I))
Core.magnitude(planes.get(0), planes.get(1), planes.get(0));// planes.get(0) = 幅度
Mat magI = planes.get(0);

Python

    cv.split(complexI, planes) # planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
    cv.magnitude(planes[0], planes[1], planes[0])# planes[0] = 幅度
magI = planes[0]

切換到對數尺度

傅立葉係數的動態範圍太大，無法在螢幕上顯示。我們有一些小值和一些變化大的值，無法這樣觀察。因此，高值將全部顯示為白點，而小值將顯示為黑點。為了使用灰度值進行視覺化，我們可以將線性尺度轉換為對數尺度

\[M_1 = \log{(1 + M)}\]

轉換為OpenCV程式碼

C++

magI += Scalar::all(1); // 切換到對數尺度
log(magI, magI);

Java

Mat matOfOnes = Mat.ones(magI.size(), magI.type());
Core.add(matOfOnes, magI, magI); // 切換到對數尺度
Core.log(magI, magI);

Python

matOfOnes = np.ones(magI.shape, dtype=magI.dtype)
    cv.add(matOfOnes, magI, magI) # 切換到對數尺度
    cv.log(magI, magI)

裁剪和重排

還記得第一步中我們擴充套件了影像嗎？現在是時候丟棄新引入的值了。為了視覺化目的，我們還可以重新排列結果的象限，以便原點（零，零）與影像中心對應。

C++

    // 如果頻譜的行數或列數為奇數，則進行裁剪
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);

Java

        // 如果頻譜的行數或列數為奇數，則進行裁剪
magI = magI.submat(new Rect(0, 0, magI.cols() & -2, magI.rows() & -2));
 
        // 重新排列傅立葉影像的象限，使原點位於影像中心
        int cx = magI.cols()/2;
        int cy = magI.rows()/2;
 
Mat q0 = new Mat(magI, new Rect(0, 0, cx, cy)); // 左上角 - 為每個象限建立一個感興趣區域
Mat q1 = new Mat(magI, new Rect(cx, 0, cx, cy)); // 右上角
Mat q2 = new Mat(magI, new Rect(0, cy, cx, cy)); // 左下角
Mat q3 = new Mat(magI, new Rect(cx, cy, cx, cy)); // 右下角
 
Mat tmp = new Mat(); // 交換象限（左上角與右下角）
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
 
q1.copyTo(tmp); // 交換象限（右上角與左下角）
q2.copyTo(q1);
tmp.copyTo(q2);

Python

magI_rows, magI_cols = magI.shape
    # 如果頻譜的行數或列數為奇數，則進行裁剪
magI = magI[0:(magI_rows & -2), 0:(magI_cols & -2)]
cx = int(magI_rows/2)
cy = int(magI_cols/2)
 
q0 = magI[0:cx, 0:cy] # 左上角 - 為每個象限建立一個感興趣區域
q1 = magI[cx:cx+cx, 0:cy] # 右上角
q2 = magI[0:cx, cy:cy+cy] # 左下角
q3 = magI[cx:cx+cx, cy:cy+cy] # 右下角
 
tmp = np.copy(q0) # 交換象限（左上角與右下角）
magI[0:cx, 0:cy] = q3
magI[cx:cx + cx, cy:cy + cy] = tmp
 
tmp = np.copy(q1) # 交換象限（右上角與左下角）
magI[cx:cx + cx, 0:cy] = q2
magI[0:cx, cy:cy + cy] = tmp

歸一化

再次這樣做是為了視覺化目的。我們現在有了幅度，但這些值仍然超出我們影像顯示範圍（0到1）。我們使用 cv::normalize() 函式將我們的值歸一化到此範圍。

C++

    normalize(magI, magI, 0, 1, NORM_MINMAX); // 將浮點值矩陣轉換為
                                            // 視覺化影像形式（浮點值介於0和1之間）。

Java

Core.normalize(magI, magI, 0, 255, Core.NORM_MINMAX, CvType.CV_8UC1); // 將浮點值矩陣
                                                                            // 轉換為視覺化影像形式（浮點值介於
                                                                            // 0和255之間）。

Python

    cv.normalize(magI, magI, 0, 1, cv.NORM_MINMAX) # 將浮點值矩陣轉換為
    

結果

一個應用想法是確定影像中存在的幾何方向。例如，我們能否判斷一段文字是水平的還是傾斜的？觀察一些文字，你會發現文字行會形成水平線，而字母會形成垂直線。文字片段的這兩個主要組成部分也可以在傅立葉變換中看到。讓我們使用這張水平文字影像和這張旋轉文字影像。

對於水平文字

對於旋轉文字

你可以看到頻率域中最具影響力的分量（幅度影像上最亮點）遵循影像上物體的幾何旋轉。由此我們可以計算偏移量並執行影像旋轉以糾正最終的未對齊。