/****************************************************/
/*This program is used to calculate separability    */
/*between training patterns.                        */
/*source:computer processing of remotely-sensed     */
/*images. pp.321-327. by P. M. Mather.              */
/*COPY RIGHT: Brandt Tso           */
/****************************************************/

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

#define band_max 22
#define class_max 13

void showmatrix(z,nn,cc)
/* FILE *file1;*/
double z[band_max][band_max];
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");
		}
      }
}



void cleanmatrix(z,nn,cc)/*to clear matrix z's elements*/
double z[band_max][band_max];
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*/


void transpose(z,nn,cc,tz)/*transpose z matrix*/
double z[band_max][band_max],tz[band_max][band_max];
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*/

void addmatrix(z1,z2,row,col,item)/*add z1 and z2 matrix*/
double z1[band_max][band_max],z2[band_max][band_max],item[band_max][band_max];
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*/

void minusmatrix(z1,z2,row,col,item)/*minus z1 and z2 matrix*/
double z1[band_max][band_max],z2[band_max][band_max],item[band_max][band_max];
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*/

void multiply(z1, z2, row,col1,col2,item)/*multiply matrices z1[row][col1] and z2[col1][col2]*/
double z1[band_max][band_max],z2[band_max][band_max],item[band_max][band_max];
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)*/
void inverse(z,nn,cc,item)
double z[band_max][band_max],item[band_max][band_max];
int nn,cc;
{
  FILE *fout;
  int i,j,k;
  double l[band_max][band_max],u[band_max][band_max];
  double sum,determ;
  double il[band_max][band_max],iu[band_max][band_max];/*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;
      }

  determ = 1;
  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;

      determ = determ*u[i][i];
    }

  if((fout = fopen("1_determ_co","w")) == NULL) /*file used to save determ*/
       {
         fprintf(stderr,"Failed to open %s\n","1_determ_co");
         exit(EXIT_FAILURE);
       }
  fprintf(fout,"%lf ",determ);
  fclose(fout);

  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;
       }
    }
    /*
   multiply(z,item,6,6,6,temp1);
   showmatrix(fout3,temp1,6,6);
   */
/*=======================================================================================*/
  /*             showmatrix(u,nn,cc);
                  showmatrix(iu,nn,cc);
                  showmatrix(l,nn,cc);
                  showmatrix(il,nn,cc);
                  showmatrix(item,nn,cc);*/
/*=======================================================================================*/

}/*end of inverse()*/



void main(void)
{
  FILE 	*fin[band_max], *fout1;
  int   i,j,num_of_class,num_of_band;
  int   count1, count2, count3,store1, num_of_pixel[band_max];
  double avg_div, count_avg;
  double s[band_max][band_max][band_max],is[band_max][band_max][band_max],mean_xi[band_max][band_max],minus_x[band_max][band_max];
  double s1[band_max][band_max],s2[band_max][band_max],is1[band_max][band_max],is2[band_max][band_max],mean1[band_max][band_max],mean2[band_max][band_max];
  double sum_of_wi=0;
  double temp1[band_max][band_max],temp2[band_max][band_max],temp3[band_max][band_max],temp4[band_max][band_max],temp5[band_max][band_max],temp6[band_max][band_max],temp7[band_max][band_max],temp8[band_max][band_max];
  double xi[band_max][band_max], p0 = 2, p1=0.5, p2 = 1.25;
  double f1,f2,diverg,diverg1,mean_xi_1[band_max][band_max];
  unsigned char store[band_max][2],in_file[50];
  fpos_t file_start[band_max];



      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 total number of image bands(maximum %d):",band_max-2);
      scanf("%d",&num_of_band);
	  for(count1=1; count1<=num_of_band; 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("The output file name for recording divergence indices:");
	   scanf("%s",in_file);

      if((fout1 = fopen(in_file,"w")) == NULL)/*output file*/
       {
         fprintf(stderr,"Failed to open %s\n",in_file);
         exit(EXIT_FAILURE);
       }

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


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

      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[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 = (int) store[count3][0];
	         xi[count3][0] = (double) store1;
			 mean_xi[count3][count1] = mean_xi[count3][count1]+xi[count3][0];
			}
	   }
       printf("\nTraining mean for class %d:\n",count1);
       for(count2 = 0; count2 < num_of_band; count2++)
		{
          mean_xi[count2][count1] = (double) mean_xi[count2][count1]/num_of_pixel[count1];
	      printf(" %lf ",mean_xi[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[count2][count1];

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

       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 = (int) 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];
	     }

       }/*End of count1*/


      for(count3 = 0; count3 < num_of_class; 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,is1);

                 for(i = 0; i < num_of_band; i++)
					for(j = 0; j < num_of_band; j++)
						is[count3][i][j] = is1[i][j];

				}



      cleanmatrix(s1,num_of_band,num_of_band);
      cleanmatrix(is1,num_of_band,num_of_band);
      cleanmatrix(s2,num_of_band,num_of_band);
      cleanmatrix(is2,num_of_band,num_of_band);
      cleanmatrix(mean1,num_of_band,num_of_band);
      cleanmatrix(mean2,num_of_band,num_of_band);

      avg_div = 0;
      count_avg = 0;

      for(count1 = 0;count1 < num_of_class;count1++)
		{

        for(i = 0; i < num_of_band; i++)
			mean1[i][count1] =  mean_xi[i][count1];


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

             }

         for(i = 0; i < num_of_band; i++)
             for(j = 0; j < num_of_band; j++)
             {
				 is1[i][j]=is[count1][i][j];

             }


         for(count2 = count1+1;count2 < num_of_class;count2++)
           {
           for(i = 0; i < num_of_band; i++)
			mean2[i][count2] =  mean_xi[i][count2];

            for(i = 0; i < num_of_band; i++)
             for(j = 0; j < num_of_band; j++)
             {
				    s2[i][j]=s[count2][i][j];

             }

            for(i = 0; i < num_of_band; i++)
             for(j = 0; j < num_of_band; j++)
             {
				    is2[i][j]=is[count2][i][j];

             }


             minusmatrix(s1,s2,num_of_band,num_of_band,temp1);
             minusmatrix(is2,is1,num_of_band,num_of_band,temp2);
             minusmatrix(mean1,mean2,num_of_band,1,temp3);
             transpose(temp3,num_of_band,1,temp4);
             addmatrix(is1,is2,num_of_band,num_of_band,temp8);
             multiply(temp1,temp2,num_of_band,num_of_band,num_of_band,temp5);
             multiply(temp8,temp3,num_of_band,num_of_band,1,temp6);
             multiply(temp6,temp4,num_of_band,1,num_of_band,temp7);
             f1 = 0; f2 = 0;
             for(i = 0; i < num_of_band; i++)
	       {
                f1 = f1 + temp5[i][i];
                f2 = f2 + temp7[i][i];
	       }
             diverg = (double) (f1+f2)/2;

             diverg1 = (double) -(f1+f2)/8;
             diverg1 = 2000*(1-exp(diverg1));
	      /*transform Divergence*/
            fprintf(fout1,"calss [%d, %d] %lf or (%lf)\n",count1,count2,diverg,diverg1);
 	        printf("calss [%d, %d] %lf or(%lf)\n",count1,count2,diverg,diverg1);
            avg_div = avg_div + diverg;
            count_avg++;

	   }
	}


      avg_div = (double) avg_div/count_avg;
      fprintf(fout1,"Average divergence %f\n",avg_div);
      printf("Average divergence %f\n",avg_div);



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

      fclose(fout1);

}