/****************************************/
/*This program is used to estimate sub- */
/*pixel landercover  percentage, based  */
/*on linear mixture model using 2 bands.*/
/*from: INT. J. Remote Sensing 1993     */
/*vol. 14, no. 6, pp 1159 - 1177.       */
/****************************************/


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


float m[16][16],m_hat[16][16],n_matrix[16][16];
main(int argc, char **argv)
{
  int   c,n,t,i,j,w, img_length,l,m;
  unsigned char  ch,ch1,*store1,*store2;
  float im[16][16],f_hat[16][16],f[16][16];
  FILE *fin1, *fin2, *fout;




/*where n is the number of spectra and c is the component number(end member)*/
/* m[][] is the end member matrix for each band f[] is the proportion matrix  */
/* x[] is the DN of observed pixels, im[][] is the inverse of m[][], tm[][] is*/
/* transpose matrix of matrix m[][], mm[][] is T(m)*m matrix, e[][] is  the   */
/* noise matrix, tf[][] is the transpose matrix of f[][]                      */

 if(argc == 6)
    {
      img_length = atoi(argv[3]);
      img_width = atoi(argv[4]);

      store1 = (unsigned char *)malloc(2);
      store2 = (unsigned char *)malloc(2);
      if((fin1 = fopen(argv[1],"rb")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n",argv[1]);
	  exit(EXIT_FAILURE);
	}
      if((fin2 = fopen(argv[2],"rb")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n",argv[2]);
	  exit(EXIT_FAILURE);
	}
      if((fout = fopen(argv[5],"w")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n",argv[4]);
	  exit(EXIT_FAILURE);
	}

  printf("Welcome! The purpose of this program is to test linear mixture model!\n");
  printf("Two variables need to be defined!\n");
  c = number_of_component(); /*number_of_component*/
  n = c;   /*Apply c = n+1,n is number of band, the last row of M[][] is f1+f2+...+fn = 1*/
  openmatrix(n,c,1);/*1:open m[n][c] square matrix*/
  openmatrix(n,c,2);/*2:training pixels*/
  menu();
  while ((ch = getchar()) != 'z')
    {
    if (ch != '\n')
      {
      switch(ch)
	{
	case 'a':
                c = number_of_component();/*To get number of component*/
                n = c;
                openmatrix(n,c,1);
                openmatrix(n,c,2);
                menu();/*main menu*/
                break;
        case 'b':
                showmatrix(m,n,c);/*This function is used to show matrix*/
                inverse(m,n,c,im);/*This function is used to calculate inverse matrix of m[][]*/
                for(l=0;l<img_length;l++)
                 for(m=0;m<img_width;m++)
		   {
                   fread(store1,1,1,fin1);
                   fread(store2,1,1,fin2);
                   proportion(fout,n,c,f,store1,store2);/*This function  estimates proporation of components in a pixel*/
		   }
                menu();
                break;
        case 'c':
                for(l=0;l<img_length;l++)
                 for(m=0;m<img_length;m++)
		   {
                   fread(store1,1,1,fin1);
                   fread(store2,1,1,fin2);
                   estimator_11(fout,n,c,f_hat,store1,store2);/*output f_hat*/
                   /*showmatrix(f_hat,n,1);*/ /*estimator (11)*/
		   }
                menu();
                break;
        default:
                printf("please enter the correct alphabet!\n");
                menu();
                break;

	}
      }
    }
    printf("\n");

    fclose(fin1);
    fclose(fin2);
    free(store1);
    free(store2);


   }
  else
      printf("\nUsage:\n\tLinear_mix   'band1 file' 'band2 file' 'rows in file' 'columns' 'output file' .\n\n");

}


menu()/*The choise of user*/
{

  printf("\n");
  printf("Please make a choice:\n");
  printf("Renew again                                                               --->a \n");
  printf("Calculate the proportions of each component using estimator (3),No Noise.  --->b \n");
  printf("Calculate the proportions of each component using estimator (11)           --->c \n");
  printf("Exit the program                                                   --->z \n");
  printf("\n");
}

/*next three function is used to calculate estimator(11)*/
estimator_11(fout1,nn,cc,part_1_2_3,s1,s2)
FILE *fout1;
int   nn,cc;
float part_1_2_3[16][16];
char s1[2],s2[2];
{
   int   i,k;
   float a,j[16][16];
   float u[16][16],tm[16][16],in[16][16],r_iu[16][16],iu[16][16],multi_uj[16][16];
   float r_multi_ujtju[16][16],multi_ujtju[16][16],minus_u_auju[16][16],x_dn[16][16];
   float r_multi_tminx[16][16],multi_tminx[16][16],part_2_3[16][16];

   /*next four functions are used to calculate U=inverse(tm*in*m);*/
   transpose(m,nn,cc,tm);
/*==================================1=====================================================
              showmatrix(tm,nn,cc);
=======================================================================================*/

   inverse(n_matrix,nn,cc,in);
/*===================================2====================================================
              showmatrix(in,nn,cc);
=======================================================================================*/

   multiply(tm,in,nn,cc,cc,r_iu);
/*===================================3====================================================
              showmatrix(r_iu,nn,cc);
=======================================================================================*/

   multiply(r_iu,m,nn,cc,cc,iu);
/*===================================4====================================================
              showmatrix(iu,nn,cc);
=======================================================================================*/

   inverse(iu,nn,cc,u);/*output U*/
/*==================================5=====================================================
              showmatrix(u,nn,cc);
=======================================================================================*/

   a = 0;
   for(i=1;i <= nn;i++)
      for(k=1;k <= cc;k++)
	{
         a = a + u[i][k];
         j[i][k] = 1;
	}
/*==================================5=====================================================
             printf("the number a is %f",a);
=======================================================================================*/

   multiply(u,j,nn,cc,1,multi_uj);/*calculate Uj then output multi_uj*/

   if(a !=0)
     for(i=1;i <= nn;i++)
       multi_uj[i][1] = multi_uj[i][1]/a;/*part1:aUj*/
   else
     for(i=1;i <= nn;i++)
       multi_uj[i][1] = multi_uj[i][1]*a;/*part1:aUj*/
/*=================================6======================================================
              showmatrix(multi_uj,nn,1);
=======================================================================================*/

   multiply(u,j,nn,cc,cc,r_multi_ujtju);/*calculate UJ*/
/*================================7=======================================================
              showmatrix(r_multi_ujtju,nn,cc);
=======================================================================================*/

   multiply(r_multi_ujtju,u,nn,cc,cc,multi_ujtju);/*calculate UJU*/
/*================================8=======================================================
              showmatrix(multi_ujtju,nn,cc);
=======================================================================================*/

   if(a != 0)
    for(i=1;i <= nn;i++)
     for(k=1;k <= cc;k++)
        multi_ujtju[i][k] = multi_ujtju[i][k]/a;/*calculate aUJU*/
   else
     for(i=1;i <= nn;i++)
      for(k=1;k <= cc;k++)
        multi_ujtju[i][k] = multi_ujtju[i][k]*a;/*calculate aUJU*/
/*================================9=======================================================
              showmatrix(multi_ujtju,nn,cc);
=======================================================================================*/

   minusmatrix(u,multi_ujtju,nn,cc,minus_u_auju);/*part2:(U - aUJU)*/
/*================================10=======================================================
              showmatrix(minus_u_auju,nn,cc);
=======================================================================================*/

   multiply(tm,in,nn,cc,cc,r_multi_tminx);
   get_x(nn,s1,s2,x_dn);
   multiply(r_multi_tminx,x_dn,nn,cc,1,multi_tminx);/*part3:tMiNx*/

   multiply(minus_u_auju,multi_tminx,nn,cc,1,part_2_3);/*multiply part2, part3*/
   addmatrix(multi_uj,part_2_3,nn,1,part_1_2_3);/*add part1 to part2*part3,finished*/
   for(i=1; i <= nn ; i++)
    fprintf(fout1,"%f ",part_1_2_3[i][1]);
}


get_x(row,ss1,ss2,xx)
int row;
char ss1[2],ss2[2];
float xx[16][16];
{
  int i;

/*  for(i=1;i <= row-1;i++)
    {
     printf("please input pixel DN in band[%d]:",i);
     scanf("%f",&xx[i][1]);
    }
    */
  xx[0][1] = ss1[0];
  xx[1][1] = ss2[0];
  xx[2][1] = 1;
/*  printf("Input finished, thank you!\n");
  printf("\n");*/

}





int number_of_component()/*To read number of component*/
{
  int cc;

  printf("please input number of components(ground objects):");
  scanf("%d",&cc);
  return(cc);
}


int number_of_training()/*To read number of band*/
{
  int tt;

  printf("please input number of training pixels:");
  scanf("%d",&tt);
  return(tt);
}

/*This function is to show matrix*/
showmatrix(z,nn,cc)
float z[16][16];
int   nn,cc;
{
  int i,j;
  printf("\n");
  printf("The matrix is:!\n");
  for(i=1; i <= nn; i++)
    for(j=1; j <= cc; j++)
      {
        printf("%15f",z[i][j]);
        if (j == cc)
           printf("\n");
      }
}

/*This function is used to open component number m[][],x[][] and e[][]*/
openmatrix(row,col,ss)/*ss = 1 for opening m[][],2 for get M_hat and N_hat*/
int row,col,ss;/*row is number of band+1, col is number of components*/
{
  int    i,j,k,tt;
  float  x_dn[16][16],record_of_x[32][32],f[16][16],tf[16][16],sum_of_etf[16][16],mf[16][16],record_of_f[32][32];
  float  e[16][16],etf[16][16],ftf[16][16],sum_of_ftf[16][16],inver_ftf[16][16],x_minus_mhatf[16][16];
  float  etf_ftf[16][16],mhatf[16][16],sum_cx_x_mhatf[16][16],cx_minus_mhatf[16][16],multi_cx_x[16][16];
  float  t_x_minus_mhatf[16][16],pro;

  switch(ss)
    {
    case 1:
          printf("The number of bands is %d , and The number of components is %d\n",row-1,col);
          printf("\n");
          for(i=1;i <= row-1; i++)
            {
              printf("Band[%d]:please input End-Member DN:\n",i);
              for(j=1; j <= col; j++)
        	{
                 printf("member[%d]: ",j);
                 scanf("%f",&m[i][j]);
                }
              printf("Band[%d] End-Member finished!\n",i);
              printf("\n");
            }
          for(j=1; j <= col; j++)
              m[row][j] = 1;/*setup last line*/
          for(i=1; i <= row; i++)
              for(j=1; j <= col; j++)
                {
                 printf(" %15f ",m[i][j]);
                 if (j == col)
                    printf("\n");
                }

          printf("Above is your End-Member matrix!\n");
          printf("\n");
          break;
    case 2:
          tt = number_of_training();/*number of training pixels*/
          cleanmatrix(sum_of_ftf,row,col);
          cleanmatrix(sum_of_etf,row,col);
          for(i=1;i <= tt;i++)/*begin training from training pixel 1 to tt*/
	    {
              for(j=1;j <= row-1;j++)
		{
                 printf("please input observed DN of training pixel[%d] in band[%d]:",i,j);
                 scanf("%f",&x_dn[j][1]);
                 record_of_x[j][i] = x_dn[j][1];
                }
              printf("please input proportion of each component in training pixel[%d]\n",i);
              pro = 0;
              while( pro != 1)
		{
                for(j=1;j <= col;j++)
	   	 {
                  printf("proportion of component[%d]:",j);
                  scanf("%f",&f[j][1]);
                  pro = pro + f[j][1];
                  record_of_f[j][i] = f[j][1];/*recording f*/
                 }
                if(pro != 1)
		  {
                  printf("Total proportion is not equal 1, please input again\n");
                  pro = 0;
		  }
		}
              x_dn[row][1] = 1;/*setup last line*/
              record_of_x[row][1] = 1;
              multiply(m,f,row,col,1,mf);/*comput m*f then output mf*/
              minusmatrix(x_dn,mf,row,1,e);/*x_dn minus mf then output e*/
              transpose(f,row,1,tf);/*transpose f[][] then output tf*/
              multiply(e,tf,row,1,col,etf);/*calculate e*tf then output etf*/
              multiply(f,tf,row,1,col,ftf);/*calculate f*tf then output ftf*/
              addmatrix(sum_of_etf,etf,row,col,sum_of_etf);/*calculate sum of etf*/
	      addmatrix(sum_of_ftf,ftf,row,col,sum_of_ftf);/*calculate sum of ftf*/
	    }
/*=======================================================================================*/
              showmatrix(x_dn,row,1);
              showmatrix(e,row,1);
              showmatrix(etf,row,col);
              showmatrix(ftf,row,col);
              showmatrix(sum_of_etf,row,col);
              showmatrix(sum_of_ftf,row,col);
/*=======================================================================================*/

          for(i=1;i <= row;i++)
             for(j = 1;j <= col;j++)
	       {
                sum_of_etf[i][j] = sum_of_etf[i][j]/tt;/*claculate average:<etf>*/
                sum_of_ftf[i][j] = sum_of_ftf[i][j]/tt;/*calculate average:<ftf>*/
	       }
/*=======================================================================================*/
              showmatrix(sum_of_etf,row,col);
              showmatrix(sum_of_ftf,row,col);
/*=======================================================================================*/

          inverse(sum_of_ftf,row,col,inver_ftf);/*inverse ftf then output inver_ftf*/

/*=======================================================================================*/
              showmatrix(inver_ftf,row,col);
/*=======================================================================================*/

          multiply(sum_of_etf,inver_ftf,row,col,col,etf_ftf);/*calculate sum_of_etf*inver_ftf then output ef_ftf*/

/*=======================================================================================*/
              showmatrix(etf_ftf,row,col);
/*=======================================================================================*/
          cleanmatrix(m_hat,row,col);
          addmatrix(m,etf_ftf,row,col,m_hat);/*get unbias estimated m_hat[][]*/


          printf("\n");
          for(i=1; i <= row; i++)
              for(j=1; j <= col; j++)
                {
                 printf(" %15f ",m_hat[i][j]);
                 if (j == col)
                    printf("\n");
                }
          printf("Above is your unbias estimated End-Member matrix!\n");
          printf("\n");

          /*============next to calculate n_matrix[][]==============*/

          cleanmatrix(sum_cx_x_mhatf,row,col);
          for(i=1;i <= tt;i++)
	    {
             multiply(m_hat,record_of_f,row,col,1,mhatf);/*calculate mhat*f then output mhatf*/
             for(j=1;j <= row;j++)
	       {
                record_of_f[j][1] = record_of_f[j][i+1];
                x_dn[j][1] = record_of_x[j][1];
                record_of_x[j][1] = record_of_x[j][i+1];
               }
             minusmatrix(x_dn,mhatf,row,1,x_minus_mhatf);/*calculate x - mhatf then output x_minus_mhatf*/
             for(j=1;j <= row;j++)
                x_dn[j][1] = col*x_dn[j][1];/*calculate cx*/
             minusmatrix(x_dn,mhatf,row,1,cx_minus_mhatf);/*calculate cx - mhatf then output cx_minus_mhatf*/
             transpose(x_minus_mhatf,row,col,t_x_minus_mhatf);/*transpose x_minus_mhatf then output t_..*/
             multiply(cx_minus_mhatf,t_x_minus_mhatf,row,1,col,multi_cx_x);
             addmatrix(sum_cx_x_mhatf,multi_cx_x,row,col,sum_cx_x_mhatf);
            }
          for(i=1;i <= row;i++)
            for(j=1;j <= col;j++)
               n_matrix[i][j] = (sum_cx_x_mhatf[i][j]/tt)*(k/(k-col));/*claculate (k/(k-c))*<(cx - mhatf)*t(x - mhatf)>*/

          for(i=1; i <= row; i++)
             for(j=1; j <= col; j++)
                {
                 printf(" %15f ",n_matrix[i][j]);
                 if (j == col)
                    printf("\n");
                }
          printf("Above is your unbias estimated N matrix!\n");
          printf("\n");
          break;
    }

  printf("Input finished! Thanks for your data!\n");
  printf("\n");

}

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

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


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

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

}/*end of transpose*/

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


  for(i=1; i <= row; i++)
     for(j=1; 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*/
float z1[16][16],z2[16][16],item[16][16];
int   row,col;
{
 int   i,j;


  for(i=1; i <= row; i++)
     for(j=1; 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]*/
float z1[16][16],z2[16][16],item[16][16];
int row,col1,col2;
{
  int i,j,k;
  float sum;

  sum = 0;
  cleanmatrix(item,row,col2);
  for(i=1; i <= row; i++)
    {
     for(j=1; j <= col2; j++)
       {
        for(k=1; 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)
float z[16][16],item[16][16];
int nn,cc;
{
  int i,j,k,p,q,r,s;
  float l[16][16],u[16][16];
  float append,determ,sum;
  float il[16][16],iu[16][16];/*il[][] and iu[][] is the inverse matrix of l[][] and u[][]*/


  for(i=1; i <= nn; i++)/*setup the elements of u[][] and l[][]*/
    for(j=1; j <= cc; j++)
      {
        if(j < i)/*lower triangle elements = 0*/
          u[i][j] = 0;
        else
	  {
            u[i][j] = z[i][j];
            for(k=1; 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=1 ; 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=1 ; i <= nn; i++)
    for(j=1; j <= cc; j++)
      {
        il[i][j] = 0;
        iu[i][j] = 0;
      }

  for(i=1; 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=1; i < cc; 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=1; 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=1; i <= nn; i++)
    {
     for(j=1; j <= cc; j++)
       {
        for(k=1; 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()*/


/*Next function is used to estimate the proporation of each element in a pixel*/
proportion(fout1,nn,cc,ff,s1,s2)
FILE *fout1;
char s1[2],s2[2];
float ff[16][16];
int nn,cc;
{
  int i,j;
  float item,x_dn[16][16],iz[16][16];

  get_x(nn,s1,s2,x_dn);/*input x DN*/
  inverse(m,nn,cc,iz);
  /*showmatrix(iz,nn,cc);*/
  item = 0;
  for(i=1; i <= nn ; i++)
    {
     for(j=1; j <= cc ; j++)
       {
        item = item + iz[i][j]*x_dn[j][1];
       }
     ff[i][1] = item;
     item = 0;
    }
  printf("The proportion of each element is:\n");
  for(i=1; i <= nn ; i++)
    fprintf(fout1,"%f ",ff[i][1]);
  /*printf("\n");*/
}