公式チュートリアルのDiscrete Fourier Transformの要約。
筆者の環境はdebian9。
##やること
入力画像を離散フーリエ変換して周波数画像を出力する
###何が嬉しいのか
入力画像の幾何的な構造についての情報が得られる。今回の例では文字列が傾いているとかがわかりやすくなる。
###理論
フーリエ変換は画像をサイン・コサインの要素に分解する。つまり空間ドメインから周波数ドメインへ変換する。フーリエ変換のアイデアは、いかなる関数も無限のsin,cos関数の和で近似できるというもの。
画像のフーリエ変換は次式で表される。
F(k,l)= \sum_{i=0}^{N−1}\sum_{j=0}^{N−1}f(i,j)e^{−i2π(\frac{ki}{N}+\frac{lj}{N})}
e^{ix}=cosx+isinx
$f$は空間ドメインの画素値、$F$はその周波数ドメイン。結果は複素数になる。実数部、虚数部、絶対値、位相それぞれで画像を作れるが、ほしい情報をもっているのは絶対値(magnitude)画像のみ。
ピクセルの値は離散値なので離散フーリエ変換(Discrete Fourier Transform (DFT))を使う。グレースケール画像で行う。
##コード
dft.cpp
#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(void)
{
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
<< "Usage:" << endl
<< "./discrete_fourier_transform [image_name -- default ../data/lena.jpg]" << endl;
}
int main(int argc, char ** argv)
{
help();
const char* filename = argc >=2 ? argv[1] : "../data/lena.jpg";
Mat I = imread(filename, IMREAD_GRAYSCALE);
if( I.empty()){
cout << "Error opening image" << endl;
return -1;
}
Mat padded;
// 入力画像をDFTに最適な大きさに広げる。画像サイズが2,3,5の倍数のときに高速になる
int m = getOptimalDFTSize( I.rows );
int n = getOptimalDFTSize( I.cols );
// 入力画像を中央に置き、周囲は0で埋める
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
//dtfの結果は複素数（１要素につき２つ値がある）であり、また値の表現範囲が広い。
//そのためfloatと複素数を保持するもので2枚作る
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
// 2枚の画像を、2チャネルを持った1枚にする
merge(planes, 2, complexI);
// complexIにdftを適用し、complexIに結果を戻す
dft(complexI, complexI);
// 絶対値を計算し、logスケールにする
// => 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);
// 結果の値は大きすぎるものと小さすぎるものが混じっているので、logを適用して抑制する
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);
// 見ることができる値(float[0,1])に変換
normalize(magI, magI, 0, 1, NORM_MINMAX);
//表示
imshow("Input Image" , I );
imshow("spectrum magnitude", magI);
waitKey();
return 0;
}
###コンパイル
g++ dft.cpp -I/usr/local/include/opencv2 -I/usr/local/include/opencv -L/usr/local/lib -lopencv_core -lopencv_imgcodecs -lopencv_highgui
##結果
次のような感じになる
入力
結果
入力
結果