/*	Semi_variogram.c
**	This program is used to calculate
**      Semi_variogram and Fractal
**      dimension used for user
**      defined transect data!
*/

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#define cell_num 500
#define pixel_list 2000
main(int argc, char **argv)
{
  FILE 		*fin, *fout;
  unsigned long int img_length, y1, y2;
  unsigned long int img_width, x1, x2, transact[pixel_list], pix_num, temp[pixel_list],temp_num;
  unsigned long int low_bound, up_bound;
  unsigned long int count1, count2, mid_l, mid_h, a, b;
  unsigned long int step, lag, N_pairs, num, x_head, y_head, x_step, y_step;
  double x_add, y_add;
  double vari_sum, semi_variogram, sdx, sdy, sumy, sumx, avgx, avgy, intercept, r;
  double lagd, covall, covx, covy, slope, x_arry[cell_num], y_arry[cell_num], fractal;
  long int z1, z2, x_temp, y_temp, x_minus, y_minus, op_num, x_t, y_t, tune1, tune2;
  unsigned char *store;
  long address1,address2;
  fpos_t *img_start, *lag_start;
  img_start = &address1;
  lag_start = &address2;



  if(argc == 11)
    {
      img_length = atoi(argv[3]);
      img_width = atoi(argv[4]);
      low_bound = atoi(argv[5]);
      up_bound = atoi(argv[6]);
      x1 = atoi(argv[7]);/*Transact upper x axis*/
      y1 = atoi(argv[8]);/*Transact upper y axis*/
      x2 = atoi(argv[9]);/*Transact lower x axis*/
      y2 = atoi(argv[10]);/*Transact lower y axis*/
      store = (unsigned char *)malloc(1);

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

       if((fout = fopen(argv[2],"w")) == NULL)
	{
	  fprintf(stderr,"Failed to open %s\n",argv[2]);
	  exit(EXIT_FAILURE);
	}
      fprintf(fout,"Geostatistic ****************SEMIVARIOGRAM ******************* Analysis\n");
      fprintf(fout,"========================================================================\n");
      fprintf(fout,"                           by Brandt Tso\n");
      fprintf(fout,"\n");


      printf("\n");
      printf("Start Calculating ......\n");
      if(x1 >= x2 && y1 >= y2)
	{
         printf("Invalid input! Please check input list!\n");
         exit();
	}
      a = x1; b= y1;
      pix_num = 0;
      if( labs(y2-y1) < labs(x2-x1))
	{
         if((y2 - y1) != 0)
	   {
            x_minus = x2 - x1;
            y_minus = y2 - y1;
            x_add = (double) x_minus/y_minus;
            printf("x_minus = %d\n",x_minus);
            printf("y_minus = %d\n",y_minus);
            printf("x_add = %lf\n",x_add);
            x_temp = x_add;

            if((x_add-x_temp) > 0.5)
               x_temp = x_temp + 1;
            if((x_add-x_temp) < -0.5)
	      {
               x_temp = x_temp - 1;
               x_temp = -x_temp;
	      }
            if(x_temp < 0)
               x_temp = -x_temp;
            printf("x_temp = %d\n",x_temp);
            x_t = x_temp;
            x_head = x_temp/2 + 1;
            x_step = y2 - y1 - 1;
            printf("x_head = %d\n",x_head);
            printf("x_step = %d\n",x_step);
            op_num = x_temp*x_step + 2*x_head;
            if(x_minus > 0)
	      {
               op_num = op_num - x_minus;
               op_num = op_num/2;
	      }
            else
	      {
               op_num = op_num + x_minus;
               op_num = op_num/2;
	      }

            if(op_num >= 0)
	      {
               mid_l = x_step*0.5 - op_num;
               mid_h = x_step*0.5 + op_num;
	      }
            else
	      {
               mid_l = x_step*0.5 + op_num;
               mid_h = x_step*0.5 - op_num;
	      }
            printf("%d\n",mid_l);
            printf("%d\n",mid_h);

            if(x_minus > 0)
	      {
               for(count1 = 0; count1 < y1; count1++)
                  fseek(fin,img_width,1);
               fseek(fin,x1,1);
               for(count1 = 0; count1 < x_head; count1++)
	         {
                  fread(store,1,1,fin);
                  transact[pix_num] = store[0];
                  pix_num++;
                  printf("(%d,%d)",a++,b);
                 }
               fseek(fin,img_width,1);
               b++;
               for(count1 = 0; count1 < x_step; count1++)
    	         {
                  if(count1 >= mid_l && count1 <= mid_h)
		    {
                    if(op_num > 0)
                       x_temp = x_t - 1;
                    if(op_num < 0)
                       x_temp = x_t + 1;
                    }
                  else
                    x_temp = x_t;
                  for(count2 = 0; count2 < x_temp; count2++)
		    {
                     fread(store,1,1,fin);
                     transact[pix_num] = store[0];
                     pix_num++;
                     printf("(%d,%d)",a++,b);
		    }
                  fseek(fin,img_width,1);
                  b++;
		 }
               for(count1 = 0; count1 < x_head; count1++)
	         {
                  fread(store,1,1,fin);
                  transact[pix_num] = store[0];
                  pix_num++;
                  printf("(%d,%d)",a++,b);
                 }
	      }
            else
	      {
               for(count1 = 0; count1 < y1; count1++)
                  fseek(fin,img_width,1);
               fseek(fin,x1-x_head+1,1);
               fgetpos(fin,lag_start);
               temp_num = 0;

               for(count1 = 0; count1 < x_head; count1++)
	         {
                  fread(store,1,1,fin);
                  temp[temp_num] = store[0];
                  temp_num++;
                 }

               for(count1 = 0; count1 < x_step; count1++)
    	         {
                  fsetpos(fin,lag_start);
                  fseek(fin,img_width-x_temp,1);
                  fgetpos(fin,lag_start);
                  if(count1 >= mid_l && count1 <= mid_h)
		    {
                    if(op_num > 0)
                       x_temp = x_t - 1;
                    if(op_num < 0)
                       x_temp = x_t + 1;
                    }
                  else
                    x_temp = x_t;

                  for(count2 = 0; count2 < x_temp; count2++)
		    {
                     fread(store,1,1,fin);
                     temp[temp_num] = store[0];
                     temp_num++;
		    }
		 }

               fsetpos(fin,lag_start);
               fseek(fin,img_width-x_head+1,1);

               for(count1 = 0; count1 < x_head; count1++)
	         {
                  fread(store,1,1,fin);
                  temp[temp_num] = store[0];
                  temp_num++;
                 }

               pix_num = temp_num;
               tune2 = pix_num - x_head;
               tune1 = 0;
               for(count1 = 0; count1 < x_head; count1++)
		 {
                  transact[tune2+count1] = temp[tune1];
                  tune1++;
		 }

               for(count1 = 0; count1 < x_step; count1++)
		 {
                  if(count1 >= mid_l && count1 <= mid_h)
		    {
                    if(op_num > 0)
                       x_temp = x_t - 1;
                    if(op_num < 0)
                       x_temp = x_t + 1;
                    }
                  else
                    x_temp = x_t;
                  tune2 = tune2 - x_temp;
                  for(count2 = 0; count2 < x_temp; count2++)
		    {
                     transact[tune2+count2]= temp[tune1];
                     tune1++;
		    }
		 }

               for(count1 = 0; count1 < x_head; count1++)
		 {
                  transact[count1] = temp[tune1];
                  tune1++;
		 }

	      }
	   }
         else
	   {
             for(count1 = 0; count1 < y1; count1++)
                  fseek(fin,img_width,1);
             fseek(fin,x1,1);
             for(count1 = 0; count1 < (x2-x1+1); count1++)
	         {
                  fread(store,1,1,fin);
                  transact[pix_num] = store[0];
                  pix_num++;
                  printf("(%d,%d)",a++,b);
                 }

	   }

	}
      else
	{
         if((x2-x1) != 0)
	   {
            y_minus = y2 - y1;
            x_minus = x2 - x1;
            y_add = y_minus/x_minus;
            if(y_add < 0)
              y_add = -y_add;
            y_temp = y_add;
            printf("x_minus = %d\n",x_minus);
            printf("y_minus = %d\n",y_minus);
            printf("y_add = %lf\n",y_add);
            if((y_add-y_temp) > 0.5)
               y_temp = y_temp + 1;

            y_head = y_temp/2 + 1;
            printf("y_temp = %d\n",y_temp);
            y_t = y_temp;
            if(x_minus > 0)
              y_step = x2 - x1 - 1;
            else
	      {
               y_step = x2 - x1 + 1;
               y_step = -y_step;
	      }
            printf("y_head = %d\n",y_head);
            printf("y_step = %d\n",y_step);
            op_num = y_temp*y_step + 2*y_head;
            op_num = op_num - y_minus;
            op_num = op_num/2;

            if(op_num >= 0)
	      {
               mid_l = y_step*0.5 - op_num;
               mid_h = y_step*0.5 + op_num;
	      }
            else
	      {
               mid_l = y_step*0.5 + op_num;
               mid_h = y_step*0.5 - op_num;
	      }
            printf("%d\n",mid_l);
            printf("%d\n",mid_h);
            if(x_minus > 0)
	      {
               for(count1 = 0; count1 < y1; count1++)
                  fseek(fin,img_width,1);
               fseek(fin,x1,1);
               for(count1 = 0; count1 < y_head; count1++)
	         {
                  fgetpos(fin,lag_start);
                  fread(store,1,1,fin);
                  fsetpos(fin,lag_start);
                  transact[pix_num] = store[0];
                  pix_num++;
                  fseek(fin,img_width,1);
                  printf("(%d,%d)",a,b++);
                 }
               fseek(fin,1,1);
               a++;
               for(count1 = 0; count1 < y_step; count1++)
    	         {
                  if(count1 >= mid_l && count1 <= mid_h)
		    {
                    if(op_num > 0)
                       y_temp = y_t - 1;
                    if(op_num < 0)
                       y_temp = y_t + 1;
                    }
                  else
                    y_temp = y_t;
                  for(count2 = 0; count2 < y_temp; count2++)
		    {
                     fgetpos(fin,lag_start);
                     fread(store,1,1,fin);
                     transact[pix_num] = store[0];
                     pix_num++;
                     fsetpos(fin,lag_start);
                     fseek(fin,img_width,1);
                     printf("(%d,%d)",a,b++);
		    }
                  fseek(fin,1,1);
                  a++;
		 }
               for(count1 = 0; count1 < y_head; count1++)
	         {
                  fgetpos(fin,lag_start);
                  fread(store,1,1,fin);
                  transact[pix_num] = store[0];
                  pix_num++;
                  fsetpos(fin,lag_start);
                  fseek(fin,img_width,1);
                  printf("(%d,%d)",a,b++);
                 }
	      }
            else
	      {
               for(count1 = 0; count1 < y1; count1++)
                  fseek(fin,img_width,1);
               fseek(fin,x1,1);
               for(count1 = 0; count1 < y_head; count1++)
	         {
                  fgetpos(fin,lag_start);
                  fread(store,1,1,fin);
                  transact[pix_num] = store[0];
                  pix_num++;
                  fsetpos(fin,lag_start);
                  if(count1 == y_head - 1)
                    fseek(fin,img_width-1,1);
                  else
                    fseek(fin,img_width,1);
                  printf("(%d,%d)",a,b++);
                 }
               a--;
               for(count1 = 0; count1 < y_step; count1++)
    	         {
                  if(count1 >= mid_l && count1 <= mid_h)
		    {
                    if(op_num > 0)
                       y_temp = y_t - 1;
                    if(op_num < 0)
                       y_temp = y_t + 1;
                    }
                  else
                    y_temp = y_t;
                  for(count2 = 0; count2 < y_temp; count2++)
		    {
                     fgetpos(fin,lag_start);
                     fread(store,1,1,fin);
                     transact[pix_num] = store[0];
                     pix_num++;
                     fsetpos(fin,lag_start);
                     if(count2 == y_temp - 1)
                       fseek(fin,img_width-1,1);
                     else
                       fseek(fin,img_width,1);
                     printf("(%d,%d)",a,b++);
		    }
                   a--;
		 }


               for(count1 = 0; count1 < y_head; count1++)
	         {
                  fgetpos(fin,lag_start);
                  fread(store,1,1,fin);
                  transact[pix_num] = store[0];
                  pix_num++;
                  fsetpos(fin,lag_start);
                  fseek(fin,img_width,1);
                  printf("(%d,%d)",a,b++);
                 }
	      }
	   }
         else
	   {
             for(count1 = 0; count1 < y1; count1++)
                  fseek(fin,img_width,1);
             fseek(fin,x1,1);
             for(count1 = 0; count1 < (y2-y1+1); count1++)
	         {
                  fgetpos(fin,lag_start);
                  fread(store,1,1,fin);
                  transact[pix_num] = store[0];
                  pix_num++;
                  fsetpos(fin,lag_start);
                  fseek(fin,img_width,1);
                  printf("(%d,%d)",a,b++);
                 }

	   }

	}
 /********Above programme is used for extracting transact pixel data*/
 /*Follows starts to calculate the semivariogram*/

      printf("\n");
      for(count1 = 0; count1 < pix_num; count1++)
         printf(" %d ",transact[count1]);
      printf("\n");

      printf("%d\n",pix_num);
      sumy = 0;
      sumx = 0;
      num = 0;
      for(lag = low_bound; lag <= up_bound; lag++)
        {

          vari_sum = 0;
          N_pairs = pix_num - lag;
          for(count1 = 0; count1 < N_pairs ; count1++)
            {
             z1 = transact[count1];
             z2 = transact[count1+lag];
             vari_sum = vari_sum + (z1-z2)*(z1-z2);
	    }

          vari_sum = vari_sum/(2*N_pairs);
          printf("Semivariogram: r(%3d) = %12lf\n",lag,vari_sum);
          fprintf(fout,"Semivariogram: r(%3d) = %12lf\n",lag,vari_sum);
          lagd = lag;
          sumy = sumy + log(vari_sum);
          sumx = sumx + log(lagd);
          x_arry[num] = log(lagd);
          y_arry[num] = log(vari_sum);
          num++;
          if(lag < up_bound)
              printf("Calculating lag *%d* ...\n",lag+1);
          fsetpos(fin,img_start);
	}

       avgx = sumx / num;
       avgy = sumy / num;
       covx = 0;
       covy = 0;
       covall = 0;
       for(count1 = 0; count1 < num; count1++)
	 {
          covx = covx + (x_arry[count1]-avgx)*(x_arry[count1]-avgx);
          covy = covy + (y_arry[count1]-avgy)*(y_arry[count1]-avgy);
          covall = covall + (x_arry[count1]-avgx)*(y_arry[count1]-avgy);
	 }
       covx = covx/num;
       covy = covy/num;
       sdx = pow(covx,0.5);
       sdy = pow(covy,0.5);
       covall = covall/num;
       r = covall/(sdx*sdy);
       slope = r*sdy/sdx;
       intercept = avgy - (avgx*slope);
       fractal = (4-slope)/2;
       printf("\n");
       fprintf(fout,"\n");
       printf("Variance of ln(Lag)------------->%12lf\n",covx);
       fprintf(fout,"Variance of ln(Lag)------------->%12lf\n",covx);
       printf("Variance of ln(Semivariogram)--->%12lf\n",covy);
       fprintf(fout,"Variance of ln(Semivariogram)--->%12lf\n",covy);
       printf("The slope of regression line---->%12lf\n",slope);
       fprintf(fout,"The slope of regression line---->%12lf\n",slope);
       printf("The intercept of regression line:%12lf\n",intercept);
       fprintf(fout,"The intercept of regression line:%12lf\n",intercept);
       printf("\n");
       fprintf(fout,"\n");
       printf("The estimated FRACTAL DIMENSION :%12lf\n",fractal);
       fprintf(fout,"The estimated FRACTAL DIMENSION :%12lf\n",fractal);
       printf("\n");
       fprintf(fout,"\n");
      fclose(fin);
      fclose(fout);
      free(store);
      printf("****************************Done**************************\n");
      printf("\n");
    }
  else
      printf("\nUsage:\n\tsemivariogram   'input file'  'output file'  'rows in file' ' cols' 'lag_low_bound' 'lag_upper_bound' 'upper x axis' 'upper y axis' 'down x axis' 'down y axis'.\n\n");
}