/*	s_MPM.c ***Maximum a Posterior Margin
**	this program is used to accomplish Markovian classification using
**      Maximise A Posterior Margin based on the paper proposed by
**      Marroquin et al. 1987, "Probabilistic solution of ill-posed
**      problem in computational vision."
**      Journal of the American Statistical Association
**      vol.82, 397, pp.76-89.
**      Note that the input file(s) is(are) the type of following 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:  BRANDT TSO
*/

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

#define num_of_run_at_T 1 /*200 define number of iterations for runing N-inner time at temporature T*/
#define win_size 3
#define k 120 /*for determining how slow the temporature will cool down*/
#define centre 1
#define band_max 12

void main(void)
{
  FILE 	 *fin[band_max];
  FILE 	 *fout[5];
  int img_row;
  int img_col;
  int count[10],i;
  double weight[band_max];
  double **temp, **image[band_max],**window;
  int    num_of_class, num_of_band, kk,**mpm,nn;
  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,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);
		}

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

      T = 1;
      kk=20; /*Stability threshold*/
      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");
      mpm=(int **)calloc(img_row,sizeof(int *));

	  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++)
		{
        mpm[i] = (int *)calloc(img_col*num_of_class,sizeof(int));
		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] = 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("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[8] = 1; count[8] <= 200; count[8]++) /*number of iteration for running under fixed temperature*/
      {
        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()%3; if(gibbs_x < 0) gibbs_x = -gibbs_x;
           gibbs_y = (int) rand()%3; if(gibbs_y < 0) gibbs_y = -gibbs_y;
           perturb_z = temp_img[count[5]-centre+gibbs_x][count[6]-centre+gibbs_y] - 1;

           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(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;
            delta = temp[site1][site2+perturb_z]+u_f - (max1+tot1[0]*beta[0]+beta[1]*tot1[1]+beta[2]*tot1[2]+beta[3]*tot1[3]);
	        tot1[0]=0; tot1[1]=0; tot1[2]=0; tot1[3]=0;
            if(delta >= 0)
               temp_img[count[5]][count[6]] = perturb_z+1;

            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;

			}

            if(count[8] > kk)
			{
               nn = temp_img[count[5]][count[6]]-1;
               mpm[count[5]][count[6]*num_of_class+nn]++;
			}

	     energy_sum = energy_sum + temp[site1][site2+ temp_img[count[5]][count[6]]-1]-u_f;
	  }/*end of count[6]*/

	 }/*end of count[5]*/
       printf("energy sum = %lf\n",energy_sum/(img_row*img_col));
      }/*end of count[8]*/


       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;
			max1 =  temp[site1][site2];
            classname = 1;
            for(count[0] = 0; count[0] < num_of_class; count[0]++)
				{
               if(temp[site1][site2+count[0]] > max1)
					{
                  max1 = temp[site1][site2+count[0]];
                  classname = count[0]+1;
					}
				}
            temp_img[count[5]][count[6]]= classname;
			}
		}

       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];
          site2 = count[6]*num_of_class;

          max1=mpm[site1][site2];
          classname1=1;
          for(count[0] = 0; count[0] < num_of_class; count[0]++)
	      {

               if(mpm[site1][site2+count[0]] > max1)
	          {
                max1 = mpm[site1][site2+count[0]];
                classname1 = count[0]+1;
	          }
	      }

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


    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++)
        {
         free(mpm[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(mpm);
      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]);


}