script day.log

大学生がなんとなく始めた、趣味やら生活のことを記録していく。

フィルタリング(周波数領域でのフィルタリング~low-pass-filter)

以前は空間領域でのフィルタリングを行ったので,
次は周波数領域かなと考えています.
やることは周波数領域でのフィルタリングについての説明(?),
パワースペクトルを出力,ローパスフィルタの実装です.

makose3p1229.hatenablog.com
makose3p1229.hatenablog.com

周波数領域でのフィルタリングとは

周波数領域でのフィルタリングは,画像に対して2次元フーリエ変換を施して空間周波数領域へ変換したデータに
フィルタ関数を乗じ,2次元逆フーリエ変換を行うことで所望の画像を得る操作である.周波数領域におけるフィルタリングでは,
特定の周波数成分の強調や除去を容易に行うことが出来る.

画像の周波数領域でのフィルタリングの基礎となるのが2次元離散フーリエ変換(DFT)および逆変換(IDFT),
そして畳み込み定理である.濃度値f(m.n)を持つM \times N画素の画像に対する2次元フーリエ変換対は
以下の式で定義される.

\begin{align}
 F(u,v) = \sum^{M-1}_{m=0}\sum^{N-1}_{n=0}f(m,n)W_{M}^{mu}W_{N}^{nv} 
\end{align}
\begin{align}
 f(m,n) = \frac{1}{MN}\sum^{M-1}_{u=0}\sum^{N-1}_{v=0}f(u,v)W_{M}^{-mu}W_{N}^{-nv}
\end{align}

ただし,j^{2} = -1としW_{M} = \exp(-2\pi j/M),W_{N} = \exp(-2\pi j/N)で,f(m,n)及びF(u,v)には
以下の周期性を仮定する.

\begin{align}
 f(m + \alpha M,n + \beta N) &= f(m,n) \\
 F(u + \alpha M,v + \beta N) &= F(u,v),\quad \alpha \beta=0,\pm1,\pm2,\cdots 
\end{align}

上式から明らかなようにF(u,v)複素数値であり,その絶対値|F(u,v)|^{2}パワースペクトルと呼ばれ,
原画像中に各周波数成分がどの程度の強さで含まれているかを示す.さて,2次元DFTにおいて,
結果を格納する2次元配列F(u,v)は右上隅が直流成分に対応し,各隅の周辺が低周波成分,中央付近が高周波成分に対応する.
中心が直流成分で,中心から離れるにしたがって高周波成分になるように並べ替える.ただし,周波数成分の配置を並べ替えた後の
データを逆フーリエ変換しても,元の画像は得られない.

一般に空間領域でのフィルタリングは画像f(m,n)とフィルタ関数h(m,n)との畳み込み演算

\begin{align}
 \label{math-4}
 g(m,n) &= h * f(m,n) \\ 
        &= \sum^{M-1}_{\mu = 0}\sum^{N-1}_{\nu = 0}h(m-\mu,n-\nu)f(\mu,\nu)
\end{align}

で実行される.ただし,h,fともに上に述べた周期性を仮定する.ここでg,h,fの2次元DFTをG,H,Fとすると,
畳み込み定理により

\begin{align}
 G(u,v) = H(u,v)F(u,v)
\end{align}

と表される.つまり,空間領域における畳み込み演算が周波数領域では積で表現される.フィルタ処理が施された
画像を得るためにはG(u,v)の逆フーリエ変換を行えばよい.M \times N画素の画像に対する畳み込み演算による
フィルタリングの計算量はO(M^{2}N^{2})であるが,W_{M},W_{N}の周期性を利用した
高速フーリエ変換FFT)を用いることで周波数領域でのフィルタリングの計算量はO(MN \log MN)となり,
計算量の節約が可能である.

パワースペクトルの出力(C言語

パワースペクトルを表示するプログラムを作成する.
ただし,低周波成分が中央になるように並べ替えを行う.
プログラムでパワースペクトルのデータファイルを作成し,\texttt{GNUPLOT}を用いて,パワースペクトルを表示する.
流れとしてはFFT2して,FFTSHIFTして,パワースペクトルを計算(実数,虚数を2乗する)する.
あとはファイルに書き込んで,GNUPLOT使ってグラフで見たりとか,正規化してあげて,
画像として出力するとかしました.
log10してるのは、直流成分の2乗とかしてるので,ものすごく大きい値になるので,
スケールしてます.
以下に実装部分を示す.

FFT2(image, fft_data);

    //FFTSHIFT
    //1 to 3
    for(int i = 0; i < HEIGHT/2; i++){
        for (int j = 0; j < WIDTH/2; j++) {
            double temp_re, temp_im;
            temp_re = fft_data[i][j].re;
            fft_data[i][j].re = fft_data[HEIGHT/2 + i][WIDTH/2 + j].re;
            fft_data[HEIGHT/2 + i][WIDTH/2 + j].re = temp_re;
            temp_im = fft_data[i][j].im;
            fft_data[i][j].im = fft_data[HEIGHT/2 + i][WIDTH/2 + j].im;
            fft_data[HEIGHT/2 + i][WIDTH/2 + j].im = temp_im;
        }
    }

    //2 to 4
    for(int i = HEIGHT/2 + 1; i < HEIGHT; i++){
        for (int j = 0; j < WIDTH/2; j++) {
            double temp_re, temp_im;
            temp_re = fft_data[i][j].re;
            fft_data[i][j].re = fft_data[i - HEIGHT/2][WIDTH/2 + j].re;
            fft_data[i - HEIGHT/2][WIDTH/2 + j].re = temp_re;
            temp_im = fft_data[i][j].im;
            fft_data[i][j].im = fft_data[i - HEIGHT/2][WIDTH/2 + j].im;
            fft_data[i - HEIGHT/2 ][WIDTH/2 + j].im = temp_im;
        }
    }

    //power-spectrum
    fp = fopen("spectrum.dat", "wb");
    double max = log10(sqrt(pow(fft_data[0][0].re,2.0) + pow(fft_data[0][0].im,2.0)));
    double min = log10(sqrt(pow(fft_data[0][0].re,2.0) + pow(fft_data[0][0].im,2.0)));
    for (int i = 0; i < 256; i++){
        for (int j = 0; j < 256; j++){
            double spectrum = pow(fft_data[i][j].re,2.0) + pow(fft_data[i][j].im,2.0);
            fprintf(fp,"%3d %3d %f\n",i,j,spectrum);
            spectrum = log10(sqrt(spectrum));
            if (max < spectrum) {
                max = spectrum;
            }
            if (min > spectrum) {
                min = spectrum;
            }
        }
    }
    fclose(fp);

    fp = fopen("spectrum_log.dat", "wb");
    for (int i = 0; i < 256; i++){
        for (int j = 0; j < 256; j++){
            double spectrum = pow(fft_data[i][j].re,2.0) + pow(fft_data[i][j].im,2.0);
            spectrum = log10(sqrt(spectrum));
            spectrum = (unsigned char)(255 * (spectrum - min)/(max -min));
            out_image[i][j] = spectrum;
            fprintf(fp,"%3d %3d %f\n",i,j,spectrum);
        }
    }
fclose(fp);

low-paas filterの実装(C言語

FFTSHIFTで低周波成分を中央付近に固めたので,中心からある半径以降はゼロにして,
高周波成分をカットしようってことです.
ただ座標軸が円の中心が原点Oになっていないので注意して.
実装部分を以下に示す.distがフィルタの半径ですね.
dist = 14.1が今回使った画像の最適値でした.
distが最適値より小さすぎるとボケすぎるし,
distが最適値より大きいとノイズがカットできないので,
パワースペクトル見ながら,試しつつ,最適値に近づけていってください.

FFT2(image, fft_data);

    //fft-shift
    //1 to 3
    for (int i = 0; i < HEIGHT / 2; i++)
    {
        for (int j = 0; j < WIDTH / 2; j++)
        {
            double temp_re, temp_im;
            temp_re = fft_data[i][j].re;
            fft_data[i][j].re = fft_data[HEIGHT / 2 + i][WIDTH / 2 + j].re;
            fft_data[HEIGHT / 2 + i][WIDTH / 2 + j].re = temp_re;
            temp_im = fft_data[i][j].im;
            fft_data[i][j].im = fft_data[HEIGHT / 2 + i][WIDTH / 2 + j].im;
            fft_data[HEIGHT / 2 + i][WIDTH / 2 + j].im = temp_im;
        }
    }

    //2 to 4
    for (int i = HEIGHT / 2 + 1; i < HEIGHT; i++)
    {
        for (int j = 0; j < WIDTH / 2; j++)
        {
            double temp_re, temp_im;
            temp_re = fft_data[i][j].re;
            fft_data[i][j].re = fft_data[i - HEIGHT / 2][WIDTH / 2 + j].re;
            fft_data[i - HEIGHT / 2][WIDTH / 2 + j].re = temp_re;
            temp_im = fft_data[i][j].im;
            fft_data[i][j].im = fft_data[i - HEIGHT / 2][WIDTH / 2 + j].im;
            fft_data[i - HEIGHT / 2][WIDTH / 2 + j].im = temp_im;
        }
    }

    //low-pass
    double dist = 14.1;
    int hw = WIDTH / 2;
    int hh = HEIGHT / 2;
    for (int i = 0; i < WIDTH; i++)
    {
        for (int j = 0; j < HEIGHT; j++)
        {
            if (sqrt(pow(hw - i, 2) + pow(hh - j, 2)) > dist)
            {
                fft_data[i][j].re = 0;
                fft_data[i][j].im = 0;
            }
        }
    }

IFFT2(fft_data, out_image);