/***************************************************************************************
 * pixel_pull.c                                                                        *
 *                                                                                     *
 * Takes S sets of align files and N sets of image files, and generates S sets of      *
 * bead files, each a binary file having 7 colums x O objects (1 per bead).  The 7     *
 * colums are (in order) frame_num, object_num, object_size, centroid_x, centroid_y,   *
 *                       border_flag, object_value                                     *
 *                                                                                     *
 * TO COMPILE: g++ pixel_pull.c -O3 -ffloat-store -foptimize-sibling-calls             * 
 *             -fno-branch-count-reg -o pixel_pull -ltiff -lm -I../../INCLUDE          * 
 *             -L../../LIB                                                             *
 *                                                                                     *
 * IMPORTANT: pixel_pull expects to find white-light images in IMAGES/RAW/999/         *
 *                                                                                     *
 * USAGE:                                                                              *
 *       ./pixel_pull xcols yrows base_dir/ set1 set2 set3 ... setN                    *
 *       WHERE:                                                                        *
 *             xcols = number of colums in the images (e.g. 1000)                      *
 *             yrows = number of rows in the images (e.g. 1000)                        *
 *             base_dir/ = the base dir (e.g. /mnt/data_seq1/B6/)                      *
 *             set is the index of a set (e.g. 4 for set 004                           *
 *                                                                                     *
 * OUTPUT goes to PROC/BEAD/                                                           *
 *                                                                                     *
 * WRITTEN BY: Greg Porreca, Jay Shendure (Church Lab) 02-20-2005                      *
 *                                                                                     *
 *                                                                                     *
 ***************************************************************************************/

/* NOTE: for objects > 65535 pixels (max_val(short unsigned int)), their reported size
will not be accurate since object_size is cast as short unsigned int for output.  object
intensity will be accurate, though, since object_size is maintained in memory as int.  No
other values in the output should exceed 65535.
*/
 
#include <stdio.h>
#include <tiffio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <dirent.h>
#include <limits.h>
#include <map>
#include <iostream>
#include <string>
#include <algorithm>
#include <iterator>
#include <fstream>
#include <fnmatch.h>
using namespace std;

#define MAX_ROOTPATH_LEN 100
#define MAX_FULLPATH_LEN 150
#define MAX_OBJECTS 64000 //max objects per frame
#define MAX_NUM_SETS 20 //maximum number of image sets to process
#define MAX_FILESIZE 2252800  // maximum size (in bytes) of image files allowed; 2200kB
#define PATH_RAW "IMAGES/RAW/"
#define PATH_OBJECT "IMAGES/OBJECT/"
#define PATH_ALIGN "PROC/ALIGN/"
#define PATH_BEADS "PROC/BEAD/"
#define PATH_BKGRND "PROC/BKGRND/"
#define PATH_BEADMEAN "PROC/BEADMEAN/"
#define OBJ_COUNT_OUTFILENAME "IMAGES/OBJECT/obj_counts.dat"
//#define DEBUG 1

unsigned short int xcols, yrows;

int scandir_sel(struct dirent*);
void load_image(char*, char*);

int main(int argc, char *argv[]){

  char raw_fn[MAX_FULLPATH_LEN], 
    object_fn[MAX_FULLPATH_LEN],
    obj_count_fn[MAX_FULLPATH_LEN], 
    align_fn[MAX_FULLPATH_LEN], 
    bead_fn[MAX_FULLPATH_LEN],
    bkgrnd_fn[MAX_FULLPATH_LEN],
    beadmean_fn[MAX_FULLPATH_LEN],
    object_dir[MAX_FULLPATH_LEN],
    ROOT_PATH[MAX_ROOTPATH_LEN],
    curr_fn_index[5],
    curr_set_index[4],
    curr_obj_count[20],
    text[500];

  char *raw_image,
    *object_image;

  int num_frames,
    num_sets,
    set_index[MAX_NUM_SETS],
    curr_set,
    curr_frame,
    numeric,
    curr_obj,
    OBJ_COUNT,
    count2,
    count4,
    curr_y_offset,
    curr_x_offset,
    curr_pixel,
    count,
    count3,
    count5,
    field_size,
    success,
    object_size_short;

  double bkgrnd_sum[MAX_NUM_SETS],
    beadmean_sum[MAX_NUM_SETS],
    bkgrnd_count[MAX_NUM_SETS];

  
  FILE *obj_count_infile,
    *infile[MAX_NUM_SETS],
    *outfile[MAX_NUM_SETS],
    *bkgrndfile[MAX_NUM_SETS],
    *beadmeanfile[MAX_NUM_SETS];

  int *align_x[MAX_NUM_SETS], 
    *align_y[MAX_NUM_SETS], 
    *object_val[MAX_NUM_SETS],
    *object_size, 
    *object_border, 
    *obj_count;
  
  float *centroids_x, 
    *centroids_y;

  short unsigned int *IP,
    *IP2;

  short unsigned int output_frame,
    output_set,
    output_size,
    output_centroid_x,
    output_centroid_y,
    output_border,
    output_value,
    bkgrnd_val,
    beadmean_val;

  struct dirent **namelist;
  
  // Set warning handler for libtiff to void
  TIFFSetWarningHandler(0);
  
  field_size = sizeof(unsigned short int);
  

  if(argc < 5){
    fprintf(stdout, "\n");
    fprintf(stdout, "To run this program, execute the following:\n");
    fprintf(stdout, " ./pixel_pull xcols yrows /base/dir/ cyclenums\n");
    fprintf(stdout, "  WHERE:\n");
    fprintf(stdout, "        xcols = the number of columns (in pixels) in the image; our camera is 1000\n");
    fprintf(stdout, "        yrows = the number of rows (in pixels) in the image; our camera is 1000\n");
    fprintf(stdout, "        /base/dir/ = the full path to the data directory; if your data directory is \"E1\", and its full\n");
    fprintf(stdout, "                     path is \"/mnt/data1/\", you would specify \"/mnt/data1/E1/\"\n");
    fprintf(stdout, "        cyclenums = numbers of each cycle directory to process; for example, if you want to process\n");
    fprintf(stdout, "                    images from cycles 001, 002, 003, 004, and 005, specify \"1 2 3 4 5\"\n");
    fprintf(stdout, "\n");
    exit(0);
  }


  //****************** Parse command-line arguments **********************
  //
  yrows = atoi(argv[1]);
  xcols = atoi(argv[2]);
  strcpy(ROOT_PATH, argv[3]);
  num_sets = argc - 4;
  //
  //**********************************************************************

  //***************** Allocate memory for object image *******************
  //
  if((object_image = (char*) malloc(MAX_FILESIZE*sizeof(char))) == NULL){
    fprintf(stdout, "Could not allocate enough memory for the image.\n");
    exit(42);
  }
  //
  //**********************************************************************

  //******************* Allocate memory for raw image ********************
  //
  if((raw_image = (char*) malloc(MAX_FILESIZE*sizeof(char))) == NULL){
    fprintf(stdout, "Could not allocate enough memory for the image.\n");
    exit(42);
  }
  //
  //**********************************************************************

  //****************** Determine number of frames ************************
  //                                 
  // Generate path to object images
  strcpy(object_dir, "");
  strcat(object_dir, ROOT_PATH);
  strcat(object_dir, PATH_OBJECT);
  // Get number of images in this dir
  num_frames = scandir(object_dir, &namelist, (int (*)(const struct dirent *))(scandir_sel), alphasort);
  free(namelist);
  //
  //**********************************************************************

  //****************** Allocate memory for align data ********************
  //

#ifdef DEBUG
  fprintf(stdout, "Allocating memory for %d(%d) x %d sets of alignment data.\n", num_sets,MAX_NUM_SETS, num_frames);
#endif

  for(curr_set=0; curr_set < num_sets; curr_set++){
    if((align_x[curr_set] = (int*) malloc(num_frames * sizeof(int))) == NULL){
      fprintf(stdout, "Could not allocate enough memory for align_x list.\n");
      exit(42);
    }
    if((align_y[curr_set] = (int*) malloc(num_frames * sizeof(int))) == NULL){
      fprintf(stdout, "Could not allocate enough memory for align_y list.\n");
      exit(42);
    }
  }
  //
  //**********************************************************************

  //**************** Load number of objects for each frame ***************
  //

#ifdef DEBUG
  fprintf(stdout, "Loading number of objects per frame from file.\n");
#endif
  
  if((obj_count = (int*) malloc(num_frames * sizeof(int))) == NULL){
    fprintf(stdout, "Could not allocate enough memory for the object count matrix.\n");
    exit(42);
  }
  // Open object count file
  strcpy(obj_count_fn, "");
  strcat(obj_count_fn, ROOT_PATH);
  strcat(obj_count_fn, OBJ_COUNT_OUTFILENAME);
  obj_count_infile = fopen(obj_count_fn, "r");
  for(curr_frame = 0; curr_frame < num_frames; curr_frame++){
    fgets(curr_obj_count, 19, obj_count_infile);
    fgets(curr_obj_count, 19, obj_count_infile);
    *(obj_count + curr_frame) = atoi(curr_obj_count);
  }
  fclose(obj_count_infile);
  //
  //**********************************************************************

  //*********** Open input align files and output bead files *************
  //

#ifdef DEBUG
  fprintf(stdout, "Opening input align and output bead files (%d sets).\n",num_sets);
#endif

  for(curr_set = 0; curr_set < num_sets; curr_set++){
    set_index[curr_set] = atoi(argv[curr_set+4]); //set indices come from cmd-line args
    sprintf(curr_set_index, "%03d", *(set_index + curr_set)); //convert to 0-padded 4-digit
    // Generate filename for current input (align) file
    strcpy(align_fn, "");
    strcat(align_fn, ROOT_PATH);
    strcat(align_fn, PATH_ALIGN);
    strcat(align_fn, curr_set_index);
    strcat(align_fn, ".align");
    // Generate filename for current output (bead) file
    strcpy(bead_fn, "");
    strcat(bead_fn, ROOT_PATH);
    strcat(bead_fn, PATH_BEADS);
    strcat(bead_fn, curr_set_index);
    strcat(bead_fn, ".beads");
    // Open files
#ifdef DEBUG
    fprintf(stdout, "Opening file %s for input; ", align_fn);
#endif
    if((*(infile + curr_set) = fopen(align_fn, "r")) == NULL){
      fprintf(stdout, "ERROR OPENING INFILE.\n");
      exit(42);
    }
#ifdef DEBUG
    fprintf(stdout, "Opening file %s for output.\n", bead_fn);
#endif
    if((*(outfile + curr_set) = fopen(bead_fn, "w")) == NULL){
      fprintf(stdout, "ERROR OPENING OUTFILE.\n");
      exit(42);
    }
  }
  //
  //**********************************************************************

  //*********************** Open output bkgrnd files *********************
  //
  for(curr_set = 0; curr_set < num_sets; curr_set++){
    set_index[curr_set] = atoi(argv[curr_set+4]); //set indices come from cmd-line args
    sprintf(curr_set_index, "%03d", *(set_index + curr_set)); //convert to 0-padded 4-digit
    strcpy(bkgrnd_fn, "");
    strcat(bkgrnd_fn, ROOT_PATH);
    strcat(bkgrnd_fn, PATH_BKGRND);
    strcat(bkgrnd_fn, curr_set_index);
    strcat(bkgrnd_fn, ".bkgrnd");
#ifdef DEBUG
    fprintf(stdout, "Opening file %s for output.\n", bkgrnd_fn);
#endif
    if((*(bkgrndfile + curr_set) = fopen(bkgrnd_fn, "w")) == NULL){
      fprintf(stdout, "ERROR OPENING OUTFILE.\n");
      exit(42);
    }
  }

  //*********************** Open output beadmean files *********************
  //
  for(curr_set = 0; curr_set < num_sets; curr_set++){
    set_index[curr_set] = atoi(argv[curr_set+4]); //set indices come from cmd-line args
    sprintf(curr_set_index, "%03d", *(set_index + curr_set)); //convert to 0-padded 4-digit
    strcpy(beadmean_fn, "");
    strcat(beadmean_fn, ROOT_PATH);
    strcat(beadmean_fn, PATH_BEADMEAN);
    strcat(beadmean_fn, curr_set_index);
    strcat(beadmean_fn, ".beadmean");
#ifdef DEBUG
    fprintf(stdout, "Opening file %s for output.\n", beadmean_fn);
#endif
    if((*(beadmeanfile + curr_set) = fopen(beadmean_fn, "w")) == NULL){
      fprintf(stdout, "ERROR OPENING OUTFILE.\n");
      exit(42);
    }
  }


  //******************* Load align files into memory *********************
  //
#ifdef DEBUG
  fprintf(stdout, "Loading align data files into memory (%d sets).\n", num_sets);
#endif
  for(curr_set = 0; curr_set < num_sets; curr_set++){
    for(curr_frame = 0; curr_frame < num_frames; curr_frame++){
      fscanf(*(infile + curr_set), "%s", text);

      fscanf(*(infile + curr_set), "%s", text);
      *(align_x[curr_set]+curr_frame) = atoi(text);

      fscanf(*(infile + curr_set), "%s", text);
      *(align_y[curr_set]+curr_frame) = atoi(text);

      fscanf(*(infile + curr_set), "%s", text);
    }
  }
  //
  //**********************************************************************

  //**************** ITERATE OVER ALL FRAMES AND ALL SETS ****************
  //
  for(curr_frame = 0; curr_frame < num_frames; curr_frame++){
    fprintf(stdout, "Processing frame %d of %d.\n", curr_frame + 1, num_frames);

    //************* Open object image for current frame ********************
    //
    //current image index, as a 4-digit string
    sprintf(curr_fn_index, "%04d", curr_frame + 1);    
    //generate complete filename for current object image
    strcpy(object_fn, "");
    strcat(object_fn, ROOT_PATH);
    strcat(object_fn, PATH_OBJECT);
    strcat(object_fn, curr_fn_index);
    strcat(object_fn, ".tif");
    //load object image for current frame
    load_image(object_fn, object_image);
    IP = (short unsigned int*) object_image;
    //
    //**********************************************************************

    //*************** Allocate memory for centroids, etc *******************
    //
#ifdef DEBUG
  fprintf(stdout, "Allocating memory for centroids, etc..\n");
#endif
  
    if((centroids_x = (float*) malloc(*(obj_count + curr_frame) *sizeof(float))) == NULL){
      fprintf(stdout, "Could not allocate enough memory for centroids_x list.\n");
      exit(42);
    }
    if((centroids_y = (float*) malloc(*(obj_count + curr_frame)*sizeof(float))) == NULL){
      fprintf(stdout, "Could not allocate enough memory for centroids_y list.\n");
      exit(42);
    }
    if((object_size = (int*) malloc(*(obj_count + curr_frame)*sizeof(int))) == NULL){
      fprintf(stdout, "Could not allocate enough memory for object_size list.\n");
      exit(42);
    }
    if((object_border = (int*) malloc(*(obj_count + curr_frame)*sizeof(int))) == NULL){
      fprintf(stdout, "Could not allocate enough memory for object_border list.\n");
      exit(42);
    }

    for(curr_set = 0; curr_set < num_sets; curr_set++){
      if((object_val[curr_set] = (int*) malloc(*(obj_count + curr_frame)*sizeof(int))) == NULL){
	fprintf(stdout, "Could not allocate enough memory for object_val list.\n");
	exit(42);
      }
      for(curr_obj = 0; curr_obj < *(obj_count + curr_frame); curr_obj++){
	object_val[curr_set][curr_obj] = 0;
      }
    }
#ifdef DEBUG
    fprintf(stdout, "Memory for centroids, etc. allocated successfully\n");
#endif
    //
    //**********************************************************************
    
    //**************** Initialize vars for current frame *******************
    //
    if(*(obj_count + curr_frame) == 0){
#ifdef DEBUG
    fprintf(stdout, "Special case: No objects in this frame.\n");
#endif
      *(centroids_x) = 0;
      *(centroids_y) = 0;
      *(object_size) = 0;
      *(object_border) = 0;
    }
  
    for (count = 0; count < *(obj_count + curr_frame); count++) {
      *(centroids_x + count) = 0;
      *(centroids_y + count) = 0;
      *(object_size + count) = 0;
      *(object_border + count) = 0;
    }
    for (curr_set = 0; curr_set < num_sets; curr_set++){
      *(bkgrnd_sum + curr_set) = 0;
      *(bkgrnd_count + curr_set) = 0;
      *(beadmean_sum + curr_set) = 0;
    }

    if(*(obj_count + curr_frame) == 0){
#ifdef DEBUG
      fprintf(stdout, "Special case: No objects in this frame, go on to next frame.\n");
#endif
      continue;
    }

    OBJ_COUNT = 0;
    //
    //**********************************************************************

    //********* Calculate centroids & obj sizes for current frame **********
    //
    for(curr_pixel = 0; curr_pixel < (xcols * yrows); curr_pixel++){
      if (*(IP+curr_pixel) > 0) {
	if (*(IP+curr_pixel) > OBJ_COUNT){
	  (OBJ_COUNT) = *(IP+curr_pixel);
	}
	*(object_size + *(IP + curr_pixel) - 1) = *(object_size + *(IP + curr_pixel) - 1) + 1;
	*(centroids_x + *(IP + curr_pixel) - 1) += (curr_pixel % xcols) + 1;
	*(centroids_y + *(IP + curr_pixel) - 1) += (int) floor(float(curr_pixel)/xcols) + 1;
      }
    }
    
    for(curr_obj = 0; curr_obj < *(obj_count + curr_frame); curr_obj++){
      *(centroids_x + curr_obj) = *(centroids_x + curr_obj) / *(object_size + curr_obj);
      *(centroids_y + curr_obj) = *(centroids_y + curr_obj) / *(object_size + curr_obj);
    }
    //
    //**********************************************************************
    

    //***** Iterate over all sets for the current frame & pull pixels ******
    //
    for(curr_set = 0; curr_set < num_sets; curr_set++){

#ifdef DEBUG
      fprintf(stdout, "Generating complete filename for current raw image file.\n");
#endif
      //generate complete filename for current raw image file
      sprintf(curr_set_index, "%03d", *(set_index + curr_set));
      strcpy(raw_fn, "");
      strcat(raw_fn, ROOT_PATH);
      strcat(raw_fn, PATH_RAW);
      strcat(raw_fn, curr_set_index);
      strcat(raw_fn, "/SC_");
      strcat(raw_fn, curr_fn_index);
      strcat(raw_fn, ".tif");
      
      //load raw image for current frame,set
      load_image(raw_fn, raw_image);
      IP2 = (short unsigned int*) raw_image;
      
      //iterate over all pixels in the object image to calculate centroids
#ifdef DEBUG
      fprintf(stdout, "Calculating centroids for current image.\n");
#endif
      for(curr_pixel = 0; curr_pixel < xcols * yrows; curr_pixel++){

	*(bkgrnd_sum + curr_set) += *(IP2 + curr_pixel);
	*(bkgrnd_count + curr_set) = *(bkgrnd_count + curr_set) + 1;

	if( *(IP + curr_pixel) > 0){
	  
	  count3 = (curr_pixel % xcols) + 1;
	  count4 = (int) floor(float(curr_pixel) / xcols) + 1;
	  count3 += align_x[curr_set][curr_frame];
	  count4 += align_y[curr_set][curr_frame];
	  if (count3 >= 1 && count3 <= 1000 && count4 >= 1 && count4 <= 1000) {
	    count5 = ((count4 - 1) * xcols) + count3 - 1;
	    object_val[curr_set][*(IP + curr_pixel) - 1] += *(IP2 + count5);
	    *(beadmean_sum + curr_set) += *(IP2 + count5);
	  } else {
	    *(object_border + *(IP + curr_pixel) - 1) = 1;
	  }
	}
      }//end for(curr_pixel)

#ifdef DEBUG
      fprintf(stdout, "Writing data (frame %d, set %d)\n", curr_frame + 1, curr_set + 1);
      fprintf(stdout, "Iterating over %d objects:\n", *(obj_count + curr_frame));
#endif
      for(curr_obj = 0; curr_obj < *(obj_count + curr_frame); curr_obj++){
	// make sure object size will not overflow unsigned short int
	/*	if (curr_frame > 91){
	  fprintf(stdout,"%d\t",curr_obj);
	  fprintf(stdout,"%d\n",*(object_size + curr_obj));
	  }*/
	if(*(object_size + curr_obj) > 65535){
	  object_size_short = 65535;
	  //fprintf(stdout, "Got here.\n");
	}
	else{
	  object_size_short = *(object_size + curr_obj);
	}

	/*	if(curr_frame > 91){
	fprintf(stdout,"%d\t",curr_frame);
	fprintf(stdout,"%d\t",curr_obj);
	fprintf(stdout,"%d\t",*(object_size_short));
	fprintf(stdout,"%d\t",*(centroids_x + curr_obj));
	fprintf(stdout,"%d\t",*(centroids_y + curr_obj));
	fprintf(stdout,"%d\t",*(object_border + curr_obj));
	fprintf(stdout,"%d\t",object_val[curr_set][curr_obj]);
	}*/
	
#ifdef DEBUG
	fprintf(stdout, "Assigning values to output variables\n");
#endif

	output_frame = (short unsigned int) curr_frame + 1;
	output_set = (short unsigned int) curr_obj + 1;
	output_size = (short unsigned int) object_size_short;
	output_centroid_x = (short unsigned int) *(centroids_x + curr_obj);
	output_centroid_y = (short unsigned int) *(centroids_y + curr_obj);
	output_border = (short unsigned int) *(object_border + curr_obj);
	output_value = (short unsigned int) (object_val[curr_set][curr_obj] / *(object_size + curr_obj));

	
#ifdef DEBUG
	fprintf(stdout, "Writing values to file\n");
#endif
	if((success=fwrite(&output_frame, field_size, 1, *(outfile + curr_set))) < 1){
	  fprintf(stdout, "ERROR writing output_frame to bead file %4.0d", curr_set); 
	  exit(42);
	}
	if((success=fwrite(&output_set, field_size, 1, *(outfile + curr_set))) < 1){
	  fprintf(stdout, "ERROR writing output_set to bead file %4.0d", curr_set); 
	  exit(42);
	}
	if((success=fwrite(&output_size, field_size, 1, *(outfile + curr_set))) < 1){
	  fprintf(stdout, "ERROR writing output_size to bead file %4.0d", curr_set); 
	  exit(42);
	}
	if((success=fwrite(&output_centroid_x, field_size, 1, *(outfile + curr_set))) < 1){
	  fprintf(stdout, "ERROR writing output_centroid_x to bead file %4.0d", curr_set); 
	  exit(42);
	}
	if((success=fwrite(&output_centroid_y, field_size, 1, *(outfile + curr_set))) < 1){
	  fprintf(stdout, "ERROR writing output_centroid_y to bead file %4.0d", curr_set); 
	  exit(42);
	}
	if((success=fwrite(&output_border, field_size, 1, *(outfile + curr_set))) < 1){
	  fprintf(stdout, "ERROR writing output_border to bead file %4.0d", curr_set); 
	  exit(42);
	}
	if((success=fwrite(&output_value, field_size, 1, *(outfile + curr_set))) < 1){
	  fprintf(stdout, "ERROR writing output_value to bead file %4.0d", curr_set); 
	  exit(42);
	}
	
	
	/*	if(curr_frame > 91){
	fprintf(stdout, "%d\t%d\t%d\t%d\t%d\t%d\t%d\n", 
			  (short unsigned int)curr_frame + 1, 
			  (short unsigned int)curr_obj + 1, 
			  (short unsigned int)*(object_size_short),
			  (short unsigned int)*(centroids_x + curr_obj),
			  (short unsigned int)*(centroids_y + curr_obj),
			  (short unsigned int)*(object_border + curr_obj),
			  (short unsigned int)(object_val[curr_set][curr_obj] / *(object_size + curr_obj)));
			  }*/
      }

#ifdef DEBUG
	fprintf(stdout, "Computing background value for current frame\n");
#endif
      bkgrnd_val = (unsigned short int) (*(bkgrnd_sum + curr_set) / *(bkgrnd_count + curr_set));
#ifdef DEBUG
	fprintf(stdout, "Writing background value for current frame to file\n");
#endif
      if((success=fwrite(&bkgrnd_val, field_size, 1, *(bkgrndfile + curr_set))) < 1){
	fprintf(stdout, "ERROR writing bkgrnd_val to bkgrnd file %4.0d", curr_set);
	exit(42);
      }

      beadmean_val = (unsigned short int) (*(beadmean_sum + curr_set) / 1000000);
      if((success=fwrite(&beadmean_val, field_size, 1, *(beadmeanfile + curr_set))) < 1){
	fprintf(stdout, "ERROR writing beadmean_val to beadmean file %4.0d", curr_set);
	exit(42);
      }

      //fprintf(*(bkgrndfile + curr_set), "%d\t%d\n", curr_frame, (short unsigned int) (*(bkgrnd_sum+curr_set) / *(bkgrnd_count + curr_set)));
    }//end for(current_set)
    
    //
    //**********************************************************************
    
    //******************* Free frame-specific variables ********************
    //
#ifdef DEBUG
    fprintf(stdout, "Freeing frame-specific variables\n");
#endif
    free(centroids_x);
    free(centroids_y);
    free(object_size);
    free(object_border);
    for(curr_set = 0; curr_set < num_sets; curr_set++){
      free(object_val[curr_set]);
    }
    //
    //**********************************************************************
    
  }//end for(current_frame)
  

#ifdef DEBUG
      fprintf(stdout, "Closing all input and output files\n");
#endif
  //***************** Close all files (input & output) *******************
  //
  for(curr_set = 0; curr_set < num_sets; curr_set++){
    fclose(*(infile + curr_set));    
    fclose(*(outfile + curr_set));
    fclose(*(bkgrndfile + curr_set));
    fclose(*(beadmeanfile + curr_set));
  }
  //
  //**********************************************************************

  free(raw_image);
  free(object_image);
}



int scandir_sel(struct dirent * ent){
  char* FILEMASK = "????.tif";
  if(1)
    return(!fnmatch(FILEMASK, ent->d_name, 0));
  return(0);
}


void load_image(char FILENAME[MAX_FULLPATH_LEN], char* img){
  TIFF *image;
  tsize_t stripSize;
  unsigned long imageOffset, result;
  int stripMax, stripCount;
  
#ifdef DEBUG
  fprintf(stdout, "Opening file %s to load image into memory.\n", FILENAME);
#endif
  
  // Open the TIFF image
  if((image = TIFFOpen(FILENAME, "r")) == NULL){
    fprintf(stdout, "Could not open image %s\n", FILENAME);
    exit(42);
  }
  
  stripSize = TIFFStripSize (image);
  stripMax = TIFFNumberOfStrips (image);
  imageOffset = 0;
  
  // Load data to buffer
  for (stripCount = 0; stripCount < stripMax; stripCount++){
    if((result = TIFFReadEncodedStrip (image, stripCount,
				       img + imageOffset,
				       stripSize)) == -1){
      fprintf(stdout, "Read error on input strip number %d\n", stripCount);
      exit(42);
    }
    imageOffset += result;
  }
  // Close image
  TIFFClose(image);   
#ifdef DEBUG
  fprintf(stdout, "File %s loaded successfully.\n", FILENAME);
#endif
  

}