離散フーリエ変換

目次

• 目的
• ソースコード
• Explanation
  • 画像を最適なサイズに拡張する
  • 複素数値と実数値の両方のための場所を確保する
  • 離散フーリエ変換を行う
  • 実数値と複素数値を振幅に変換する
  • 対数スケールに切り替える
  • 切り抜きと並べ替え
  • 正規化
• 結果

前のチュートリアル: 画像のコントラストと明るさを変える!
次のチュートリアル: 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
        <<  "This program demonstrated the use of the discrete Fourier transform (DFT). " << endl
        <<  "The dft of an image is taken and it's power spectrum is displayed."  << endl << endl
        <<  "Usage:"                                                                      << endl
        << argv[0] << " [image_name -- default 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 << "Error opening image" << endl;
 return EXIT_FAILURE;
    }
 
 Mat padded;                            //expand input image to optimal size
 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));
 
 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
 
 // 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))
    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
    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, NORM_MINMAX); // Transform the matrix with float values into a
 // viewable image form (float between values 0 and 1).
 
    imshow("Input Image"       , I   );    // Show the result
    imshow("spectrum magnitude", magI);
    waitKey();
 
 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("" +
 "This program demonstrated the use of the discrete Fourier transform (DFT). \n" +
 "The dft of an image is taken and it's power spectrum is displayed.\n" +
 "Usage:\n" +
 "./DiscreteFourierTransform [image_name -- default ../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("Error opening image");
            System.exit(-1);
        }
 
        Mat padded = new Mat();                     //expand input image to optimal size
 int m = Core.getOptimalDFTSize( I.rows() );
 int n = Core.getOptimalDFTSize( I.cols() ); // on the border add zero values
        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);         // Add to the expanded another plane with zeros
 
        Core.dft(complexI, complexI);         // this way the result may fit in the source matrix
 
 // compute the magnitude and switch to logarithmic scale
 // => 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) = magnitude
        Mat magI = planes.get(0);
 
        Mat matOfOnes = Mat.ones(magI.size(), magI.type());
        Core.add(matOfOnes, magI, magI);         // switch to logarithmic scale
        Core.log(magI, magI);
 
 // crop the spectrum, if it has an odd number of rows or columns
        magI = magI.submat(new 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 = new Mat(magI, new Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
        Mat q1 = new Mat(magI, new Rect(cx, 0, cx, cy));  // Top-Right
        Mat q2 = new Mat(magI, new Rect(0, cy, cx, cy));  // Bottom-Left
        Mat q3 = new Mat(magI, new Rect(cx, cy, cx, cy)); // Bottom-Right
 
        Mat tmp = new Mat();               // 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);
 
        magI.convertTo(magI, CvType.CV_8UC1);
        Core.normalize(magI, magI, 0, 255, Core.NORM_MINMAX, CvType.CV_8UC1); // Transform the matrix with float values
 // into a viewable image form (float between
 // values 0 and 255).
 
        HighGui.imshow("Input Image"       , I   );    // Show the result
        HighGui.imshow("Spectrum Magnitude", magI);
        HighGui.waitKey();
 
        System.exit(0);
    }
}
 
 
public class DiscreteFourierTransform {
 public static void main(String[] args) {
 // Load the native library.
        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('''
    This program demonstrated the use of the discrete Fourier transform (DFT).
    The dft of an image is taken and it's power spectrum is displayed.
    Usage:
    discrete_fourier_transform.py [image_name -- default 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('Error opening image')
 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)         # Add to the expanded another plane with zeros
 
 cv.dft(complexI, complexI)         # this way the result may fit in the source matrix
 
 cv.split(complexI, planes)                   # planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
 cv.magnitude(planes[0], planes[1], planes[0])# planes[0] = magnitude
    magI = planes[0]
 
    matOfOnes = np.ones(magI.shape, dtype=magI.dtype)
 cv.add(matOfOnes, magI, magI) #  switch to logarithmic scale
 cv.log(magI, magI)
 
    magI_rows, magI_cols = magI.shape
 # crop the spectrum, if it has an odd number of rows or columns
    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]         # Top-Left - Create a ROI per quadrant
    q1 = magI[cx:cx+cx, 0:cy]     # Top-Right
    q2 = magI[0:cx, cy:cy+cy]     # Bottom-Left
    q3 = magI[cx:cx+cx, cy:cy+cy] # Bottom-Right
 
    tmp = np.copy(q0)               # swap quadrants (Top-Left with Bottom-Right)
    magI[0:cx, 0:cy] = q3
    magI[cx:cx + cx, cy:cy + cy] = tmp
 
    tmp = np.copy(q1)               # swap quadrant (Top-Right with Bottom-Left)
    magI[cx:cx + cx, 0:cy] = q2
    magI[0:cx, cy:cy + cy] = tmp
 
 cv.normalize(magI, magI, 0, 1, cv.NORM_MINMAX) # Transform the matrix with float values into a
 
 cv.imshow("Input Image"       , I   )    # Show the result
 cv.imshow("spectrum magnitude", magI)
 cv.waitKey()
 
if __name__ == "__main__":
 main(sys.argv[1:])

解説

フーリエ変換は画像を正弦成分と余弦成分に分解する。言い換えると、画像を空間領域から周波数領域に変換する。その考え方は、任意の関数は無限個の正弦関数と余弦関数の和によって厳密に近似できるというものである。フーリエ変換はそれを実現する方法の一つである。数学的に、2次元画像のフーリエ変換は次のとおりである:

\[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 は周波数領域における値である。変換の結果は複素数になる。これを表示するには、実部 (real) 画像と 虚部 (complex) 画像で表す方法と、振幅 (magnitude) 画像と 位相 (phase) 画像で表す方法がある。ただし画像処理アルゴリズム全般において重要なのは 振幅 (magnitude) 画像だけであり、これは画像の幾何学的構造について必要な情報をすべて含んでいる。とはいえ、これらの形式で画像に何らかの修正を加え、その後に逆変換する必要がある場合は、両方を保持しておく必要がある。

このサンプルでは、フーリエ変換の 振幅 (magnitude) 画像を計算して表示する方法を示す。デジタル画像の場合、画像は離散的である。これは、与えられた領域内の値を取りうることを意味する。例えば、基本的なグレースケール画像では値は通常0から255の間である。したがってフーリエ変換も離散型である必要があり、その結果が離散フーリエ変換 (DFT) となる。画像の構造を幾何学的観点から判定したいときには、これを使うとよい。以下に従うべき手順を示す(グレースケール入力画像 I の場合):

画像を最適なサイズに拡張する

DFTの性能は画像サイズに依存する。2、3、5の倍数となるサイズで最速になる傾向がある。したがって、最大の性能を得るためには、こうした性質を持つサイズになるよう画像に境界値をパディングするのが一般によい考えである。getOptimalDFTSize() はこの最適なサイズを返し、copyMakeBorder() 関数を使って画像の境界を拡張できる(追加されたピクセルはゼロで初期化される):

C++

 Mat padded;                            //expand input image to optimal size
 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));

Java

        Mat padded = new Mat();                     //expand input image to optimal size
 int m = Core.getOptimalDFTSize( I.rows() );
 int n = Core.getOptimalDFTSize( I.cols() ); // on the border add zero values
        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])

複素数値と実数値の両方を格納する場所を確保する

フーリエ変換の結果は複素数である。これは、各画像値に対して結果が2つの画像値(成分ごとに1つ)になることを意味する。さらに、周波数領域の値の範囲は空間領域に比べてはるかに大きい。したがって、これらは通常、少なくとも float 形式で格納する。そのため、入力画像をこの型に変換し、複素数値を保持するためにもう1つチャンネルを追加して拡張する:

C++

 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

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);         // Add to the expanded another plane with zeros

Python

    planes = [np.float32(padded), np.zeros(padded.shape, np.float32)]
    complexI = cv.merge(planes)         # Add to the expanded another plane with zeros

離散フーリエ変換を行う

インプレース計算(入力と出力が同じ)も可能である:

C++

 dft(complexI, complexI);            // this way the result may fit in the source matrix

Java

        Core.dft(complexI, complexI);         // this way the result may fit in the source matrix

Python

 cv.dft(complexI, complexI)         # this way the result may fit in the source matrix

実数値と複素数値を振幅に変換する

複素数は実部 (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] = magnitude
 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) = magnitude
        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] = magnitude
    magI = planes[0]

対数スケールに切り替える

フーリエ係数のダイナミックレンジは大きすぎて画面に表示できないことがわかる。このままでは観察できないような、小さい値と大きく変化する値が混在している。そのため、大きい値はすべて白い点として表れ、小さい値は黒として表れる。可視化のためにグレースケール値を使うには、線形スケールを対数スケールに変換すればよい:

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

OpenCVのコードに翻訳すると:

C++

    magI += Scalar::all(1);                    // switch to logarithmic scale
    log(magI, magI);

Java

        Mat matOfOnes = Mat.ones(magI.size(), magI.type());
        Core.add(matOfOnes, magI, magI);         // switch to logarithmic scale
        Core.log(magI, magI);

Python

    matOfOnes = np.ones(magI.shape, dtype=magI.dtype)
 cv.add(matOfOnes, magI, magI) #  switch to logarithmic scale
 cv.log(magI, magI)

切り出しと並べ替え

最初のステップで画像を拡張したことを思い出してほしい。ここで、新たに追加した値を破棄する時が来た。可視化のために、結果の象限を並べ替えて、原点(ゼロ, ゼロ)が画像の中心に対応するようにしてもよい。

C++

 // crop the spectrum, if it has an odd number of rows or columns
    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);

Java

 // crop the spectrum, if it has an odd number of rows or columns
        magI = magI.submat(new 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 = new Mat(magI, new Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
        Mat q1 = new Mat(magI, new Rect(cx, 0, cx, cy));  // Top-Right
        Mat q2 = new Mat(magI, new Rect(0, cy, cx, cy));  // Bottom-Left
        Mat q3 = new Mat(magI, new Rect(cx, cy, cx, cy)); // Bottom-Right
 
        Mat tmp = new Mat();               // 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);

Python

    magI_rows, magI_cols = magI.shape
 # crop the spectrum, if it has an odd number of rows or columns
    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]         # Top-Left - Create a ROI per quadrant
    q1 = magI[cx:cx+cx, 0:cy]     # Top-Right
    q2 = magI[0:cx, cy:cy+cy]     # Bottom-Left
    q3 = magI[cx:cx+cx, cy:cy+cy] # Bottom-Right
 
    tmp = np.copy(q0)               # swap quadrants (Top-Left with Bottom-Right)
    magI[0:cx, 0:cy] = q3
    magI[cx:cx + cx, cy:cy + cy] = tmp
 
    tmp = np.copy(q1)               # swap quadrant (Top-Right with Bottom-Left)
    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); // Transform the matrix with float values into a
 // viewable image form (float between values 0 and 1).

Java

        Core.normalize(magI, magI, 0, 255, Core.NORM_MINMAX, CvType.CV_8UC1); // Transform the matrix with float values
 // into a viewable image form (float between
 // values 0 and 255).

Python

 cv.normalize(magI, magI, 0, 1, cv.NORM_MINMAX) # Transform the matrix with float values into a
 

結果

応用例の一つとして、画像中に存在する幾何学的な向きを判定することが考えられる。例えば、あるテキストが水平かどうかを調べてみよう。テキストを見ると、テキスト行がいわば水平の線を形成し、文字がいわば垂直の線を形成していることに気づくだろう。テキストの断片におけるこれら2つの主要な成分は、フーリエ変換の場合にも見ることができる。テキストについての この水平な画像と この回転した画像を使ってみよう。

水平なテキストの場合:

回転したテキストの場合:

周波数領域で最も影響力のある成分(振幅画像上で最も明るい点)が、画像中のオブジェクトの幾何学的な回転に従っていることがわかる。ここからオフセットを計算し、画像を回転させて最終的な位置ずれを補正できる。