/*	Co-occurrence.c
**	This program is used to calculate
**      Textural feature files based on
**      oc-currence matrix, the measure is
**      from user difined window sizes and
**      pixel grey level!
**
**      Output: five feature files called
**      1)Angular Second Moment,
**      2)Entropy,
**      3)Inverse Difference Moment,
**      4)Contrast,
**      5)Correlation.
**      related to each pixel.

**======IEEE Transactions on Systems, man==*
**======,and Cybernetics. Vol. SMC-3, No.6=*
**======1973, pp.610-621.Textural Features=*
**======for Image Classification===========*
**COPY RIGHT: Brandt Tso
*/

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

#define window_size_r 256
#define gray_num 256
#define image_size 1024

unsigned char **image;
unsigned char window[window_size_r][window_size_r];
int **ph,**pv;
int **pld,**prd;
unsigned int N_r,gray_lev;


void start_initial(void)
  {
unsigned int i,j;

      for(i=0;i<= gray_lev ;i++)
         for(j=0;j <= gray_lev;j++)
               {
               ph[i][j] = 0;
               pv[i][j] = 0;
               pld[i][j] = 0;
               prd[i][j] = 0;
	       }
      for(i=0;i<= window_size_r ;i++)
         for(j=0;j <= window_size_r;j++)
              window[i][j] = 0;
  }

void subtract_ph(w_size,dd)
unsigned int w_size,dd;
  {
unsigned int i,gray1,gray2;

  for(i = 0;i < w_size;i++)
    {
      gray1 = window[i][0];
      gray2 = window[i][0+dd];
      ph[gray1][gray2]--;
      ph[gray2][gray1]--;
     }
  }

void add_ph(w_size,dd)
unsigned int w_size,dd;
  {
unsigned int i,gray1,gray2;

   for(i = 0;i < w_size;i++)
    {
      gray1 = window[i][w_size-1];
      gray2 = window[i][w_size-1-dd];
      ph[gray1][gray2]++;
      ph[gray2][gray1]++;
     }
  }

void subtract_pv(w_size,dd)
unsigned int w_size,dd;
  {
long int i,gray1,gray2;

   i = 0;
   while((i+dd) < w_size)
     {
      gray1 = window[i][0];
      gray2 = window[i+dd][0];
      pv[gray1][gray2]--;
      pv[gray2][gray1]--;
      i++;
     }
  }

void add_pv(w_size,dd)
unsigned int w_size,dd;
  {
long int i,gray1,gray2;

   i = 0;
   while((i+dd) < w_size)
     {
      gray1 = window[i][w_size-1];
      gray2 = window[i+dd][w_size-1];
      pv[gray1][gray2]++;
      pv[gray2][gray1]++;
      i++;
     }
  }

void subtract_pld(w_size,dd)
unsigned int w_size,dd;
  {
long int i,gray1,gray2;

  i = 0;
   while((i+dd) < w_size)
     {
      gray1 = window[i][0];
      gray2 = window[i+dd][0+dd];
      pld[gray1][gray2]--;
      pld[gray2][gray1]--;
      i++;
     }
  }

void add_pld(w_size,dd)
unsigned int w_size,dd;
  {
long int i,gray1,gray2;

   i = 0;
   while((i+dd) < w_size)
     {
      gray1 = window[i][w_size-1-dd];
      gray2 = window[i+dd][w_size-1];
      pld[gray1][gray2]++;
      pld[gray2][gray1]++;
      i++;
     }
  }

void subtract_prd(w_size,dd)
unsigned int w_size,dd;
  {
long int i,gray1,gray2;

  i = 0;
   while((i+dd) < w_size)
     {
      gray1 = window[i][0+dd];
      gray2 = window[i+dd][0];
      prd[gray1][gray2]--;
      prd[gray2][gray1]--;
      i++;
     }
  }

void add_prd(w_size,dd)
unsigned int w_size,dd;
  {
long int i,gray1,gray2;

   i = 0;
   while((i+dd) < w_size)
     {
      gray1 = window[i][w_size-1];
      gray2 = window[i+dd][w_size-1-dd];
      prd[gray1][gray2]++;
      prd[gray2][gray1]++;
      i++;
     }
  }

void cal_ph(w_size,dd)
unsigned int w_size,dd;
{
unsigned int i,j;
int gray_tone1, gray_tone2;


    for(i=0;i<w_size;i++)
      for(j=0;j<w_size;j++)
	{
         if(j == 0)
	   {
            gray_tone1 = window[i][j];
            if((j+dd) < w_size)
	      {
               gray_tone2 = window[i][j+dd];
               ph[gray_tone1][gray_tone2]++;

	      }
           }
         else if(j == w_size - 1)
	   {
            gray_tone1 = window[i][j];
            if(((j-dd) < w_size)&&((j-dd)>=0))
	      {
               gray_tone2 = window[i][j-dd];
               ph[gray_tone1][gray_tone2]++;

	      }
           }
         else
	   {
            gray_tone1 = window[i][j];
            if((j+dd) < w_size)
	      {
               gray_tone2 = window[i][j+dd];
               ph[gray_tone1][gray_tone2]++;

	      }
            if(((j-dd) < w_size)&&((j-dd)>=0))
	      {
               gray_tone2 = window[i][j-dd];
               ph[gray_tone1][gray_tone2]++;

	      }
	   }
	}
/*
   for(i=0;i<3;i++)
      {
     for(j=0;j<3;j++)
         printf("%d",window[i][j]);
     printf("\n");
      }
*/
/*    for(i=0;i<=5;i++)
      {
      for(j=0;j<=5;j++)
        {
         zz = (double) ph[i][j]/(2*127*128);
         printf(" %d ",ph[i][j]);
          printf(" %lf ", zz);
        }
      printf("\n");
      }*/
}

void cal_pv(w_size,dd)
unsigned int w_size,dd;
{
unsigned int i,j;
int gray_tone1, gray_tone2;

    for(i=0;i<w_size;i++)
      for(j=0;j<w_size;j++)
	{
         if(j == 0)
	   {
            gray_tone1 = window[j][i];
            if((j+dd) < w_size)
	      {
               gray_tone2 = window[j+dd][i];
               pv[gray_tone1][gray_tone2]++;
	      }
       }
         else if(j == w_size - 1)
	   {
            gray_tone1 = window[j][i];
            if(((j-dd) < w_size)&&((j-dd)>=0))
	      {
               gray_tone2 = window[j-dd][i];
               pv[gray_tone1][gray_tone2]++;
              }
           }
         else
	   {
            gray_tone1 = window[j][i];
            if((j+dd) < w_size)
	      {
               gray_tone2 = window[j+dd][i];
               pv[gray_tone1][gray_tone2]++;

	      }
             if(((j-dd) < w_size)&&((j-dd)>=0))
	      {
               gray_tone2 = window[j-dd][i];
               pv[gray_tone1][gray_tone2]++;

	      }
	   }
	}
/*     for(i=1;i<=4;i++)
      {
      for(j=1;j<=4;j++)
         printf("%d",pv[i][j]);
      printf("\n");
      }*/
}

void cal_pld(w_size,dd)
unsigned int w_size,dd;
{
unsigned int i,j;
int gray_tone1, gray_tone2;

    for(i=0;i<w_size;i++)
      for(j=0;j<w_size;j++)
	{
         if((j == 0)||(i==0))
	   {
            gray_tone1 = window[i][j];
            if(((j+dd) < w_size)&&((i+dd) < w_size))
	      {
               gray_tone2 = window[i+dd][j+dd];
               pld[gray_tone1][gray_tone2]++;

	      }
           }
         else if((j == w_size - 1)||(i == w_size - 1))
	   {
            gray_tone1 = window[i][j];
            if(((j-dd) < w_size)&&((j-dd)>=0)&&((i-dd) < w_size)&&((i-dd)>=0))
	      {
               gray_tone2 = window[i-dd][j-dd];
               pld[gray_tone1][gray_tone2]++;

	      }
           }
         else
	   {
            gray_tone1 = window[i][j];
            if(((j+dd) < w_size)&&((i+dd) < w_size))
	      {
               gray_tone2 = window[i+dd][j+dd];
               pld[gray_tone1][gray_tone2]++;

	      }
            if((((j-dd) < w_size)&&((j-dd)>=0))&&(((i-dd) < w_size)&&((i-dd)>=0)))
	      {
               gray_tone2 = window[i-dd][j-dd];
               pld[gray_tone1][gray_tone2]++;

	      }
	    }
	}
/*     for(i=1;i<=4;i++)
      {
      for(j=1;j<=4;j++)
         printf("%d",pld[i][j]);
      printf("\n");
      }*/
}


void cal_prd(w_size,dd)
unsigned int w_size,dd;
{
unsigned int i,j;
int gray_tone1, gray_tone2;

    for(i=0;i<w_size;i++)
      for(j=0;j<w_size;j++)
	{
         if((i == 0)||(j == w_size-1))
	   {
            gray_tone1 = window[i][j];
            if((((j-dd) < w_size)&&((j-dd) >= 0))&&((i+dd) < w_size))
	      {
               gray_tone2 = window[i+dd][j-dd];
               prd[gray_tone1][gray_tone2]++;

	      }
           }
         else if((i == w_size - 1)||(j == 0))
	   {
            gray_tone1 = window[i][j];
            if(((j+dd) < w_size)&&(((i-dd) < w_size)&&((i-dd)>=0)))
	      {
               gray_tone2 = window[i-dd][j+dd];
               prd[gray_tone1][gray_tone2]++;

	      }
           }
         else
	   {
            gray_tone1 = window[i][j];
            if((((j-dd) < w_size)&&((j-dd) >= 0))&&((i+dd) < w_size))
	      {
               gray_tone2 = window[i+dd][j-dd];
               prd[gray_tone1][gray_tone2]++;

	      }
             if(((j+dd) < w_size)&&(((i-dd) < w_size)&&((i-dd)>=0)))
	      {
               gray_tone2 = window[i-dd][j+dd];
               prd[gray_tone1][gray_tone2]++;
	      }
	   }
	}
/*    for(i=1;i<=4;i++)
      {
      for(j=1;j<=4;j++)
         printf("%d",prd[i][j]);
      printf("\n");
      }*/
}


double cal_asm(w_size)
unsigned int w_size;
{
unsigned int i,j;
double final_asm,asm_h,asm_v,asm_ld,asm_rd;
double temp1,temp2,temp3,temp4,pow_num;


    pow_num = 2;

    final_asm=0;
    asm_h=0;
    asm_v=0;
    asm_ld=0;
    asm_rd=0;

    for(i = 0;i < gray_lev;i++)
     for(j = 0;j < gray_lev;j++)
       {
        if(ph[i][j] != 0)
	  {
           temp1 = (double) ph[i][j]/(2*w_size*(w_size-1));
           asm_h = asm_h + pow(temp1,pow_num);
	  }
        if(pv[i][j] != 0)
	  {
           temp2 = (double) pv[i][j]/(2*w_size*(w_size-1));
           asm_v = asm_v + pow(temp2,pow_num);
	  }
        if(pld[i][j] != 0)
	  {
           temp3 = (double) pld[i][j]/(2*(w_size-1)*(w_size-1));
           asm_ld = asm_ld + pow(temp3,pow_num);
	  }
        if(prd[i][j] != 0)
	  {
           temp4 = (double) prd[i][j]/(2*(w_size-1)*(w_size-1));
           asm_rd = asm_rd + pow(temp4,pow_num);
	  }
       }
 /*printf("\n asm %lf %lf %lf %lf\n",asm_h,asm_v,asm_ld,asm_rd);*/

     final_asm = asm_h+asm_v+asm_ld+asm_rd;
     return(final_asm);
}

double cal_ent(w_size)
unsigned int w_size;
{
unsigned int i,j;
double final_ent,ent_h,ent_v,ent_ld,ent_rd;
double temp1,temp2,temp3,temp4;


    final_ent=0;
    ent_h=0;
    ent_v=0;
    ent_ld=0;
    ent_rd=0;

    for(i = 0;i < gray_lev;i++)
     for(j = 0;j < gray_lev;j++)
       {
        if(ph[i][j] != 0)
	  {
           temp1 = (double) ph[i][j]/(2*w_size*(w_size-1));
           ent_h = ent_h - temp1*log(temp1);
	  }
        if(pv[i][j] != 0)
	  {
           temp2 = (double) pv[i][j]/(2*w_size*(w_size-1));
           ent_v = ent_v - temp2*log(temp2);
	  }
        if(pld[i][j] != 0)
	  {
           temp3 = (double) pld[i][j]/(2*(w_size-1)*(w_size-1));
           ent_ld = ent_ld - temp3*log(temp3);
	  }
        if(prd[i][j] != 0)
	  {
           temp4 = (double) prd[i][j]/(2*(w_size-1)*(w_size-1));
           ent_rd = ent_rd - temp4*log(temp4);
	  }
       }
 /*printf("ent %lf %lf %lf %lf\n",ent_h,ent_v,ent_ld,ent_rd);*/

     final_ent = ent_h+ent_v+ent_ld+ent_rd;
     return(final_ent);
}

double cal_idm(w_size)
unsigned int w_size;
{
unsigned int i,j;
double final_idm,idm_h,idm_v,idm_ld,idm_rd;
double temp1,temp2,temp3,temp4;
double ii,jj,pow_num,t0,t1;

    pow_num = 2;

    final_idm=0;
    idm_h=0;
    idm_v=0;
    idm_ld=0;
    idm_rd=0;

    for(i = 0;i < gray_lev;i++)
     for(j = 0;j < gray_lev;j++)
       {
        ii = i; jj = j;
        if(ii > jj)
           t0 = ii - jj;
        else
           t0 = jj - ii;

        if(t0 != 0)
           t1 = (double) 1/(1+pow(t0,pow_num));
        else
           t1 = 1;

        if(ph[i][j] != 0)
	  {
           temp1 = (double) ph[i][j]/(2*w_size*(w_size-1));
           idm_h = idm_h + t1*temp1;
	  }
        if(pv[i][j] != 0)
	  {
           temp2 = (double) pv[i][j]/(2*w_size*(w_size-1));
           idm_v = idm_v + t1*temp2;
	  }
        if(pld[i][j] != 0)
	  {
           temp3 = (double) pld[i][j]/(2*(w_size-1)*(w_size-1));
           idm_ld = idm_ld + t1*temp3;
	  }
        if(prd[i][j] != 0)
	  {
           temp4 = (double) prd[i][j]/(2*(w_size-1)*(w_size-1));
           idm_rd = idm_rd + t1*temp4;
	  }
       }
	/*printf("idm %lf %lf %lf %lf\n",idm_h,idm_v,idm_ld,idm_rd); */

     final_idm = idm_h+idm_v+idm_ld+idm_rd;
     return(final_idm);
}

double cal_con(w_size)
unsigned int w_size;
{
unsigned int i,j;
unsigned int n;
double final_con,con_h,con_v,con_ld,con_rd;
double temp1,temp2,temp3,temp4;
double t1,t2,t3,t4;


    final_con=0;
    con_h=0;
    con_v=0;
    con_ld=0;
    con_rd=0;
    t1=0; t2=0; t3=0; t4=0;

    for(n = 1;n < gray_lev;n++)
      {
      for(i = 0;i < gray_lev;i++)
       for(j = 0;j < gray_lev;j++)
        {
        if((i-j)==n || (j-i) == n)
	  {
           if(ph[i][j] != 0)
	     {
              temp1 = (double) ph[i][j]/(2*w_size*(w_size-1));
              t1 = t1 + temp1;
	     }
           if(pv[i][j] != 0)
	     {
              temp2 = (double) pv[i][j]/(2*w_size*(w_size-1));
              t2 = t2 + temp2;
	     }
           if(pld[i][j] != 0)
	     {
              temp3 = (double) pld[i][j]/(2*(w_size-1)*(w_size-1));
              t3 = t3 + temp3;
	     }
           if(prd[i][j] != 0)
	     {
              temp4 = (double) prd[i][j]/(2*(w_size-1)*(w_size-1));
              t4 = t4 + temp4;
	     }
	  }
        }
       con_h = con_h + n*n*t1;
       con_v = con_v + n*n*t2;
       con_ld = con_ld + n*n*t3;
       con_rd = con_rd + n*n*t4;
       t1=0; t2=0; t3=0; t4=0;
      }
      /*printf("con %lf %lf %lf %lf\n",con_h,con_v,con_ld,con_rd); */

     final_con = con_h+con_v+con_ld+con_rd;
     return(final_con);
}


double cal_crr(w_size)
unsigned int w_size;
{
unsigned int i,j;
double pow1 = 2, pow2 = 0.5, temp1, temp2, temp3, temp4, temp5, temp6, temp7, temp8;
double sumhx,sumhy,meanhx,meanhy,sdvhx,sdvhy;
double sumvx,sumvy,meanvx,meanvy,sdvvx,sdvvy;
double sumldx,sumldy,meanldx,meanldy,sdvldx,sdvldy;
double sumrdx,sumrdy,meanrdx,meanrdy,sdvrdx,sdvrdy;
double final_cor,cor_h,cor_v,cor_ld,cor_rd;

    sumhx = 0;
    sumvx = 0;
    sumldx = 0;
    sumrdx = 0;
    sumhy = 0;
    sumvy = 0;
    sumldy = 0;
    sumrdy = 0;
    sdvhx = 0;
    sdvvx = 0;
    sdvldx = 0;
    sdvrdx = 0;
    sdvhy = 0;
    sdvvy = 0;
    sdvldy = 0;
    sdvrdy = 0;
    final_cor=0;
    cor_h=0;
    cor_v=0;
    cor_ld=0;
    cor_rd=0;


    for(i = 0;i < gray_lev;i++)
     for(j = 0;j < gray_lev;j++)
       {
        temp1 = (double) ph[i][j]/(2*(w_size-1)*w_size);
        sumhx = sumhx + temp1;
        temp2 = (double) pv[i][j]/(2*(w_size-1)*w_size);
        sumvx = sumvx + temp2;
        temp3 = (double) pld[i][j]/(2*(w_size-1)*(w_size-1));
        sumldx = sumldx + temp3;
        temp4 = (double) prd[i][j]/(2*(w_size-1)*(w_size-1));
        sumrdx = sumrdx + temp4;
        temp5 = (double) ph[j][i]/(2*(w_size-1)*w_size);
        sumhy = sumhy + temp5;
        temp6 = (double) pv[j][i]/(2*(w_size-1)*w_size);
        sumvy = sumvy + temp6;
        temp7 = (double) pld[j][i]/(2*(w_size-1)*(w_size-1));
        sumldy = sumldy + temp7;
        temp8 = (double) prd[j][i]/(2*(w_size-1)*(w_size-1));
        sumrdy = sumrdy + temp8;
       }
   /*printf("\n sum %lf %lf %lf %lf %lf %lf %lf %lf\n", sumhx, sumvx, sumldx, sumrdx, sumhy, sumvy, sumldy, sumrdy);*/
    meanhx = (double) sumhx/(gray_lev*gray_lev);
    meanvx = (double) sumvx/(gray_lev*gray_lev);
    meanldx = (double) sumldx/(gray_lev*gray_lev);
    meanrdx = (double) sumrdx/(gray_lev*gray_lev);
    meanhy = (double) sumhy/(gray_lev*gray_lev);
    meanvy = (double) sumvy/(gray_lev*gray_lev);
    meanldy = (double) sumldy/(gray_lev*gray_lev);
    meanrdy = (double) sumrdy/(gray_lev*gray_lev);

    /*printf("\n mean %lf %lf %lf %lf %lf %lf %lf %lf\n", meanhx, meanvx, meanldx, meanrdx, meanhy, meanvy, meanldy, meanrdy);*/

    for(i = 0;i < gray_lev;i++)
     for(j = 0;j < gray_lev;j++)
       {
        temp1 = (double) ph[i][j]/(2*(w_size-1)*w_size);
        sdvhx = sdvhx + pow(temp1-meanhx,pow1);
        temp2 = (double) pv[i][j]/(2*(w_size-1)*w_size);
        sdvvx = sdvvx + pow(temp2-meanhx,pow1);
        temp3 = (double) pld[i][j]/(2*(w_size-1)*(w_size-1));
        sdvldx = sdvldx + pow(temp3-meanhx,pow1);
        temp4 = (double) prd[i][j]/(2*(w_size-1)*(w_size-1));
        sdvrdx = sdvrdx + pow(temp4-meanhx,pow1);
        temp5 = (double) ph[j][i]/(2*(w_size-1)*w_size);
        sdvhy = sdvhy + pow(temp5-meanhx,pow1);
        temp6 = (double) pv[j][i]/(2*(w_size-1)*w_size);
        sdvvy = sdvvy + pow(temp6-meanhx,pow1);
        temp7 = (double) pld[j][i]/(2*(w_size-1)*(w_size-1));
        sdvldy = sdvldy + pow(temp7-meanhx,pow1);
        temp8 = (double) prd[j][i]/(2*(w_size-1)*(w_size-1));
        sdvrdy = sdvrdy + pow(temp8-meanhx,pow1);
       }
       /*printf("\n sdv %lf %lf %lf %lf %lf %lf %lf %lf\n", sdvhx, sdvvx, sdvldx, sdvrdx, sdvhy, sdvvy, sdvldy, sdvrdy);   */
        sdvhx = (double) sdvhx/(gray_lev*gray_lev);
        sdvvx = (double) sdvvx/(gray_lev*gray_lev);
        sdvldx = (double) sdvldx/(gray_lev*gray_lev);
        sdvrdx = (double) sdvrdx/(gray_lev*gray_lev);
        sdvhy = (double) sdvhy/(gray_lev*gray_lev);
        sdvvy = (double) sdvvy/(gray_lev*gray_lev);
        sdvldy = (double) sdvldy/(gray_lev*gray_lev);
        sdvrdy = (double) sdvrdy/(gray_lev*gray_lev);

        sdvhx = pow(sdvhx,pow2);
        sdvvx = pow(sdvvx,pow2);
        sdvldx = pow(sdvldx,pow2);
        sdvrdx = pow(sdvrdx,pow2);
        sdvhy = pow(sdvhy,pow2);
        sdvvy = pow(sdvvy,pow2);
        sdvldy = pow(sdvldy,pow2);
        sdvrdy = pow(sdvrdy,pow2);


    for(i = 0;i < gray_lev;i++)
     for(j = 0;j < gray_lev;j++)
       {
         temp1 = (double) ph[i][j]/(2*(w_size-1)*w_size);
	 cor_h = cor_h + i*j*temp1;
         temp2 = (double)  pv[i][j]/(2*(w_size-1)*w_size);
	 cor_v = cor_v + (double) temp2;
        temp3 = (double) pld[i][j]/(2*(w_size-1)*(w_size-1));
	 cor_ld = cor_ld +(double)i*j*temp3;
         temp4 = (double) prd[i][j]/(2*(w_size-1)*(w_size-1));
	 cor_rd = cor_rd + i*j*temp4;
       }

    cor_h = (double) (cor_h-meanhx*meanhy)/(sdvhx*sdvhy);
    cor_v = (double) (cor_v-meanhx*meanhy)/(sdvvx*sdvvy);
    cor_ld = (double) (cor_ld-meanhx*meanhy)/(sdvldx*sdvldy);
    cor_rd = (double) (cor_rd-meanhx*meanhy)/(sdvrdx*sdvrdy);

    /*printf(" cor %lf \n", cor_h); */
     final_cor = cor_h+cor_v+cor_ld+cor_rd;
     return(final_cor);
}

/*
double cal_f10(w_size)
unsigned int w_size;
{
unsigned int i,j;
double final_f10,f10_h,f10_v,f10_ld,f10_rd;
double temp1,temp2,temp3,temp4;
double t1,t2,t3,t4,m1,m2,m3,m4,pow1 = 2;


    final_f10=0;
    f10_h=0;
    f10_v=0;
    f10_ld=0;
    f10_rd=0;
    t1=0; t2=0; t3=0; t4=0;


      for(i = 0;i < gray_lev;i++)
       for(j = 0;j < gray_lev;j++)
        {
        if((i-j)==1 || (i-j) == -1)
	  {
           if(ph[i][j] != 0)
	     {
              temp1 = (double) ph[i][j]/(2*w_size*(w_size-1));
              t1 = t1 + temp1;
	     }
           if(pv[i][j] != 0)
	     {
              temp2 = (double) pv[i][j]/(2*w_size*(w_size-1));
              t2 = t2 + temp2;
	     }
           if(pld[i][j] != 0)
	     {
              temp3 = (double) pld[i][j]/(2*(w_size-1)*(w_size-1));
              t3 = t3 + temp3;
	     }
           if(prd[i][j] != 0)
	     {
              temp4 = (double) prd[i][j]/(2*(w_size-1)*(w_size-1));
              t4 = t4 + temp4;
	     }
	  }
        }

       m1 = (double) t1/(gray_lev*gray_lev);
       m2 = (double) t2/(gray_lev*gray_lev);
       m3 = (double) t3/(gray_lev*gray_lev);
       m4 = (double) t4/(gray_lev*gray_lev);

      for(i = 0;i < gray_lev;i++)
       for(j = 0;j < gray_lev;j++)
        {
        if((i-j)==1 || (i-j) == -1)
	  {
           f10_h = f10_h + pow(ph[i][j]/(2*w_size*(w_size-1))-m1,pow1);
           f10_v = f10_v + pow(pv[i][j]/(2*w_size*(w_size-1))-m2,pow1);
           f10_ld = f10_ld + pow(pld[i][j]/(2*(w_size-1)*(w_size-1))-m3,pow1) ;
           f10_rd = f10_rd + pow(prd[i][j]/(2*(w_size-1)*(w_size-1))-m4,pow1);
	  }
	}

      f10_h = (double) f10_h/(gray_lev*gray_lev);
      f10_v = (double) f10_v/(gray_lev*gray_lev);
      f10_ld = (double) f10_ld/(gray_lev*gray_lev);
      f10_rd = (double) f10_rd/(gray_lev*gray_lev);


     final_f10 = f10_h+f10_v+f10_ld+f10_rd;
     return(final_f10);
}
*/

void main(void)
{
  FILE 		*fin, *fout[10];
  unsigned int img_row;
  unsigned int img_col;
  unsigned int d,cc;
  unsigned int n_w, n_l, count1, count2, count3, count4;
  unsigned int step,total;
  double asm,ent,idm,con,crr;
  unsigned int win_size;
  unsigned char *store,*win_inside, in_file[50];
  fpos_t img_start;




      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 the value for d:");
	  scanf("%d",&d);

      printf("Please input the value of gray levels:");
	  scanf("%d",&gray_lev);

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


      image = (unsigned char **)calloc(img_row,sizeof(unsigned char *));
      for(count1 = 0 ; count1 < img_row; count1++)
          image[count1] = (unsigned char *)calloc(img_col,sizeof(unsigned char));

      pv = (int **)calloc(gray_lev+2,sizeof(int *));
      ph = (int **)calloc(gray_lev+2,sizeof(int *));
	  pld = (int **)calloc(gray_lev+2,sizeof(int *));
	  prd = (int **)calloc(gray_lev+2,sizeof(int *));
      for(count1 = 0 ; count1 < gray_lev+2; count1++)
	  {
          pv[count1] = (int *)calloc(gray_lev+2,sizeof(int));
		  ph[count1] = (int *)calloc(gray_lev+2,sizeof(int));
		  pld[count1] = (int *)calloc(gray_lev+2,sizeof(int));
		  prd[count1] = (int *)calloc(gray_lev+2,sizeof(int));
	  }




      store = (unsigned char *)calloc(img_col,sizeof(unsigned char));
      win_inside = (unsigned char *)calloc(window_size_r,sizeof(unsigned char));

      if((fin = fopen(in_file,"rb")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n",in_file);
	  exit(EXIT_FAILURE);
	}
      if((fout[1] = fopen("1asm_file","w")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n","2_asm");
	  exit(EXIT_FAILURE);
	}
      if((fout[2] = fopen("1entropy_file","w")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n","2_ent");
	  exit(EXIT_FAILURE);
	}
      if((fout[3] = fopen("1idm_file","w")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n","2_idm");
	  exit(EXIT_FAILURE);
	}
      if((fout[4] = fopen("1contrast_file","w")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n","2_con");
	  exit(EXIT_FAILURE);
	}
      if((fout[5]= fopen("1correlation_file","w")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n","2_crr");
	  exit(EXIT_FAILURE);
	}
 /*     if((fout[6] = fopen("f10","w")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n","2_f10");
	  exit(EXIT_FAILURE);
	}
	*/


      fgetpos(fin,&img_start);
      fsetpos(fin,&img_start);
      printf("\n");
      printf("Start Calculating ......\n");
      total = 0;
      n_w = img_col - win_size + 1;  /*printf("nw=%d\n",n_w);*/
      n_l = img_row - win_size + 1; /*printf("nl=%d\n",n_l);*/
      step = img_col - win_size;     /*printf("step=%d\n",step);*/
      N_r = n_w*n_l;                   /*printf("nr=%d\n",N_r);*/
      cc = 0;

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

      for(count1 = 0; count1 < n_l ; count1++)
	{
        start_initial();     /*initial window content and co-occurrence matrices*/
        for(count2 = 0;count2 < win_size;count2++)
          for(count3 = 0;count3 < win_size;count3++)
	    {
	     window[count2][count3] = image[count1+count2][count3];
	    }
        printf("processed line %d\n",count1);
        cal_ph(win_size,d);
        cal_pv(win_size,d);
        cal_pld(win_size,d);
        cal_prd(win_size,d);
        asm = cal_asm(win_size);
		ent = cal_ent(win_size);
        idm = cal_idm(win_size);
        con = cal_con(win_size);
        crr = cal_crr(win_size);
        /*f10 = cal_f10(win_size);*/
        fprintf(fout[1],"%lf ", asm);
        fprintf(fout[2],"%lf ", ent);
        fprintf(fout[3],"%lf ", idm);
        fprintf(fout[4],"%lf ", con);
        fprintf(fout[5],"%lf ", crr);
        /*fprintf(fout[6],"%lf ", f10);*/

        total++;

        for(count2 = 0; count2 < n_w-1; count2++)
	  {
          total = total + 1;

          subtract_ph(win_size,d);
          subtract_pv(win_size,d);
          subtract_pld(win_size,d);
          subtract_prd(win_size,d);

          for(count3 = 0;count3 < win_size;count3++)
           for(count4 = 0;count4 < win_size-1;count4++)
	    {
	     window[count3][count4] = window[count3][count4+1];
	    }

          for(count3 = 0;count3 < win_size;count3++)
	    {
	    window[count3][win_size-1] = image[count1+count3][count2+win_size];
	    }

          add_ph(win_size,d);
          add_pv(win_size,d);
          add_pld(win_size,d);
          add_prd(win_size,d);
          asm = cal_asm(win_size);
     	  ent = cal_ent(win_size);
          idm = cal_idm(win_size);
          con = cal_con(win_size);
          crr = cal_crr(win_size);
          /*f10 = cal_f10(win_size);*/
          fprintf(fout[1],"%lf ", asm);
          fprintf(fout[2],"%lf ", ent);
          fprintf(fout[3],"%lf ", idm);
          fprintf(fout[4],"%lf ", con);
          fprintf(fout[5],"%lf ", crr);
          /*fprintf(fout[6],"%lf ", f10);*/
	  }
	}

      printf("total pixels have be processed:%d\n",total);
      fclose(fin);
      fclose(fout[1]);
      fclose(fout[2]);
      fclose(fout[3]);
      fclose(fout[4]);
      fclose(fout[5]);
      /*fclose(fout[6]);*/


      for(count1 = 0 ; count1 < img_row; count1++)
           free(image[count1]);
	  free(image);

      for(count1 = 0 ; count1 < gray_lev; count1++)
	  {
		  free(pv[count1]);
		  free(ph[count1]);
		  free(pld[count1]);
		  free(prd[count1]);
	  }
	  free(pv);
	  free(ph);
	  free(pld);
	  free(prd);


      free(win_inside);
      free(store);
      printf("********************Done******************\n");
      printf("\n");

}