/*	s_MAP.c
**	this program is used to accomplish annealing using
**      Maximise A Posterior -MAP based on Simulating Annealing.
**      German and German 1984, IEEE Trans PAMI;721-741
**      Note that the input file(s) is(are) the type of following matrix:
**      User can use the programme Maximum_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
*******************************************
COPYRIGHT: BRANDT TSO
*/

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

#define band_max 12 /*maximal image bands*/
#define win_size 3
#define centre 1

void main(void)
{
  FILE 	 *fin[band_max];
  FILE 	 *fout[5];
  int img_row;
  int img_col;
  int count[10],i,k,num_of_run_at_T;
  double weight[band_max],substitute;
  double **temp, **image[band_max],**window;
  int    num_of_class, num_of_band;
  int    site1, site2,classname,classname1,mark,tot[10],tot1[10],perturb_z,gibbs_x,gibbs_y;
  unsigned int **temp_img;
  unsigned char in_file[50];
  double max1,normal_threshold,energy_sum;
  double beta[10],delta,u_f,u_f1, T,alfa;



      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);
		}

      printf("Ouput energy variation name:");
	  scanf("%s",in_file);
      if((fout[3] = 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("Please input initial temperture for annealing:(suggested value 1 ~ 4)");
	  scanf("%lf",&T);
	  printf("Please input the number of iterations for cooling down schedule:(suggested value 10 ~ 120)");
	  scanf("%d",&k);/*for determining how slow the temporature will cool down*/
	  printf("Please input the number of iterations running under Temperutre T:");
	  scanf("%d",&num_of_run_at_T); /*define number of iterations for running N-inner time at temporature T*/


      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));
	  }


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

     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] = 0; 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]];

               if(temp[site1][site2+count[0]] > max1)
				  {
                max1 = temp[site1][site2+count[0]];
                classname1 = count[0]+1;
				}
			}
	  /* temp_img[count[5]][count[6]]= (int) random()%8+1;*/
	  temp_img[count[5]][count[6]] = classname1;/*class nnumber start from 1 to 8*/
	  }
	}

    printf("Reading finished, Start calculating!!!\n");

    tot[0]=0;tot[1]=0; tot[2]=0; tot[3]=0; tot1[0]=0;tot1[1]=0; tot1[2]=0; tot1[3]=0;
    for(count[7] = 1; count[7] <= k; count[7]++) /*number of iteration for cooling down temperature*/
      {
      fprintf(fout[3],"\ntemperature start %d \n",count[7]);
      for(count[8]=1; count[8] < num_of_run_at_T; count[8]++)/*num of iterations running under temperture T*/
       {
       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*/
           mark = temp_img[count[5]][count[6]]-1;/*find xt*/
	       max1=temp[site1][site2+mark];
	       gibbs_x = (int) rand()%win_size; if(gibbs_x < 0) gibbs_x = -gibbs_x;
           gibbs_y = (int) rand()%win_size; if(gibbs_y < 0) gibbs_y = -gibbs_y;
           perturb_z = temp_img[count[5]-centre+gibbs_x][count[6]-centre+gibbs_y] - 1;


	     if(perturb_z<0) perturb_z = -perturb_z;
         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 = 0;
	        if(((site1-centre+count[3])>=0)&&((site1-centre+count[3])<(img_row-1))&&((site2-centre*num_of_class+count[4]*num_of_class)>=0)&&((site2-centre*num_of_class+count[4]*num_of_class)<(img_col-1)*num_of_class))
		  {
                  max2=temp[site1-centre+count[3]][site2-centre*num_of_class+count[4]*num_of_class];
                  for(count[0] = 0; count[0] < num_of_class; count[0]++)
                    if(temp[site1-centre+count[3]][site2-centre*num_of_class+count[4]*num_of_class+count[0]] > max2)
	 	    {
                      max2 = temp[site1-centre+count[3]][site2-centre*num_of_class+count[4]*num_of_class+count[0]];
                      classname = count[0];
 		    }
		 }*/
               classname = temp_img[count[5]-centre+count[3]][count[6]-centre+count[4]]-1;
	       if(perturb_z==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]++;
			}
           if(mark==classname)
            {
			 if((count[3]==0 && count[4]==0) || (count[3]==2 && count[4]==2)) tot1[0]++;
	         if((count[3]==0 && count[4]==1) || (count[3]==2 && count[4]==1)) tot1[1]++;
			 if((count[3]==0 && count[4]==2) || (count[3]==2 && count[4]==0)) tot1[2]++;
			 if((count[3]==1 && count[4]==0) || (count[3]==1 && count[4]==2)) tot1[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;
            u_f1 = tot1[0]*beta[0]+beta[1]*tot1[1]+beta[2]*tot1[2]+beta[3]*tot1[3];
            delta = temp[site1][site2+perturb_z]+u_f - (max1+u_f1);
			tot1[0]=0; tot1[1]=0; tot1[2]=0; tot1[3]=0;
            if(delta >= 0)
              {
               temp_img[count[5]][count[6]] = perturb_z+1;
			   energy_sum = energy_sum +temp[site1][site2+ temp_img[count[5]][count[6]]-1]- u_f;
			}
            else
              {
              alfa = (double) delta/T;
              alfa = exp(alfa);
              normal_threshold = rand();
              if(normal_threshold<0)
                   normal_threshold=-normal_threshold;

              while (normal_threshold > 1)
                   normal_threshold =normal_threshold/10;

              if(alfa > normal_threshold)
				{
                 temp_img[count[5]][count[6]] = perturb_z+1;
                 energy_sum = energy_sum +temp[site1][site2+ temp_img[count[5]][count[6]]-1]- u_f;
				}
              else
                 energy_sum = energy_sum +temp[site1][site2+ temp_img[count[5]][count[6]]-1]- u_f1;
	      }

	  }/*end of count[6]*/

	 }/*end of count[5]*/
          fprintf(fout[3],"energy = %lf \n",energy_sum/(img_row*img_col));
       }/*end of count[8] for number iteratin running under temper T*/
        substitute = count[7]+1;
        T = (double) T*(log(substitute)/log(substitute+1)); /*cooling*/
        /*	T = (double) C/log(substitute);*/
        printf("T=%lf\n",T);
      }/*end of k 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]++)
		{
         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]);
	  fclose(fout[3]);


}