Q.4. 大津の二値化

大津の二値化を実装せよ。大津の二値化とは判別分析法と呼ばれ、二値化における分離の閾値を自動決定する手法である。これはクラス内分散とクラス間分散の比から計算される。

・閾値t未満をクラス0, t以上をクラス1とする。
・w0, w1 ... 閾値tにより分離された各クラスの画素数の割合 (w0 + w1 = 1を満たす)
・S0^2, S1^2 ... 各クラスの画素値の分散
・M0, M1 ... 各クラスの画素値の平均値
とすると、

クラス内分散 Sw^2 = w0 * S0^2 + w1 * S1^2
クラス間分散 Sb^2 = w0 * (M0 - Mt)^2 + w1 * (M1 - Mt)^2 = w0 * w1 * (M0 - M1) ^2
画像全体の画素の分散 St^2 = Sw^2 + Sb^2 = (const)
以上より、分離度は次式で定義される。
分離度 X = Sb^2 / Sw^2 = Sb^2 / (St^2 - Sb^2)

となるので、

argmax_{t} X = argmax_{t} Sb^2

となる。すなわち、Sb^2 = w0 * w1 * (M0 - M1) ^2 が最大となる、閾値tを二値化の閾値とすれば良い。

入力	出力

回答

大津の二値化する部分のコードです。

ex4.c

void find_t_num(int image_gry[128][128]){
    double sb2_max = 0;
    for(int t_num = 0; t_num < 256; t_num ++){
        double sb2 = 0;
        double w0 = 0, w1 = 0, m0 = 0, m1 = 0;
        int class0 = 0, sum0 = 0, ave0 = 0;
        int class1 = 0, sum1 = 0, ave1 = 0;
        for(int i = 0; i < 128; i ++){
            for(int j = 0; j < 128; j ++){
                if(image_gry[i][j] < t_num){
                    sum0 = sum0 + image_gry[i][j];
                    class0 ++;
                }else{
                    sum1 = sum1 + image_gry[i][j];
                    class1 ++;
                }
            }
        }
        w0 = class0/(double)16384;
        w1 = class1/(double)16384;
        if(class0 == 0){
            m0 = 0;
        }else{
            m0 = sum0/class0;
        }
        if(class1 == 0){
            m1 = 0;
        }else{
            m1 = sum1/class1;
        }
        sb2 = w0 * w1 * (m0 - m1) * (m0 - m1);
        if(sb2 > sb2_max){
            sb2_max = sb2;
            t = t_num;
        }
    }
}

解説

上記のプログラムは閾値tを求める部分のプログラムです。
基本的には問題文にしたがって必要な値を求め、閾値を算出しました。
tの値を0から255まで変化させながらfor文を回していますが、最初のfor文は画素数の割合や平均を出すために各クラスの合計や要素数を求めています。
その後w0,w1,m0,m1を求めていくのですが、t=0の時とt=255の時はclass0、class1の値が片方0になってしまいますのでそのまま計算させるとFloating point exceptionが出てしまいます。したがって分子分母に0が来ないように条件分岐で処理しています。
常にsb^2の最大値とtの値を更新し、sb^2が最大になるようなtの値を最後算出して、その値を閾値としています。

完成したプログラム（全文）

1 ~ 3でのコメントでいろんなアドバイスを頂いたので、それを参考に全体的にプログラムを少し改良しました。
それぞれの機能は関数化しました。

ex4.c

# include <stdio.h>
# define N 256

FILE *infp;
FILE *outfp;
char magic[64];
char str[256];
int width;
int height;
int max = 0;
char readline[N] = {'\0'};
double image[128][384];
double image_gry_i[128][128];
int image_gry[128][128];
void convert_ppm_to_pgm(double image[128][384]);
void readImage();
void printImage(int image_gry[128][128], int width, int height, int max);
void find_t_num(int image_gry[128][128]);
int t_num = 0;
int t = 0;



int main(void) {
    readImage();
    convert_ppm_to_pgm(image);
    find_t_num(image_gry);
    printImage(image_gry, width, height, max);
    return 0;
}


void convert_ppm_to_pgm(double image[128][384]){
    for(int i = 0; i < height; i ++){
        int k = 0;
        for(int j = 0; j < width*3-2; j = j+3){
            image[i][j] = image[i][j] * 0.2126;
            image[i][j+1] = image[i][j+1] * 0.7152;
            image[i][j+2] = image[i][j+2] * 0.0722;
            image_gry_i[i][k] = image[i][j] + image[i][j+1] + image[i][j+2];
            image_gry[i][k] = (int)image_gry_i[i][k];
            k ++;
        }
    }
}

void readImage(){
    infp = fopen("../imori.ppm", "r");
    fgets(magic, N, infp);
    int num[4];
    for(int i = 0; i < 2; i ++){
        fscanf(infp, "%d", &num[i]);
    }
    width = num[0];
    height = num[1];
    fscanf(infp, "%d", &max);
    for(int i = 0; i < height; i ++){
        for(int j = 0; j < width*3; j ++){
            fscanf(infp, "%lf", &image[i][j]);
        }
    }
}

void printImage(int image_gry[128][128], int width, int height, int max){
    outfp = fopen("imori.pgm", "w");
    fprintf(outfp, "P2\n");
    fprintf(outfp, "%d ", width);
    fprintf(outfp, "%d\n", height);
    fprintf(outfp, "%d\n", max);
    for(int i = 0; i < height; i ++){
        for(int j = 0; j < width; j ++){
            if(image_gry[i][j] < t){
                fprintf(outfp, "0 ");
            }else{
                fprintf(outfp, "255 ");
            }
        }
        fprintf(outfp, "\n");
    }
    printf("%d\n", t);
}

void find_t_num(int image_gry[128][128]){
    double sb2_max = 0;
    for(int t_num = 0; t_num < 256; t_num ++){
        double sb2 = 0;
        double w0 = 0, w1 = 0, m0 = 0, m1 = 0;
        int class0 = 0, sum0 = 0, ave0 = 0;
        int class1 = 0, sum1 = 0, ave1 = 0;
        for(int i = 0; i < 128; i ++){
            for(int j = 0; j < 128; j ++){
                if(image_gry[i][j] < t_num){
                    sum0 = sum0 + image_gry[i][j];
                    class0 ++;
                }else{
                    sum1 = sum1 + image_gry[i][j];
                    class1 ++;
                }
            }
        }
        w0 = class0/(double)16384;
        w1 = class1/(double)16384;
        if(class0 == 0){
            m0 = 0;
        }else{
            m0 = sum0/class0;
        }
        if(class1 == 0){
            m1 = 0;
        }else{
            m1 = sum1/class1;
        }
        sb2 = w0 * w1 * (m0 - m1) * (m0 - m1);
        if(sb2 > sb2_max){
            sb2_max = sb2;
            t = t_num;
        }
    }
}

謝辞

コメント、アドバイスをくださった @fujitanozomu 様、ありがとうございました。

リンク

画像処理100本ノックを作ったった
 画像処理100本ノックC言語で挑戦してみる　イントロダクション