// グレイスケール画像を判別分析法(Discriminant analysis method)を使って２値化するプログラム
// プログラム名：binary_da.cpp

#include "eps-header.hh"

// 元画像(src)を判別分析法により２値化した画像を処理済画像(dst)に入れる関数
void binary_discriminant(
	class image& src,		// 元画像クラス(入力)
	class image& dst)		// 処理済画像クラス(出力)
{
	if (src.p == NULL) {
		printf("元画像がありません。");
		return;
	}

	printf("\nグレイスケール画像を判別分析法を使って２値化します。\n");

	int nr_lines  = src.get_nr_lines();		// 元画像のライン数   (縦の画素数)
	int nr_pixels = src.get_nr_pixels();	// 元画像のピクセル数 (横の画素数)
	dst.gmake(nr_lines, nr_pixels);			// 処理済画像の配列を確保

	// 判別分析法を使って最適な閾値を決める。
	double avef = 0.0;						// 全画素の輝度値の平均値
	for (int iy = 0; iy < nr_lines; iy ++) {
		for (int ix = 0; ix < nr_pixels; ix ++) {
			avef += src.p[iy][ix];
		}
	}
	avef /= (nr_lines * nr_pixels);

	// クラス１：閾値未満の画素領域
	// クラス２：閾値以上の画素領域
	int    n1;		int    n2;		// クラス１，２の画素数
	double ave1;	double ave2;	// クラス１，２の輝度値の平均値
	double var1;	double var2;	// クラス１，２の輝度値の分散

	double var_w;		// クラス内分散
	double var_b;		// クラス間分散

	int th_opt = MIN_GRAY_LEVEL + 1;	// 最適な閾値
	double max_var_ratio = 0.0;			// 分散比（クラス間分散／クラス内分散）の最大値

	for (int th = MIN_GRAY_LEVEL + 1; th <= MAX_GRAY_LEVEL; th ++) {
		n1   = n2   = 0;
		ave1 = ave2 = 0.0;
		var1 = var2 = 0.0;
		for (int iy = 0; iy < nr_lines; iy ++) {
			for (int ix = 0; ix < nr_pixels; ix ++) {
				int p = src.p[iy][ix];
				if (p < th) {
					n1 ++;		// クラス１
					ave1 += p;
					var1 += p * p;
				} else {
					n2 ++;		// クラス２
					ave2 += p;
					var2 += p * p;
				}
			}
		}
		if (n1 != 0) {
			ave1 /= n1;
			var1 = var1 / n1 - ave1 * ave1;
		}
		if (n2 != 0) {
			ave2 /= n2;
			var2 = var2 / n2 - ave2 * ave2;
		}

		var_w = (n1 * var1 + n2 * var2) / (n1 + n2);
		var_b = (n1 * (ave1 - avef) * (ave1 - avef) + n2 * (ave2 - avef) * (ave2 - avef)) / (n1 + n2);
		if (var_w > 1.0e-8) {
			if (max_var_ratio < var_b / var_w) {
				th_opt = th;
				max_var_ratio = var_b / var_w;
			}
		}
	}

	printf("判別分析法による最適な閾値: %4d\n", th_opt);

	for (int iy = 0; iy < nr_lines; iy ++) {
		for (int ix = 0; ix < nr_pixels; ix ++) {
			if (src.p[iy][ix] < th_opt) {
				dst.p[iy][ix] = MIN_GRAY_LEVEL;
			} else {
				dst.p[iy][ix] = MAX_GRAY_LEVEL;
			}
		}
	}
}

// プログラムの開始位置
int main()
{
	class image src;	// 元画像クラスを宣言
	class image dst;	// 処理済画像クラスを宣言

	src.gload();		// 元画像をBMPファイルから読み込む
						// BMPファイルを指定したい場合は，
						//   src.gload("パス付きBMPファイル名");
						// とする。

	// 元画像(src)を判別分析法により２値化した画像を処理済画像(dst)に入れる
	binary_discriminant(
		src,		// 元画像クラス(入力)
		dst);		// 処理済画像クラス(出力)

	dst.gsave();		// 処理済画像をBMPファイルに保存する
						// BMPファイルを指定したい場合は，
						//   dst.gsave("パス付きBMPファイル名");
						// とする。
}