script day.log

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

TeXに苦労したこと(pdf作成)

さあ印刷しようと思ってtex->pdf変換して、
作成したDebian上のpdfをクラウドに一旦上げて、
Windowsでpdfを確認したところ、
日本語が埋め込まれてないって事態に遭遇しました。

とりあえずDebian上のpdfに対して、フォントが埋め込まれているか確認します。

//pdffonts インストールのため
apt install xpdf

//font 確認(emd:埋め込み)
pdffonts hoge.pdf 

私の場合は日本語系のフォントが埋め込まれてなかったですね。
ということで埋め込むのですが、mapを作成するのは面倒なので、
Ghostscriptにおまかせしました。

//font埋め込み
gs -q -dNOPAUSE -dBATCH -dPDFSETTINGS=/prepress -sDEVICE=pdfwrite -dEmbedAllFonts=true -sOutputFile=$YOUR_OUTPUT_FILE.pdf -f $YOUR_INPUT_FILE.pdf

一応確認
pdffonts hoge.pdf 

こんな感じの流れで一応解決しました。
しかし、クラウド上のviewerでは日本語表示出来なかったので、そちらも解決したいです。

フィルタリング(周波数領域でのフィルタリング~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);

フィルタリング(平滑化フィルタ)

差分フィルタは前回やりました.
makose3p1229.hatenablog.com

今回は平滑化フィルタ.

平滑化フィルタとは

画素の濃度値の積分を求めることにより画像中の雑音を除去し,濃度地が変化している部分をぼかす機能を
持たせることができる.なお,注目画素の処理結果に与える影響を2倍にすることでエッジの保存性を高めることも
可能である.このようなフィルタを平滑化フィルタと呼ぶ.
以下に平滑化フィルタの加重マトリックスを示す.

1/9 1/9 1/9
1/9 1/9 1/9
1/9 1/9 1/9

注目画素の強調なし

1/10 1/10 1/10
1/10 2/10 1/10
1/10 1/10 1/10

注目画素の強調あり

平滑化フィルタの実装(C言語

まずは注目画素を強調しない場合.
単純に足して,9で割るだけ.

unsigned char smoothing(unsigned char px0, unsigned char px1,unsigned char px2,unsigned char px3,unsigned char px4,unsigned char px5,unsigned char px6,unsigned char px7,unsigned char px8){
    int px[9] = {px0,px1,px2, px3,px4,px5,px6,px7,px8};

    int px_out = 0;
    for(int i = 0; i < 9; i++){
        px_out += px[i];
    }
    px_out /= 9;
    return (unsigned char)px_out;
}

続いて注目画素を強調する場合.
真ん中の画素値(i=4)だけ2倍して,単純に足して,10で割ってあげるだけ.

unsigned char smoothing(unsigned char px0, unsigned char px1,unsigned char px2,unsigned char px3,unsigned char px4,unsigned char px5,unsigned char px6,unsigned char px7,unsigned char px8){
    double px[9] = {px0,px1,px2, px3,px4,px5,px6,px7,px8};

    double px_out = 0.0;
    for(int i = 0; i < 9; i++){
        if (i == 4) {
            px[i] *= 2;
        }
        px_out += px[i];
    }
    px_out /= 10;
    return (unsigned char)px_out;
}

フィルタリング(Sobelフィルタ,ラプラシアンフィルタ)

さて、今回は空間フィルタについてです。

フィルタリングとは

画像のフィルタリングとは,画素値を用いて演算することによって,画像中のエッジの強調やノイズの除去,
特定の周波数成分のカットや強調などを行う処理である.画像のフィルタリングは,画像空間において
フィルタリングを行う空間領域でのフィルタリングと,画像にフーリエ変換を施して周波数領域に変換した後に
フィルタリングを行う周波数領域でのフィルタリングに大別出来る.

今回お話するのは空間領域でのフィルタリングです.

画像のフィルタリングとは,画素値を用いて演算することによって,画像中のエッジの強調やノイズの除去,
特定の周波数成分のカットや強調などを行う処理である.画像のフィルタリングは,画像空間において
フィルタリングを行う空間領域でのフィルタリングと,画像にフーリエ変換を施して周波数領域に変換した後に
フィルタリングを行う周波数領域でのフィルタリングに大別出来る.

空間領域でのフィルタリング

空間領域でのフィルタリングとは原画像の注目画素(m,n)の濃度値f(m,n)と注目画素の近傍画素の濃度値に対する
数値演算処理によって,新たな濃度値g(m,n)を求める操作である.ただし,近傍画素は画像全体に対して十分小さく,
具体的には3\times3画素,5\times5画素などである.変換後の濃度値g(m,n)

\begin{align}
\label{math-1}
 g(m,n) = \sum^{1}_{k=-1}\sum^{1}_{l = -1}C_{k,l}f(m + k, n + l)
\end{align}

で与えられるフィルタは線形フィルタと呼ばれ,C_{k,l}の値を変えることにより様々なフィルタを実現することが
できる.ただし,近傍画素は3\times3画素としC_{k.l}(-1\leq k,l \leq 1)は近傍の各画素に掛ける係数とする.


  • 差分フィルタ

領域の境界では濃度値の急激な変化が見られる.従って,その境界すなわち,エッジを抽出するには
画像に対して微分操作を行えばよい.いま,ディジタル画像を考えていることから,水平および垂直方向の
微分操作はそれぞれ次のように表され,差分操作によって求められる.

\begin{align}
 \Delta_{x}f(m,n) &= \frac{f(m,n+1)-f(m,n)}{n+1-n} = f(m,n+1) - f(m,n) \nonumber \\
 \Delta_{y}f(m,n) &= \frac{f(m+1,n)-f(m,n)}{m+1-m} = f(m+1,n) - f(m,n) \nonumber
\end{align}

Sobel(ソーベル)フィルタとは

微分操作を用いると雑音に弱いという欠点がある.そこで,近傍画素の濃度値も利用して安定した
1次微分が計算できるように拡張したSobelフィルタがある.
以下にSobelフィルタの加重マトリックスを示す.

-1 0 1
-2 0 2
-1 0 1

水平方向のエッジを強調

-1 -2 -1
0 0 0
1 2 1

垂直方向のエッジを強調

Sobelフィルタの実装(C言語

以下にSobelフィルタの適用部分を示す.
内容としては,各画素の濃度値を関数に渡して,returnしてるだけ.
ただ,水平方向と垂直方向の両方のエッジを強調するため,
出力画素で工夫(垂直方向と水平方向の2乗を足し、√をとって,絶対値を出力)してます.

unsigned char sobel_filtering(unsigned char px0, unsigned char px1,unsigned char px2,unsigned char px3,unsigned char px4,unsigned char px5,unsigned char px6,unsigned char px7,unsigned char px8){
    int sobel_x[9] = {-1,0,1,
                    -2,0,2,
                    -1,0,1};
    int sobel_y[9] = {-1,-2,-1,
                      0,0,0,
                      1,2,1};

    double px_x, px_y, px_out;
    px_x = px_y = px_out = 0;

    px_x = sobel_x[0] * px0 + sobel_x[1] * px1 + sobel_x[2] * px2 + sobel_x[3] * px3 + sobel_x[4] * px4 + sobel_x[5] * px5 + sobel_x[6] * px6 + sobel_x[7] * px7 + sobel_x[8] * px8;
    px_y = sobel_y[0] * px0 + sobel_y[1] * px1 + sobel_y[2] * px2 + sobel_y[3] * px3 + sobel_y[4] * px4 + sobel_y[5] * px5 + sobel_y[6] * px6 + sobel_y[7] * px7 + sobel_y[8] * px8;

    px_out = abs(sqrt(pow(px_x,2.0) + pow(px_y,2.0)));

    if(px_out > 255){
        return 255;
    }
    else{
        return (unsigned char)px_out;
    }
}

ラプラシアンフィルタとは

Sobelフィルタでは1次微分が計算されたわけだが,
ラプラシアンフィルタでは2次微分を計算する.
すると緩やかな濃度変化でも優れたエッジの検出を実行することが出来る.以下にラプラシアンフィルタの加重マトリックスを示す.

0 1 0
1 −4 1
0 1 0

4近傍

1 1 1
1 -8 1
1 1 1

8近傍

ラプラシアンフィルタの実装(C言語

以下にラプラシアンフィルタの適用部分を示す.
内容としては,各画素の濃度値を関数に渡して,returnしてるだけ.
真ん中の画素値(i=4)だけ−8倍してます.関数に渡すときや配列に格納するときに-8倍してもいいと思います.

unsigned char laplacian(unsigned char px0, unsigned char px1,unsigned char px2,unsigned char px3,unsigned char px4,unsigned char px5,unsigned char px6,unsigned char px7,unsigned char px8){
    int laplacian[9] = {px0,px1,px2, px3,px4,px5,px6,px7,px8};

    int px_out = 0;
    for (int i = 0; i < 9; i++) {
        if(i == 4){
            laplacian[i] *= -8;
        }
         px_out += laplacian[i];   
    }

    px_out = abs(px_out);
    if(px_out > 255){
        return 255;
    }
    else{
        return (unsigned char)px_out;
    }
}

次は平滑化について書きます.

2値化手法について(判別分析法)

モード法,P-タイル法をやれば、次は判別分析法(大津の2値化)ですよね.
makose3p1229.hatenablog.com

makose3p1229.hatenablog.com

判別分析法の考え方

モード法はヒストグラムが明確なピークを持たなければ適切な閾値を求めることが困難である.
しかし,判別分析法では自動的かつ統計的処理により最適な閾値を求めることが可能である.
判別分析法とは,画像の濃度値のヒストグラム閾値\thetaで2つのクラスに分割したとき,
2つのクラス間の分離が最も良くなるように閾値を定める方法である.
実際には2つのクラスのクラス間分散とクラス内分散の比を最大にするように閾値\thetaを決める.
与えられた画像データをk以下の濃度値を持つ画素と,
kより大きい濃度値を持つ画素の2つのグループに分け,
それぞれクラス1,クラス2とする.クラスiの画素数を全画素数で正規化した値を
\omega_i(k),平均濃度値をM_i(k),分散を\sigma_i^{2}(k)とする(i=1,2)
全画素の平均濃度値をM_Tとおくと,クラス間分散\sigma_B^{2}


\begin{align}
 \sigma_B^{2} &= \omega_1(M_1 - M_T)^{2} + \omega_2(M_2 - M_T)^{2} \nonumber \\
              &= \omega_1\omega_2(M_1 - M_2)^{2} \nonumber
\end{align}

で,クラス内分散\sigma_W^{2}


\begin{align}
 \sigma_W^{2} = \omega_1\sigma_1^{2} + \omega_2\sigma_2^{2} \nonumber
\end{align}

で与えられる.ここで,\sigma_B^{2}/\sigma_W^{2}を最大にするには
\sigma_B^{2}を最大にすればよく,したがって\sigma_B^{2}
最大にするkを閾値\thetaとすれば良い.
\sigma_B^{2}(k)を求めるには以下の計算をすれば良い.
画像のヒストグラムh(k)とすると,
濃度レベルkまでの濃度分布における0次モーメント\bar{\omega}(k)
及び1次モーメント\bar{M}(k)はそれぞれ次の漸化式


\begin{align}
 \bar{\omega}(k) &= \bar{\omega}(k-1) + \frac{h(k)}{MN} \nonumber \\
 \bar{M}(k) &= \bar{M}(k-1) + k\frac{h(k)}{MN} \nonumber
\end{align}

で表され,クラス間分散\sigma_B^{2}

\sigma_B^{2}(k) = \frac{[M_T\bar{\omega}(k) - \bar{M}(k)]^{2}}{\bar{\omega}(k)[1-\bar{\omega}(k)]}

となる.ただし,\bar{\omega}(0)=h(0)/MN,\bar{M}(0)=0である.
そして,\sigma_B^{2}(k)を最大にするkを閾値\thetaとする.

判別分析法の実装

私が実装した判別分析法による閾値の決定部分を以下に示す.

 int histgram[256] = {0};
 int sum = 0;
        for (int i = 0; i < HEIGHT; i++){
            for (int j = 0; j < WIDTH; j++){
                histgram[image[i][j]]++;
                sum += image[i][j];
            }
        }

        double ave = sum / (HEIGHT * WIDTH);
        double ave1,ave2;
        double max = 0.0;
    
        for(int i=0;i<256;i++){
            int cnt1 = 0;
            int cnt2 = 0;
            long data = 0;
            double brkup1 = 0.0;
            double brkup2 = 0.0;	
            double tmp = 0;
		
            for(int j=0;j<i;j++){
                cnt1 += histgram[j];
                data += histgram[j] * j;
            }

            if(cnt1 != 0){
                ave1 = (double)data / (double)cnt1;

                for(int j=0;j<i;j++){
                    brkup1 += pow( (j- ave1), 2 ) * histgram[j];
                }
                brkup1 /= (double)cnt1; 
            }
		
            data = 0;
		
            for(int j=i;j<256;j++){
                cnt2 += histgram[j];
                data += histgram[j] * j;
            }
            if(cnt2 != 0){
                ave2 = (double)data / (double)cnt2;

                for(int j=i;j<256;j++){
                    brkup2 += pow( (j - ave2), 2 ) * histgram[j];
                }
		
                brkup2 /= (double)cnt2;
            }
		
            double class1 = (cnt1 * brkup1 + cnt2 * brkup2);
            double class2 = cnt1 * pow( (ave1 - ave), 2 ) 
                          + cnt2 * pow( (ave2 - ave), 2 );
		
            tmp = class2 / class1;
		
            if(max < tmp){
                max = tmp;
                mode = i;
            }
        }

判別分析法の問題点

判別分析法ではどのようなヒストグラムでも一意に閾値を決定出来る.
しかし判別分析法も画像中のパターン領域と背景領域の画素数に大きな差がある場合,
パターン領域と背景領域が明確に分かれない閾値が設定される場合がある.
ラプラシアンフィルタを適用することで,隣接する画素との濃度値の変化分が大きい画素を
抽出することができる.抽出された画素とその周辺の画素を利用して濃度ヒストグラム
作成すると,フィルタを適用していない濃度ヒストグラムに比べて濃度ヒストグラム
双峰化する.ヒストグラムが双峰化されていると,判別分析法で精度の高い
最適なしきい値設定が可能となり、この問題を回避する事ができる.



2値化には飽きてきたので,次はフィルタリングしようかな.

2値化手法について(P-タイル法)

以前はモード法についてでしたので、今回はP-タイル法です.
makose3p1229.hatenablog.com

P-タイル法の考え方

P-タイル法では,画像の中から切り出したい対象物の面積をS_0,画像全体の面積をSとし,
その比率をp = S_0 / Sとする.そして,濃度ヒストグラムの濃度の低い方(あるいは高い方)からの画素数と,
画像全体の画素数の割合がpとなる濃度値を閾値とする.
この手法は,切り出すべき図形の面積がある程度分かっている場合によく用いられる.

P-タイル法の実装

私が作成したC言語でのP-タイル法の処理部分を以下に示します.
rateが切り出したい面積の割合です.
thが閾値となります.

        double rate = 0.5;
        int sum = 0;
        for (int i = 0; i < HEIGHT; i++){
            for (int j = 0; j < WIDTH; j++){
                histgram[image[i][j]]++;
            }
        }

        for (int i = 0; i < 256; i++){
            sum += histgram[i];
        }
    
        int th = sum * rate;
        sum = 0;
        for (int i = 0; i < 256; i++){
            sum += histgram[i];
            if(sum > th){
                th = i;
                break;
            }
        }

次は判別分析法ですかね.
微分ヒストグラム法や動的閾値法はどうしましょう…

2値化手法について(モード法~閾値の自動決定)

濃淡の表現できない媒体に画像を表示するときや,画像の特徴を解析するときなどには,
濃淡画像を2値(白・黒)の濃度値を持つ画像に変換する.M \times N画素で構成される画像の2値化は
次式の閾値処理によって行われる.

\begin{align}
 g(m,n)=\left\{
\begin{array}{l}
 255 \qquad f(m,n)\geq\theta \\
 0 \qquad\quad f(m,n)<\theta
\end{array}
\right.
\end{align}

ただし,f(m,n),g(m,n)はそれぞれ入力画像データ,作成される2値画像データとする.
ここで閾値\thetaをいかに決めるかが問題となり,その代表的な手法としてモード法がある.

モード法の考え方

与えられた画像の濃度ヒストグラムが2つのピークを持つ場合,
この2つの山の間の谷のところに閾値\thetaを決めればよい.
しかし、閾値を決めるには、濃度ヒストグラムを確認して,谷を見つける必要がある.

閾値の自動決定

Pさん、谷を見つければ良いんですよ!
今回私が谷を見つけるために使用したのは極値です.
最急降下法を実装するのは面倒だったので…
簡単に言うと、極大値が濃度ヒストグラムの山にあたり、
極小値が濃度ヒストグラムの谷に当たるわけですね…
そして極大値と極小値の差を求め、ある程度の高さがあれば、
その時の極小値を谷と考え,閾値に設定するという考え方です.
以下が私が作成したCプログラムです.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

#define HEIGHT 256
#define WIDTH 256

int maxNo(int left, int right, int data[])
{
    int maxNum = left;
    int max = data[left];
    for (int i = left; i < left + right && i < 256; i++)
    {
        if (data[i] > max)
        {
            max = data[i];
            maxNum = i;
        }
    }
    return maxNum;
}

int minNo(int left, int right, int data[])
{
    int minNum = left;
    int min = data[left];
    for (int i = left; i < right; i++)
    {
        if (min > data[i])
        {
            min = data[i];
            minNum = i;
        }
    }
    return minNum;
}

int main(int argc, char *argv[])
{
    unsigned char image[HEIGHT][WIDTH];
    unsigned char image_w[HEIGHT][WIDTH];
    int histgram[256] = {0};
    int mode;
    FILE *fp;
    int maxRNo[256];
    int minRNo[256];
    int subRNo[256];
    int range = 16;
    int thr = 150;

    if (argc < 2)
    {
        printf("don't read filename\n");
        return 0;
    }
    else if (argc > 2)
    {
        printf("many arguments\n");
        return 0;
    }
    else
    {
        fp = fopen(argv[1], "rb");
        fread(image, sizeof(unsigned char), HEIGHT * WIDTH, fp);
        fclose(fp);

        for (int i = 0; i < HEIGHT; i++)
        {
            for (int j = 0; j < WIDTH; j++)
            {
                histgram[image[i][j]]++;
            }
        }

        int cnt = 0;
        for (int i = 0; i < 256 - range; i++)
        {
            if (maxNo(i, range, histgram) == maxNo(i + 1, range, histgram))
            {
                int check = 0;
                for (int j = 1; j < range; j++)
                {
                    if (maxNo(i, range, histgram) == maxNo(i + j, range, histgram))
                    {
                        check++;
                    }
                }
                if (check == range - 1)
                {
                    maxRNo[cnt] = maxNo(i, range, histgram);
                    cnt++;
                }
            }
        }

        minRNo[0] = minNo(0, maxRNo[0], histgram);
        for (int i = 0; i < 256; i++)
        {
            if (maxRNo[i + 1] == 0)
            {
                break;
            }
            minRNo[i + 1] = minNo(maxRNo[i], maxRNo[i + 1], histgram);
        }

        int subCnt = 0;
        for (int i = 0; i < 256; i++)
        {
            if (histgram[maxRNo[i]] - histgram[minRNo[i]] >= thr)
            {
                subRNo[subCnt] = minRNo[i];
                subCnt++;
            }
        }

        if (subRNo[0] > 0)
        {
            if (subRNo[1] > 0)
            {
                printf("many peaks\n");
            }
            else
            {
                mode = subRNo[0];
                printf("閾値 = %d\n", mode);
            }
        }
        else
        {
            printf("don't appear peak\n");
            return 0;
        }

        for (int i = 0; i < HEIGHT; i++)
        {
            for (int j = 0; j < WIDTH; j++)
            {
                if (image[i][j] > mode)
                {
                    image_w[i][j] = 255;
                }
                else
                {
                    image_w[i][j] = 0;
                }
            }
        }


        return 0;
    }
}

次は判別分析法かP-タイル法を書くつもり…