/*   
**  
**  
**		PROGRAMME :som_lvq.c       
**		ORIGINAL AUTHOR: Alessandra Chiuderi (J.R.C. Italy)       
**		MODIFIED BY: Brandt Tso
**
**		This programme uses the Kohonen Self Organizing Map     
**          (Euclidean distance) and the LVQ1 algorithm 
**           to achieve image classification (unsupervised & supervised)		   
**		The program may generate 4 different results:       
**                                1) Clusters after Kohonen's training          
**				2) Clusters after LVQ training      
**				3) Classification (clusters are labelled)       
**				4) Confusion matrix       
**	
*/ 
 

#include <string.h> 
#include <stddef.h> 
#include <ctype.h>  
#include <stdlib.h> 
#include <stdio.h> 
#include <math.h> 
#include <time.h> 
 
#define MAX 200000000 
#define BIG 300000000 
#define TRUE 1 
#define FALSE 0 
 
 
FILE *fp_r ;  /* Image file pointer */ 
FILE *fp_w_k ;  /* Kohonen's map clustered image file pointer*/ 
FILE *fp_w_g ;  /*LVQ clustered image file pointer*/ 
FILE *fp_w_cl;  /* Classified image file pointer*/ 
FILE *fp_c;   /* pixel counts file pointer  */ 
 
/*cluster mean and cluster unit number pointer*/
FILE *fp_cls_mean;
FILE *fp_cls_num;


/*output cluster sequrntial file for each band*/
FILE *f1, *f2, *f3, *f4, *f5, *f6, *f7;



/* Training data file pointers */ 
FILE *fg ;   
FILE *fpt, *fpcx; 
 
unsigned char *buf_i ; /* Array containing 1 image line */ 
unsigned char *y ; /* array containing 1 output clustered image line */ 
unsigned char *g ; /* array containing 1 output classified image line */ 
 
double *x ; /* input vector for the SOM */ 
double **w ; /* weights matrix  (SOM)*/ 
double **v; /* weights matrix (LVQ) */ 

double **cls_mean; /*clsters' mean*/

double   *bb;  
int *label, correct =0, incorrect=0; 
int r,num, dimnk ,l_righe, s_time; 
long ltime; 


char store[10],s0[2],s1[2],s2[2],s3[2],s4[2],s5[2],s6[2],s7[2]; 
char testname[50]; 
char cname[50];  
char koutname[50]; 
char goutname[50]; 
char clname[80]; 
char nomefile[80] ; 
char gname[70]; /* training set file name*/ 
char savename[50]; /*  file containing network's weights */ 
char fname[50]  ; /* File containing the the following network's fundamental  
		     parameters: */ 
 
int lines;   /* number of image lines   */ 
int pixels;  /* number of image  pixel  */ 
int nbande;  /* number of input nodes*/ 
int ncols; /* number of units per row in the SOM */ 
int nrows; /* number of rows in the SOM */ 
int nuscite; /* number of classes */ 
 
int *cl_kcount;  /* Number of pixel per cluster (SOM) */ 
int *cl_gcount;  /* Numero of pixel per cluster (LVQ) */ 
int *cl_cl_count;  /* Number of pixel per class */ 
 
int nkiter=0 ;  /* Number of iter durineg previous training sessions */  
int ngiter =0; 
int np_kiter;  
int np_giter; 
 
double alfa;         /* Learning parameters */  
double gamma1 ; 
int nb_r=0; 

int num_train_cell = 0; /*number of training blocks in training set (for LVQ)*/ 
/*Header I_param, *O_param;*/ 
  
/* Procedures employed */  
 
int BEST_UNIT_CHOICE () ; 
void WEIGHTS_COMPUTE () ; 
void K_WEIGHTS_INITIALIZE () ; 
void MEMORY_ALLOCATION () ; 
void FREE_ALL(); 
void TRAIN_NET1 () ; 
void SAVE_TRAINING () ; 
void READ_TRAINING () ; 
void READ_IMAGE(); 
void PARAM_COMPUTE () ; 
void START_AGAIN(); 
int  PIXEL_COUNT (); 
double GETRANDOM(); 
int ROOF(); 
void SHOW_MAP(); 
void LABEL_UNITS(); 
double DISTANCE(); 
void LVQ(); 
void LVQ_WEIGHTS(); 
void CONFUSION_MATRIX(); 
 

main() 
{ 
  int p, q,ntot ; 
  int i,j, count1; 
  int choice,goon ; 
  int banda, dif ; 
  /***********************/
  int BestUnit ; 
  double  dist, BestDist ; 
  BestUnit = 0 ; 
  BestDist = BIG; 
  /***************************/
 
 
  printf ("                   - COUNTERPROPAGATION NEURAL NETWORK -\n") ; 
  printf ("                   *************************************\n\n\n") ; 
  printf ("\n") ; 
  printf ("    1) Start a new training session (Kohonen layer)\n") ; 
  printf ("    2) Improve previous training (Kohonen layer)\n"); 
  printf ("    3) Run the supervised LVQ1 algorithm for the first time\n") ; 
  printf ("    4) Improve previous LVQ1 training\n"); 
  printf ("    5) Test the network (Kohonen layer)\n") ; 
  printf ("    6) Test the network (Supervised layer)\n") ; 
  printf ("Enter your choice :  ") ; 
  scanf ("%d",&choice) ; 
 
  ltime=time(NULL); 
  s_time = (unsigned) ltime/100; 
 
  switch (choice) 
		  { 
		    case 1 :  
	       		    TRAIN_NET1 (1) ; 
					fclose(fp_r);
		     break ; 
 
		    case 2: 
					printf("- Enter filename of previous training session: "); 
					scanf("%s",fname); 
					READ_TRAINING (1); 
					TRAIN_NET1(2); 
					fclose(fp_r);
        	 break; 
 
		    case 3:   
                     printf("- Enter filename containing Kohonen's layer weights: "); 
                     scanf("%s", fname); 
                     READ_TRAINING(2); 
                     TRAIN_NET1(3); 
		     break; 
 
		    case 4: 
                     printf("- Enter filename of previous training session: "); 
                     scanf("%s", fname); 
                     READ_TRAINING(3); 
                     TRAIN_NET1(4); 
		     break; 
 
		    case 5:   
					 printf("- Enter filename of trained network : "); 
					scanf("%s",fname); 
					READ_TRAINING(1); 
					goon=0; 
					goto l2; 
		     break;
 
		    case 6:
                     printf("- Enter filename of trained network : ");
					 scanf("%s",fname);
                     READ_TRAINING(3);
		     break;

	    default:
		     printf("Bad choice !\n");
		     exit(1);
 
		  }  /* end switch  */
  
 
  printf ("\n\n                            START OF TEST\n") ; 
  printf ("\n\n") ; 
  printf ("      0) Kohonen map\n") ; 
  printf ("      1) Both layers(Kohonen and LVQ no labelling)\n") ; 
  printf ("      2) Supervised layer only\n") ; 
  printf ("      3) Test set (confusion matrix generation)\n") ;
  printf ("      4) No testing\n"); 
  printf ("Enter your choice :  ") ;
  scanf ("%d",&goon) ; 
 
 
l2: if (goon!=4) 
    { 
      if ((goon!=3)) 
       {  
		READ_IMAGE(); 
        /*O_param = &I_param;            
		O_param->nbands=1;
		O_param->nclass=dimnk;*/ 
        for(i=0;i<dimnk;i++) 
         { 
          cl_kcount[i]=0; 
          cl_gcount[i]=0; 
         } 
        for (i=0;i<nuscite;i++) 
          cl_cl_count[i]=0; 
        if ((goon==0)||(goon==1)) 
          {
           printf ("\n- Enter Kohonen output file name: ") ; 
           scanf ("%s",koutname) ;
           if ((fp_w_k=fopen(koutname,"wb")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", koutname) ;
                exit (1) ;
               } 
           /* 
	   fp_w_k=ErdasOpenToWrite(koutname); 
	   VaxErdasWriteHeader(fp_w_k,O_param);
           */ 
			}  /* chiude if goon==0 ||goon==1  */ 
        if (goon == 1)
          {
            printf ("- Enter LVQ output file name: ") ;
            scanf ("%s",goutname) ;
            if ((fp_w_g=fopen(goutname,"wb")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", goutname) ;
                exit (1) ;
               } 
            /*
            fp_w_g=ErdasOpenToWrite(goutname);
            VaxErdasWriteHeader(fp_w_g,O_param);
	    */
           } /*chiude if (goon==1)  */
         if (goon==2)
           {
            printf("- Enter file name for classified image: ");
            scanf("%s",clname);
            printf("- Enter training file name for labelling units: ");
            scanf("%s",gname);
            if((fg=fopen(gname,"r"))==NULL)
              {
               printf("Error in opening training file %s\n",gname);
               exit(1);
              } 
            LABEL_UNITS(v);

 /*******************Added by Brandt**********************************************/
   

			for (i=0; i<dimnk; i++)
			{
				if(label[i]+1 == 0)
				{
				BestDist = BIG;
				BestUnit = 0;
				for (count1=0; count1<dimnk; count1++) 
					{
					if((label[count1]+1 != 0) && (i != count1))
						{
						dist = 0;
						for(j=0;j<nbande;j++) 
							{
							dist = dist+(v[i][j]-v[count1][j])*(v[i][j]-v[count1][j]);
							}
						if(dist < BestDist)
							{ 
							BestDist = dist;
							BestUnit = count1;
							}
						}
					}
				/*printf("%d\n",BestUnit);*/
				label[i] = label[BestUnit];
				} 
       
			}/*END OF i*/

			for(i=0;i<dimnk;i++) 
			printf("Label of cluster %d  %d\n",i+1,label[i]+1); 
   
/******************************************************************************/

            if ((fp_w_cl=fopen(clname,"wb")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", clname) ;
                exit (1) ;
               }
            /*fp_w_cl=ErdasOpenToWrite(clname);
            O_param->nclass=nuscite; 
            VaxErdasWriteHeader(fp_w_cl,O_param);*/
           } /* end if goon==2 */ 
        
		 dif = nbande -1; 

	

          for (i=0;i<lines;i++)
			{ 
               
       	       /*dasReadImageLine(fp_r,buf_i,l_righe);*/
               fread(buf_i,nbande*pixels,1,fp_r);                
                 
	       for (j=0; j<pixels; j++) 
	          { 
	            for (banda=0; banda<nbande; banda++) 
		           x[banda]=(double) buf_i[j+banda*pixels]; 

		        p = BEST_UNIT_CHOICE (w) ; 
                   if ((goon==1)||(goon==2)) 
                     { 
                      q = BEST_UNIT_CHOICE(v); 
                      if (goon==1) 
                        { 
                          g[j] = 1 + q; 
                          cl_gcount[q]++; 
                        } 
                      if (goon==2)
                        { 
                         g[j]=1+label[q]; 
                         cl_cl_count[label[q]]++; 
                        } 
                     } 
		   y[j] = p + 1 ; 
		   cl_kcount[p]++; 

      /*******************************************added by Brandt********************
                   for (banda=0; banda<nbande; banda++)
		     {
                     if(cls_mean[p][banda] == 0)
                       cls_mean[p][banda] = cls_mean[p][banda] + x[banda];
                     else
                       cls_mean[p][banda] = (double) (cls_mean[p][banda] + x[banda])/2;
		     }
       ************************************************************************/

	          } /*end for j  */ 

               if ((goon==0)||(goon==1)) 
	         /*ErdasWriteImageLine(fp_w_k,y,pixels);*/ 
                 fwrite(y,pixels,1,fp_w_k);
               
               if (goon == 1) 
                 /*ErdasWriteImageLine(fp_w_g,g,pixels);*/
                 fwrite(g,pixels,1,fp_w_g); 
               if (goon == 2) 
                 /*ErdasWriteImageLine(fp_w_cl,g,pixels);*/
                 fwrite(g,pixels,1,fp_w_cl);
 
	      } /* chiude for i  */
          
                if ((goon==0)||(goon==1)) 
	                fclose(fp_w_k);
                if (goon == 1) 
                    fclose(fp_w_g); 
                if (goon == 2) 
                    fclose(fp_w_cl); 
				fclose(fp_r);


      /*******************************************added by Brandt output cluster mean and nuber files****************
          if ((fp_cls_mean=fopen("cluster_mean.dat","w")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", "cluster_mean.dat") ;
                exit (1) ;
               }
 
          if ((fp_cls_num=fopen("cluster_num.dat","w")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", "cluster_num.dat") ;
                exit (1) ;
               }

          for (i=0;i<dimnk;i++) 
	    {
            for (banda=0; banda<nbande; banda++)
	      {
               if(banda == nbande - 1)
                fprintf(fp_cls_mean,"%f\n",cls_mean[i][banda]);
               else
                fprintf(fp_cls_mean,"%f,",cls_mean[i][banda]);
	      }
            fprintf(fp_cls_num,"%d\n",cl_kcount[i]);
	    }
       *********************************************************************************/     
       /*******************************************added by Brandt output cluster sequential file**************
          if ((f1=fopen("cluster_b1","w")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", "cluster_b1") ;
                exit (1) ;
               }
          if ((f2=fopen("cluster_b2","w")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", "cluster_b2") ;
                exit (1) ;
               }
          if ((f3=fopen("cluster_b3","w")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", "cluster_b3") ;
                exit (1) ;
               }
          if ((f4=fopen("cluster_b4","w")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", "cluster_b4") ;
                exit (1) ;
               }
          if ((f5=fopen("cluster_b5","w")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", "cluster_b5") ;
                exit (1) ;
               }
          if ((f6=fopen("cluster_b6","w")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", "cluster_b6") ;
                exit (1) ;
               }
	      

          
           fclose(fp_w_k);
           if ((fp_w_k=fopen(koutname,"r+")) == NULL)
               {
                printf ("\nERROR IN OPENING FILE %s\n", koutname) ;
                exit (1) ;
               } 
	    
	   for(count1=1;count1<=dimnk;count1++)
	    {
            fseek(fp_r,0,0);
            fseek(fp_w_k,0,0);
            for(i=0;i<lines;i++) 
	     { 
               fread(buf_i,nbande*pixels,1,fp_r);                
	       for (j=0; j<pixels; j++) 
	          { 
                   fread(store,1,1,fp_w_k);
		   if(store[0] == count1)
		     {
                      s0[0]=buf_i[j+0*pixels];s1[0]=buf_i[j+1*pixels];s2[0]=buf_i[j+2*pixels];s3[0]=buf_i[j+3*pixels];
                      s4[0]=buf_i[j+4*pixels];s5[0]=buf_i[j+5*pixels];
		      fwrite(s0,1,1,f1); 
                      fwrite(s1,1,1,f2); 
                      fwrite(s2,1,1,f3); 
                      fwrite(s3,1,1,f4); 
	              fwrite(s4,1,1,f5); 
                      fwrite(s5,1,1,f6); 
                      
		     }
		  }
	     }
	   }
       *********************************************************************************/  

          ntot=pixels*lines; 
          fp_c=fopen("cluster.dat","w"); printf("\nOPEN Cluster.dat FILE\n");
	  fprintf(fp_c,"Number of pixel per cluster unsupervised layer:\n\n"); 
          for (i=0;i<dimnk;i++) 
             fprintf(fp_c,"Cluster %2d: %6d pixels  %6.4f\n",i+1,cl_kcount[i],(double) cl_kcount[i]/ntot);
          if(goon==1) 
           { 
	    fprintf(fp_c,"Number of pixel per cluster after LVQ1 algorithm:\n\n"); 
            for (i=0;i<dimnk;i++) 
               fprintf(fp_c,"Cluster %2d: %6d pixels  %6.4f\n",i+1,cl_gcount[i],(double) cl_gcount[i]/ntot); 
           } 
          if (goon == 2) 
            { 
             fprintf(fp_c,"Number of pixel per class:\n\n"); 
             for (i=0;i<nuscite;i++)
               fprintf(fp_c,"Class %2d: %6d pixels  %6.4f\n",i+1,cl_cl_count[i],(double) cl_cl_count[i]/ntot); 
             } 
          fclose(fp_c); 
        } /* chiude if (goon!=3)  */ 
      else if (goon==3) 
         { 
          printf("Enter file name for labelling units: "); 
          scanf("%s",gname); 
          if ((fg=fopen(gname,"r"))==NULL) 
            {
              printf("Error in opening file %s\n",gname); 
              exit(1); 
             }  
          fclose(fg); 
          printf("\n Enter file name of test set: "); 
          scanf("%s", testname); 
          if ((fpt=fopen(testname,"r")) == NULL) 
            { 
             printf ("\nERROR IN OPENING FILE %s \n", testname) ; 
             exit (1) ;
            } 
          printf("- Enter file name for confusion matrix: "); 
          scanf("%s",cname); 
          if ((fpcx=fopen(cname,"w")) == NULL) 
            { 
             printf ("\nERROR IN OPENING FILE %s \n", cname) ; 
             exit (1) ; 
            } 
          LABEL_UNITS(v); 
          CONFUSION_MATRIX();
        }   /* end if goon==3  */ 
    } /* Chiude (if goon!=4) */ 
  FREE_ALL(); 
 
} /* Fine del main  */
 
 
void CONFUSION_MATRIX() 
 
{
 int confusion[20][20], field1, field2, i, j, vai; 
 int temp, count, diagonal;/* **int_matrix();*/ 
 int sumofrows, sumofcolumns, rowtotal, columntotal; 
 double total_accuracy,  omit, commit, commitavg, omitavg; 
 
 /*confusion = int_matrix(nuscite+1, nuscite+1);*/ 
 
  for(i = 0; i <= nuscite; i++) 
    for(j = 0; j <= nuscite; j++) 
      confusion[i][j] = 0;
 
  count = 0; 
  diagonal = 0; 
  vai =1; 
  while(vai) 
   { 
    for(i=0;i<nbande;i++) 
      fscanf(fpt,"%lf",&x[i]); 
    vai = ((fscanf(fpt,"%d",&field2))!=EOF); 
    if(!vai)
      break; 
    temp = BEST_UNIT_CHOICE(v); 
    field1=1+label[temp]; 
    count++; 
    confusion[field1][field2]++; 
    if(field1 == field2) 
      diagonal++; 
   } 
  total_accuracy = (double) diagonal * 100.0/count; 
 fprintf(fpcx,"Network weights from file : %s\n",fname);
 fprintf(fpcx,"Labelling from file: %s \n \n ",gname);  
  for(i = 1; i <= nuscite; i++) 
  { 
    for(j = 1; j <= nuscite; j++) 
      fprintf(fpcx, "%4d ",confusion[i][j]); 
    fprintf(fpcx,"\n"); 
  } 
 
  fprintf(fpcx, "\n\n\n"); 

  fprintf(fpcx, "                            Accuracy of                       Accuracy of\n"); 
  fprintf(fpcx, "Row Total    Commission (%%)     Column Total    Omission (%%)"); 
  fprintf(fpcx, "\n--------------------------------------------------------------------------"); 
 
  sumofcolumns = sumofrows = 0; 
  commitavg = omitavg = 0; 
  for(i = 1; i <= nuscite; i++) 
  { 
    rowtotal = columntotal = 0; 
    for(j = 1; j <= nuscite; j++)
    { 
      rowtotal += confusion[i][j]; 
      columntotal += confusion[j][i]; 
    } 
    sumofrows += rowtotal; 
    sumofcolumns += columntotal; 
    if( rowtotal == 0 ) 
      commit = 0.0; 
    else 
      commit =  confusion[i][i] * 100.0 / rowtotal;
    if( columntotal == 0 ) 
      omit = 0.0; 
    else 
       omit =  confusion[i][i] * 100.0 / columntotal; 
    commitavg += commit; 
    omitavg += omit; 
    fprintf(fpcx, "\n%6d         %7.2f          %8d          %6.2f", rowtotal, commit, columntotal, 
omit); 
  } 
  fprintf(fpcx,"\n--------------------------------------------------------------------------");
  fprintf(fpcx, "\n%6d                          %8d", sumofrows, sumofcolumns); 
  fprintf(fpcx,"\n\nAverage Omit Precision per Class   : %6.2f%%",omitavg/nuscite); 
  fprintf(fpcx,"\nAverage Commit Precision per Class : %6.2f%%",commitavg/nuscite); 
  fprintf(fpcx,"\nOverall Classification Accuracy    : %6.2f%%\n",total_accuracy); 
  return; 
} 
 
 
double GETRANDOM(min,max) 
int min, max;
{ 
 long numero; 
 double nf, nmax; 
 
  srand(s_time);   
  numero = rand(); 
  nmax=(double)MAX; 
  nf = (double) ((numero/nmax)*(max - min)); 
  s_time = (unsigned) numero/100;   
  return(nf);
} 
 
int BEST_UNIT_CHOICE (vett) 
double **vett; 
 
/* Chooses the winning unit  */ 
 { 
  int BestUnit ; 
  int i; 
  /*double*/double  dist, BestDist ;
  BestUnit = 0 ; 
  BestDist = BIG; 
 
  for (i=0; i<dimnk; i++) 
     { 
       dist=DISTANCE(x,vett[i]); 
       if (dist < BestDist) 
	 { 
	   BestUnit = i ; 
	   BestDist = dist ;
	 } 
     }
 
  return (BestUnit) ; 
 } 
 
double DISTANCE(u,t) 
double *u; 
double *t; 
{
 int i; 
 double temp=0.0; 
 
 for(i=0;i<nbande;i++) 
   {
   temp +=(u[i] - t[i])*(u[i]-t[i]); 
   }
 return(temp) ; 
} 

void WEIGHTS_COMPUTE (u) 
int u; 
 
 /* Updates  Kohonen's weights        */ 
 
 { 
 
  int i, j,k,prod, starti, stopi, startj, stopj ; 
   
  starti=((u/ncols)<nb_r)?(-(u/ncols)):(-nb_r);
  stopi=(((u/ncols)+nb_r)>=nrows)?(nrows-1-(u/ncols)):nb_r; 
 
  startj=((u%ncols)<nb_r)?(-(u%ncols)):(-nb_r); 
  stopj=((u%ncols)+nb_r>=ncols)?(ncols-1-(u%ncols)): nb_r; 
 
  for(i=starti;i<=stopi;i++) 
   { 
    prod = u+i*ncols ; 
    for (j=startj; j<=stopj; j++) 
     {
      for (k=0;k<nbande;k++) 
        w[prod+j][k] = w[prod+j][k] + alfa*(x[k] - w[prod+j][k]) ; 
     }  
  } 
 } 
 
 
void K_WEIGHTS_INITIALIZE () 
 

 { 
  int i,j ;  
  /*double ran3();*/  
  int idum,t1; 
 
  printf("please input random number seed\n");
  scanf("%d",&idum);
  for (j=0; j<dimnk; j++) 
   { 
    for(i=0;i<nbande;i++)
     { 
      t1 = rand()*idum;
           if(t1 < 0)
            t1 = -t1;
           if(t1 > 255)
            t1 = t1%255;
      w[j][i] =(double) t1;
      
/*      w[j][i] = ran3(&idum)*255.0; 
       noise = GETRANDOM(0,255) ;
       w[j][i] = noise; 
*/ 
     } 
   } 
/* 
   printf("Initial values \n"); 
   for(i=0;i<dimnk;i++) 
    { 
     for(j=0;j<nbande;j++) 
       printf("%.2f\t",w[i][j]);
     printf("\n"); 
    } 
*/ 
 } 
 
 
 
 
void SAVE_TRAINING (s,c) 
int s,c;
 
 
 { 
  FILE *fw ; 
  int i,j ; 
  char num[28]; 
 
  if(c) 
   { 
    sprintf(num,"%s%d%s","TEMP",c,".tmp");
    strcpy(savename,num); 
   } 
 
  if ((fw=fopen(savename,"w")) == NULL) 
    { 
      printf ("\nERROR IN OPENING FILE %s \n", savename) ; 
      exit (1) ; 
    } 
 
    fprintf (fw,"nbands            %d\n",nbande) ;
    fprintf (fw,"ncolons           %d\n",ncols) ; 
    fprintf (fw,"nrows             %d\n",nrows) ; 
    fprintf (fw,"nclass            %d\n",nuscite) ; 
    fprintf (fw,"Number_of_lines   %d\n",lines) ; 
    fprintf (fw,"Number_of_pixels  %d\n",pixels) ; 
    fprintf (fw,"Number_of_K_iter  %d\n",nkiter) ; 
 
    if (c!=0) 
      { 
       fprintf(fw,"Number_of_iterm   %d\n",c);
       fprintf(fw,"Image file: %s\n",nomefile); 
      } 
    printf("- Saving Kohonen layer's weights ... \n"); 
 
    for (i=0; i<dimnk; i++) 
      {  
       for (j=0;j<nbande;j++) 
         fprintf(fw,"%lf   ",w[i][j]); 
       fprintf(fw,"\n"); 
      }
    if (s == 2) 
      { 
        printf("- Saving LVQ1's layer weights ..."); 
        fprintf(fw,"LVQ1_weights \n"); 
        for (i=0; i<dimnk; i++) 
         { 
           for(j=0;j<nbande;j++) 
             fprintf (fw,"%lf   ",v[i][j]) ; 
           fprintf(fw,"\n"); 
          }
       fprintf (fw,"Number_of_LVQ1_iter  %d\n",ngiter) ; 
      } 
                                                                                                    
    printf ("\n      The network training is saved in file %s\n\n", savename); 
    fclose (fw) ; 
 } 
 
 
 
void MEMORY_ALLOCATION ()
 
 
 { 
  int i,j ; 
 
  x = (double *)calloc(nbande,sizeof(double)) ; 
  y = (unsigned char *)calloc(pixels,sizeof(unsigned char)) ; 
  g = (unsigned char *)calloc(pixels,sizeof(unsigned char)) ; 
   

  cl_kcount = (int *)calloc(dimnk,sizeof(int)); 
  cl_gcount = (int *)calloc(dimnk,sizeof(int)); 
  cl_cl_count = (int *)calloc(nuscite,sizeof(int)); 
 
  w = (double **)calloc(dimnk,sizeof(double*)) ; 
  for (i=0;i<dimnk;i++) 
    w[i]=(double *)calloc(nbande, sizeof(double)); 
   
  v = (double **)calloc(dimnk, sizeof(double*)); 
  for (i=0;i<dimnk;i++)
    v[i]=(double *)calloc(nbande, sizeof(double)); 
  
 /********************added by Brandt*********************/

 cls_mean = (double **)calloc(dimnk, sizeof(double*)); 
  for (i=0;i<dimnk;i++) 
    cls_mean[i]=(double *)calloc(nbande, sizeof(double)); 

 for (i=0;i<dimnk;i++)
   for(j=0;j<nbande;j++)
     cls_mean[i][j] = 0;  
  
 /********************added by Brandt*********************/
 
  bb = (double *)calloc(nbande,sizeof(double)); 
 
  label = (int *)calloc(dimnk,sizeof(int)); 
  for(i=0;i<dimnk;i++) 
     label[i]=0; 

  /*O_param=(Header *)calloc(1,sizeof(Header)); */
} 
 
 
void FREE_ALL() 
 
 { 
  int i; 
 
  free(x);
  free(y); 
  free(g);

  free(cl_kcount);
  free(cl_gcount);
  free(cl_cl_count);
 
  for (i=0;i<dimnk;i++) 
    free(w[i]); 
  free(w);

  for (i=0;i<dimnk;i++) 
    free(v[i]); 
  free(v);

  for (i=0;i<dimnk;i++) 
    free(cls_mean[i]); 
  free(cls_mean); 
 
  free(label);
  free(bb); 
  free(buf_i); 
 
 } 
 
 
 
void  READ_IMAGE() 
/* Reads the image in ERDAS format and puts the data into the structure 
I_param */
 
 { 
  printf("- Enter image file name (format:band interleaved by line) "); 
  scanf("%s",nomefile); 
  /*fp_r=ErdasOpenToRead(nomefile);*/  
  /*l_righe=VaxErdasReadHeader(fp_r,&I_param) ;*/
   if ((fp_r=fopen(nomefile,"rb")) == NULL)
    {
     printf ("\nERROR IN OPENING FILE %s\n", nomefile) ;
     exit (1) ;
    } 
 printf("input Number of lines in each band image : %d\n",lines); 
 scanf("%d",&lines); 
 printf("input Number of pixels perline of the image: %d\n",pixels); 
 scanf("%d",&pixels); 
 printf("input Number of bands in the image : %d\n",nbande);
 scanf("%d",&nbande);
 l_righe = nbande*pixels+1;  
 buf_i = (unsigned char *)calloc(l_righe,sizeof(unsigned char )) ; 
 if (buf_i==NULL)
   { 
     printf("- Not enough memory! \n"); 
     exit(1); 
   } 
/* 
 lines=I_param.irows; 
 pixels= I_param.icols; 
 nbande=I_param.nbands; 
 */

} 
	   
 
void READ_TRAINING (r) 
int r; 
 
/* Reads from the file the results of previous training sessions  */ 
 
 { 
  FILE *fr ;
  int i,j ; 
  char stringa[40] ; 
 
 
  if ((fr=fopen(fname,"r")) == NULL) 
    { 
     printf ("\nERROR IN OPENING FILE %s\n", fname) ; 
     exit (1) ; 
    } 

  fscanf (fr,"%s  %d",stringa,&nbande) ; 
  fscanf (fr,"%s  %d",stringa,&ncols) ; 
  fscanf (fr,"%s  %d",stringa,&nrows) ; 
  fscanf (fr,"%s  %d",stringa,&nuscite) ; 
  fscanf (fr,"%s  %d",stringa,&lines) ; 
  fscanf (fr,"%s  %d",stringa,&pixels) ; 
  fscanf (fr,"%s  %d",stringa,&np_kiter) ; 
  printf("- Number of Kohonen iterations performed up to now: %d\n",np_kiter); 
 
  dimnk = ncols*nrows ;

  MEMORY_ALLOCATION () ;
	for (i=0; i<dimnk; i++)
		{
		printf("\nweight%d\n",i);
		for (j=0;j<nbande;j++)
			{
			fscanf (fr,"%lf",&w[i][j]) ;
			printf("%lf ",w[i][j]) ;
			}
		}

  if (r==3)
     {
      fscanf(fr,"%s",stringa);
      for (i=0; i<dimnk; i++)
	  {
		printf("\nweight%d\n",i);
        for(j=0;j<nbande;j++) 
		{
          fscanf (fr,"%lf",&v[i][j]) ; 
          printf("%lf ",v[i][j]) ;
		}
	  }
      fscanf (fr,"%s  %d",stringa,&np_giter) ; 
    } 
  printf ("\n\nTHE NETWORK IS READY !\n\n") ; 
  fclose (fr) ; 
 
 } 
 
void TRAIN_NET1 (flag) 
int flag; 
 
 /* Network's training procedure */ 
 
 { 
  int p, nc_count ; 
  int i,j,banda; 
  int opt,dif,max_no_change=0 ; 
  int iter,n_save, cont, ngb;  
 
 
 
  if ((flag == 1)||(flag==2)) 
    { 
      READ_IMAGE(); 
      printf("- Training Kohonen layer ... \n"); 
      if (flag==1) 
	{ 
	  printf ("\n\n- Enter number of iterations: ") ; 
	  scanf ("%d",&nkiter) ; 
 
	  printf ("- Enter number of Kohonen units per row :  ") ; 
	  scanf ("%d",&ncols) ; 
 
	  printf ("- Enter number of rows       :  ") ; 
	  scanf ("%d",&nrows) ; 
          
          printf("- Enter number of classes     : "); 
          scanf("%d",&nuscite); 
 
	  dimnk = ncols*nrows ; 
 
	  MEMORY_ALLOCATION () ; 
	  K_WEIGHTS_INITIALIZE () ; 
 
	} /* chiude if flag==1  */ 
 
      else if (flag==2) 
	{ 
	  printf("\n\n"); 
	  printf ("\n- Enter number of new iterations: ") ; 
	  scanf ("%d",&nkiter) ; 
	 } 
      printf ("\n\n") ; 
      printf("Enter <n_save>: "); /*period of times for saving feature map weight*/
      scanf("%d", &n_save);   
 
/* SOM training   */

      ngb =1; 
      dif = nbande-1 ;
      cont=np_kiter;
      nc_count=0;
      for (iter=0; iter<nkiter;iter++)
	 {
	   printf("\nIter=%d",iter);
	   PARAM_COMPUTE (0,ngb,cont) ;
	   for (i=0; i<lines; i++)
	     {
	       /*ErdasReadImageLine(fp_r,buf_i,l_righe);*/
              fread(buf_i,nbande*pixels,1,fp_r);
	      for (j=0; j<pixels; j++)
		 {
                  for (banda=0; banda<nbande; banda++)
                     x[banda]=(double) buf_i[j+banda*pixels];
		  p = BEST_UNIT_CHOICE (w) ;
		  WEIGHTS_COMPUTE (p) ;
	         } /*  end for j  */
	       }  /*  end for i  */
	    cont++;
	    if ((cont%n_save) == 0) 
	       SAVE_TRAINING(1,cont);
	  START_AGAIN(fp_r);
         }      /* chiude iter */

      printf ("\nKOHONEN LAYER IS TRAINED\n\n") ;
      printf ("    1) To save Kohonen layer and stop the training\n") ;
      printf ("    2) To skip\n");
      printf ("Enter your choice :  ") ;
      scanf ("%d",&opt) ;

      switch (opt)
	    {
	      case 1 : 
		nkiter=nkiter+np_kiter;
	        printf("- Enter filename for storing weights: ");
	        scanf("%s",savename);
		SAVE_TRAINING(1,0) ;                 
                if (nc_count>max_no_change)
                  max_no_change=nc_count;
		printf("Max n. of iterations without changes: %d\n",max_no_change);
		goto l1 ;
		break ;

	      case 2 : 
		break ;
	    }

    } /* Chiude (if flag==1 || flag==2) */



  nkiter=nkiter+np_kiter;
  printf ("- Enter supervised training file name : ");
  scanf ("%s",gname) ;
  /* printf("- Enter number of training block in training set :");
  scanf ("%d",&num_train_cell);
  printf("- Enter number of band");
  scanf ("%d",&nbande);*/

  if ((fg=fopen(gname,"r")) == NULL)
    {
      printf ("\nERROR IN SUPERVISED TRAINING FILE\n") ;
      exit (1) ;
    }

  printf ("- Enter number of iterations : ") ;
  scanf ("%d",&ngiter) ;


/* Supervised LVQ1 training    */



  printf("\n\nLearning Vector Quantization....\n");
  cont=np_giter;
  if(flag==3)
   {
     for(i=0;i<dimnk;i++)
       for(j=0;j<nbande;j++)
         v[i][j]=w[i][j];     /* Copia i pesi dello strato unsupervised */
    }
  LABEL_UNITS(v);
  for (iter=0; iter<ngiter; iter++)
    {
     printf("\n Iter= %d",iter);
     PARAM_COMPUTE(1,0,cont);
     LVQ();
     cont++;
     LABEL_UNITS(v);
    }      /* chiude iter */
  
  fclose(fg); 
 
  printf ("\nNETWORK TRAINING IS TERMINATED\n\n") ; 
  printf ("     1) To save the network training\n" ); 
  printf ("     0) To skip\n") ; 
  printf ("Enter your choice :  ") ; 
  scanf ("%d",&opt) ; 
  if (opt == 1) 
    { 
      printf("- Enter filename for storing weights: "); 
      scanf("%s",savename); 
      ngiter=ngiter+np_giter; 
      SAVE_TRAINING (2,0) ; 
    } 
 
l1:  printf ("\nTHE NETWORK IS READY !\n\n") ; 
 
 
 }
 
void LABEL_UNITS(vett) 
double **vett; 
{ 
 int i,p, j,vai, best, high, n_cl,ww=0; 
 int **majority; 
 
 majority = (int **) calloc(dimnk, sizeof(int *)); 
 for (i=0;i<dimnk;i++) 
   { 
    majority[i]=(int *)calloc(nuscite, sizeof(int)); 
    for (j=0;j<nuscite;j++) 
      majority[i][j]=0; 
   }  
 
 vai =1; 
 while (vai)
   /*for(count1=0;count1 < num_train_cell;count1++)*/
  { 
  
   for(i=0;i<nbande;i++)                /*training set format:[double] -|               */
     fscanf(fg,"%lf", &bb[i]);          /*                    [double] -number of bands */
   vai=((fscanf(fg,"%d",&n_cl))!=EOF);  /*                    [double] -|               */
			 	        /*                    [integer]-class type      */
   if (!vai) 
     break; 
   for (i=0;i<nbande;i++) 
     x[i]=(double) bb[i]; 
   p = BEST_UNIT_CHOICE(vett); 
   majority[p][n_cl-1]++; 

  } 
 rewind(fg); 
 for(i=0;i<dimnk;i++) 
  { 
   best=0; 
   high=0; 
   for (j=0;j<nuscite;j++) 
      { 
       if(majority[i][j]>high) 
         { 
          high=majority[i][j]; 
          best=j;  
         } 
       } 
   if(high)  
    label[i]=best; 
   else 
    label[i]=-1; 
  } 
 for(i=0;i<dimnk;i++) 
   free(majority[i]); 
 free(majority); 
/* for(i=0;i<dimnk;i++) 
     printf("Label of cluster %d  %d\n",i+1,label[i]+1); 
     */   
  
}
 
void LVQ () 
 /*  Updates the weights according to the LVQ1 algorithm calling the 
     procedure LVQ_WEIGHTS */   
 
 { 
  int p, dif ,n_class, vai; 
  int i,banda ; 

  dif = nbande -1; 
  vai=1;  

  while (vai) 
    { 
     for(i=0;i<nbande;i++) 
       fscanf(fg,"%lf",&bb[i]); 
     vai=((fscanf(fg,"%d",&n_class))!=EOF); 
     if (!vai) 
       break; 
     for (banda=0;banda<nbande;banda++) 
        x[banda] = (double) bb[banda] ; 

     p = BEST_UNIT_CHOICE (v) ; 
     LVQ_WEIGHTS (p,n_class) ; 
   
       
    }  /* end while  */ 
  printf("  Correct = %d   Incorrect = %d\n",correct, incorrect); 
  correct=incorrect=0; 
  rewind(fg); 
 }
 
 
void LVQ_WEIGHTS (g, n_c) 
int g,n_c; 
 
 
 { 
  int i,j;

  j = 0; 
   
  if((n_c-1)==label[g]) 
   { 
    for (i=0;i<nbande;i++) 
      v[g][i]=v[g][i] + gamma1*(x[i]-v[g][i]);  
    correct++;  
   } 
  else 
   { 
    for (i=0;i<nbande;i++) 
      v[g][i]=v[g][i] - gamma1*(x[i]-v[g][i]); 
    incorrect++;  
  } 
 }
 
 
void START_AGAIN(ff) 
FILE *ff; 
{ 
 /*Header I_dummy;*/ 
 
 rewind(ff); 
 /*mmy = VaxErdasReadHeader(ff,&I_dummy);*/
} 
 
void PARAM_COMPUTE (sel,ng,ter)
int sel, ng,ter;
 
 
 { 
  double provv1; 
  double prv1, prv2; 
  double k2,x_1; 
 
 
  k2 = (double) 0.02;  
  x_1 = (double) 0.0023;  
  switch (sel) 
	{ 
	  case 0 : 
           if (ter<1000)  
	     alfa =(double) 0.9*(1-ter/1000.) ; 
           else 
             alfa =(double) 0.01*(1-ter/10000.); 
           if (ng) 
             { 
              if(ter < 1000) 
               { 
                provv1=(double) (ncols/3.)*(1-ter/1000.); 
                nb_r=ROOF(provv1); 
               } 
              else if (ter<2000) 
                nb_r=1; 
              else 
               nb_r=0; 
              printf("  Raggio = %d\n", nb_r); /*I think is neighborhood size parameter*/
             } 
           else  
            nb_r=0; 
           break; 
 
	  case 1 : 
	   prv1 = -(x_1*ter); 
           prv2=k2*exp (prv1); 
	   gamma1 = (double) prv2; 
           nb_r=0;	 
         break; 
	} 
 
} 
 
int ROOF(num) 
double num;
 
{ 
 int temp1; 
 double temp2; 
 
 temp1=(int) num; 
 temp2=(double) (num-temp1); 
 if (temp2 == 0.0) 
   return(temp1); 
 else 
   return(temp1+1); 
}