/*	mar.c
**	    This program is used to calculate
**      texture image parameters based on  the
**      'multiplicative autoregressive rendom field' Technique.
**======IEEE Trans.Gs&RS,vol.GE-25,no.2====*
**======pp.195-207, 1987===================*
**======Lognormal Random-Field Models and==*
**======Their Applications to Radar Image==*
**====== Synthesis.========================*
**
**      This programme outputs three files storing
**      Sida, covarience, and mean value (all in floating format),
**      respectively.
**      User may use boundary_f_img.c to convert each file
**      into an texture image (mapped to 0-255).
**      COPYRIGHT: Brandt Tso
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define window_size_max 30
#define num_of_n 3 /*Three model supports defined as (0,-1),(-1,-1),(-1,0)*/
#define image_size 1024

unsigned char image[window_size_max][image_size];
double window[window_size_max][window_size_max];




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

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


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

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

}/*end of transpose*/


void multiply(z1,z2,row,col1,col2,item)/*multiply matrices z1[row][col1] and z2[col1][col2]*/
double z1[window_size_max][window_size_max],z2[window_size_max][window_size_max],item[window_size_max][window_size_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*/


void inverse(magic,nn,cc,item)
double magic[window_size_max][window_size_max],item[window_size_max][window_size_max];
int nn,cc;
{
  int i,j,k;
  double l[window_size_max][window_size_max],u[window_size_max][window_size_max];
  double sum;
  double il[window_size_max][window_size_max],iu[window_size_max][window_size_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] = magic[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] = magic[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=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; 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 is to show matrix*/
void showmatrix(magic,nn,cc)
double magic[window_size_max][window_size_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("%lf ",magic[i][j]);
        if ((j+1) == cc)
           printf("\n");
      }
}




void main(void)
{
  FILE 		*fin, *fout[12];
  unsigned long int img_row, img_col;
  unsigned long int count[10];
  unsigned long int site;
  double z[window_size_max][window_size_max],sida[window_size_max][window_size_max],temp_z[window_size_max][window_size_max];
  double inv_temp_z[window_size_max][window_size_max],temp_sida[window_size_max][window_size_max], covar,p1=0.5, p2=2.0;
  double trans_sida[window_size_max][window_size_max],lnsum,lnmean, avg_sida;
  double sida_z[window_size_max][window_size_max], luck;
  unsigned long int win_size;
  unsigned char in_file[50];
  fpos_t img_start;



      printf("This programme processes an image and output three files\n");
	  printf("known as MAR average sida, covarience, and mean, all in floating values\n");
	  printf("If user want to view how these values look like, you can use the programme\n");
	  printf("called boundary_f_img.c to translate these floating files into images\n\n");

	  printf("Please input the number of image rows:");
      scanf("%d",&img_row);
	  printf("Please input the number of image columns:");
      scanf("%d",&img_col);

	  printf("Please input the window size (3, or 5, or 7,...etc.):");
	  scanf("%d",&win_size);

      printf("Please input image for proccessing:");
	  scanf("%s",in_file);

      if((fin = fopen(in_file,"rb")) == NULL)
		{
		  fprintf(stderr,"Failed to open %s\n",in_file);
		  exit(EXIT_FAILURE);
		}

      printf("Please input file name for storing average sida values:");
	  scanf("%s",in_file);
      if((fout[1] = fopen(in_file,"w")) == NULL)
		{
			fprintf(stderr,"Failed to open %s\n",in_file);
			exit(EXIT_FAILURE);
		}
	  printf("Please input file name for storing covarience values:");
	  scanf("%s",in_file);
      if((fout[2] = fopen(in_file,"w")) == NULL)
		{
			fprintf(stderr,"Failed to open %s\n",in_file);
			exit(EXIT_FAILURE);
		}
	  printf("Please input file name for storing mean values:");
	  scanf("%s",in_file);
      if((fout[3] = fopen(in_file,"w")) == NULL)
		{
			fprintf(stderr,"Failed to open %s\n",in_file);
			exit(EXIT_FAILURE);
		}


	  fgetpos(fin,&img_start);
      for(count[5] = 0; count[5] < img_row-win_size+1; count[5]++)
	  {
        fsetpos(fin,&img_start);
        for(count[3] = 0; count[3] < win_size; count[3]++)
			{
              fread(image[count[3]],1,img_col,fin);
			  if(count[3] == 0)
				  fgetpos(fin,&img_start);
			}

        for(count[6] = 0; count[6] < img_col-win_size+1; count[6]++)
		{
          site = count[6];
          lnsum = 0;
		  for(count[3] = 0; count[3] < win_size; count[3]++)
			{
			for(count[4] = 0; count[4] < win_size; count[4]++)
				{
				luck = image[count[3]][site+count[4]];
				if(luck == 0) luck = 1;
				window[count[3]][count[4]] = log(luck);
				lnsum = lnsum + window[count[3]][count[4]];
				}
			}

          lnmean = (double) lnsum/(win_size*win_size);
          covar = 0;
          for(count[3] = 0; count[3] < num_of_n; count[3]++)
			{
            for(count[4] = 0; count[4] < num_of_n; count[4]++)
				{
				 temp_sida[count[3]][count[4]] = 0;
				 sida[count[3]][count[4]] = 0;
				 temp_z[count[3]][count[4]] = 0;
				 inv_temp_z[count[3]][count[4]] = 0;
				}
            }

          for(count[3] = 1; count[3] < win_size; count[3]++)
            for(count[4] = 1; count[4] < win_size; count[4]++)
			{
               z[0][0] = window[count[3]][count[4]-1]-lnmean;
               z[1][0] = window[count[3]-1][count[4]-1]-lnmean;
               z[2][0] = window[count[3]-1][count[4]]-lnmean;
               for(count[1]= 0; count[1] < num_of_n; count[1]++)
                 for(count[2] = 0; count[2] < num_of_n; count[2]++)
                    temp_z[count[1]][count[2]] = temp_z[count[1]][count[2]]+z[count[1]][0]*z[count[2]][0];


               for(count[1] = 0; count[1] < num_of_n; count[1]++)
                 temp_sida[count[1]][0] = temp_sida[count[1]][0]+z[count[1]][0]*(window[count[3]][count[4]]-lnmean);
			}



           inverse(temp_z,num_of_n,num_of_n,inv_temp_z);
           multiply(inv_temp_z,temp_sida,num_of_n,num_of_n,1,sida);
           avg_sida = (double)(sida[0][0]+sida[1][0]+sida[2][0])/3;

           for(count[3] = 0; count[3] < win_size; count[3]++)
             for(count[4] = 0; count[4] < win_size; count[4]++)
			   {
               transpose(sida,3,1,trans_sida);
               multiply(trans_sida,z,1,3,1,sida_z);
               covar = covar+(window[count[3]][count[4]]-lnmean-sida_z[0][0]);
				}
           covar = (double) covar/(win_size*win_size);

           fprintf(fout[1],"%lf ",avg_sida);
           fprintf(fout[2],"%lf ",covar);
           fprintf(fout[3],"%lf ",lnmean);


		}/*end of count[6]*/
		printf("Processing row %d\n",count[5]+1);
	  }/*end of count[5]*/




      fclose(fin);
      fclose(fout[1]);
      fclose(fout[2]);
      fclose(fout[3]);


      printf("********************Done******************\n");
      printf("\n");

}