DFTを使って、ガウシアンフィルタ（Gaussian Filter）を高速化する方法を説明します。基本的なアイデアは、対象画像を一旦フーリエ変換し、フーリエ領域（Fourie domain）でマスク処理を行うことで、計算量を減らそうというものです。
ガウシアンフィルタとは？
ガウシアンフィルタでは、対象画素に近い画素に大きな重みをつけ、遠い画素には小さい重みを付けた加重平均を取ることで、自然なぼかしを実現しています。この重みづけを下のガウス関数を用いて行っているのでガウシアンフィルタと呼ばれます。
この式からも分かるようにσの値が小さいとぼかし効果は小さく、σの値が大きくなるとそれに従ってぼかし効果も大きくなります。
ガウシアンフィルタはボックスフィルタに比べて非常にきれいなぼかしが可能です。ボックスフィルタの場合、ぼかしを強くしたとき画質劣化がひどいのに対して、ガウシアンフィルタでは自然なぼかしを実現しています。
高速化のアイデア
ガウシアンフィルタを高速化する方基本的なアイデアは、上にも書いたように、時間のかかる畳み込み演算の部分をフーリエ領域に持ち込んで、計算量をへらそうというものです。図示するとこんな感じ。
空間領域(spatial domain)での畳み込み演算がフーリエ領域(fourie domain)では乗算に帰着できる理論はこの辺の論文を参考にしてみてください。
- P. M. Morse and H. Feshbach, “Fourier transforms,” in Methods of Theoretical Physics (Part I, pp. 453–471), New York: McGraw-Hill, 1953.
- R. Bracewell, “Convolution” and “Two-dimensional convolution,” in The Fourier Transform and Its Applications (pp. 25–50 and 243–244), New York: McGraw-Hill, 1965.
ガウシアンフィルタの高速化
1. たたみ込み積分での処理
畳み込み演算をした場合には、下の図のようにガウシアンフィルタのカーネルと、対象画像の画素との畳み込み演算( convolution )を行うことで処理を行います。
ガウシアンフィルタのカーネルサイズをMxM、対象の画像サイズをNxNとした場合、たたみ込み積分を行うと計算量は O(N2M2)になります。このように、畳み込み演算では、ガウシアンフィルタのカーネルサイズが大きくなればなるほど処理時間がかかるという問題があります。
2. フーリエ変換を用いた処理
一方で、DFTを用いて処理した場合の処理の流れは下のようになります。対象画像をDFTでフーリエ変換し、その結果に対して高周波部分を全て０にするようなマスクをかけます（乗算処理）。こうすることで高周波部分が除去され、ローパスフィルタをかけたのと同じ効果が得られます。
DFTを用いた場合の計算量はO（N2logN）に比例します。このように、DFTを用いることでガウシアンフィルタのカーネルサイズにかかわらず一定時間で処理できるようになるのです。
高速ガウシアンフィルタのプログラム
上に書いたように、DFTを用いた高速ガウシアンフィルタのプログラムのコア部分を載せます。このプログラムはこの本を参考させて頂きました。
Gary Bradski Adrian Kaehler
オライリージャパン
売り上げランキング: 15137
この関数の引数、srcとfilterがそれぞれ元画像とガウシアンフィルタのカーネルに相当します。まずは、それぞれの行列をcvDFT()関数を用いてフーリエ変換し、その結果どおしをcvMulSpectrums()関数で乗算し、最後に逆フーリエ変換することで、ぼかしのかかった画像を作成しています。
void speedy_conv (
const CvMat* src,
const CvMat* filter,
CvMat* dst
){
int dft_M = cvGetOptimalDFTSize( src->rows+dst->rows-1 );
int dft_N = cvGetOptimalDFTSize( src->cols+dst->cols-1 );
CvMat* dft_A = cvCreateMat( dft_M, dft_N, CV_32F );
CvMat* dft_B = cvCreateMat( dft_M, dft_N, CV_32F);
CvMat tmp;
// copy A to dft_A and pad dft_A with zeros
//
cvGetSubRect( dft_A, &tmp, cvRect(0,0,src->cols,src->rows));
cvCopy( src, &tmp );
cvGetSubRect( dft_A, &tmp, cvRect( src->cols, 0, dft_A->cols-src->cols, src->rows));
cvZero( &tmp );
// no need to pad bottom part of dft_A with zeros because of
// use nonzero_rows parameter in cvDFT() call below
//
cvDFT( dft_A, dft_A, CV_DXT_FORWARD, src->rows );
// repeat the same with the second array
//
cvGetSubRect( dft_B, &tmp, cvRect(0,0,filter->cols,filter->rows) );
cvCopy( filter, &tmp );
cvGetSubRect(dft_B, &tmp, cvRect( filter->cols, 0, dft_B->cols-filter->cols, filter->rows ));
cvZero( &tmp );
// no need to pad bottom part of dft_B with zeros because of
// use nonzero_rows parameter in cvDFT() call below
//
cvDFT( dft_B, dft_B, CV_DXT_FORWARD, filter->rows );
// or CV_DXT_MUL_CONJ to get correlation rather than convolution
//
cvMulSpectrums( dft_A, dft_B, dft_A, 0 );
// calculate only the top part
//
cvDFT( dft_A, dft_A, CV_DXT_INV_SCALE, dst->rows );
cvGetSubRect( dft_A, &tmp, cvRect(0,0,dst->cols,dst->rows) );
cvCopy( &tmp, dst );
cvReleaseMat( &dft_A );
cvReleaseMat( &dft_B );
}
ドライバプログラムも含めたソースファイルをここに置いておきます。参考にしてみてください。
高速ガウシアンフィルタプログラム
