/* -------------------------------
1999年05月30日 13時35分10秒
最大値の途中経過を表示

1999年05月29日 11時24分20秒
対応表を３種類出力する

ay1.cls: Average最大
ay2.cls: Overall最大
ay3.cls: Kappa最大
（ay4.cls: 特例措置）

comp6.c
1999年05月25日
Type指定付加
0=Average, 1=Overall, 2=Kappa
Typeごとの判断基準で、最も高い分類精度の項目をクラスタ順に
調査して割り当てる


植生項目自動割当て
Ｏｖｅｒａｌｌ分類精度評価による最大値の対応表を使用

 ---------------------------------

植生図と分類結果の画素ごとの対応を調査
草地の割り当てを２位順で行った場合ＤＮ値判定をしない
変数名の整理した
1999年04月20日
1999年04月27日 10時44分02秒
広田湾付近の調査で、針葉樹と広葉樹の入れ替えでKHATが向上した
念のため、僅少差であった項目のメッセージを出力する
1位と2位の優先順位で
どちらも３割以上の場合、とりあえず１位を選ぶが、その事実があったことを表示

---------------------------------- */
#define Max_array=10000000  /* 植生図データ読み込み配列サイズ限界 */
#define Max_cluster=100     /* 教師無し分類のクラスタ最大数 */
#define Max_vegetation=100  /* 植生図項目の最大数 */
#define Max_band=6          /* 使用する衛星画像の最大バンド数 */

#include <stdio.h>
#include <math.h>
main(argv,argc)
int argv;
char *argc[];

{
	static int  ic,ii,total_data,dmax;
	static int  ix,iy,pnum,i,m,mm,n,j,k;
	static int  vegetation_max,cluster_max,flag;
	static int  pixel,line,veg_freq[10];
	static float f_x,f_y,f_z;
	/* 分類精度評価値 a_x;average o_x;overall k_x;kappa */
	static float a_x,o_x,k_x;
	/* コードブックの微分係数　１次、２次 */
	static int  deriv[Max_cluster][Max_band];
	static int  deriv1[Max_cluster][Max_band];
	static int  deriv2[Max_cluster][Max_band];
	/* コードブックの文字読みとり用 */
	static char f1[128]="",f2[32]="",f3[32]="";
	/* coinは一致数行列、priortyは優先順位割当表 rsは、画素数および割当番号 */
	static int  coin[Max_cluster][Max_vegetation];
	static int  priorty[Max_cluster][Max_vegetation];
	static int  priortys[Max_cluster][Max_vegetation]; /* save用 */
	/* 割当て結果の分類精度を格納 */
	static float  max_a,max_o,max_k;
	static float  ma,mo,mk;
	/* cls_aはAverage割当表、oはOvaerall,kはKappa格納用 */
	static int  next_p_num;
	/* 0=Average, 1=Ovaerall, 2=Kappa */
	static int  type; /* 0=Average, 1=Ovaerall, 2=Kappa */
	static int  cls_a[Max_cluster],cls_o[Max_cluster],cls_k[Max_cluster];
	/* clsは割当表、assign1は対角要素格納用、assgin2は行方向要素格納用 */
	static int  assign1[Max_cluster];
	static int  assign2[Max_cluster];
	static int  rs[Max_cluster][2];
	/* 植生項目ごとの総画素数 */
	static int  vg[Max_cluster];
	/* 植生データ、分類クラスタデータを一時格納する配列 */
	static int  vegetation[Max_array];
	static int  cluster[Max_array];
	static char fn[1024]="";
	FILE *fp1,*fp2,*fp3,*fp4;

	/* argument parameter read  */
	if (argv!=6){
		printf("--- Automatic Category Assignment Ver.7 ------\n");
		printf("Usage : comp7 PIXEL LINE veg_image class_image flag\n");
		printf("type=0 Average Type=1 Overall Type=2 Kappa \n");
		printf("flag=0 ; Output to 1 column  flag=1 ; table output\n");
		exit(1);
	}

	if (( fp1=fopen(argc[3],"r") )==NULL) {
		printf("Can't open Vegetation Inage\n");
		exit(1);
	}

	if (( fp2=fopen(argc[4],"r") )==NULL) {
		printf("Can't open Classification Image\n");
		exit(1);
	}

	sscanf(argc[1],"%d",&pixel);
	sscanf(argc[2],"%d",&line);
	sscanf(argc[5],"%d",&flag);

	total_data = pixel*line;

	/* Get Image */
	vegetation_max=0;    /* vegetation_max 植生図における最大数 */
	cluster_max=0;    /* cluster_max 分類結果における最大数 */
	for (i=0; i<total_data; i++) {
		m=getc(fp1);
		vegetation[i]=m;
		veg_freq[m]++;    /* 植生図のデータ数チェック */
		vegetation_max=(m>vegetation_max)? m:vegetation_max;
		m=getc(fp2);
		cluster[i]=m;
		cluster_max=(m>cluster_max)? m:cluster_max;
	}
	/* 全座標での対応値を調査 */
	for (i=0; i<total_data; i++) {
		/* coin[a][b]にはクラスタ分析結果と植生図データが合致した総数が入る */
		coin[cluster[i]][vegetation[i]]+=1;
		/* vg[a]:植生図項目ごとの総画素数 */
		vg[vegetation[i]]+=1;
	}

	/* 各クラスタ値と植生図の一致数の表示とソート */
	/* ここからは表示出力のための文 */
	if ( flag !=0 ) { 
		printf("Cluster |"); 
	}
	for (i=0; i<=vegetation_max; i++) {
		if ( flag !=0 ) { 
			printf("  Veg%1d ",i); 
		}
	}
	if ( flag !=0 ) {
		printf("  | Total | Assign Order\n");
		for (j=0; j<90; j++){
			printf("-"); 
		}
		printf("\n");
	}
	/* 表示出力文の終了 */

	/* クラスタ値と植生図の一致数の調査と優先付け */
	/* 対応表の作成と優先順位表の作成 */
	for (i=0; i<=cluster_max; i++)
	{
		if ( flag !=0 ) { 
			printf("%2d      | ",i); 
		}
		ii=0;
		for (j=0; j<=vegetation_max; j++)
		{
			m=coin[i][j];
			/* カテゴリー対応表deriv[a][0]は、植生番号が入る */
			/* カテゴリー対応表deriv[a][1]は、植生図とクラスタ図の一致画素数が入る */
			deriv[j][0]=j;
			deriv[j][1]=m;
			if ( flag !=0 ) { 
				printf("%5d  ",m); 
			}
			ii+=m;
		}
		if ( flag !=0 ) { 
			printf(" | %5d |",ii); 
		}
		/* カテゴリー対応表deriv[a][b]のソート deriv[0][0]が最も優先度が高い */
		for (k=0; k<vegetation_max; k++)
		{
			for (j=0; j<vegetation_max; j++)
			{
				if( deriv[j][1]<deriv[j+1][1] )
				{
					m=deriv[j][0];
					deriv[j][0]=deriv[j+1][0];
					deriv[j+1][0]=m;
					m=deriv[j][1];
					deriv[j][1]=deriv[j+1][1];
					deriv[j+1][1]=m;
				}
			}
		}
		/* rs[i][0]には、クラスタ番号ごとの画素総数が入っている */
		/* rs[i][1]には、最も大きい植生凡例番号（カテゴリー対応番号）が入っている */
		/* もしflag=1なら、カテゴリー対応優先順の表が表示される */
		rs[i][0]=ii;
		rs[i][1]=deriv[0][0];

		/* 優先順位表へコピー */
		for (j=0; j<=vegetation_max; j++)
		{
			priorty[i][j]=deriv[j][0];
			priortys[i][j]=deriv[j][0];
		}

		/* ここからは表示出力のための文 */
		if ( flag !=0 ) {
			for (j=0; j<=vegetation_max; j++)
			{
				printf("%2d ",deriv[j][0]);
			}
			printf("\n");
		}

	}


	/*  -------------------------------
　カテゴリー対応

　当面、特例措置は行わない

---------------------------------- */
	/* 最適割当のための繰り返しの開始 */
	for (type=0; type<=2; type++)
	{
        max_a=0;
	max_o=0;
	max_k=0;

	if ( flag !=0 ) {
		printf("-------- Adjust type = %d ---------\n",type);
	}

	for (iy=0; iy<=cluster_max; iy++)	
	{
        ma=0;
	mo=0;
	mk=0;

	for (ix=0; ix<=4; ix++)	
	{
	priorty[iy][0]=priortys[iy][ix];

	/* 対角要素の決定 */
		/* 対角要素行列クリア */
		for (i=0; i<=vegetation_max; i++)
		{
		assign1[i]=0;
		}
		/* 優先表に基づき要素の加算格納 */
		for (i=0; i<=cluster_max; i++)
		{
		m=priorty[i][0]; /* 選択要素位置 */
		assign1[m]+=coin[i][m];
		}

	/* 行方向要素数の決定 */
		/* 行列要素行列クリア */
		for (i=0; i<=vegetation_max; i++)
		{
		assign2[i]=0;
		}
		/* 優先表に基づき要素の加算格納 */
		for (i=0; i<=cluster_max; i++)
		{
		m=priorty[i][0]; /* 選択要素位置 */
		assign2[m]+=rs[i][0]; /* クラスタごとの総画素数を加算 */
		}

	/*
	if ( flag !=0 ) {
		printf("---------- Chack Message from comp6 Assign 2 ------------\n");
		for (i=0; i<=vegetation_max; i++)
		{
		printf("%d ",assign2[i]);
		}
		printf("\n");
	}
	*/

	/* Average誤差 a_x の算出 */
	a_x=0;
	 for (i=0; i<=vegetation_max; i++)
	 {
		 if (vg[i] != 0 )
		 {
		 a_x += (float)assign1[i]/(float)vg[i];
		}
	 }
	 a_x= a_x/(float)(1+vegetation_max);

	/* Overall誤差 b_x の算出 */
	o_x=0;
	 for (i=0; i<=vegetation_max; i++)
	 {
		 if (vg[i] != 0 )
		 {
		 o_x += (float)assign1[i]/(float)total_data;
		}
	 }
	 /*
	 if ( flag !=0 ) {
		 printf("Overall = %f\n",o_x);
	 }
	*/

	/* Kappa誤差 c_x の算出 */
	f_x=0;
	f_y=0;
	 for (i=0; i<=vegetation_max; i++)
	 {
		f_x += (float)assign1[i];
		f_y += (float)assign2[i]*(float)vg[i];
	 }
	 k_x = ((float)total_data*f_x - f_y)/((float)total_data*(float)total_data - f_y);

	/* 同クラスタ内での最大値判断 */
	if ( type == 0 ){
	if (ma < a_x) {
		ma=a_x;
		next_p_num= priorty[iy][0];
	}
	}

	if ( type == 1 ){
	if (mo < o_x) {
		mo=o_x;
		next_p_num= priorty[iy][0];
	}
	}

	if ( type == 2 ){
	if (mk < k_x) {
		mk=k_x;
		next_p_num= priorty[iy][0];
	}
	}


	/* 精度の最大値判断 */
	if (max_a < a_x) {
		max_a=a_x;
		for (i=0; i<=cluster_max; i++)
		{
		cls_a[i]=priorty[i][0];
		}
	}

	if (max_o < o_x) {
		max_o=o_x;
		for (i=0; i<=cluster_max; i++)
		{
		cls_o[i]=priorty[i][0];
		}
	}

	if (max_k < k_x) {
		max_k=k_x;
		for (i=0; i<=cluster_max; i++)
		{
		cls_k[i]=priorty[i][0];
		}
	}


	/* 最適割当ての終了検出 */

	 if ( flag !=0 )
	 {
	  for (i=0; i<=cluster_max; i++)
	  {
	  printf("%d ",priorty[i][0]);
	  }
 	  printf(" |  %f %f %f max=%d\n", a_x,o_x,k_x,next_p_num);
	 }


	} /* iy end */

	priorty[iy][0]=next_p_num;

	} /* ix end */



	 if ( flag !=0 ) {
 	  printf("Max_average = %f , overall = %f , kappa = %f\n", max_a,max_o,max_k);
	 }

	/* 各typeでの結果の表示  Reclass File 作成用table */
        if ( flag == 0 )
        {

                if (type == 0) {
		if (( fp3=fopen("ay0.cls","w") )==NULL) {
			printf("Can't open ay0.cls file\n");
			exit(1); }
                fprintf(fp3,"* new number of class : %d\n",vegetation_max+1);
		for (i=0; i<=cluster_max; i++) {
                        fprintf(fp3,"%2d %d\n",i,cls_a[i]); }
		fclose(fp3);
                }

                if (type == 1) {
		if (( fp3=fopen("ay1.cls","w") )==NULL) {
			printf("Can't open ay1.cls file\n");
			exit(1); }
                fprintf(fp3,"* new number of class : %d\n",vegetation_max+1);
		for (i=0; i<=cluster_max; i++) {
                        fprintf(fp3,"%2d %d\n",i,cls_o[i]); }
		fclose(fp3);
                }

                if (type == 2) {
		if (( fp3=fopen("ay2.cls","w") )==NULL) {
			printf("Can't open ay2.cls file\n");
			exit(1); }
                fprintf(fp3,"* new number of class : %d\n",vegetation_max+1);
		for (i=0; i<=cluster_max; i++) {
                        fprintf(fp3,"%2d %d\n",i,cls_k[i]); }
		fclose(fp3);
                }

        }

	} /* type end */

	fclose(fp1);
	fclose(fp2);

	exit(0);
}