OpenCV基礎篇之圖像的DFT頻域變換

程序及分析

/*
 * FileName : fft2.cpp
 * Author   : xiahouzuoxin @163.com
 * Version  : v1.0
 * Date     : Wed 30 Jul 2014 09:42:12 PM CST
 * Brief    : 
 * 
 * Copyright (C) MICL,USTB
 */

#include <iostream>
#include <cv.h>
#include <highgui.h>
#include "imgproc/imgproc.hpp"

using namespace std;
using namespace cv;

int main(int argc, char *argv[])
{
    if (argc < 2) {
        cout<<"Usage:./fft2 [image name]"<<endl;
        return -1;
    }

    // Read as grayscale image
    Mat image = imread(argv[1], CV_LOAD_IMAGE_GRAYSCALE);
    if (!image.data) {
        cout << "Read image error"<<endl;
        return -1;
    }

    Mat padded;
    int m = getOptimalDFTSize(image.rows);  // Return size of 2^x that suite for FFT
    int n = getOptimalDFTSize(image.cols);
    // Padding 0, result is @padded
    copyMakeBorder(image, padded, 0, m-image.rows, 0, n-image.cols, BORDER_CONSTANT, Scalar::all(0));

    // Create planes to storage REAL part and IMAGE part, IMAGE part init are 0
    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
    Mat complexI;
    merge(planes, 2, complexI);

    dft(complexI, complexI);

    // compute the magnitude and switch to logarithmic scale
    split(complexI, planes);
    magnitude(planes[0], planes[0], planes[1]);
    Mat magI = planes[0];

    // => log(1+sqrt(Re(DFT(I))^2+Im(DFT(I))^2))
    magI += Scalar::all(1);
    log(magI, magI);

    // crop the spectrum
    magI = magI(Rect(0, 0, magI.cols & (-2), magI.rows & (-2)));
    Mat _magI = magI.clone();
    normalize(_magI, _magI, 0, 1, CV_MINMAX);

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

    // exchange Top-Left and Bottom-Right
    Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    // exchange Top-Right and Bottom-Left
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);

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

    imshow("Input image", image);
    imshow("Spectrum magnitude before shift frequency", _magI);
    imshow("Spectrum magnitude after shift frequency", magI);
    waitKey();

    return 0;
}

本程序的作用是：將圖像從空間域轉換到頻率域，並繪制頻域圖像。

1. 二維圖像的DFT（離散傅裡葉變換），

   圖像的頻域表示的是什麼含義呢？又有什麼用途呢？圖像的頻率是表征圖像中灰度變化劇烈程度的指標，是灰度在平面空間上的梯度。圖像的邊緣部分是突變部分，變化較快，因此反應在頻域上是高頻分量；圖像的噪聲大部分情況下是高頻部分；圖像大部分平緩的灰度變化部分則為低頻分量。也就是說，傅立葉變換提供另外一個角度來觀察圖像，可以將圖像從灰度分布轉化到頻率分布上來觀察圖像的特征。

   頻域在圖像處理中，就我所知的用途主要在兩方面：圖像壓縮和圖像去噪。關於這兩點將在下面給出圖片DFT的變換結果後說明。

   有關DFT的更多性質請參考胡廣書教授的《數字信號處理》教材。

2. 請注意讀圖片的函數與之前有所不同：

   Mat image = imread(argv[1], CV_LOAD_IMAGE_GRAYSCALE);

   CV_LOAD_IMAGE_GRAYSCALE參數表示將原圖像轉換為灰度圖後讀入，這是因為後面的DFT變換都是基於二維信號的，而彩色圖像是三維信號。當然，也可以對RGB每一通道都進行DFT運算。

3. DFT算法的原理要求輸入信號的長度最好為2^n，這樣可以使用快速傅裡葉變換算法（FFT算法）進行加速。所以程序中使用

   copyMakeBorder(image, padded, 0, m-image.rows, 0, n-image.cols, BORDER_CONSTANT, Scalar::all(0));

   填充0使橫縱長度都為2^n。

   對於一維信號，原DFT直接運算的復雜度是O(N^2)，而快速傅裡葉變換的復雜度降低到O(Nlog2(N))，假設N為512，足足提高了512/9≈57倍。

4. 由DFT的性質知，輸入為實信號（圖像）的時候，頻域輸出為復數，因此將頻域信息分為幅值和相位。頻域的幅值高的代表高頻分量，幅值低的地方代表低頻分量，因此程序中使用

   // => log(1+sqrt(Re(DFT(I))^2+Im(DFT(I))^2))
   magI += Scalar::all(1);
   log(magI, magI);
   
   // crop the spectrum
   magI = magI(Rect(0, 0, magI.cols & (-2), magI.rows & (-2)));
   Mat _magI = magI.clone();
   normalize(_magI, _magI, 0, 1, CV_MINMAX);

   進行log幅值計算及歸一化幅值（歸一化目的主要是方便將頻域通過圖像的形式進行顯示）。

5. 關於頻域中心平移：將圖像的高頻分量平移到圖像的中心，便於觀測。

   int cx = magI.cols/2;
   int cy = magI.rows/2;
   
   Mat q0(magI, Rect(0,0,cx,cy));    // Top-Left
   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
   
   // exchange Top-Left and Bottom-Right
   Mat tmp;
   q0.copyTo(tmp);
   q3.copyTo(q0);
   tmp.copyTo(q3);
   
   // exchange Top-Right and Bottom-Left
   q1.copyTo(tmp);
   q2.copyTo(q1);
   tmp.copyTo(q2);

   其原理就是將左上角的頻域和右下角的互換，右上角和左下角互換。

   請注意：頻域點和空域點的坐標沒有一一對應的關系，兩者的關系只是上面的DFT公式所見到的。

6. 本程序因為使用到圖像處理相關的函數，所以包含了頭文件imgproc/imgproc.hpp，該文件位於opencv安裝目錄的include/opencv2/目錄下，在編寫Makefile時也要增加相關的頭文件路徑和庫，本程序使用的Makefile如下：

   TARG=fft2
   SRC=fft2.cpp
   LIB=-L/usr/local/lib/
   INC=-I/usr/local/include/opencv/ -I/usr/local/include/opencv2
   CFLAGS=
   
   $(TARG):$(SRC)
     g++ -g -o $@ ${CFLAGS} $(LIB) $(INC) \
         -lopencv_core -lopencv_highgui -lopencv_imgproc \
         $^
   
   .PHONY:clean
   
   clean:
     -rm $(TARG) tags -f

   其中Makefile中的\表示換行（反斜槓後不能再有任何字符，包括空格），如上庫增加了-lopencv_imgproc，頭文件路徑增加了-I/usr/local/include/opencv2。

效果

1. 上圖從左到右分別是：原始灰度圖（我大愛的楊過啊）、頻域平移前的頻域圖像、頻域中心平移後的頻域圖像。

2. 提到圖像頻域變換的用途：壓縮和去噪。壓縮的原理就是在頻域中，大部分頻域的值為0（或接近0，可以進行有損壓縮，如jpeg圖像），只要壓縮頻域中的少數非0值即可達到圖片壓縮的目的。去噪則是通過頻域的濾波實現，因為噪聲大部分情況下體現為高頻信號，使用低通濾波器即可濾除高頻噪聲（當然，也會帶來損失，那就是邊緣會變得模糊（之前說過，邊緣也是高頻信號））。

--------------------------------------分割線 --------------------------------------

Ubuntu Linux下安裝OpenCV2.4.1所需包 http://www.linuxidc.com/Linux/2012-08/68184.htm

Ubuntu 12.04 安裝 OpenCV2.4.2 http://www.linuxidc.com/Linux/2012-09/70158.htm

CentOS下OpenCV無法讀取視頻文件 http://www.linuxidc.com/Linux/2011-07/39295.htm

Ubuntu 12.04下安裝OpenCV 2.4.5總結 http://www.linuxidc.com/Linux/2013-06/86704.htm

Ubuntu 10.04中安裝OpenCv2.1九步曲 http://www.linuxidc.com/Linux/2010-09/28678.htm

基於QT和OpenCV的人臉識別系統 http://www.linuxidc.com/Linux/2011-11/47806.htm

--------------------------------------分割線 --------------------------------------

OpenCV的詳細介紹：請點這裡
OpenCV的下載地址：請點這裡