/*	s_icm.c
**	this program is used for classification using
**      Iterative Conditional Modes - ICM proposed by
**      Besag, J. (1986) On the statistical analysis of dirty pictures"
**      Journal of the Royal Statistical Society, 48(3), 259-302.
**      Note that the inputs are in following format (probability matrix).
**      User can use the programme Maximul_likelihood.c
**      to generate such matrix.
*******************************************
U Mtrix:
       pixel number 1   2   3   4   5 ...
class1             0.1 0.3  .   .   .
class2             0.5  .
.                   .   .
.                   .
.                   .
                   sum = 1
*******************************************
COPY RIGHT: PAUL M. MATHER & BRANDT TSO
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define num_of_converge 12 /*define number of iterations for runing ICM to converged*/
#define band_max 12
#define centre 1
#define win_size 3

void main(void)
{
  FILE 	 *fin[band_max];
  FILE 	 *fout[5];
  unsigned char in_file[50];
  int img_row;
  int img_col;
  int count[10],i;
  double weight[20];
  double good[16],num_units[16], energy_sum;
  double sum_row[16],confusion[16][16];
  double **temp, **image[band_max],**window;
  int    num_of_class, num_of_band;
  int    site1, site2,classname,classname1,mark,tot[10],**temp_img;
  double max1,max2, max3;
  double beta[10],u_f;




      printf("Please input the total number of input matrices files(maximum %d):",band_max-2);
      scanf("%d",&num_of_band);
	  for(count[1]=1; count[1]<=num_of_band; count[1]++)
		{
	     printf("Please input the file name for file %d:", count[1]);
	     scanf("%s",in_file);
      	 if((fin[count[1]] = fopen(in_file,"r")) == NULL)
			{
	         printf("Failed to open %s\n",in_file);
	         exit(EXIT_FAILURE);
			}
		}


	  printf("Ouput probability file name:");
	  scanf("%s",in_file);
      if((fout[1] = fopen(in_file,"w")) == NULL)
		{
	      printf("Failed to open %s\n",in_file);
	      exit(EXIT_FAILURE);
		}

	  printf("Ouput classified image name:");
	  scanf("%s",in_file);
      if((fout[2] = fopen(in_file,"wb")) == NULL)
		{
	      printf("Failed to open %s\n",in_file);
	      exit(EXIT_FAILURE);
		}

      for(count[1]=1; count[1]<=num_of_band; count[1]++)
		{
		 printf("please input weighting factor for file %d:",count[1]);
		 scanf("%lf", &weight[count[1]-1]);
		}

      printf("Please input the total of classes:");
	  scanf("%d",&num_of_class);


      printf("Please input the beta values (Markovian neighbourhood coefficients)!\n");
	  for(count[1]=0; count[1]<4; count[1]++)
	  {
	    printf("Please input beta value %d:",count[1]+1);
		scanf("%lf",&beta[count[1]]);
	  }

	  printf("Input image row:");
	  scanf("%d",&img_row);
	  printf("Input image col:");
	  scanf("%d",&img_col);

      printf("If you are making multiple sources classification\n");
	  printf("Please make sure the class SEQUENCE from your input files are the SAME\n");
      printf("***********************************************************************\n");

      printf("Assigning Memory, Please Wait!!!\n");


	  for(count[1]=1; count[1]<=num_of_band; count[1]++)
          image[count[1]] = (double **)calloc(img_row,sizeof(double *));
      temp = (double **)calloc(img_row,sizeof(double *));
      temp_img =(unsigned int **)calloc(img_row,sizeof(unsigned int *));
      window = (double **)calloc(img_row,sizeof(double *));

      for(i = 0 ; i < img_row; i++)
		{

		for(count[1]=1; count[1]<=num_of_band; count[1]++)
           image[count[1]][i] = (double *)calloc(img_col*num_of_class,sizeof(double));

        temp[i] = (double *)calloc(img_col*num_of_class,sizeof(double));
        temp_img[i] = (unsigned int *)calloc(img_col,sizeof(unsigned int));
        window[i] = (double *)calloc(img_col*num_of_class,sizeof(double));
	  }






	   for(count[1] = 0; count[1] < num_of_class; count[1]++)
		   {
		   for(count[2] = 0; count[2] < num_of_class; count[2]++)
		    {
		      confusion[count[1]][count[2]] = 0;
	            }
		good[count[1]] = 0;  num_units[count[1]] = 0;
        	 sum_row[count[1]] = 0;
              	 }


      printf("Reading probabilities, Please wait!!\n");

/*     for(count[1] =0;count[1]<img_row;count[1]++)
       for(count[2]=0;count[2]<img_col;count[2]++)
	  {
	   fread(pix,1,1,fin5);
	   truth[count[1]][count[2]]=pix[0];
	   }
*/


      for(count[1] = 0;count[1] < img_row;count[1]++)/*read probability into buffer*/
	{
         for(count[2] = 0;count[2] < img_col*num_of_class;count[2]++)
	   {
			 for(count[3]=1; count[3]<=num_of_band; count[3]++)
	            fscanf(fin[count[3]],"%lf ",&image[count[3]][count[1]][count[2]]);

            temp[count[1]][count[2]] = 0;
	   }
	}


       printf("Initialising Image\n");



       for(count[5] = 0; count[5] < img_row; count[5]++)
	 	{
         for(count[6] = 0; count[6] < img_col; count[6]++)
 		  {

          site1 = count[5];
          site2 = count[6]*num_of_class;
          temp[site1][site2] = 0;
		  for(count[1]=1; count[1]<=num_of_band; count[1]++)
			temp[site1][site2] = temp[site1][site2]+weight[count[1]-1]*image[count[1]][site1][site2];

          max1=temp[site1][site2];
          classname1=1;
          for(count[0] = 1; count[0] < num_of_class; count[0]++)
	      {
			  temp[site1][site2+count[0]] = 0;
			  for(count[1]=1; count[1]<=num_of_band; count[1]++)
                   temp[site1][site2+count[0]] = temp[site1][site2+count[0]]+weight[count[1]-1]*image[count[1]][site1][site2+count[0]];
                     /* weight[0]*image1[site1][site2+count[0]]+weight[1]*image2[site1][site2+count[0]]
                     +weight[2]*image3[site1][site2+count[0]]+weight[3]*image4[site1][site2+count[0]]
                     +weight[4]*image5[site1][site2+count[0]];*/
               if(temp[site1][site2+count[0]] > max1)
	          {
                max1 = temp[site1][site2+count[0]];
                classname1 = count[0]+1;
	          }
	      }

           temp_img[count[5]][count[6]] = classname1;/*class nnumber start from 1 to 8*/
		}
	}


      printf("start calculating\n");
      tot[0]=0; tot[1]=0; tot[2]=0; tot[3]=0;
      for(count[1]=0;count[1]<10;count[1]++) /*run 5 or 6 times to converge*/
	{
        printf("Iteration %d\n",count[1]+1);
        energy_sum =0;
       for(count[5] = 1; count[5] < img_row-1; count[5]++)
	 {
        for(count[6] = 1; count[6] < img_col-1; count[6]++)
	  {
          site1 = count[5];/*set pixel t*/
          site2 = count[6]*num_of_class;/*set pixel t*/
          max2=temp[site1][site2];
          mark=0;
	  for(count[2] = 0; count[2] < num_of_class; count[2]++)
            {
	   for(count[3] = 0; count[3] < win_size; count[3]++)
	     {
            for(count[4] = 0; count[4] < win_size; count[4]++)
	      {
	       if((count[3]==centre)&&(count[4] == centre)) count[4]++;

	       classname = temp_img[count[5]-centre+count[3]][count[6]-centre+count[4]]-1;
	       /* if(count[1]==1&&count[2]==0) printf(" %d-%d ",count[5]-centre+count[3],count[6]-centre+count[4]);*/
	       if(count[2]==classname)
                 {
                  if((count[3]==0 && count[4]==0) || (count[3]==2 && count[4]==2)) tot[0]++;
                  if((count[3]==0 && count[4]==1) || (count[3]==2 && count[4]==1)) tot[1]++;
                  if((count[3]==0 && count[4]==2) || (count[3]==2 && count[4]==0)) tot[2]++;
                  if((count[3]==1 && count[4]==0) || (count[3]==1 && count[4]==2)) tot[3]++;
                 }
	      }
	     }
            u_f = beta[0]*tot[0]+beta[1]*tot[1]+beta[2]*tot[2]+beta[3]*tot[3];
            tot[0]=0; tot[1]=0; tot[2]=0; tot[3]=0;


            if((temp[site1][site2+count[2]]+u_f)>max2)
	      {
               max2 = temp[site1][site2+count[2]]+u_f;
               temp_img[count[5]][count[6]] = count[2]+1;
               max3 = temp[site1][site2+count[2]]-u_f;
	      }
	    }
           energy_sum = energy_sum + max3;

	  }

	 }
         printf("energy sum = %lf\n",energy_sum/(img_row*img_col));
	}/*end of converge running*/


      printf("Writing probabilities, Please wait!!\n");
      for(count[1] = 0;count[1] < img_row;count[1]++)/*writing probability*/
	{
         for(count[2] = 0;count[2] < img_col*num_of_class;count[2]++)
	   {
           fprintf(fout[1],"%lf ",temp[count[1]][count[2]]);
	   }
	}


	    printf("Writing image, Please wait!!\n");
        for(count[1] = 0;count[1] < img_row;count[1]++)/*writing image*/
	{
         for(count[2] = 0;count[2] < img_col;count[2]++)
	   {
           fprintf(fout[2],"%c",temp_img[count[1]][count[2]]);
	   }
	}
      printf("\nDone!!\n");


      for(i = 0 ; i < img_row; i++)
        {

		 for(count[1]=1;count[1]<=num_of_band;count[1]++)
			free(image[count[1]][i]);

         free(window[i]);
         free(temp[i]);
         free(temp_img[i]);
        }


      free(window);
      free(temp);
      free(temp_img);
	  for(count[1]=1;count[1]<=num_of_band;count[1]++)
			free(image[count[1]]);

	  for(count[1]=1;count[1]<=num_of_band;count[1]++)
			fclose(fin[count[1]]);

      fclose(fout[1]);
      fclose(fout[2]);



}