/*=============================================================================

Program name: image_get_int.c
              ===============

$Date: 2009/05/27 15:13:48 $
$Author: schlepuetz $
$Source: /cvs/G/SPEC/local/X04SA/ES3/src/image_get_int.c,v $
$Revision: 1.4 $
$Tag: $

Description:  This program extracts various intensity information from binary
              image files.

Note:         This macro is based on an earlier version by O. Bunk,
              November 2004.

Authors(s):           O. Bunk (OB), C. M. Schlepuetz (CS)
Co-author(s):         R. Herger (RH), P. R. Willmott (PW)
Address:              Surface Diffraction Station
                      Materials Science Beamline X04SA
                      Swiss Light Source (SLS)
                      Paul Scherrer Institute
                      CH - 5232 Villigen PSI
Created:              2005/07/13

Compile:      g++ -o image_get_int image_get_int.c

Change Log:
-----------

2005/07/13 (CS):
- added information and comments, changed formatting to comply with X04SA
  conventions.

2005/12/12 (CS):
- implemented flexibility to read different binary file formats (different
  image sizes, header lengths, color depths).
- implemented a fourth threshold level.
- information on the image which is to be evaluated (name, header length,
  dimensions, box coordinates, etc.) is now stored in a file which is passed as
  a command line argument.

2006/01/27 (CS):
- fixed wrong order of arguments for median3x3

2006/09/21 (CS,PW):
- created mode for selectively median filtering of random hot and dead pixels:
  - additional input argument <filter_nsigma>: specifies by how many standard
    deviations the value of a pixel is allowed to differ from the average of
    its surrounding pixels.

2006/11/29 (PW,CS):
- added the wait() function.
- included wait loop when reading the image file to check that the file is
  complete. If file is not complete after 10 seconds (1000 * 0.01 seconds), an
  error message is issued and the program aborted.

=============================================================================*/

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <time.h>

/* uncomment the next lines to print debug info during execution */
//#define print_debug_info
//#define print_debug_info_2
//#define print_debug_info_3

#define PIXEL_TOTAL (image_xdim * image_ydim)
#define I(x,y) (img[(x) + image_xdim*(y)])
#define FALSE 0

/*---------------------------------------------------------------------------*/
/* wait function */
/*---------------------------------------------------------------------------*/

void wait ( double seconds )
{
  clock_t endwait;
  endwait = clock () + int(floor(seconds * CLOCKS_PER_SEC)) ;
  while (clock() < endwait) {}
}

/*---------------------------------------------------------------------------*/
/* functions to perform a 3x3 median filtering */
/*---------------------------------------------------------------------------*/

#define MED_BOX 3
#define MED_BOX_HALF (MED_BOX/2)

/*---------------------------------------------------------------------------*/
/* sub-function for the median filter below */
/* needed to sort the array elements */
/* refer to the qsort documentation in cstdlib.h */

int qsort_comp_func(const void *I1,const void *I2)
{
  if (*(int *)I1 < *(int *)I2)
    return -1;
  if (*(int* )I1 > *(int *)I2)
    return 1;
  return 0;
}

/*---------------------------------------------------------------------------*/
/* 3x3 median filter to the image data *img */

void median_3x3(unsigned long int *img,unsigned long int *filt_img,int image_xdim,int image_ydim,int x1,int y1,int x2,int y2)
{
  int x,y;

  if (x1 < MED_BOX_HALF)
    x1 = MED_BOX_HALF;
  if (y1 < MED_BOX_HALF)
    y1 = MED_BOX_HALF;
  if (x2 > image_xdim -1 -MED_BOX_HALF)
    x2 = image_xdim -1 -MED_BOX_HALF;
  if (y2 > image_ydim -1 -MED_BOX_HALF)
    y2 = image_ydim -1 -MED_BOX_HALF;
  #ifdef print_debug_info
  printf("median: (x1,y1)=(%d %d) (x2,y2)=(%d %d) box-size=%d half-box %d\n",
	 x1,y1,x2,y2,MED_BOX,MED_BOX_HALF);
  #endif
  
  for (x=x1;x <= x2;x++)
    for (y=y1;y <= y2;y++)
      {
	int I_values[MED_BOX*MED_BOX];
	int xm,ym,i = 0;
	for (xm=x-MED_BOX_HALF;xm <= x+MED_BOX_HALF;xm++)
	  for (ym=y-MED_BOX_HALF;ym <= y+MED_BOX_HALF;ym++)
	    {
	      int img_index =  xm + ym*image_xdim;

#ifdef print_debug_info
  printf("median box %5d %5d %5d\n",xm,ym,img[img_index]);
#endif

	      I_values[i++] = img[img_index];
	    }
	qsort(I_values,i,sizeof(*I_values),qsort_comp_func);

#ifdef print_debug_info
  for (i=0;i < MED_BOX*MED_BOX;i++)
    printf("qsort %d %5d\n",i,I_values[i]);
  printf("%d %d: %d\n",x,y,I_values[(MED_BOX*MED_BOX)/2]);
#endif	  
	filt_img[x + y*image_xdim] = I_values[(MED_BOX*MED_BOX)/2];
      }

  for (x=x1;x <= x2;x++)
    for (y=y1;y <= y2;y++)
      {
	int img_index =  x + y*image_xdim;
	img[img_index] = filt_img[img_index];
      }
  return;
}

/*---------------------------------------------------------------------------*/
/* main function image_get_int */
/*               ============= */
/*---------------------------------------------------------------------------*/

int main(const int no_of_params,const char **params)
{
  FILE *file_img_info;
  char filename_img[255];
  FILE *file_img;
  double image_timestamp;
  int image_header_length;
  int image_xdim,image_ydim,image_color_depth;
  int x1,y1,x2,y2;
  int bx1,by1,bx2,by2;
  int x,y,i;
  long int I_sum = 0,I_sum_bgr = 0,area_I,area_bgr;
  int thresh1_count = 0,thresh2_count = 0,thresh3_count = 0,thresh4_count = 0;
  int filter_median = 0;
  int parameter_okay = 1;
  int thresh1,thresh2,thresh3,thresh4;
  int retry;
  float filter_nsigma;

  /* check input parameter */
  if ((no_of_params < 6) || 
      (no_of_params > 8) ||
      (sscanf(params[ 2],"%d",&thresh1) != 1) ||
      (sscanf(params[ 3],"%d",&thresh2) != 1) ||
      (sscanf(params[ 4],"%d",&thresh3) != 1) ||
      (sscanf(params[ 5],"%d",&thresh4) != 1) ||
      (thresh1 < 0) || (thresh2 < 0) || (thresh3 < 0) || (thresh4 < 0)){
    parameter_okay = FALSE;
  }

#ifdef print_debug_info
  printf("No. of command-line parameters: %d\n",no_of_params);
#endif

  if (no_of_params > 6)
    {
      if (strcmp(params[ 6],"-median") == 0)
	filter_median = !FALSE;
      else
	parameter_okay = FALSE;
    }

  if (no_of_params > 7)
    {

#ifdef print_debug_info_2
  printf("params(7) = %s\n",params[7]);
#endif

      if (sscanf(params[ 7],"%g",&filter_nsigma) != 1)
        {
	  parameter_okay = FALSE;
        }
      if (filter_nsigma < 0)
        {
	  parameter_okay = FALSE;
        }

#ifdef print_debug_info_2
  printf("filter_nsigma = %f\n",filter_nsigma);
#endif

    }
  else
    {
      filter_nsigma = 0;
    }

  /* open image info text-file for read-only access */
  if ((file_img_info = fopen(params[1],"r")) == NULL)
    {
      fprintf(stderr,"Can't open image info file %s\n",params[1]);
      return 1;
    }
  else
    {

#ifdef print_debug_info
    printf("read image info file %s\n",params[1]);
#endif

      int read_count = fscanf (file_img_info,
                   "%s %f %d %d %d %d %d %d %d %d %d %d %d %d",
                   filename_img,&image_timestamp,&image_header_length,
                   &image_xdim,&image_ydim,&image_color_depth,
                   &x1,&y1,&x2,&y2,
                   &bx1,&by1,&bx2,&by2);

      if (read_count != 14)
        {
          fprintf(stderr,"Can't read image info file %s\n",params[1]);
          return 1;
        }
      fclose(file_img_info);

      /* check consistency of image information */
      if ((x1 < 0) || (y1 < 0) || (x2 < 0) || (y2 < 0) ||
          (x1 >= image_xdim) || (y1 >= image_ydim) ||
          (x2 >= image_xdim) || (y2 >= image_ydim) || 
          (x2 < x1) || (y2 < y1) || 
          (bx1 < 0) || (by1 < 0) || (bx2 < 0) || (by2 < 0) ||
          (bx1 >= image_xdim) || (by1 >= image_ydim) ||
          (bx2 >= image_xdim) || (by2 >= image_ydim) || 
          (bx2 < bx1) || (by2 < by1))
        {
          parameter_okay = FALSE;
        }
    }

  /* print help text in case of wrong parameters */
  if (!parameter_okay)
    {
      fprintf(stderr,"Get the number of pixel within a box (x1,y1,x2,y2)\n");
      fprintf(stderr,"with counts above four threshold values.\n");
      fprintf(stderr,"(0,0) is the upper left corner of the image.\n");
      fprintf(stderr,"(%d,%d) is the lower right corner.\n",
	      image_xdim-1,image_ydim-1);
      fprintf(stderr,"Median filtering can be used to remove dead and hot pixels.\n\n");
      fprintf(stderr,"A second integration box can be used to determine the background\n");
      fprintf(stderr,"(coordinates (bx1,by1) (bx2,by2)).\n");
      fprintf(stderr,"Usage:\n");
      fprintf(stderr,"%s <image_info_filename> <thresh1> <thresh2> <thresh3> <thresh4> [<-median> <filter_nsigma>]\n\n",
	      params[0]);
      fprintf(stderr,"All parameters are mandatory except the last -median switch and filter_nsigma value.\n");
      fprintf(stderr,"If the filter_nsigma value is set to zero or is unspecified, all the image will be automatically median filtered.\n");
      fprintf(stderr,"Positive values of filter_nsigma will lead to selective median filtering of individual pixels that differ by more than (filter_nsigma x st. dev.) of the eight neighbouring pixels.\n");
      fprintf(stderr,"The order of the parameters is fixed.\n");
      fprintf(stderr,"The image information file is automatically created by SPEC.\n");
      return 1;
    }

  unsigned long int img[PIXEL_TOTAL];

  area_I   = (x2- x1+1)*( y2- y1+1);
  area_bgr = (bx2-bx1+1)*(by2-by1+1);

#ifdef print_debug_info
  printf("( x1, y1)=(%3d %3d) ( x2, y2)=(%3d %3d) area=%ld\n",
         x1,y1,x2,y2,area_I);
  printf("(bx1,by1)=(%3d %3d) (bx2,by2)=(%3d %3d) area=%ld\n",
         bx1,by1,bx2,by2,area_bgr);
  if (filter_median)
    printf("with median filtering\n");
  else
    printf("no median filtering\n");
#endif

  /* open the actual image as text-file for read-only access */
  retry = 0;
  int pixel_count = 0;
  while ((pixel_count != PIXEL_TOTAL) & (retry < 1000)){

    /*fprintf(stderr,"Waiting for file %s\tretry = %d\n",filename_img,retry);*/

    if ((file_img = fopen(filename_img,"r")) == NULL)
      {
        fprintf(stderr,"Can't open %s\n",filename_img);
        return 1;
      }

    if(image_header_length >0)
      {
        fseek (file_img,image_header_length,SEEK_SET);
      }

    if(image_color_depth == 16)
      {
        unsigned short int img16[PIXEL_TOTAL];
        pixel_count = fread(img16,sizeof(unsigned short int),
                            PIXEL_TOTAL,file_img); 
        for (i=0; i<PIXEL_TOTAL; i++);
          {
             img[i] = (unsigned long int) img16[i];
          }
      }
    else if(image_color_depth == 32)
      {
        pixel_count = fread(img,sizeof(unsigned long int),
                            PIXEL_TOTAL,file_img);
      }
    else
      {
        fprintf(stderr,"Color-depth %d not supported.\n",image_color_depth);
        return 1;  
      }

    fclose(file_img);
    retry++;
    wait(0.01);
  }

  if (pixel_count != PIXEL_TOTAL)
    {
      fprintf(stderr,"Can't read pixel data.\n");
      fprintf(stderr,"read %d pixels of %d\n",pixel_count,PIXEL_TOTAL);
      return 1;
    }


  /*-------------------------------------*/
  /* optionally perform median filtering */ 

  if (filter_nsigma > 0)
    {
      unsigned long int filt_img[PIXEL_TOTAL];

      /* filter only bounding box of ROIs */
      int fx1 = (bx1<x1) ? bx1 : x1;
      int fy1 = (by1<y1) ? by1 : y1;
      int fx2 = (bx2>x2) ? bx2 : x2;
      int fy2 = (by2>y2) ? by2 : y2;

      /* make sure we do not operate on the image border */
      fx1 = (fx1>0) ? fx1 : 1;
      fy1 = (fy1>0) ? fy1 : 1;
      fx2 = (fx2<(image_xdim-1)) ? fx2 : (image_xdim-2);
      fy2 = (fy2<(image_ydim-1)) ? fy2 : (image_ydim-2);
        
      for (x=fx1; x<=fx2; x++){
        for (y=fy1; y<=fy2; y++){

          /* calculate average of eight surrounding pixels */
          int xb,yb,i = 0;
          double box_mean = 0;
          for (xb=x-MED_BOX_HALF;xb <= x+MED_BOX_HALF;xb++){
	    for (yb=y-MED_BOX_HALF;yb <= y+MED_BOX_HALF;yb++){
	      int img_index =  xb + yb*image_xdim;
                
#ifdef print_debug_info
  printf("median box %5d %5d %5d\n",xb,yb,img[img_index]);
#endif
              if (!((xb == x) && (yb == y)))
                {
    	          box_mean += img[img_index];
	        }
              }
            }
          box_mean = box_mean / (MED_BOX*MED_BOX-1);

#ifdef print_debug_info_2
  printf("box_mean = %f\n",box_mean);
#endif

          /* calculate standard deviation of eight surrounding pixels */
          xb = 0; yb = 0; i = 0;
          double box_sigma = 0;
          for (xb=x-MED_BOX_HALF;xb <= x+MED_BOX_HALF;xb++){
	    for (yb=y-MED_BOX_HALF;yb <= y+MED_BOX_HALF;yb++){
	      int img_index =  xb + yb*image_xdim;
                
#ifdef print_debug_info
  printf("median box %5d %5d %5d\n",xb,yb,img[img_index]);
#endif
              if (!((xb == x) && (yb == y)))
                {
    	          box_sigma += pow((img[img_index]-box_mean),2);
	        }
              }
            }
          box_sigma = sqrt(box_sigma / (MED_BOX*MED_BOX-1));

#ifdef print_debug_info_2
  printf("box_sigma = %f\n",box_sigma);
  printf("filter_nsigma = %f\n",filter_nsigma);
#endif

          if (fabs(I(x,y)-box_mean) > (filter_nsigma * box_sigma))
            {

#ifdef print_debug_info_3
  printf("filtered (%d,%d): %d --> ",x,y,I(x,y));
#endif

              median_3x3(img,filt_img,image_xdim,image_ydim,x,y,x,y);

#ifdef print_debug_info_3
  printf("%d\n",I(x,y));
#endif

            }
        }
      }
    }
  else if (filter_median)
    {

      unsigned long int filt_img[PIXEL_TOTAL];

      /* filter only bounding box of ROIs */
      int fx1 = (bx1<x1) ? bx1 : x1;
      int fy1 = (by1<y1) ? by1 : y1;
      int fx2 = (bx2>x2) ? bx2 : x2;
      int fy2 = (by2>y2) ? by2 : y2;

      /* make sure we do not operate on the image border */
      fx1 = (fx1>0) ? fx1 : 1;
      fy1 = (fy1>0) ? fy1 : 1;
      fx2 = (fx2<(image_xdim-1)) ? fx2 : (image_xdim-2);
      fy2 = (fy2<(image_ydim-1)) ? fy2 : (image_ydim-2);


      /* perform median filtering on the whole image */
      median_3x3(img,filt_img,image_xdim,image_ydim,fx1,fy1,fx2,fy2);
    }

  /* calculate the results (1 of 2): intensity integration area */
  for (x=x1;x <= x2;x++)
    {
      for (y=y1;y <= y2;y++)
	{
#ifdef print_debug_info
	  printf("(%5d,%5d) / index %5d: %5u\n",
		 x,y,x+y*image_xdim,img[x+y*image_xdim]);
#endif
	  I_sum += I(x,y);
	  if (I(x,y) > thresh1)
	    thresh1_count++;
	  if (I(x,y) > thresh2)
	    thresh2_count++;
	  if (I(x,y) > thresh3)
	    thresh3_count++;
	  if (I(x,y) > thresh4)
	    thresh4_count++;
	}
    }

  /* calculate the results (2 of 2): background integration area */
  for (x=bx1;x <= bx2;x++)
    {
      for (y=by1;y <= by2;y++)
	{
          /* only sum up background outside the intensity integration box */
          if ((x < x1) || (x > x2) || (y < y1) || (y > y2))
            {
#ifdef print_debug_info
              printf("(%5d,%5d) / index %5d: %5u\n",
                     x,y,x+y*image_xdim,img[x+y*image_xdim]);
#endif
              I_sum_bgr += I(x,y);
            }
          /* if pixel is inside the intensity integration box subtract it */
          /* from the background box area */
          else
            {
              area_bgr += -1;
            }
	}
    }

  /* print results */
  printf("I_sum %ld   area_I %ld  Threshold1 %d  Threshold2 %d  Threshold3 %d  Threshold4 %d  I_bgr %ld  area_bgr %ld\n",
	 I_sum,area_I,thresh1_count,thresh2_count,thresh3_count,thresh4_count,
         I_sum_bgr,area_bgr);
  return 0;
}

/*============================================================================

####################################################
 emacs setup:  force text mode to get no help with #
               indentation and force use of spaces #
               when tabbing.                       #
 Local Variables:                                  #
 mode:text                                         #
 indent-tabs-mode:nil                              #
 End:                                              #
####################################################

$Log: image_get_int.c,v $
Revision 1.4  2009/05/27 15:13:48  schlepuetz
commented out status message when waiting for file on server

Revision 1.3  2006/12/01 08:31:00  schlepuetz
added the wait() function; included wait loop when reading the image file to check that the file is complete

Revision 1.2  2006/09/21 11:42:39  schlepuetz
added mode for selectively filtering random hot and dead pixels

Revision 1.1  2006/03/02 13:24:27  schlepuetz
Program to extract intensities from pixel detector data - first release, verified to work on PILATUS I and PILATUS II in conjunction with .img files.


============================= End of $RCSfile: image_get_int.c,v $ ===*/