/*	
 *	clipimg.c
 *
 *	ユーティリティ：画像処理
 *
 *	CLiP - Common Library for P/ECE
 *	Copyright (C) 2001-2015 Naoyuki Sawa
 *
 *	* Sun Sep 06 23:32:49 JST 2015 Naoyuki Sawa
 *	- 1st リリース。
 *	- discriminant_analysis()を作成しました。
 *	* Mon Sep 07 22:07:58 JST 2015 Naoyuki Sawa
 *	- image_labeling()を追加しました。
 *	* Sun Sep 20 18:19:03 JST 2015 Naoyuki Sawa
 *	- discriminant_analysis()の処理の一部を、32ビット整数処理⇒実数処理に変更しました。
 */
#include "clip.h"
/*****************************************************************************
 *	判別分析法(大津の二値化)
 *****************************************************************************/
//□参考資料
//│①「画像処理ソリューション」(http://imagingsolution.blog107.fc2.com/)さんの、判別分析法(大津の二値化)」(http://imagingsolution.blog107.fc2.com/blog-entry-113.html)の記事より引用:
//│＞判別分析法【discriminant analysis method】は大津の二値化とも言われ、分離度(separation metrics)という値が最大となるしきい値を求め、自動的に二値化を行う手法です。
//│＞分離度はクラス間分散(between-class variance)とクラス内分散(within-class variance)との比で求める事ができ、以下の様に求める。
//│＞しきい値ｔで二値化したとき、しきい値よりも輝度値が小さい側(黒クラス)の画素数をω1、平均をｍ1、分散をσ1、輝度値が大きい側(白クラス)の画素数をω2、平均をｍ2、分散をσ2、
//│＞画像全体の画素数をωt、平均をｍt、分散をσtとしたときクラス内分散σw2は
//│＞　　σw^2 = (ω1・σ1^2 + ω2・σ2^2) / (ω1 + ω2)
//│＞クラス間分散σb^2は
//│＞　　σb^2 = (ω1・(ｍ1 - ｍt)^2 + ω2・(ｍ2 - ｍt)^2) / (ω1 + ω2) = (ω1・ω2・(ｍ1 - ｍ2)^2) / ((ω1 + ω2)^2)
//│＞としてあらわす事ができる。ここで、全分散(total variance)σtは
//│＞　　σt^2 = σb^2 + σw^2
//│＞としてあらわす事ができることから、求めるクラス間分散とクラス内分散との比である分離度は
//│＞　　(σb^2) / (σw^2) = (σb^2) / (σt^2 - σb^2)
//│＞となり、この分離度が最大となるしきい値ｔを求めればよい。ここで、全分散σtはしきい値に関係なく一定なので、クラス間分散σb^2が最大となるしきい値を求めればよい事が分かる。
//│＞さらにクラス間分散の式の分母もしきい値に関係なく一定なので、クラス間分散の分子
//│＞　　ω1・ω2・(ｍ1 - ｍ2)^2
//│＞が最大となるしきい値ｔを求めればよい。
//│＞結局、標準偏差とか関係なく、黒、白それぞれの領域のヒストグラムから、画素数ωと輝度値の平均値ｍから上記の値が最大となるしきい値ｔをしらみつぶしに求めればいいので、以外と簡単…
//│②「CodeDaily」(http://ayb52895.hatenablog.com/)さんの、「画像の2値化処理～判別分析法～」(http://ayb52895.hatenablog.com/entry/2013/05/24/091131)の記事
//記事①は理論の説明のみで，記事②は実装例のみなので、両方を合わせて参考にさせて頂きました。
//記事①の式変形はなぜそうなるのかまだ良く理解できていない部分が有るのですが、自分で実装してみたコードと記事②の実装例がほぼ同じになったので、動作は多分合っていると思います。
/*---------------------------------------------------------------------------*/
//画像の2値化処理を行うための、しきい値を求める。
//[in]
//	p		画像のアドレス
//	n		画像のバイト数
//[out]
//	戻り値		しきい値
//[note]
//	- 'しきい値以下'と'しきい値超過'に分けるか，'しきい値未満'と'しきい値以上'に分けるかは、(理論上は)どちらでも構いません。
//	  当関数の実装は、{0,1,2,...,14,15}に対して'7'を返すので、'しきい値以下'と'しきい値超過'に分ける方が良いかも知れません。
//<使用例>
//	│render_object(&render,0,0,0,0,DISP_X,DISP_Y,DRW_NOMAL);
//	│t=discriminant_analysis(vbuff,sizeof vbuff);
//	│for(i=0;i<sizeof vbuff;i++){
//	│  vbuff[i]=(vbuff[i]<=t)?0:15;//'高速切替・16階調表示'の使用を想定
//	│}
//{{2015/09/20変更:discriminant_analysis()の処理の一部を、32ビット整数処理⇒実数処理に変更しました。
//int discriminant_analysis(const uint8_t* p, int n) {
//	int t = 0;
//	if(n) {
//		int w[256];//スタックを1キロバイト使う事に注意せよ。
//		int w1, n1, m1, w2, n2, m2, i, x, max;
//		memset(w, 0, sizeof w);
//		w1 = n;
//		n1 = 0;
//		do {
//			x = *p++;
//			n1 += x;
//			w[x]++;
//		} while(--n);
//		w2 = 0;
//		n2 = 0;
//		max = 0;
//		i = 255;
//		do {
//			if(w1 && w2) {
//				m1 = n1 / w1;
//				m2 = n2 / w2;
//				x = w1 * w2 * (m1 - m2) * (m1 - m2);
//				if(x > max) {
//					t   = i;
//					max = x;
//				}
//			}
//			w1 -= w[i];
//			w2 += w[i];
//			n1 -= w[i] * i;
//			n2 += w[i] * i;
//		} while(--i >= 0);
//	}
//	return t;
//}
//↓2015/09/20変更:discriminant_analysis()の処理の一部を、32ビット整数処理⇒実数処理に変更しました。
int discriminant_analysis(const uint8_t* p, int n) {
	int sel = 0;
	if(n) {
		int w[256];//スタックを1キロバイト使う事に注意せよ。
		int w1, n1, w2, n2, i;
		double max;
		memset(w, 0, sizeof w);
		w1 = n;
		n1 = 0;
		do {
			int c = *p++;
			n1 += c;
			w[c]++;
		} while(--n);
		w2 = 0;
		n2 = 0;
		max = 0.0;
		i   = 255;
		do {
			if(w1 && w2) {
				double m1 = (double)n1 / (double)w1;	//┐
				double m2 = (double)n2 / (double)w2;	//│32ビット整数で取り扱うと、約200x200ドット以上の画像に適用するとオーバーフローする。
				double d  = m1 - m2;			//│200x200ドット程度の画像には対応出来ないと実用的でないので、実数で取り扱う事にした。
				d = d * d * (double)w1 * (double)w2;	//┘
				if(d > max) {
					max = d;
					sel = i;
				}
			}
			w1 -= w[i];
			w2 += w[i];
			n1 -= w[i] * i;
			n2 += w[i] * i;
		} while(--i >= 0);
	}
	return sel;
}
//}}2015/09/20変更:discriminant_analysis()の処理の一部を、32ビット整数処理⇒実数処理に変更しました。
/****************************************************************************
 *	ラベリング処理
 ****************************************************************************/
//aの属すグループの代表に向かって経路圧縮(代表を返す)
static int image_labeling_compress(uint16_t* parents, int a) {
	//aがグループの代表でなければ…
	while(parents[a] != a) {
		//aの属すグループの代表のラベルに置き換える。
		a = parents[a] = parents[parents[a]];
	}
	return a;
}
/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
//aの属すグループとbの属すグループを併合(併合後の代表を返す)
static int image_labeling_link(uint16_t* parents, int a, int b) {
	a = image_labeling_compress(parents, a);
	b = image_labeling_compress(parents, b);
	return (a <= b) ? (parents[b] = a) : (parents[a] = b);
}
/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
//番号とびとびなラベルを0,1,2,...に貼り替え
static int image_labeling_relabel(uint16_t* parents, int regions) {
	int index = 0, a;
	for(a = 0; a < regions; a++) {
		//aがグループの代表ならば…
		if(parents[a] == a) {
			//aに新しいラベルを割り当てる。
			parents[a] = index++;
		//aがグループの代表でなければ…
		} else {
			//aの属すグループの代表と同じラベルを割り当てる。
			parents[a] = parents[parents[a]];
		}
	}
	return index;
}
/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
//画像のラベリング処理を行う。
//[in]
//	image[w*h]		画像のメモリアドレス。
//	label[w*h]		当関数が成功したら、画像の各画素に対応するラベルを格納する。
//	diag			0を指定すると4近傍、0以外の値を指定すると8近傍で処理を行う。
//	parents[parents_max]	当関数が処理の途中で使用するワークエリア。最終的に意味の有る値は格納されません。
//[out]
//	戻り値			成功したら、領域の数(=ラベルの最大値+1)を返す。
//				失敗したら、0を返す。
//				当関数が失敗する要因は、ワークエリア不足です。parents[parents_max]を増やして、再度呼び出してください。
//				最悪な状況(市松模様)では、parents[w*h]のワークエリア容量が必要となります。
//[note]
//	- □参考資料
//	  │「A440Hz ブログ」(http://d.hatena.ne.jp/wosugi/)さんの、「ほぼワンパスなラベリング処理」(http://d.hatena.ne.jp/wosugi/20110628)の記事を参考にさせて頂きました。
//	- 上記資料の実装はparents[]を固定サイズで確保していますが、実際に適切な容量は個々のアプリケーションによって異なるので、アプリケーション側で確保するようにしました。
int image_labeling(const uint8_t* image/*[w*h]*/, int w, int h, uint16_t* label/*[w*h]*/, int diag/*0=4近傍,1=8近傍*/, uint16_t* parents/*[parents_max]*/, int parents_max) {
	int regions = 0, x, y, c, flagA, flagB, flagC, flagD, index;
	for(y = 0; y < h; y++) {
		for(x = 0; x < w; x++) {
			//隣接画素(4近傍,又は,8近傍)との連結チェック
			c = image[(w * y) + x];
			flagA = 0;            if( x                ) { flagA = (image[(w *  y     ) + (x - 1)] == c); }		//左　
			flagB = 0;            if(                 y) { flagB = (image[(w * (y - 1)) +  x     ] == c); }		//　上
			flagC = 0; if(diag) { if( x            && y) { flagC = (image[(w * (y - 1)) + (x - 1)] == c); } }	//左上
			flagD = 0; if(diag) { if((x < (w - 1)) && y) { flagD = (image[(w * (y - 1)) + (x + 1)] == c); } }	//右上
			//着目画素と連結画素を併合
			if(regions >= parents_max) { return 0; }	//ワークエリア不足。parents[parents_max]を増やして再試行せよ。
			parents[regions] = regions;
			index = regions++;
			if(flagA | flagB | flagC | flagD) {
				regions--;
				if(flagA) { index = image_labeling_link(parents, index, label[(w *  y     ) + (x - 1)]); }
				if(flagB) { index = image_labeling_link(parents, index, label[(w * (y - 1)) +  x     ]); }
				if(flagC) { index = image_labeling_link(parents, index, label[(w * (y - 1)) + (x - 1)]); }
				if(flagD) { index = image_labeling_link(parents, index, label[(w * (y - 1)) + (x + 1)]); }
			}
			label[(w * y) + x] = index;
		}
	}
	//再ラベリング
	regions = image_labeling_relabel(parents, regions);
	for(y = 0; y < h; y++) {
		for(x = 0; x < w; x++) {
			label[(w * y) + x] = parents[label[(w * y) + x]];
		}
	}
	return regions;
}
/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
//{{使用例
#if 0
//「画像処理ソリューション」(http://imagingsolution.blog107.fc2.com/)さんの、
//「ラベリング処理アルゴリズム」(http://imagingsolution.blog107.fc2.com/blog-entry-193.html)の記事で提示されている、
//ラベリング処理をすると失敗しやすい特徴的な画像の例から二例(◇Level.1、◇Level.2)と、
//必要なワークエリアの容量が最大となる最悪ケースである、一例(◇市松模様)を使って、検証を行う事にしました。
#define _ 111
#define X 222
//◇Level.1　⇒　ok
static uint8_t image1[9][9]={
  {_,_,_,_,_,_,_,_,_},
  {_,X,_,_,_,_,_,X,_},
  {_,_,X,_,_,_,X,_,_},
  {_,_,_,X,_,X,_,_,_},
  {_,_,_,_,X,_,_,_,_},
  {_,_,_,X,_,X,_,_,_},
  {_,_,X,_,_,_,X,_,_},
  {_,X,_,_,_,_,_,X,_},
  {_,_,_,_,_,_,_,_,_}};
//◇Level.2　⇒　ok
static uint8_t image2[9][9]={
  {_,_,_,_,_,_,_,_,_},
  {_,_,X,X,_,X,X,_,_},
  {_,X,_,_,X,_,_,X,_},
  {_,_,_,X,_,X,_,_,_},
  {_,_,X,_,_,_,X,_,_},
  {_,_,_,X,_,X,_,_,_},
  {_,X,_,_,X,_,_,X,_},
  {_,_,X,X,_,X,X,_,_},
  {_,_,_,_,_,_,_,_,_}};
//◇市松模様
static uint8_t image3[9][9]={
  {X,_,X,_,X,_,X,_,X},
  {_,X,_,X,_,X,_,X,_},
  {X,_,X,_,X,_,X,_,X},
  {_,X,_,X,_,X,_,X,_},
  {X,_,X,_,X,_,X,_,X},
  {_,X,_,X,_,X,_,X,_},
  {X,_,X,_,X,_,X,_,X},
  {_,X,_,X,_,X,_,X,_},
  {X,_,X,_,X,_,X,_,X}};
#undef _
#undef X
static uint16_t label[9][9];
static uint16_t parents[(9*9)/*調整可*/];//Level.1と2は[6]で足りる｡市松模様が最悪ｹｰｽで(w*h)必要｡一般的には(max(予想されるﾗﾍﾞﾙ数*2,32))程度で良いと思う｡
static void test() {
  int regions,x,y;
  //◇Level.1　⇒　ok
  regions=image_labeling(&image1[0][0],9,9,&label[0][0],1/*8近傍*/,parents,ARRAY_SIZE(parents));
  printf("regions=%d\n",regions);
  for(y=0;y<9;y++){for(x=0;x<9;x++){printf("%3d",label[y][x]);}putchar('\n');}
  //◇Level.2　⇒　ok
  regions=image_labeling(&image2[0][0],9,9,&label[0][0],1/*8近傍*/,parents,ARRAY_SIZE(parents));
  printf("regions=%d\n",regions);
  for(y=0;y<9;y++){for(x=0;x<9;x++){printf("%3d",label[y][x]);}putchar('\n');}
  //◇市松模様　⇒　ok
  regions=image_labeling(&image3[0][0],9,9,&label[0][0],0/*4近傍*/,parents,ARRAY_SIZE(parents));
  printf("regions=%d\n",regions);
  for(y=0;y<9;y++){for(x=0;x<9;x++){printf("%3d",label[y][x]);}putchar('\n');}
}
#endif
//}}使用例
