/*	Lacunarity.c
**	This program is used to calculate
**      Lacunarity Index, can be used from
**      user difined window sizes and 
**      different pixel ranges of insterest!
**======Landscape Ecology vol.8, no.3======* 
**======pp.201-211, 1993===================*
**======Lacunarity Index===================* 
*/

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#define window_size_r 800
#define frequency_S 640000

 long int n[frequency_S];
 long float sq[frequency_S], q[frequency_S];
 long float s2q[frequency_S];
 long int N_r;

main(int argc, char **argv)
{
  FILE 		*fin, *fout;
  unsigned long int count, img_length, topat;
  unsigned long int max_single, starcount, startot, img_width;
  unsigned long int mass1, mass2, low_bound, up_bound;
  unsigned long int n_w, n_l, c1,count1, count2, count3, count4, S, r;
  unsigned long int step, win_sum; 
  long float lacu_index, cal_lacu();
  long int win_size; 
  unsigned char *store,*win_inside ;
  long address1,address2;
  fpos_t *img_start, *win_start;
  img_start = &address1;
  win_start = &address2;



  if(argc == 9)
    {
      img_length = atoi(argv[3]);
      img_width = atoi(argv[4]);
      mass1 = atoi(argv[5]);
      mass2 = atoi(argv[6]); 
      low_bound = atoi(argv[7]); 
      up_bound = atoi(argv[8]);
      store = (unsigned char *)malloc(img_width);
      win_inside = (unsigned char *)malloc(window_size_r);
      
      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 ***************** LACUNARITY ********************* Analysis\n");
      fprintf(fout,"========================================================================\n");
      fprintf(fout,"                           by Brandt Tso\n");
      fprintf(fout,"\n");
      
      fgetpos(fin,img_start);
      /*
      for(count1=0;count1<img_length;count1++)
	{
         fread(store,1,img_width,fin);
           for(count2=0;count2<img_length;count2++)
             printf("%d ",store[count2]);
           printf("\n");
	}
	*/
      fsetpos(fin,img_start);
      printf("\n");
      printf("Start Calculating ......\n");
      for(win_size = low_bound; win_size <= up_bound; win_size++)
        {
          /*printf("calculating window size %d , Please wait..",win_size);*/
          start_initial(win_size);
          n_w = img_width - win_size + 1;
          n_l = img_length - win_size + 1;
          step = img_width - win_size;
          N_r = n_w*n_l;
          
          for(count1 = 0; count1 < n_l ; count1++)
            {
              for(count2 = 0; count2 < n_w; count2++)
                {           
                  fgetpos(fin,win_start);  
                  win_sum = 0;    
	          for(count3 = 0; count3 < win_size; count3++)
	            {
                      fread(win_inside,1,win_size,fin);
                      for(count4 = 0; count4 < win_size; count4++)
                        {
	                  if( win_inside[count4] >= mass1 && win_inside[count4]<= mass2 )
	                    win_sum++;
                        }
                      if(count3 < win_size - 1)
                        fseek(fin,step,1);
                      /*
                      if(count2 == 1)            
                      for (count4 = 0;count4 < win_size;count4++)
                        printf("%d ",win_inside[count4]);
                      printf("\n");
		      */
                                        
		    }                 
                  n[win_sum]++; 
                  fsetpos(fin,win_start);
		  if(count2 < n_w-1)
                    fseek(fin,1,1);
		}                  
              if(count1 < n_l-1)
                fseek(fin,win_size,1);
	    }
          fprintf(fout,"Window size %d\n",win_size);
	  lacu_index = cal_lacu(win_size,fout);
           
          fprintf(fout,"window size %d--> lacunarity = %f\n",win_size,lacu_index);
          if(win_size < up_bound)
              printf("Calculating window size *%d* ...\n",win_size+1);          
          fsetpos(fin,img_start);
	}
	
      fclose(fin);
      fclose(fout);   
      free(win_inside);
      free(store); 
      printf("********************Done******************\n");
      printf("\n");
    }
  else
      printf("\nUsage:\n\tlacunarity   'input file'  'output file'  'rows in file' ' cols' 'value_insterest_low_bound' 'value_insterest_upper_bound' 'win_low_bound' 'win_upper_bound' .\n\n");
}

start_initial(input_win_size)
long int input_win_size;
  {
long int i,j;

      j = input_win_size*input_win_size;
      for(i=0;i<=j;i++)
               {
               n[i] = 0;
               q[i] = 0;
               sq[i] = 0;
               s2q[i] = 0;
	       }
  }

long float cal_lacu(in_win_size,output)
long int in_win_size;
FILE *output;
  {
long int i,j;
long float sum_sq,sum_s2q,lacu; 
       
      sum_sq = 0;
      sum_s2q = 0; 
      j = in_win_size*in_win_size;
      fprintf(output,"         S    n(S,r)         Q(S,r)          SQ(S,r)         S2Q(Sr)\n");
      for(i=0;i<=j;i++)
	if(n[i] != 0)
               {
		fprintf(output,"%10d",i);
                q[i] = (float) n[i]/N_r;
		fprintf(output,"%10d",n[i]);
                fprintf(output,"%16f",q[i]);
                
	        sq[i] = i*q[i];
                fprintf(output,"%16f",sq[i]);
                s2q[i] = i*i*q[i];
                fprintf(output,"%16f\n",s2q[i]);
                sum_sq = sum_sq + sq[i];
                sum_s2q = sum_s2q + s2q[i];
	       }
       if(sum_sq == 0)
	 {
          printf("No pixel value of interest, Please Check!\n");
          exit();
	 }
        fprintf(output,"Z[1] =     %f\n",sum_sq);
        fprintf(output,"Z[2] =     %f\n",sum_s2q);
/* printf("sum ====   %f\n",sum_sq);*/
       sum_sq = sum_sq*sum_sq;
       lacu = (float) sum_s2q/sum_sq;
       printf("window-%d: Lacunarity==>  %f\n",in_win_size,lacu);
       return(lacu);
  }