/*	
 *	clipimg.c
 *
 *	ユーティリティ：画像処理
 *
 *	CLiP - Common Library for P/ECE
 *	Copyright (C) 2001-2016 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ビット整数処理⇒実数処理に変更しました。
 *	* Sat Dec 17 21:25:03 JST 2016 Naoyuki Sawa
 *	- k_means_clustering_init(),k_means_clustering_init_r(),k_means_clustering_step(),k_means_clustering_step_r()を追加しました。
 *	* Sun Dec 25 21:14:35 JST 2016 Naoyuki Sawa
 *	- 動的メモリ割り当てを使って、速度を改善したバージョンの、k_means_clustering_init_r()を実装しました。
 */
#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
//}}使用例
/*****************************************************************************
 *	k-means法,及び,k-means++法によるクラスタリング処理
 *****************************************************************************/
#undef  K_MEANS_USE_DYNAMIC_MEMORY_ALLOCATION
#define K_MEANS_USE_DYNAMIC_MEMORY_ALLOCATION	//このシンボルを定義すると、動的メモリ割り当てを使って、速度を改善したバージョンを使用します。このシンボルを定義しなければ、動的メモリ割り当てを行わない、低速版を使用します。どちらを使用しても、結果は同じになります。
/*--------------------------------------------------------------------------*/
//k-means++法の初期化処理を実行する。
//[in]
//	vector		ベクトルの配列						(in)
//	_n		ベクトルの数						(in)
//	center		クラスタ中心の配列					(unused)	この引数は、分かり易さのために、k_means_clustering_init()とk_means_clustering_step()の引数を揃えるための、ダミー引数です。当関数はこの引数を参照しないので、実際の引数の値は何でも構いません。
//	_k		クラスタ中心の数					(in)
//	label		ラベルの配列						(out)		ラベルとは、各ベクトルが属するクラスタのインデクス(0～(k-1))の事です。
//	width		ベクトル,及び,クラスタ中心の配列要素のバイト数		(in)
//	distance	ベクトルとクラスタ中心の距離を求める関数
//	arg		アプリケーション定義の任意の引数
//[out]
//	戻り値		引数_nをそのまま返します。
//[note]
//	- □参考資料
//	  │「Wikipedia - k-means++法」(http://ja.wikipedia.org/wiki/k-means++法)
//	- 今回は、動的メモリ割り当てを行わない事を優先して実装したので、実行速度はかなり遅いです。
//	  今後、動的メモリ割り当てを使って、速度を改善したバージョンも実装する予定です。	⇒{{2016/12/25コメント追記:動的メモリ割り当てを使って、速度を改善したバージョンの、k_means_clustering_init_r()を実装しました。}}
int k_means_clustering_init(  const void* vector/*[width*n]*/, size_t _n, const void* center/*[width*k]*/, size_t _k, int label[/*n*/], size_t width, double (*distance)(const void* vector, const void* center           )           ) {
	return k_means_clustering_init_r(vector, _n, center, _k, label, width, (double (*)(const void*, const void*, void*))distance, NULL);	//k_means_clustering_init()とk_means_clustering_init_r()は、distance関数の形式が異なりますが、'cdecl呼出規約'ならば、k_means_clustering_init()形式のdistance関数において、余分な引数が無視されるので問題ありません。
}
#ifndef K_MEANS_USE_DYNAMIC_MEMORY_ALLOCATION
//□低速版(動的メモリ割り当てを使用しない)
int k_means_clustering_init_r(const void* vector/*[width*n]*/, size_t _n, const void* center/*[width*k]*/, size_t _k, int label[/*n*/], size_t width, double (*distance)(const void* vector, const void* center, void* arg), void* arg) {
	int n = _n, k = _k;
	int i, l;
	//ベクトル数が0以下ならば何もしない。そうでなければ…
	if(n > 0) {
		//クラスタ数をベクトル数以下に切り詰める。
		if(k > n) { k = n; }
		//無条件で、最初のベクトルを、最初のクラスタ(0)中心として選択する。
		label[0] = 0;
		//二番目以降の全てのベクトルに、ラベル(k-1)を設定しておく。
		// - ラベル(k-1)は、このベクトルがまだ、クラスタ中心として選択されていない事を示す。
		//   尚、最後まで選択されなかったベクトルは、全てクラスタ(k-1)に属すると見なすので、ここで設定したラベル(k-1)そのままで良い。
		for(i = 1; i < n; i++) { label[i] = (k - 1); }
		//二番目のクラスタ(1)から、最後の一つ前のクラスタ(k-2)までの、クラスタ中心となるベクトルを選択する。
		for(l = 1; l < (k - 1); l++) {
			//最近傍のクラスタ中心までの距離が最大であるベクトルを検索する。
			// - 「Wikipedia - k-means++法」(http://ja.wikipedia.org/wiki/k-means++法)では、最近傍のクラスタ中心までの距離を重みと見なし、重みつき確率分布を用いてランダムで選ぶ事になっている。
			//   しかし、最近傍のクラスタ中心までの距離が最大であるベクトルを、単純に選択しても良いと思うので、そうする事にした。
			double d_nearest_max = 0.0;	//警告抑制
			int    j_nearest_max =  -1;
			for(i = 1; i < n; i++) {
				//このベクトルがまだ、クラスタ中心として選択されていなければ…
				if(label[i] == (k - 1)) {
					const void* v = (const char*)vector + (width * i);
					//最近傍のクラスタ中心までの距離を求める。
					double d_nearest = 0.0;	//警告抑制
					int j, j_nearest =  -1;
					for(j = 0; j < n; j++) {
						//このベクトルが、クラスタ中心として選択されていたら…
						if(label[j] != (k - 1)) {
							const void* c = (const char*)vector + (width * j);
							double d = (*distance)(v, c, arg);
							if(d < 0.0) { DIE(); }
							//初回,又は,これまでで最小の距離ならば…
							if((j_nearest == -1) || (d < d_nearest)) {
								d_nearest = d;
								j_nearest = j;
							}
						}
					}
					if(j_nearest == -1) { DIE(); }	//バグ
					//初回,又は,これまでで最大の距離ならば…
					if((j_nearest_max == -1) || (d_nearest > d_nearest_max)) {
						d_nearest_max = d_nearest;
						j_nearest_max = i;	//'j_nearest_max = j_nearest;'としないよう注意せよ。
					}
				}
			}
			if(j_nearest_max == -1) { DIE(); }	//バグ
			//最近傍のクラスタ中心までの距離が最大であるベクトルを、l番目のクラスタ中心に選択する。
			label[j_nearest_max] = l;
		}
	}
	//全てのベクトルのクラスタを付け替えたと見なして、ベクトル数を返す。
	return n;
}
#else //K_MEANS_USE_DYNAMIC_MEMORY_ALLOCATION
//□高速版(動的メモリ割り当てを使用する)
int k_means_clustering_init_r(const void* vector/*[width*n]*/, size_t _n, const void* center/*[width*k]*/, size_t _k, int label[/*n*/], size_t width, double (*distance)(const void* vector, const void* center, void* arg), void* arg) {
	int n = _n, k = _k;
	//ベクトル数が0以下ならば何もしない。そうでなければ…
	if(n > 0) {
		double d, d_nearest_max, *TBL_d_nearest/*[n]*/;
		int i, l, j_nearest_max;
		const void *v, *c;
		//最近傍のクラスタ中心までの距離を格納する配列のメモリを確保する。
		TBL_d_nearest/*[n]*/ = malloc(sizeof(double) * n);
		if(!TBL_d_nearest) { DIE(); }
		//クラスタ数をベクトル数以下に切り詰める。
		if(k > n) { k = n; }
		//無条件で、最初のベクトルを、最初のクラスタ(0)中心として選択する。
		      j_nearest_max  = 0;
		label[j_nearest_max] = 0;
		//二番目以降の全てのベクトルに、ラベル(k-1)を設定しておく。
		// - ラベル(k-1)は、このベクトルがまだ、クラスタ中心として選択されていない事を示す。
		//   尚、最後まで選択されなかったベクトルは、全てクラスタ(k-1)に属すると見なすので、ここで設定したラベル(k-1)そのままで良い。
		for(i = 1; i < n; i++) { label[i] = (k - 1); }
		//二番目のクラスタ(1)から、最後の一つ前のクラスタ(k-2)までの、クラスタ中心となるベクトルを選択する。
		for(l = 1; l < (k - 1); l++) {
			c = (const char*)vector + (width * j_nearest_max);	//直前に選択されたクラスタ中心
			//二番目以降の全てのベクトルに対して…
			for(i = 1; i < n; i++) {
				//このベクトルがまだ、クラスタ中心として選択されていなければ…
				if(label[i] == (k - 1)) {
					//直前に選択されたクラスタ中心までの距離を求め、最近傍のクラスタ中心までの距離を更新する。
					v = (const char*)vector + (width * i);
					d = (*distance)(v, c, arg);
					if(d < 0.0) { DIE(); }	//バグ
					if((l == 1) || (d < TBL_d_nearest[i])) { TBL_d_nearest[i] = d; }	//初回,又は,これまでで最小の距離ならば…
				}
			}
			//最近傍のクラスタ中心までの距離が最大であるベクトルを検索する。
			// - 「Wikipedia - k-means++法」(http://ja.wikipedia.org/wiki/k-means++法)では、最近傍のクラスタ中心までの距離を重みと見なし、重みつき確率分布を用いてランダムで選ぶ事になっている。
			//   しかし、最近傍のクラスタ中心までの距離が最大であるベクトルを、単純に選択しても良いと思うので、そうする事にした。
			d_nearest_max = 0.0;	//警告抑制
			j_nearest_max =  -1;
			for(i = 1; i < n; i++) {
				//このベクトルがまだ、クラスタ中心として選択されていなければ…
				if(label[i] == (k - 1)) {
					//最近傍のクラスタ中心までの距離を更新する。
					d = TBL_d_nearest[i];
					if(d < 0.0) { DIE(); }	//バグ
					if((j_nearest_max == -1) || (d > d_nearest_max)) {			//初回,又は,これまでで最大の距離ならば…
						d_nearest_max = d;
						j_nearest_max = i;
					}
				}
			}
			if(j_nearest_max == -1) { DIE(); }	//バグ
			//最近傍のクラスタ中心までの距離が最大であるベクトルを、l番目のクラスタ中心に選択する。
			label[j_nearest_max] = l;
		}
		//最近傍のクラスタ中心までの距離を格納する配列のメモリを解放する。
		free(TBL_d_nearest);
	}
	//全てのベクトルのクラスタを付け替えたと見なして、ベクトル数を返す。
	return n;
}
#endif//K_MEANS_USE_DYNAMIC_MEMORY_ALLOCATION
/*--------------------------------------------------------------------------*/
//k-means法(,又は,k-means++法)の繰り返し処理を1ステップ実行する。
//[in]
//	vector		ベクトルの配列						(in)
//	_n		ベクトルの数						(in)
//	center		クラスタ中心の配列					(in)
//	_k		クラスタ中心の数					(in)
//	label		ラベルの配列						(in,out)	ラベルとは、各ベクトルが属するクラスタのインデクス(0～(k-1))の事です。
//	width		ベクトル,及び,クラスタ中心の配列要素のバイト数		(in)
//	distance	ベクトルとクラスタ中心の距離を求める関数
//	arg		アプリケーション定義の任意の引数
//[out]
//	戻り値		最も近いクラスタの付け替えが発生した回数
//[note]
//	- □参考資料
//	  │「Wikipedia - k平均法」(http://ja.wikipedia.org/wiki/k平均法)
//	- 今回は、動的メモリ割り当てを行わない事を優先して実装したので、実行速度はかなり遅いです。
//	  今後、動的メモリ割り当てを使って、速度を改善したバージョンも実装する予定です。	⇒{{2016/12/25コメント追記:k_means_clustering_step_r()は、元々そんなに遅くはないので、現在のままでも構わない気がして来ました。}}
int k_means_clustering_step(  const void* vector/*[width*n]*/, size_t _n, const void* center/*[width*k]*/, size_t _k, int label[/*n*/], size_t width, double (*distance)(const void* vector, const void* center           )           ) {
	return k_means_clustering_step_r(vector, _n, center, _k, label, width, (double (*)(const void*, const void*, void*))distance, NULL);	//k_means_clustering_step()とk_means_clustering_step_r()は、distance関数の形式が異なりますが、'cdecl呼出規約'ならば、k_means_clustering_step()形式のdistance関数において、余分な引数が無視されるので問題ありません。
}
int k_means_clustering_step_r(const void* vector/*[width*n]*/, size_t _n, const void* center/*[width*k]*/, size_t _k, int label[/*n*/], size_t width, double (*distance)(const void* vector, const void* center, void* arg), void* arg) {
	int n = _n, k = _k;
	int i, relabel = 0;
	//各ベクトルについて…
	for(i = 0; i < n; i++) {
		const void* v = (const char*)vector + (width * i);
		double d_nearest = 0.0;	//警告抑制
		int j, j_nearest =  -1;
		//各クラスタ中心について…
		for(j = 0; j < k; j++) {
			const void* c = (const char*)center + (width * j);
			//ベクトルとクラスタ中心の距離を求める。
			double d = (*distance)(v, c, arg);
			if(d < 0.0) { DIE(); }
			//初回,又は,これまでで最小の距離ならば…
			if((j_nearest == -1) || (d < d_nearest)) {
				d_nearest = d;
				j_nearest = j;
			}
		}
		if(j_nearest == -1) { DIE(); }
		//最も近いクラスタの付け替えが発生したら…
		if(j_nearest != label[i]) {
			label[i] = j_nearest;
			relabel++;
		}
	}
	//最も近いクラスタの付け替えが発生した回数を返す。
	return relabel;
}
/*--------------------------------------------------------------------------*/
//{{使用例
#if 0
static void test(int n, int k);  //前方参照
int app_main(int argc, char* argv[]) {
  test(1000/*調整可*/, 50/*調整可*/);
  return 0;
}
static double distance(const void* vector, const void* center) {
  double v = *(double*)vector;
  double c = *(double*)center;
  return (v - c) * (v - c);
}
static void test(int n, int k) {
  double* vector = malloc(sizeof(double) * n);
  double* center = malloc(sizeof(double) * k);
  int* label = malloc(sizeof(int) * n);
  int i, j, try;
  //ベクトルを作成する。
  for(i = 0; i < n; i++) {
    vector[i] = drand48();
  }
  //k-means++法の初期化処理を実行する。
  k_means_clustering_init(vector, n, center/*ダミー*/, k, label, sizeof vector[0], distance);
  //k-means法(,又は,k-means++法)の繰り返し処理を実行する。
  for(try = 0; try < 100/*調整可*/; try++) {  //極端に収束が悪い場合のために、繰り返し回数を制限する。
    //クラスタ中心を求める。
    for(j = 0; j < k; j++) {
      double sum = 0.0;
      int    cnt = 0;
      for(i = 0; i < n; i++) {
        if(label[i] == j) {
          sum += vector[i];
          cnt++;
        }
      }
      if(cnt) { sum /= cnt; }
      center[j] = sum;
    }
    //k-means法(,又は,k-means++法)の繰り返し処理を1ステップ実行する。
    //最も近いクラスタの付け替えが発生しなければ(又は少なければ)、終了する。
    if(k_means_clustering_step(vector, n, center, k, label, sizeof vector[0], distance) <= 1/*調整可*/) { break; }
  }
  //結果を表示する。
  for(i = 0; i < n; i++) {
    printf("%f\t%d\n", vector[i], label[i]);
  }
  free(vector);
  free(center);
  free(label);
}
#endif
//}}使用例
