/***************************************************************************************
 * find_objects.c                                                                      *
 *                                                                                     *
 * Takes a raw transmitted light image, finds all objects, and outputs a new image     *
 * file where every object from the original image is numbered sequentially (same      *
 * as bwlabel in Matlab) in the 'object' file.                                         *
 *                                                                                     *
 * TO COMPILE: g++ find_objects.c -O3 -ffloat-store -foptimize-sibling-calls           * 
 *             -fno-branch-count-reg -o find_objects -ltiff -lm -I../../INCLUDE        * 
 *             -L../../LIB                                                             *
 *                                                                                     *
 * USAGE:                                                                              *
 *       ./find_objects xcols yrows base_dir/ cutoff {-1,1}                            *
 *       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/)                      *
 *             cutoff = the threshold value (e.g. 9000)                                *
 *             {-1, 1} = whether bead val is 0 or 65535 (black beads = -1, white = 1)  *
 *                                                                                     *
 * OUTPUT goes to base_dir/IMAGES/OBJECT/obj_counts.dat and to base_dir/IMAGES/OBJECT/ *
 *                                                                                     *
 * WRITTEN BY: Jay Shendure, Greg Porreca (Church Lab) 01-30-2005                      *
 *                                                                                     *
 *                                                                                     *
 ***************************************************************************************/

#include <stdio.h>
#include <tiffio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <dirent.h>
#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 PATH_FULL_RAW "IMAGES/RAW/999/"
#define PATH_FULL_CONV "IMAGES/CONV/"
#define PATH_RAW "IMAGES/RAW/"
#define PATH_CONV "IMAGES/CONV/"
#define MAX_FILESIZE 2252800 //2200KB

//function prototypes
void load_image(char*, char*);
int save_image(char*, char*);
int scandir_sel(struct dirent * ent);

//global variables
unsigned short int yrows, xcols;
char ROOT_PATH[MAX_ROOTPATH_LEN];
int STRIPSIZE = 8000;
int STRIPMAX = 250;  


int main(int argc, char *argv[]){
  
  char *raw_image;//, *conv_image;
  short unsigned int *conv_image_short;

  long int *conv_image;

  /*int kernel[][] = {{0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0},
		    {0,  0,  0,  0, -5, -5, -5,  0,  0,  0,  0},
		    {0,  0, -3, -5, -3, -3, -3, -5, -3,  0,  0},
		    {0, -3, -5, -3, -1,  1, -1, -3, -5, -3,  0},
		    {0, -5, -5, -1,  4,  8,  4, -1, -5, -5,  0},
		    {0, -5, -5,  1,  8,  8,  8,  1, -5, -5,  0},
		    {0, -5, -5, -1,  4,  8,  4, -1, -5, -5,  0},
		    {0, -3, -5, -3, -1,  1, -1, -3, -5, -3,  0},
		    {0,  0, -3, -5, -3, -3, -3, -5, -3,  0,  0},
		    {0,  0,  0,  0, -5, -5, -5,  0,  0,  0,  0},
		    {0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0}};
  */

  int kernel[] = {0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
		  0,  0,  0,  0, -5, -5, -5,  0,  0,  0,  0,
		  0,  0, -3, -5, -3, -3, -3, -5, -3,  0,  0,
		  0, -3, -5, -3, -1,  1, -1, -3, -5, -3,  0,
		  0, -5, -5, -1,  4,  8,  4, -1, -5, -5,  0,
		  0, -5, -5,  1,  8,  8,  8,  1, -5, -5,  0,
		  0, -5, -5, -1,  4,  8,  4, -1, -5, -5,  0,
		  0, -3, -5, -3, -1,  1, -1, -3, -5, -3,  0,
		  0,  0, -3, -5, -3, -3, -3, -5, -3,  0,  0,
		  0,  0,  0,  0, -5, -5, -5,  0,  0,  0,  0,
		  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0};

  int kernel_len = 121;
  int kernel_size = 11;
  int kernel_precalc_indx[] = {0, -5, -3, -1, 1, 4, 8};
  int kernel_precalc_size = 7;
  long int* kernelimg_precalc[kernel_precalc_size];
  int curr_x, curr_y, curr_kernel_x, curr_kernel_y;

  int *kernel_precalc;

  long int min_val, max_val;
  double curr_pixel_val;
  
  //--------------pathname generation
  char raw_image_dir[MAX_FULLPATH_LEN];
  unsigned short int current_frame;
  char curr_raw_fn[MAX_FULLPATH_LEN];
  char curr_conv_fn[MAX_FULLPATH_LEN];
  char curr_fn_index[5];
  int num_frames = 0;
  struct dirent **namelist;
  int i, j;
  int FLAG_edge;
  int curr_pixel;


  yrows = atoi(argv[1]);
  xcols = atoi(argv[2]);
  strcpy(ROOT_PATH, argv[3]);

  
  
  if((kernel_precalc = (int*) malloc(kernel_len * sizeof(int))) == NULL){
    fprintf(stderr, "Could not allocate enough memory for the kernel.\n");
    exit(42);
  }

  for(i=0; i< kernel_precalc_size; i++){
    fprintf(stderr, "allocating for kernelimg_precalc[%d] %ld bytes\n", i,xcols * yrows * sizeof(long int)); 
    if((kernelimg_precalc[i] = (long int*) malloc(xcols * yrows * sizeof(long int)))==NULL){
      fprintf(stderr, "Could not allocate memory for kernel precalc image %d.\n", i);
      exit(42);
    }
  }

  for(curr_pixel = 0; curr_pixel < xcols*yrows; curr_pixel++){
    for(i=0; i< kernel_precalc_size; i++){
      *(kernelimg_precalc[i] + curr_pixel) = 0;
    }
  }

  for(i=0; i<kernel_len; i++){
    for(j=0; j< kernel_precalc_size; j++){
      if(*(kernel_precalc_indx + j) == *(kernel + i)){ 
	*(kernel_precalc + i) = j;
      }
    }
  }

  //turn annoying TIFF warnings off
  TIFFSetWarningHandler(0);
  
  //determine number of RAW files to process
  sprintf(raw_image_dir, "%s%s", ROOT_PATH, PATH_FULL_RAW);
  num_frames = scandir(raw_image_dir, &namelist, (int (*)(const struct dirent *))(scandir_sel), alphasort);
  free(namelist);
  
  //allocate memory for raw image
  if((raw_image = (char *) malloc(MAX_FILESIZE)) == NULL){
    fprintf(stderr, "Could not allocate enough memory for the raw image.\n");
    exit(42);
  }
  
  //allocate memory for convolved image
  if((conv_image = (long int *) malloc(xcols * yrows * sizeof(long int))) == NULL){
    fprintf(stderr, "Could not allocate enough memory for the convolved image.\n");
    exit(42);
  }

  //allocate memory for convolved image
  if((conv_image_short = (short unsigned int *) malloc(xcols * yrows * sizeof(short unsigned int))) == NULL){
    fprintf(stderr, "Could not allocate enough memory for the convolved image.\n");
    exit(42);
  }
  
  //  for(current_frame = 1; current_frame < num_frames + 1; current_frame++){
  for(current_frame = 1; current_frame < 2; current_frame++){
    
    //current image index, as a 4-digit string
    sprintf(curr_fn_index, "%04d", current_frame);
    
    //generate complete filename for current raw image
    strcpy(curr_raw_fn, "");
    strcat(curr_raw_fn, ROOT_PATH);
    strcat(curr_raw_fn, PATH_FULL_RAW);
    strcat(curr_raw_fn, "WL_");
    strcat(curr_raw_fn, curr_fn_index);
    strcat(curr_raw_fn, ".tif");
    
    //generate complete filename for current convolved image
    strcpy(curr_conv_fn, "");
    strcat(curr_conv_fn, ROOT_PATH);
    strcat(curr_conv_fn, PATH_FULL_CONV);
    strcat(curr_conv_fn, curr_fn_index);
    strcat(curr_conv_fn, ".tif");
    
    //allocate memory and load tiff
    load_image(curr_raw_fn, raw_image);
    
    //generate n kernel precalc images (n = # of distinct vals in the kernel)
    for(curr_pixel = 0; curr_pixel < xcols * yrows; curr_pixel++){

      for(j=0; j< kernel_precalc_size; j++){
	*(kernelimg_precalc[j] + curr_pixel) = (long int)*(raw_image + curr_pixel) * (long int)(*(kernel_precalc_indx + j)); 
      }
    }

    min_val = *(raw_image);
    max_val = *(raw_image);

    //generate the convolved image using the kernel precalc images (only need to take sums now)
    for(curr_pixel = 0; curr_pixel < xcols * yrows; curr_pixel++){

      curr_x = curr_pixel % xcols;
      curr_y = (int)floor(float(curr_pixel / xcols));

      if((curr_x < ((kernel_size/2)+1)) || (curr_x > (xcols - ((kernel_size/2)+1))) || (curr_y < ((kernel_size/2)+1)) || (curr_y > (yrows - ((kernel_size/2)+1)))){
	FLAG_edge = 1;
      }
      else{
	FLAG_edge = 0;
      }


      if(!FLAG_edge){
	for(i=0; i< kernel_len; i++){
	  curr_kernel_x = (curr_pixel % xcols) + (i % kernel_size) - 5;
	  curr_kernel_y = (int) floor(float(curr_pixel / xcols)) + (int) floor(float(i / kernel_size)) - 5;
	  //fprintf(stderr, "curr_x%d curr_y%d curr_kernel_x%d, curr_kernel_y%d, kernelimg%ld\n", curr_x, curr_y, curr_kernel_x, curr_kernel_y,  *(kernelimg_precalc[*(kernel_precalc + i)] + (curr_kernel_y*xcols) + curr_kernel_x));
	  *(conv_image + curr_pixel) = *(conv_image + curr_pixel) + *(kernelimg_precalc[*(kernel_precalc + i)] + (curr_kernel_y*xcols) + curr_kernel_x);
	  if( *(conv_image + curr_pixel) > max_val){
	    max_val = *(conv_image + curr_pixel);
	  }
	  else if( *(conv_image + curr_pixel) < min_val){
	    min_val = *(conv_image + curr_pixel);
	  }
	}
      }
      else{
	*(conv_image + curr_pixel) = *(raw_image + curr_pixel);

	if( *(conv_image + curr_pixel) > max_val){
	  max_val = *(conv_image + curr_pixel);
	}
	else if( *(conv_image + curr_pixel) < min_val){
	  min_val = *(conv_image + curr_pixel);

	}
      }
      //      fprintf(stdout, "%ld\n", (long int) *(conv_image + curr_pixel));
    }
    
    fprintf(stderr, "Frame %d convolved.\n", current_frame);
    for(i=0; i< 11; i++){
      for(j=0; j<11; j++){
	fprintf(stderr, "%d\t",*(kernel_precalc + (i * 11) + j));
      }
      fprintf(stderr, "\n");
    }

    for(curr_pixel = 0; curr_pixel < xcols * yrows; curr_pixel++){
      curr_pixel_val = *(conv_image + curr_pixel) - min_val;
      curr_pixel_val = curr_pixel_val / (max_val - min_val);
      curr_pixel_val = curr_pixel_val * 16383;
      *(conv_image_short + curr_pixel) = (short unsigned int) curr_pixel_val;
    }

    //save convolved image just created
        save_image(curr_conv_fn, (char*)conv_image_short);
    
    //free memory
  }
  free(raw_image);
  free(conv_image);
  exit(0);
}  


//make scandir only return files matching pattern FILEMASK
int scandir_sel(struct dirent * ent){
  char* FILEMASK = "*.tif";
  if(1)
    return(!fnmatch(FILEMASK, ent->d_name, 0));
  return(0);
}



int save_image(char *save_image_fn, char *save_image_data){
  TIFF *image2;
  tsize_t stripSize;
  unsigned long imageOffset, result;
  int stripMax, stripCount;
  unsigned long count;
  
  // Open the output TIFF image
  if((image2 = TIFFOpen(save_image_fn, "w")) == NULL){
    fprintf(stderr, "Could not open outgoing image\n");
    exit(42);
  }
  
  // We need to set some values for basic tags before we can add any data
  // Hard-coding this for now, ideally it would be flexible to different image sizes
  TIFFSetField(image2, TIFFTAG_IMAGEWIDTH, xcols);
  TIFFSetField(image2, TIFFTAG_IMAGELENGTH, yrows);
  TIFFSetField(image2, TIFFTAG_BITSPERSAMPLE, 16);
  TIFFSetField(image2, TIFFTAG_ROWSPERSTRIP, 4);
  TIFFSetField(image2, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK);
  TIFFSetField(image2, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
  TIFFSetField(image2, TIFFTAG_XRESOLUTION, 72.0);
  TIFFSetField(image2, TIFFTAG_YRESOLUTION, 72.0);
  TIFFSetField(image2, TIFFTAG_RESOLUTIONUNIT, RESUNIT_INCH);
  
  // Read in the number and size of strips from the other file
  stripSize = STRIPSIZE;
  stripMax = STRIPMAX;
  //  stripSize = TIFFStripSize (image1);
  //stripMax = TIFFNumberOfStrips (image1);
  imageOffset = 0;
 
  // Write the data
  for (stripCount = 0; stripCount < stripMax; stripCount++){
    if((result = TIFFWriteEncodedStrip (image2, stripCount, save_image_data + imageOffset, stripSize)) == -1){
      fprintf(stderr, "Write error on input strip number %d\n", stripCount);
      exit(42);
    }
    imageOffset += result;
  }
  
  //TIFFClose(image1);
  TIFFClose(image2);
  return 0;
}


void load_image(char FILENAME[MAX_FULLPATH_LEN], char* img_data){
  TIFF *image;
  tsize_t stripSize;
  unsigned long imageOffset, result;
  int stripMax, stripCount;
  unsigned long count;
    
  // Open the TIFF image
  if((image = TIFFOpen(FILENAME, "r")) == NULL){
    fprintf(stderr, "Could not open image %s\n", FILENAME);
    exit(42);
  }
  
  stripSize = TIFFStripSize (image);
  stripMax = TIFFNumberOfStrips (image);
  //STRIPSIZE = stripSize;
  //STRIPMAX = stripMax;
  imageOffset = 0;
  
  // Load data to buffer
  for (stripCount = 0; stripCount < stripMax; stripCount++){
    if((result = TIFFReadEncodedStrip (image, stripCount,
				       img_data + imageOffset,
				       stripSize)) == -1){
      fprintf(stderr, "Read error on input strip number %d\n", stripCount);
      exit(42);
    }
    
    imageOffset += result;
  }
  
  // Close image
  TIFFClose(image);   
}