日記/2010-04-04/octave で画像処理 - 2 次元 FFT

octave で画像処理 - 2 次元 FFT †

octave での画像の読み書きができたところで，カッ飛ばしていきます．

２次元 FFT †

ここで 2 次元 FFT について説明します．

まず，元の画像を横方向に１ラインずつ抜きだし，（１次元の）FFT をかけます． それらを寄せ集め，今度は縦方向に１ラインずつ抜きだし，FFT をかけ，寄せ集めます． この結果が 2 次元 FFT となります．

なお，横→縦と１次元 FFT を実行しても縦→横と１次元 FFT を実行しても同じ結果になることが知られています．

2 次元逆 FFT も同様に定義されます． こちらでも「横→縦」の逆変換でも「縦→横」の逆変換でも同じ結果になることが知られています．

octave では †

octave では fft2(), ifft2() という関数で２次元 FFT・２次元逆 FFT を行うことができます．

octave:119> [ x, m, a ] = imread ( 'DSCF0135.JPG' );  # 画像ファイルを読み込み
octave:120> y = fft2 ( x );                           # ２次元 FFT の実行
octave:121> xx = ifft2 ( y );                         # ２次元逆 FFT の実行

実際に 2 次元 FFT をやってみる †

理屈をグダグダ言っててもつまらないので，実際にやってみます． 画像はこんなやつを使ってみました．

まずは画像ファイルを読み込んで

octave:26> [ x, m, a ] = imread ( 'DSCF0135.JPG' );

2 次元 FFT を実行します．

octave:27> y = fft2(x);

FFT では複素数で波の位相を表していますが，そのままでは表示に無理があるので，絶対値 ( √(Re^2 + Im^2) ) で画像にして表示してみます．

octave:28> imshow ( abs(y));

む，画像が真っ白だ．値が大きすぎてサチってるようです． 値を適当にスケーリングします． いろいろ試行錯誤したら，100000 で割ると，何となくいい感じのようです．

octave:29> imshow ( abs(y)/100000);

こんな感じの画像が表示されます．

左上 1/4 が周波数成分（の大きさ）で，右上，左下，右下には横方向・縦方向のエイリアジングが現れています． 1 次元の FFT のときのように，これらは負の周波数の成分，という見かたもできます．

というわけで，十文字に４等分して，上下左右を入れ替えてみます． が，そういう関数が見つからなかったので，適当に作ってみました． swap2.m という名前で保存すると，octave 上から swap2 (x) のように呼び出すことができます．

function y = swap2 ( x )
  r = rows(x);
  c = columns(x);
  r2 = floor(r/2);
  c2 = floor(c/2);

  y = [x(r2+1:r, c2+1:c, : ) x(r2+1:r, 1:c2, : );
       x(1:r2,   c2+1:c, : ) x(1:r2,   1:c2, : ) ];
endfunction

使ってみます．

octave:33> imshow ( abs(swap2(y))/100000 );

こんな感じの画像が表示されます． 画像の中央が低い周波数成分で，端に行くほど高い周波数成分を表しています．