[C/C++]平均誤差最小法による二値化の実装[コード付]

C++

どうもー
今回も前回までの続きで二値化を行なっていきます!

今回は条件付き閾値決定法の一種である平均誤差最小法を用いて二値化を行なっていきます!

条件付き閾値決定法とは?

ここは前回も説明しましたが今回もやっておきます。

独立閾値決定法では、対象の画素を白にするか黒にするかをその画素の濃度値のみで決定していました。

それに対して条件付き閾値決定法その画素だけでなく、周辺の画素の濃度値も参照しながら閾値を決定します。

すなわち濃度値を持つ各画素を二値化するという操作をある範囲で行うことによって、各画素を単独で二値化するよりも画質が良くなるんじゃね?っていう考え方です。
代表的な手法として

  • 平均値制限法
  • 平均誤差最小法
  • 誤差拡散法

があげられます。
その中でも今回の記事では平均誤差最小の実装をしていきます!

平均誤差最小法について

平均誤差最小法はその名の通り、近傍画素を二値化した際の誤差を最小にすることを目標とした方法です。つまり閾値も誤差が小さくなるようなものに決定します。

対象の画像をグレースケール化して、その上で濃度値を正規化したものf(x,y)を考えます。ここでいう正規化とは階調数-1で割る、つまりグレースケール化した画像の全ての画素の濃度値を255で割った値を保持しているのがf(x,y)ってことですね。x,yはその画素の座標のことです。
ここで二値化した画像をG(x,y)とすれば(黒なら1、白なら0)誤差err(x,y)は

err(x,y)=f(x,y)-G(x,y)

と表すことができます。
つまりこの誤差err(x,y)が小さくなるような閾値を設定していくわけです。

でも誤差を計算するにはそもそも二値化が完了している必要があるよね?だとしたらその後に閾値を求めても無駄なような、、

そうなんです。
なので今回のアルゴリズムは、
ある画素の閾値を求めるために、既に二値化が完了した画素の誤差を用いる
といった表現の方がわかりやすいですかね。

より具体的に説明すると、
このアルゴリズムでは毎度のよう画像の各画素をラスタ走査によって注目していきます。
二値化の順番もラスタ走査の通りにやっていきます。つまりある画素を注目するとき、その画素よりy座標が小さい画素と、y座標が同じでもx座標が小さい画素は既に二値化が完了しているわけです。

これらの画素の誤差を用いようぜ!ってことです。
それでは閾値の式をどうぞ

ここでのprevは先ほど説明したようにラスタ走査の順番において、その画素より手前の画素のことです。
またαは重みの変数であり、

上の画像のように、注目画素に近ければ重く、遠ければ軽くなるような数値にします。
上の重み付けが一般的なものであり、具体的に式に起こすと

こんな感じですね。注目画素からのマンハッタン距離によって重みが変わる感じです。
0以下の値は0とみなします。
今回の実装ではこの重みの付け方を採用してます。

コードの実装

読み書きの関数

#include <iostream>
#include <opencv4/opencv2/opencv.hpp>
#define WIDTH 640
#define HEIGHT 480
using namespace std;


void readImage_gray( const char *fname, unsigned char img[HEIGHT][WIDTH])
{
    cv::Mat ReadImage;

    ReadImage = cv::imread(fname,0);

    for( int w=0; w<WIDTH; w++ ){
        for( int h=0; h<HEIGHT; h++ ){
            img[h][w] = ReadImage.data[w+h*WIDTH];
        }
    }
}

void writeImage_binary( const char *fname, bool binary_img[HEIGHT][WIDTH])
{
    cv::Mat writeImage( HEIGHT, WIDTH, CV_8U );

    for( int w=0; w<WIDTH; w++ ){
        for( int h=0; h<HEIGHT; h++ ){
            if(binary_img[h][w])
                writeImage.data[w+h*WIDTH]=0;
            else
                writeImage.data[w+h*WIDTH]=255;
        }
    }

    cv::imwrite( fname , writeImage );
}

もし上の関数でわからないことがあれば下の記事を参考にしてください!

閾値を計算する関数

まず閾値を計算して、その値を返す関数を作成します。
重みを乗算するあたりがややこしいんです、、

float getThreshold(int y,int x,float err[HEIGHT][WIDTH]){
    float threshold=0.5;
    
    if(y==0 && x==0)
        return threshold;
    
    float numerator=0;//分子
    float denominator=0;//分母
    
    for(int h=0; h<=y; h++){
        if(h!=y){//最後の行ではない時
            for(int w=0; w<WIDTH; w++){
                if((abs(h-y)+abs(w-x))<=4){
                    int weight=(4-abs(h-y)-abs(w-x))*2+1;
                    denominator += weight;
                    numerator += weight*err[h][w];
                }
            }
        }else{//最後の行の時
            for(int w=0; w<x; w++){
                if((abs(h-y)+abs(w-x))<=4){
                    int weight=(4-abs(h-y)-abs(w-x))*2+1;
                    denominator += weight;
                    numerator += weight*err[h][w];
                }
            }
        }
    }
    
    threshold += numerator/denominator;
    return threshold;
}

とりあえず関数の型はfloatです。今回は画素の濃度値を正規化したものを扱う関係で、閾値の範囲も0から1になります。てことで少数も扱えるfloatにしています。

そんで引数ですがint x,yは対象となる画素の座標ですね。
float err[HEIGHT][WIDTH]は各画素の誤差が格納されている配列です。各画素といってもx,yより前の画素の誤差しか入ってないですけどね。

とりあえず最初の文について

    float threshold=0.5;
    
    if(y==0 && x==0)
        return threshold;

    float numerator=0;//分子
    float denominator=0;//分母

一旦のちの返り値となるthresholfdを0.5で定義します。
そんでx,yが共に0のとき、つまりラスタ走査の最初ですね。この時はそれ以前に二値化した画素がないので0.5のまま閾値を返します。そのためのif文てことです。

その次は閾値計算の右の項の分数のために分母、分子を格納する変数を0で定義してます。

次は二重ループについて

    for(int h=0; h<=y; h++){
        if(h!=y){//最後の行ではない時
            for(int w=0; w<WIDTH; w++){
                if((abs(h-y)+abs(w-x))<=4){
                    int weight=(4-abs(h-y)-abs(w-x))*2+1;
                    denominator += weight;
                    numerator += weight*err[h][w];
                }
            }
        }else{//最後の行の時
            for(int w=0; w<x; w++){
                if((abs(h-y)+abs(w-x))<=4){
                    int weight=(4-abs(h-y)-abs(w-x))*2+1;
                    denominator += weight;
                    numerator += weight*err[h][w];
         }
            }
        }
    }

この部分で、画像の最初の画素から注目画素までラスタ走査していきます。

外側のfor文はy座標の走査のためのループです。
そんでいきなりif文なんですが、これは現在走査している画素のy座標が、注目画素のy座標と同じかどうかを判定しています。
なんでこんな判定をしているかというと、それによってx座標の走査する範囲が変わってくるからですね。
基本的にいつも通りw<WIDTHの範囲で走査すれば良いですが、注目画素と同じy座標の行を走査するときは注目画素を超えないようにしなくてはなりません。
なのでif文とelse文の後のx座標に関するfor文の継続条件が異なっています。

if(h!=y){//最後の行ではない時
            for(int w=0; w<WIDTH; w++){
                ///////
            }
        }else{//最後の行の時
            for(int w=0; w<x; w++){
               ////////
            }
        }

そして次は二重ループの中身についてです。
上で説明した通りx座標に関するループは二つありますが、その中身は同じなので片方だけを説明していきます。

                if((abs(h-y)+abs(w-x))<=4){
                    int weight=(4-abs(h-y)-abs(w-x))*2+1;
                    denominator += weight;
                    numerator += weight*err[h][w];
                }

まずif文で重みが0で無い範囲にある画素なのかを判定します。
上でも説明した通り、重みは注目画素からのマンハッタン距離で決まるので絶対値関数abs()を使って上のように定義しています。

if文の中ではまず重みを計算してweightに格納します。
その後最初に定義した分母と分子の変数にそれぞれ重みや、重みと誤差の積を足し合わせています。

最後です。

    threshold += numerator/denominator;
    return threshold;

thresholdにループ内で計算した右の項を足してそれをreturnしています。

平均誤差最小法に基づく二値化を行う関数

void MinimizedAverageError(unsigned char img[HEIGHT][WIDTH], bool binary_img[HEIGHT][WIDTH]){
    float err[HEIGHT][WIDTH]={0};
    float f[HEIGHT][WIDTH]={0};  
    float G[HEIGHT][WIDTH]={0};
    
    for(int h=0; h<HEIGHT; h++){
        for(int w=0; w<WIDTH; w++){
            f[h][w]=(float)img[h][w]/255;
        }
    }
    
    for(int h=0; h<HEIGHT; h++){
        for(int w=0; w<WIDTH; w++){
            
            if(f[h][w]>getThreshold(h, w, err)){
                G[h][w]=1;//黒
                binary_img[h][w]=false;
            }else{
                G[h][w]=0;//白
                binary_img[h][w]=true;
            }
            err[h][w]=G[h][w]-f[h][w];
        }
    }
    
}

まず引数ですがグレースケール画像が格納されている二次元配列imgと二値化した画像を入れるbool型配列binary_imgです。

本文で最初に定義しているerr、f、Gはそれぞれ
画像の各画素の誤差、濃度値を正規化した値、二値化した値を格納する配列です。説明の時と同じですね。

そして最初のループ

    for(int h=0; h<HEIGHT; h++){
        for(int w=0; w<WIDTH; w++){
            f[h][w]=(float)img[h][w]/255;
        }
    }

ここは正規化を行なっておます。キャストを忘れずに。

for(int h=0; h<HEIGHT; h++){
        for(int w=0; w<WIDTH; w++){
            
            if(f[h][w]>getThreshold(h, w, err)){
                G[h][w]=1;//黒
                binary_img[h][w]=false;
            }else{
                G[h][w]=0;//白
                binary_img[h][w]=true;
            }
            err[h][w]=G[h][w]-f[h][w];
        }
    }

次の二重ループは毎度おなじみの二値化を行う部分ですね。
if文で閾値を先ほどの関数の返り値にしてます。
また一つの画素の二値化を終えたら誤差を計算するようにしています。

実際に使ってみた

それでは実行してみましょう。
今回も用いるのはこちらの画像

そして上のプログラムを実行して生成されたのがこちら

おおー、なんか独特な雰囲気がありますね。
ちなみに前回作成した平均値制限法による二値化の結果はこんな感じでした。

どちらも構造は残ってますが印象が違いますね。

全体のコード

#include <iostream>
#include <opencv4/opencv2/opencv.hpp>
#define WIDTH 640
#define HEIGHT 480
using namespace std;
using namespace cv;

float getThreshold(int y,int x,float err[HEIGHT][WIDTH]){
    float threshold=0.5;
    
    if(y==0 && x==0)
        return threshold;
    
    float numerator=0;//分子
    float denominator=0;//分母
    
    for(int h=0; h<=y; h++){
        if(h!=y){//最後の行ではない時
            for(int w=0; w<WIDTH; w++){
                if((abs(h-y)+abs(w-x))<=4){
                    int weight=(4-abs(h-y)-abs(w-x))*2+1;
                    denominator += weight;
                    numerator += weight*err[h][w];
                }
            }
        }else{//最後の行の時
            for(int w=0; w<x; w++){
                if((abs(h-y)+abs(w-x))<=4){
                    int weight=(4-abs(h-y)-abs(w-x))*2+1;
                    denominator += weight;
                    numerator += weight*err[h][w];
                }
            }
        }
    }
    
    threshold += numerator/denominator;
    return threshold;
}


void MinimizedAverageError(unsigned char img[HEIGHT][WIDTH], bool binary_img[HEIGHT][WIDTH]){
    float err[HEIGHT][WIDTH]={0};
    float f[HEIGHT][WIDTH]={0};
    
    
    
    float G[HEIGHT][WIDTH]={0};
    
    for(int h=0; h<HEIGHT; h++){
        for(int w=0; w<WIDTH; w++){
            f[h][w]=(float)img[h][w]/255;
        }
    }
    
    for(int h=0; h<HEIGHT; h++){
        for(int w=0; w<WIDTH; w++){
            
            if(f[h][w]>getThreshold(h, w, err)){
                G[h][w]=1;//黒
                binary_img[h][w]=false;
            }else{
                G[h][w]=0;//白
                binary_img[h][w]=true;
            }
            err[h][w]=G[h][w]-f[h][w];
        }
    }
    
}

int main(int argc, const char * argv[]) {
    unsigned char img[HEIGHT][WIDTH]={0};
    bool binary_img[HEIGHT][WIDTH] = {0};
    cv::Mat showImage(HEIGHT, WIDTH, CV_8U );
    readImage_gray("hana.png", img);
    MinimizedAverageError(img, binary_img);

    writeImage_binary("MinimizedAverageError.png", binary_img);
    
    return 0;
}


コメント

タイトルとURLをコピーしました