/*	Parameter_es.c
**	This program is used to estimate
**      MRF model-associated parameters beta[1] to beta[4] based on
**      the 'least squares error' technique.
**======IEEE Trans. Pattern Analysis and============*
**======Machine Intelligence, Vol PAMI-9============*
**======pp.39-55, 1987  ============================*
**======Modelling and Segmentation of noisy and=====*
**======Textured Images Using Gibbs Random Fields===*
**Note that the class number must start from 1
**COPYRIGHT: Brandt Tso
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define win_size 3
#define image_size 512
#define class_max 12
#define equation_max 2000

double det[10];
double p,x1[equation_max][10], x2[10][equation_max], x0[equation_max][10],x3[10][10],x4[10][10], x5[10][equation_max],x6[10][10],yy1[equation_max][10], beta[20][10];




/*This function is to show matrix*/
void showmatrix(z,nn,cc)
/* FILE *file1;*/
double z[700][10];
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++)
      {
        /*fprintf(file1,"%15f",z[i][j]); */
        printf("%15lf",z[i][j]);
        if (j == cc-1)
	  {
           printf("\n");

	  }
      }
}



void cleanmatrix(z,nn,cc)/*to clear matrix z's elements*/
double z[700][10];
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[700][10],tz[10][700];
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[16][16],z2[16][16],item[16][16];
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[16][16],z2[16][16],item[16][16];
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 multiply1(z1,z2,row,col1,col2,item)/*multiply matrices z1[row][col1] and z2[col1][col2]*/
double z1[10][700],z2[700][10],item[10][10];
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 multiply2(z1,z2,row,col1,col2,item)/*multiply matrices z1[row][col1] and z2[col1][col2]*/
double z1[10][10],z2[10][700],item[10][700];
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[10][10],item[10][10];
int nn,cc;
{
  int i,j,k;
  double l[10][10],u[10][10];
  double sum;
  double il[10][10],iu[10][10];/*il[][] and iu[][] is the inverse matrix of l[][] and u[][]*/




  for(i=0; i < nn; i++)
    for(j=0; j < cc; j++)
      {
        if(j < i)
          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;
        else if(j == i)
           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++)
    {
      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++)
    {
      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];
	    }

        }
    }

  for(i=0; i < nn; i++)
    {
      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];
	}
    }



  /*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*/

double determ(z,nn,cc,class_num)
double z[16][16];
int nn,cc,class_num;
{
  int i,j,k;
  double l[16][16],u[16][16];
  double item=0;
//  double il[16][16],iu[16][16];/*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];
	  }
      }
   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);
}




void main(void)
{
  FILE 		*fin, *fout;
  long int img_row, img_col;
  long int site1,site2, num_of_class;
  long int count1,count2, count3, count4, count5, count6,count7, i1, j1, k1;
  long int centre,window[win_size][win_size],v[10],v1[10],w[10][100];
  long int set_num,mark[10][100],p_sort1[class_max][100],p_sort2[class_max][100],max, sort_num1, sort_num2;
  unsigned char store[2],in_file[50],**image;
  double par1, par2, par3;


      printf("Parameter estimations, please note the class number must start from 1\n");
	  printf("Please input the total number for texture classes:");
      scanf("%d",&num_of_class);
	  printf("Please input the number of image rows:");
      scanf("%d",&img_row);
	  printf("Please input the number of image columns:");
      scanf("%d",&img_col);


      image =(unsigned char **)calloc(img_row+5,sizeof(unsigned char *));

      for(i1 = 0 ; i1 < img_row+5; i1++)
		{
        image[i1] = (unsigned char *)calloc(img_col+5,sizeof(unsigned char));
		}




      printf("Please enter how many items you want to choose to form equations:");
      scanf("%d",&j1);/*number of items been chosen to form equations*/


	  printf("Please input image file name:");
	  scanf("%s",in_file);
      if((fin = fopen(in_file,"rb")) == NULL)
		{
	     printf("Failed to open %s\n",in_file);
	     exit(EXIT_FAILURE);
		}

      printf("Output file name:");
	  scanf("%s",in_file);
      if((fout = fopen(in_file,"w")) == NULL)
	   {
		fprintf(stderr,"Failed to open %s\n",in_file);
		exit(EXIT_FAILURE);
		}


      for(count1 = 0;count1 < img_row;count1++)
        for(count2 = 0;count2 < img_col;count2++)
		{
          fread(store,1,1,fin);
          image[count1][count2] = store[0];
		}

      centre = (win_size+1)/2-1;
	  for(count2=0;count2<=num_of_class;count2++)
       {
        for(count1=0;count1<81;count1++)
			{
			w[count2][count1]=0;
			mark[count2][count1]=0;
			}
       }

      printf("Start Calculating ......\n");


    for(count7 = 1; count7<= num_of_class; count7++)
      {
      for(count5 = 1; count5 < img_row; count5++)
       {
        for(count6 = 1; count6 < img_col; count6++)
		{
          site1 = count5;
          site2 = count6;

		  for(count3 = 0; count3 < win_size; count3++)
		   for(count4 = 0; count4 < win_size; count4++)
				{
                window[count3][count4] =  image[site1-centre+count3][site2-centre+count4];
				}
          v[0]=0; v[1]=0; v[2]=0; v[3]=0;

          if(count7 == window[centre][centre])
			{
            for(count3 = 0; count3 < win_size; count3++)
              for(count4 = 0; count4 < win_size; count4++)
				{
                 if((count3==centre)&&(count4==centre)) count4++;

				 if(window[centre][centre]== window[count3][count4])
					{
                   if((count3==0 && count4==0) || (count3==2 && count4==2)) v[2]++;
                   if((count3==0 && count4==1) || (count3==2 && count4==1)) v[1]++;
                   if((count3==0 && count4==2) || (count3==2 && count4==0)) v[3]++;
                   if((count3==1 && count4==0) || (count3==1 && count4==2)) v[0]++;
					}
                 else if(window[centre][centre]!= window[count3][count4])
					{
                   if((count3==0 && count4==0) || (count3==2 && count4==2)) v[2]--;
                   if((count3==0 && count4==1) || (count3==2 && count4==1)) v[1]--;
                   if((count3==0 && count4==2) || (count3==2 && count4==0)) v[3]--;
                   if((count3==1 && count4==0) || (count3==1 && count4==2)) v[0]--;
					}
				}

              for(count3 = 0; count3 < 4; count3++)
				{
                 if(v[count3] == 2)
                   v1[count3]=2;
                 else if(v[count3] == 0)
                   v1[count3]= 1;
                 else
                   v1[count3]=0;
				}

              set_num = 3*3*3*v1[0] + 3*3*v1[1] + 3*v1[2] + v1[3];
              x0[set_num][0]= v[0]; x0[set_num][1]= v[1]; x0[set_num][2]= v[2]; x0[set_num][3]= v[3];/*configuration*/
              w[count7][set_num]++;
			}/*end if*/

          }/*end of count6*/

       }/*end of count5*/

      }/*end of count7*/


      for(count7=1;count7<= num_of_class;count7++)
		{
         fprintf(fout,"\n class %d\n",count7);
		 for(count1=0;count1<81;count1++)
			{
			max = 0;
			for(count2=0;count2<81;count2++)
				{
				if((mark[count7][count2]==0)&&(w[count7][count2]>=max))
					{
					max =w[count7][count2];
					count3=count2;
					}
				}

			mark[count7][count3]=1;
			p_sort1[count7][count1]=count3;/*configuration pointer, i.e. 2222, 2221, 0r 0000 etc.*/
			p_sort2[count7][count1]=max; /*number of occurances*/
			p=max;
			if(p !=0)
				{
				p = (double) log(p);
				/*printf(" %d-%d-%lf ",p_sort1[count7][count1],p_sort2[count7][count1],p);*/
				fprintf(fout," %d-%d-%lf ",p_sort1[count7][count1],p_sort2[count7][count1],p);
				}
			}
		}

     for(count1=0; count1 < equation_max; count1++)
           x1[count1][0]=1;

     i1 = 0;

     printf("\nForming equations and starting regression process!\n");

     for(count1=1; count1 <= num_of_class; count1++)
       {
      for(count2=0; count2 < j1; count2++)
		{
        sort_num1 = p_sort1[count1][count2];
		   for(count3=count2+1; count3 < j1; count3++)
			{
          	sort_num2 = p_sort1[count1][count3];
            x1[i1][1]= x0[sort_num1][0]-x0[sort_num2][0];
			x1[i1][2]= x0[sort_num1][1]-x0[sort_num2][1];
			x1[i1][3]= x0[sort_num1][2]-x0[sort_num2][2];
			x1[i1][4]= x0[sort_num1][3]-x0[sort_num2][3];
			par1 = (double) p_sort2[count1][count3];
			if(p_sort2[count1][count2]!=0)
				{
				  par2 = (double) p_sort2[count1][count2];
				  par3 = (double) log(par1)-log(par2);
				}
			else
				  par3 =  (double) log(par1);
			yy1[i1][0] = par3;
			i1++;
			}
	  }
	 }


     for(count1=1; count1 <= num_of_class-1; count1++)
       {
      for(count2=0; count2 < j1; count2++)
		{
        sort_num1 = p_sort1[count1][count2];
        for(count3 = 1; count3 <= num_of_class-count1; count3++)
			{
          for(count4=0; count4 < j1; count4++)
				{
				sort_num2 = p_sort1[count1+count3][count4];

				if(sort_num1 != sort_num2)
					{
					x1[i1][1]= x0[sort_num1][0]-x0[sort_num2][0];
					x1[i1][2]= x0[sort_num1][1]-x0[sort_num2][1];
					x1[i1][3]= x0[sort_num1][2]-x0[sort_num2][2];
					x1[i1][4]= x0[sort_num1][3]-x0[sort_num2][3];
					par1 = (double) p_sort2[count1+count3][count4];
					if(p_sort2[count1][count2]!=0)
					{
					  par2 = (double) p_sort2[count1][count2];
					  par3 = (double) log(par1)-log(par2);
					}
					else
					  par3 =  (double) log(par1);

					yy1[i1][0] = par3;
					i1++;
					}
				}
			}
		}
       }


      printf("\nTotal number of equations = %ld \n",i1);
      fprintf(fout,"\nTotal number of equations = %ld \n",i1);

      k1 = 5;
      transpose(x1,i1,k1,x2);
      multiply1(x2, x1, k1, i1, k1, x3);
      inverse(x3,k1,k1,x4);
      multiply2(x4, x2, k1, k1, i1,x5);
      j1 = 1;
      multiply1(x5,yy1, k1, i1,j1,beta);
      printf("\n");
      showmatrix(beta,k1,j1);
	  printf("The estimated coefficients are:\n");
	  printf("%lf %lf %lf %lf\n",beta[1][0],beta[2][0],beta[3][0],beta[4][0]);
      fprintf(fout,"The estimated coefficients are:\n");
	  fprintf(fout,"%lf %lf %lf %lf\n",beta[1][0],beta[2][0],beta[3][0],beta[4][0]);

	  for(i1 = 0 ; i1 < img_row; i1++)
        {
         free(image[i1]);
        }
      free(image);


	  fclose(fin);
      fclose(fout);
      printf("********************Done******************\n");
      printf("\n");
}