/****************************************************/
/*This program implements the maximum likelihood    */
/*classification procedure.                         */
/*ouput:1.classified image, and 2. probability file */
/*Note: For constructong variance-                  */
/*covariance matrix must be generic binary file.    */
/****************************************************/

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

#define img_maximum 1024

double det[22];
main()
{
  FILE 		*fin[22], *fout[5];
  int   i,j, img_row, img_col, num_of_class, num_of_band;
  int   count1, count2, count3, count4, store1;
  int   num_of_pixel[22],clas;
  double pai=6.28,mini_f;
  double num_band, t_2,t_3;
  double s[22][22][22],s1[22][22],is[22][22],mean_xi_1[22][22];
  double minus_x[22][22],sum_of_x[22][22],mean_xi_a[22][200];
  double m[22][22], tx[22][22], sum_of_wi=0,temp[2][2];
  double xi[22][22], p0 = 2, p1=0.5, p2 = 1.25;
  unsigned char store[22][1],read_pixel[22][img_maximum],in_file[50], out_file[50];
  long address[22];
  fpos_t *file_start[22];


      printf("\n====Maximal likelihood classification====\n");
	  printf("Please input the file name recording the training pixel number for each class:");
      scanf("%s",in_file);
      if((fin[0] = fopen(in_file,"r")) == NULL)/*pixel number in each training group*/
       {
         printf("Failed to open %s\n",in_file);
         exit(EXIT_FAILURE);
       }

	  printf("Please input the file name for the image to be classified (band interleaved by line):");
      scanf("%s",in_file);
      if((fin[21] = fopen(in_file,"rb")) == NULL)/*being classified image: band interleaved by line*/
       {
         printf("Failed to open %s\n",in_file);
         exit(EXIT_FAILURE);
       }

	  printf("Please input the total number of image bands(maximum 20):");
      scanf("%d",&num_of_band);
	  for(count1=1; count1<=num_of_band; count1++)
		{
		 file_start[count1]=&address[count1];
         printf("Please input training file name for band %d:", count1);
	     scanf("%s",in_file);
      	 if((fin[count1] = fopen(in_file,"rb")) == NULL)
			{
	         printf("Failed to open %s\n",in_file);
	         exit(EXIT_FAILURE);
			}
		}

      printf("\nInput images completed!\n");
      printf("Output probability record file:");
	  scanf("%s",out_file);
      if((fout[1] = fopen(out_file,"w")) == NULL)/*output probability file*/
       {
         printf("Failed to open %s\n",out_file);
         exit(EXIT_FAILURE);
       }

      printf("Output classified image:");
	  scanf("%s",out_file);
      if((fout[2] = fopen(out_file,"wb")) == NULL)/*output classified file*/
       {
         printf("Failed to open %s\n",out_file);
         exit(EXIT_FAILURE);
       }

	  printf("Please input the number for total classes:");
      scanf("%d",&num_of_class);
	  printf("Please input the number of image rows in each band:");
      scanf("%d",&img_row);
	  printf("Please input the number of image columns in each band:");
      scanf("%d",&img_col);


      cleanmatrix(xi,22,1);
      cleanmatrix(sum_of_x,22,1);

	  for(count1 = 1; count1<=num_of_band; count1++)
		  fgetpos(fin[count1],file_start[count1]);

    printf("/nComputing mean vlaues for each class\n");
    for(count1=0;count1 < num_of_class;count1++)/*mean value of each class*/
       {
       fscanf(fin[0],"%d",&num_of_pixel[count1]);/*read training pixel numbers*/
	   printf(" %d ",num_of_pixel[count1]);
       for(count2 = 0; count2 < num_of_band; count2++)
	          mean_xi_a[count2][count1] = 0;

       for(count2=0;count2 < num_of_pixel[count1]; count2++)
		{
         for(count3=0;count3 < num_of_band; count3++)
			{
		     fread(store[count3],1,1,fin[count3+1]);
			 store1 = store[count3][0];
	         xi[count3][0] = (double) store1;
			 mean_xi_a[count3][count1] = mean_xi_a[count3][count1]+xi[count3][0];
			}
	   }

       printf("\nTraining mean for class %d:\n",count1);
       for(count2 = 0; count2 < num_of_band; count2++)
		{
          mean_xi_a[count2][count1] = (double) mean_xi_a[count2][count1]/num_of_pixel[count1];
	      printf(" %f ",mean_xi_a[count2][count1]);
		 }
	   printf("\n");

       }/*End of count1*/




      for(count1 =1; count1<=num_of_band; count1++)
	      fsetpos(fin[count1],file_start[count1]);


     for(count1=0;count1 < num_of_class;count1++)/*covariance matrix and determine*/
       {
       for(count2 = 0; count2 < num_of_band; count2++)
			mean_xi_1[count2][0] = mean_xi_a[count2][count1];

       cleanmatrix(s1,22,22);
       cleanmatrix(s[count1],22,22);

       for(count2=0;count2 < num_of_pixel[count1]; count2++)
		{
		 for(count3 = 0; count3<num_of_band; count3++)
			{
			 fread(store[count3],1,1,fin[count3+1]);
			 store1 = store[count3][0];
             xi[count3][0] = (double) store1;
			}

          minusmatrix(xi,mean_xi_1,num_of_band,1,minus_x);

          for(i = 0; i < num_of_band; i++)
           for(j = 0; j < num_of_band; j++)
			 {
               s[count1][i][j] = s[count1][i][j] + minus_x[i][0]*minus_x[j][0];
               s1[i][j] = s1[i][j] + minus_x[i][0]*minus_x[j][0];
			}

         }/*End of count2*/

       for(i = 0; i < num_of_band; i++)
         for(j = 0; j < num_of_band; j++)
	     {
              s[count1][i][j] = (double)  s[count1][i][j]/num_of_pixel[count1];
              s1[i][j] = (double) s1[i][j]/num_of_pixel[count1];
	     }


       determ(s1,num_of_band,num_of_band,count1);
       printf("determ of covariance matrices:");
       printf(" %f ",det[count1]);
       }/*End of count1*/



      for(count1 = 0; count1 < img_row; count1++)    /*classifier*/
		{
        printf("processing row %d \n", count1+1);

        for(count2 =0;count2 < num_of_band;count2++)/*read tatal band pixels row by row*/
			 fread(read_pixel[count2],1,img_col,fin[21]);

        for(count2 = 0;count2 < img_col ; count2++)
			{
		    for(i = 0; i < num_of_band; i++)
			   xi[i][0] = read_pixel[i][count2];

			mini_f = 9999999999;

			for(count3 = 0; count3 < num_of_class; count3++)
				{
				for(count4 = 0; count4 < num_of_band; count4++)
					mean_xi_1[count4][0] = mean_xi_a[count4][count3];

				for(i = 0; i < num_of_band; i++)
					for(j = 0; j < num_of_band; j++)
						s1[i][j] = s[count3][i][j];

				inverse(s1,num_of_band,num_of_band,is);
				num_band = num_of_band;/*from int to double*/
  				minusmatrix(xi,mean_xi_1,num_of_band,1,minus_x);
				transpose(minus_x,num_of_band,1,tx);
				multiply(tx,is,1,num_of_band,num_of_band,m);
				multiply(m,minus_x,1,num_of_band,1,temp);
				t_2 = 0.5*(temp[0][0]+log(det[count3])+num_band*log(pai));


            if(t_2 < mini_f)
			{
              mini_f = t_2;
              clas = count3;
			}
            t_3 = exp(-t_2);
            fprintf(fout[1],"%lf ",-t_2);
		  }   /*End of count3 num_of_class*/
	    fprintf(fout[2],"%c",clas+1);
		}/*End of count2 img_col*/
	}/*End of count1  img_row*/

     for(count1 = 0; count1<=num_of_band; count1++)
	 {
        fclose(fin[count1]);
	 }
	 fclose(fin[21]);
     fclose(fout[1]);
     fclose(fout[2]);

}


/*This function is to show matrix*/
showmatrix(z,nn,cc)
/* FILE *file1;*/
double z[22][22];
int   nn,cc;
{
  int i,j;
  printf("\n");
  printf("The matrix is:!\n");

  for(i=0; i < nn; i++)
    for(j=0; j < cc; j++)
      {
        printf("%15f",z[i][j]);
        if (j == cc-1)
		{
           printf("\n");
		}
      }
}



cleanmatrix(z,nn,cc)/*to clear matrix z's elements*/
double z[22][22];
int nn,cc;
{
  int i,j;

  for(i=0; i < nn ; i++)
    for(j=0; j < cc ; j++)
        z[i][j] = 0;
}/*end of cleanmatrix*/


transpose(z,nn,cc,tz)/*transpose z matrix*/
double z[22][22],tz[22][22];
int nn,cc;
{
  int i,j;

  for(i=0; i < nn ; i++)
    for(j=0; j < cc ; j++)
         tz[j][i] = z[i][j];

}/*end of transpose*/

addmatrix(z1,z2,row,col,item)/*add z1 and z2 matrix*/
double z1[22][22],z2[22][22],item[22][22];
int   row,col;
{
 int   i,j;


  for(i=0; i < row; i++)
     for(j=0; j < col; j++)
        item[i][j] = z1[i][j]+z2[i][j];

}/*end of add*/

minusmatrix(z1,z2,row,col,item)/*minus z1 and z2 matrix*/
double z1[22][22],z2[22][22],item[22][22];
int   row,col;
{
 int   i,j;


  for(i=0; i < row; i++)
     for(j=0; j < col; j++)
        item[i][j] = z1[i][j]-z2[i][j];

}/*end of add*/

multiply(z1, z2, row,col1,col2,item)/*multiply matrices z1[row][col1] and z2[col1][col2]*/
double z1[22][22],z2[22][22],item[22][22];
int row,col1,col2;
{
  int i,j,k;
  double sum;

  sum = 0;
  cleanmatrix(item,row,col2);
  for(i=0; i < row; i++)
    {
     for(j=0; j < col2; j++)
       {
        for(k=0; k < col1; k++)
	  {
            sum = sum + z1[i][k]*z2[k][j];
	  }
        item[i][j] = sum;
        sum = 0;
       }
    }

}/*end of multiply*/



/*This function is used to calculate inverse matrix of m[][]*/
/*The method is to divide original matrix into Lower and Upper matrices respectively*/
/*from the theorem: A=LU -> Inver(A)=inver(U)inver(L)*/
/*Hence, based on the theorem, the matrix A must be a square.(n = c)*/
inverse(z,nn,cc,item)
double z[22][22],item[22][22];
int nn,cc;
{
  int i,j,k;
  double l[22][22],u[22][22];
  double sum;
  double il[22][22],iu[22][22];/*il[][] and iu[][] is the inverse matrix of l[][] and u[][]*/


  for(i=0; i < nn; i++)/*setup the elements of u[][] and l[][]*/
    for(j=0; j < cc; j++)
      {
        if(j < i)/*lower triangle elements = 0*/
          u[i][j] = 0;
        else
	  {
            u[i][j] = z[i][j];
            for(k=0; k < i ; k++)
                u[i][j] = u[i][j] - l[i][k]*u[k][j];
          }

        if(j < i)
           l[j][i] = 0;/*upper triangle elements = 0*/
        else if(j == i)/*diagonal elements = 1*/
           l[j][i] = 1;
        else
	  {
            l[j][i] = z[j][i];
            for(k=0 ; k < i ; k++)
                l[j][i] = l[j][i] - l[j][k]*u[k][i];
            l[j][i] = l[j][i]/u[i][i];
	  }
      }


  /*next to calculate the inverse of u[][] and l[][]*/
  for(i=0 ; i < nn; i++)
    for(j=0; j < cc; j++)
      {
        il[i][j] = 0;
        iu[i][j] = 0;
      }

  for(i=0; i < nn; i++)/*setup inverse matrices diagonal elements*/
    {
      il[i][i] = 1;
      if(u[i][i] != 0)
         iu[i][i] = 1/u[i][i];
      else
         iu[i][i] = 0;
    }

  for(i=0; i < cc-1; i++)/*setup L matrix rest elements*/
    {
      for(j=i+1; j < nn; j++)
        {
          il[j][i] = -l[j][i];
          for(k=i+1;k < j;k++)
	    {
              il[j][i] = il[j][i] - l[j][k]*il[k][i];
	    }

        }
    }/*L matrix end setup*/

  for(i=0; i < nn; i++)/*setup U matrix rest elements*/
    {
      for(j=i+1; j < cc; j++)
        {
          iu[i][j] = 0;
          for(k=i;k < j;k++)
	    {
              iu[i][j] = iu[i][j] - iu[i][k]*u[k][j];
	    }
          if(u[j][j] != 0)
            iu[i][j] = iu[i][j]/u[j][j];
	}
    }/*U matrix end setup*/



  /*next calculate the inverse of matrix A where A=LU,A(-1)=U(-1)L(-1)*/
  sum = 0;
  for(i=0; i < nn; i++)
    {
     for(j=0; j < cc; j++)
       {
        for(k=0; k < nn; k++)
	  {
            sum = sum +  iu[i][k]*il[k][j];
	  }
        item[i][j] = sum;
        sum = 0;
       }
    }

/*=======================================================================================*/
  /*             showmatrix(u,nn,cc);
                  showmatrix(iu,nn,cc);
                  showmatrix(l,nn,cc);
                  showmatrix(il,nn,cc);
                  showmatrix(item,nn,cc);*/
/*=======================================================================================*/

}/*end of inverse()*/

/*This function calculates determination of matrix*/

determ(z,nn,cc,class_num)
double z[22][22];
int nn,cc,class_num;
{
  int i,j,k;
  double l[22][22],u[22][22];
  double item=0;

  for(i=0; i < nn; i++)/*setup the elements of u[][] and l[][]*/
    for(j=0; j < cc; j++)
      {
        if(j < i)/*lower triangle elements = 0*/
          u[i][j] = 0;
        else
	  {
            u[i][j] = z[i][j];
            for(k=0; k < i ; k++)
                u[i][j] = u[i][j] - l[i][k]*u[k][j];
          }

        if(j < i)
           l[j][i] = 0;/*upper triangle elements = 0*/
        else if(j == i)/*diagonal elements = 1*/
           l[j][i] = 1;
        else
	  {
            l[j][i] = z[j][i];
            for(k=0 ; k < i ; k++)
                l[j][i] = l[j][i] - l[j][k]*u[k][i];
            l[j][i] = l[j][i]/u[i][i];
	  }
      }
   item = u[0][0];
   for(k=1; k < nn ; k++)
         item = item + item * u[k][k];

   printf("( %lf) ",item);

   if(item < 0.0)
      item = -item;

   det[class_num] = item;

   //return(item);
}