/* ------------------------------------------------------------------------- */
/* EZ-SCROLL: A Code for Simulating Scroll Waves Version 1.5                 */
/* Copyright (C) 1996 - 1998, 2006, 2007 Dwight Barkley and Matthew Dowle    */
/* D.Barkley@warwick.ac.uk                                                   */
/*                                                                           */
/* EZView: 2010, Vadim Biktashev vnb@liv.ac.uk                               */
/* ------------------------------------------------------------------------- */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <tiffio.h>
#include <time.h>
#include <math.h>
#include "ezscroll.h"

/* 
 * Global variables used throughout the EZ-Scroll Code (see ezscroll.h)
 * -------------------------------------------------------------------- */

Real *fields; 
Real  umin, umax, uc, vmin, vmax, vc, wmin, wmax, wc, alphamin, alphamax, grid_h, dt;
int   m=-1, mmin, mmax, mstep, nx, ny, nz, norm_res, rot_res, 
  color_mode, remove_backs, clipping, show_filament, show_surface,
  write_filament, save_images, write_history, hist_x, hist_y, hist_z, 
  verbose, field_size; 
int layer[3]={0,1,2};
float theta, phi;
char *fmt, *template;

/* Global variables for this file only
 * ----------------------------------- */

static  FILE *history_file, *filament_file; 

/* Private functions 
 * ----------------- */

static void  Initialize      (char *task);
static void  Allocate_memory (void);
static int   Read_ez         (char *template,int m);
static int   Read_tiff       (char *template,int m);
static int   Read_ppm        (char *template,int m);
/* static void  Write_history   (int m); */
static void  Finish_up       (void);

#define NEXT_LINE(fp) while(getc(fp)!='\n');  /* macro used to skip to end 
						 of input line */

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

int main(int argc, char **argv) {
  char *ptr, *argv0, *task; 

  /* Check args if any and define the task file names */
  argv0 = (ptr=strrchr(argv[0],'/'))?(ptr+1):(argv[0]);
  if ( ((argc>1) &&  (0!=strchr(argv[1],'?')))
       || ((argc>1) && (0==strncmp(argv[1],"-h",2)))
       || ((argc>1) && (0==strncmp(argv[1],"--h",3)))
       || (argc>2)
       ) { 
    fprintf(stderr,"Usage: %s [task_file (task.dat)]\n",argv[0]);
    exit(1);
  }
  task	= (argc>1) ? argv[1] : "task.dat";
  
  Initialize(task);

  /* The main loop, modelled on the simulation loop of EZSCROLL */
  /* NB m, mmin, mmax, mstep all can be changed in dialogue     */
  /* Uncomment the version you prefer */
  /* #define ENDRANGE mstep *= -1; write_history=write_filament=save_images=0 */
  #define ENDRANGE Pause()
  for (m=mmin;;m+=mstep) {
    if (verbose) printf("m=%d\n",m);
    if (m>mmax) {m=mmax; ENDRANGE; }
    if (m<mmin) {m=mmin; ENDRANGE; }
    if (0==Read_fmt(fmt,template,m)) break;
    Draw();
    if (write_history) Write_history(m);
    if (0!=Event_check()) break;
  }

  Finish_up();
  return 0;
}
/* ========================================================================= */

void Initialize (char *task) {
  /* Read task file, open output files, and initialize graphics and time
   * stepping.  */

  FILE *fp;
  #define BUFLEN 4096
  char buf[BUFLEN], *p;

  /* ---------------------------------- 
   * Read parameters from the task file
   * ---------------------------------- */

  if((fp=fopen(task,"r"))==NULL) {
    fprintf(stderr,"Cannot open task file %s\n",task);
    exit(1);
  } else {
    fgets(buf,BUFLEN,fp);p=strtok(buf," \t\r\b\n"); fmt=strdup(p);
    fgets(buf,BUFLEN,fp);p=strtok(buf," \t\r\b\n"); template=strdup(p);
    fscanf(fp,"%d,%d,%d",	&mmin,&mstep,&mmax);		NEXT_LINE(fp);
    fscanf(fp,"%d,%d,%d",	&nx,&ny,&nz);			NEXT_LINE(fp);
    fscanf(fp,"%d,%d",		&norm_res,&rot_res);		NEXT_LINE(fp);
    fscanf(fp,"%d,%f,%f,%f",	&ulayer,&umin,&umax,&uc);	NEXT_LINE(fp);
    fscanf(fp,"%d,%f,%f,%f",	&vlayer,&vmin,&vmax,&vc);	NEXT_LINE(fp);
    fscanf(fp,"%d,%f,%f,%f",	&wlayer,&wmin,&wmax,&wc);	NEXT_LINE(fp);
    fscanf(fp,"%d,%f,%f",	&color_mode,&alphamin,&alphamax);NEXT_LINE(fp);
    fscanf(fp,"%f,%f",		&theta,&phi);			NEXT_LINE(fp);
                                                                NEXT_LINE(fp);
    fscanf(fp,"%d",		&remove_backs);   		NEXT_LINE(fp);
    fscanf(fp,"%d",		&clipping);   			NEXT_LINE(fp);
    fscanf(fp,"%d",		&show_surface);   		NEXT_LINE(fp);
    fscanf(fp,"%d",		&show_filament);   		NEXT_LINE(fp);
    fscanf(fp,"%d",		&save_images);   		NEXT_LINE(fp);
    fscanf(fp,"%d",		&write_filament);   		NEXT_LINE(fp);
    fscanf(fp,"%d",		&write_history);	        NEXT_LINE(fp);
    fscanf(fp,"%d,%d,%d",	&hist_x,&hist_y,&hist_z);	NEXT_LINE(fp);
                                                                NEXT_LINE(fp);
    fscanf(fp,"%f,%f",		&grid_h,&dt);    		NEXT_LINE(fp);
                                                                NEXT_LINE(fp);
    fscanf(fp,"%d",		&verbose);
    fclose(fp);
  }

  /* ------------ 
   * Final things
   * ------------ */

  Allocate_memory();
  Read_fmt(fmt,template,mmin); 
  if(write_filament) filament_file = fopen("filament.dat", "w"); 
  Draw_ini();
}
/* ========================================================================= */

void Allocate_memory (void)
{
  /* -----------------------------------------------------------------------
   * There are NX, NY and NZ actual grid points in each direction.
   * field_size = (NX+2)*(NY+2)*(NZ+2) because of the way boundaries are
   * treated.  See ezscroll.h for definition of NX, NY, NZ.
   *                                                                         
   * I try to allocate one big block of memory for everything (u, v, w).
   * ----------------------------------------------------------------------- */

  field_size = (NZ+2)*(NY+2)*(NX+2);
  fields = (Real *) calloc(3*field_size, sizeof(Real));

  if (fields != NULL) {
    /* --------------------------------------------------- 
     * We have one big block of memory.  Set some pointers 
     * --------------------------------------------------- */
    if(verbose>2) printf("Memory allocated in one block.\n");
  }
  else {
    /* Could not find enough memory.
     * ----------------------------- */
    fprintf(stderr, "\n ****Error: Not enough memory available\n");
    exit(-1);
  }
}
/* ========================================================================= */

void Finish_up (void)
{
  QuitX();
  if(write_history) fclose(history_file);
  if(write_filament) fclose(filament_file); 
}
/* ========================================================================= */

void Write_filament_data (float x0, float y0, float z0, 
			  float x1, float y1, float z1) 
{
  /* This routine for outputting the filament data is called from
   * filament_in_cube() in ezmarching.c where the filament is computed.  */

  if (!filament_file) filament_file = fopen("filament.dat", "w");
  if (!filament_file) return;
  if (m<0) return;
  fprintf(filament_file,
	  "%d %.5f %.5f %.5f %.5f %.5f %.5f\n", 
	  m, x0, y0, z0, x1, y1, z1);
}
/* ========================================================================= */

int Read_fmt (char *fmt,char *template,int m) {
  if (0==strcasecmp(fmt,"ez"))   return Read_ez (template,m);
  if (0==strcasecmp(fmt,"tiff")) return Read_tiff (template,m);
  if (0==strcasecmp(fmt,"ppm"))  return Read_ppm (template,m);
  fprintf(stderr,"Unknown format %s\n", fmt);
  return 0;
}
/* ========================================================================= */

void Write_history (int wrt_step)
{
  if (!history_file) history_file = fopen("history.dat", "w");
  fprintf(history_file, "%.5f %.5f %.5f\n", dt*(Real)wrt_step, 
	  Fields(U_FIELD,hist_x,hist_y,hist_z), Fields(V_FIELD,hist_x,hist_y,hist_z) ); 

  if(verbose>2) {
    printf("istep = %d", wrt_step);
    printf("; (u,v) = %.5f, %.5f", 
	   Fields(V_FIELD,hist_x,hist_y,hist_z), Fields(V_FIELD,hist_x,hist_y,hist_z) ); 
    printf("\n");
  }
}
/* ========================================================================= */

/* The standard EZSCROLL format for ic.dat and fc.dat files */
int Read_ez (char *template,int m) {
  double u_in, v_in;
  Real   *u_tmp, *v_tmp, rx, ry, rz;
  int    nx_ic, ny_ic, nz_ic, npts_ic2, npts_ic3, index, i_tmp, j_tmp, k_tmp, 
         i, j, k;
  char   f_type, dummy;
  FILE  *fp;
  char filename[4096];

  sprintf(filename,template,m);
  fp=fopen(filename,"r");
  if (fp==NULL) {
    fprintf(stderr,"could not open %s\n",filename);
    return 0;
  }

  /* Read nx_ic etc following = sign on second line of file */
  NEXT_LINE(fp);                                     
  while( (dummy=getc(fp)) != '=');                   
  fscanf(fp, "%d, %d, %d", &nx_ic, &ny_ic, &nz_ic);  

  /* Skip to 10th line and read first character to determine type 
     B(inary) or A(scii) */
  for(i=0;i<8;i++) NEXT_LINE(fp); 
  f_type = getc(fp); NEXT_LINE(fp); 
  
  if ( (f_type !='B') && (f_type !='A') ) {
    fprintf(stderr,"\n ic.dat exists but of unrecognized type Binary or Ascii \n"); 
    exit(1);
  }

  if(verbose>1) printf("\nReading ic.dat with nx, ny, nz = %d, %d, %d... \n", 
		     nx_ic, ny_ic, nz_ic);

  npts_ic2 = nx_ic * ny_ic;  
  npts_ic3 = npts_ic2 * nz_ic;

  /* Allocate temporary memory and read from file */

  u_tmp =(Real *) malloc((unsigned)(npts_ic3)*sizeof(Real));
  v_tmp =(Real *) malloc((unsigned)(npts_ic3)*sizeof(Real));

  if(f_type =='B') {    
    /* Binary data file */
    fread(u_tmp, sizeof(Real), npts_ic3, fp);
    fread(v_tmp, sizeof(Real), npts_ic3, fp);
  }
  else {          
    /* Ascii data file */
    for(index=0;index<npts_ic3;index++) {
      fscanf (fp, "%lg %lg\n", &u_in, &v_in);
      u_tmp[index] = u_in; v_tmp[index] = v_in;
    }
  }

  /* Copy into u and v */
  rx = (nx_ic-1.0)/(NX-1.0);
  ry = (ny_ic-1.0)/(NY-1.0);
  rz = (nz_ic-1.0)/(NZ-1.0);

  for(k=1;k<=NZ;k++) {
    k_tmp = npts_ic2 * (int)(rz*(k-1));
    for(j=1;j<=NY;j++) {
      j_tmp = nx_ic * (int)(ry*(j-1));
      for(i=1;i<=NX;i++) {
	i_tmp = (int)(rx*(i-1));
	Fields(U_FIELD,i,j,k) = u_tmp[i_tmp + j_tmp + k_tmp];
	Fields(V_FIELD,i,j,k) = 
	Fields(W_FIELD,i,j,k) = v_tmp[i_tmp + j_tmp + k_tmp];
      }
    }
  }

  free(u_tmp);
  free(v_tmp);  
  return 1;
}
/* ========================================================================= */

/* Series of tiff files, one per each section, in one directory */
int Read_tiff (char *template,int m) {
  TIFF *tif;
  char dirname[4096];
  char filename[4096];
  uint32 wdh,hgt;
  size_t npixels;
  uint32 *raster;
  unsigned char *bytes;
  double quant=1.0/255.0, rx, ry;
  int x,y,z,i,j;
  
  sprintf(dirname,"%s/%%02d.tif",template);
  raster=NULL;
  for (z=1;z<=nz;z++) {
    sprintf(filename,dirname,m,z-1);
    tif=TIFFOpen(filename,"r");
    if (tif==NULL) {
      fprintf(stderr,"cannot open %s\n",filename); 
      break;
    }
    if (raster==NULL || bytes==NULL) {
      TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &wdh);
      TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &hgt);
      npixels = wdh * hgt;
      raster = (uint32*) _TIFFmalloc(npixels * sizeof (uint32));
      if (raster==NULL) {
	fprintf(stderr,"not enough memory in Read_tiff\b"); 
	break;
      }
      bytes = (unsigned char *)raster;
    }
    if (0==TIFFReadRGBAImage(tif, wdh, hgt, raster, 0)) {
      fprintf(stderr,"could not read from %s\n",filename);
      break;
    }
    rx = (wdh-1.0)/(NX-1.0);
    ry = (hgt-1.0)/(NY-1.0);
    for (y=1;y<=NY;y++) {
      j = (int)(ry*(y-1));
      for (x=1;x<=NX;x++) {
	i = (int)(rx*(x-1));
	Fields(0,x,y,z) = quant*bytes[0+4*(i+wdh*j)];
	Fields(1,x,y,z) = quant*bytes[1+4*(i+wdh*j)];
	Fields(2,x,y,z) = quant*bytes[2+4*(i+wdh*j)];
      }
    }
    TIFFClose(tif);
  }
  _TIFFfree(raster);
  return (z>nz);
}
/* ========================================================================= */

static unsigned char *bytes;
static int wdh, hgt, dph, max;
static inline unsigned char corrbyte(int i, int j, int k,  int l) {
  unsigned char b=bytes[l+3*(i+wdh*(j+hgt*k))];
  if (l==wlayer) return 255-b;
  if (b==255) return 0;
  return b;
}
static inline long int shift(int i, int j,  int k) {
  return (i+wdh*(j+hgt*k));
}

/* A binary ppm header followed by nz rasters */
/* Blue channel is assumed to contain tissue description and is subject to dilation */
int Read_ppm (char *template,int m) {
  FILE *ppmin;
  char filename[4096], tag[16], nl;
  Real quant, rx, ry, rz;
  long int firstbyte, lastbyte, nvoxels, npixels;
  int i,j,k,x,y,z;

  sprintf(filename,template,m);
  ppmin=fopen(filename,"r");
  if (ppmin==NULL) {
    fprintf(stderr,"could not open %s\n",filename);
    return 0;
  }
  fscanf(ppmin,"%2s ",tag);
  if (0!=strcmp(tag,"P6")) {
    fprintf(stderr,"%s is not a ppm file\n",filename);
    return 0;
  }
  if (!fscanf(ppmin,"%u ",&wdh)
   || !fscanf(ppmin,"%u ",&hgt)
   || !fscanf(ppmin,"%u",&max)
      ) {
    fprintf(stderr,"could not read dimensions from %s\n",filename);
    return 0;
  }
  quant=1.0/max;
  fscanf(ppmin,"%c",&nl);
  firstbyte=ftell(ppmin);
  fseek(ppmin,0L,SEEK_END);
  lastbyte=ftell(ppmin);
  nvoxels=lastbyte-firstbyte;
  fseek(ppmin,firstbyte,SEEK_SET);
  npixels = 3*wdh*hgt;
  if (0!=(nvoxels%npixels)) {
    fprintf(stderr,"%s nvoxels=%ld not divisible by npixels=%ld\b",filename,nvoxels,npixels);
    return 0;
  }
  dph=nvoxels/npixels;
  bytes=(unsigned char *)calloc(nvoxels,1);
  if (bytes==NULL) {
    fprintf(stderr,"not enough memory in Read_ppm\n");
    return 0;
  }
  if (nvoxels != fread(bytes,1,nvoxels,ppmin)) {
    fprintf(stderr,"error reading from %s\n",filename);
    return 0;
  }

#if 1
  /* Now the arrayes are in place and we shall dilate and invert */
  /* the tissue channel */
  {
    int arysize=(wdh*hgt*dph);
    int dk=(wdh*hgt);
    int dj=wdh;
    int di=1;
    int theshift;
    unsigned char thebyte;
    unsigned char *ary=(unsigned char *)calloc(arysize,1);
    if (ary==NULL) {
      fprintf(stderr,"not enough memory for aux array; dilation not done\n");
    } else {
      for (k=0;k<dph;k++) {
	for (j=0;j<hgt;j++) {
	  for (i=0;i<wdh;i++) {
	    theshift=shift(i,j,k);
	    ary[theshift]=bytes[wlayer+3*theshift];
	  }
	}
      }
      for (k=1;k<dph-1;k++) {
	for (j=1;j<hgt-1;j++) {
	  for (i=1;i<wdh-1;i++) {
	    theshift=shift(i,j,k);
	    thebyte=ary[theshift];
	    if (255==ary[theshift+di]) thebyte=255;
	    if (255==ary[theshift-di]) thebyte=255;
	    if (255==ary[theshift+dj]) thebyte=255;
	    if (255==ary[theshift-dj]) thebyte=255;
	    if (255==ary[theshift+dk]) thebyte=255;
	    if (255==ary[theshift-dk]) thebyte=255;

	    if (255==ary[theshift+di+dj]) thebyte=255;
	    if (255==ary[theshift-di+dj]) thebyte=255;
	    if (255==ary[theshift+di-dj]) thebyte=255;
	    if (255==ary[theshift-di-dj]) thebyte=255;

	    if (255==ary[theshift+di+dk]) thebyte=255;
	    if (255==ary[theshift-di+dk]) thebyte=255;
	    if (255==ary[theshift+di-dk]) thebyte=255;
	    if (255==ary[theshift-di-dk]) thebyte=255;

	    if (255==ary[theshift+dj+dk]) thebyte=255;
	    if (255==ary[theshift-dj+dk]) thebyte=255;
	    if (255==ary[theshift+dj-dk]) thebyte=255;
	    if (255==ary[theshift-dj-dk]) thebyte=255;

	    if (255==ary[theshift+di+dj+dk]) thebyte=255;
	    if (255==ary[theshift-di+dj+dk]) thebyte=255;
	    if (255==ary[theshift+di-dj+dk]) thebyte=255;
	    if (255==ary[theshift-di-dj+dk]) thebyte=255;
	    if (255==ary[theshift+di+dj-dk]) thebyte=255;
	    if (255==ary[theshift-di+dj-dk]) thebyte=255;
	    if (255==ary[theshift+di-dj-dk]) thebyte=255;
	    if (255==ary[theshift-di-dj-dk]) thebyte=255;

	    bytes[wlayer+3*theshift]=thebyte;
	  }
	}
      }
    }
  }
#endif

  rx = (wdh-1.0)/(NX-1.0);
  ry = (hgt-1.0)/(NY-1.0);
  rz = (dph-1.0)/(NZ-1.0);
  for (z=1;z<=NZ;z++) {
    k = (int)(rz*(z-1));
    for (y=1;y<=NY;y++) {
      j = (int)(ry*(y-1));
      for (x=1;x<=NX;x++) {
	i = (int)(rx*(x-1));
	Fields(0,x,y,z) = quant*corrbyte(i,j,k,0);
	Fields(1,x,y,z) = quant*corrbyte(i,j,k,1);
	Fields(2,x,y,z) = quant*corrbyte(i,j,k,2);
      }
    }
  }
  fclose(ppmin);
  free(bytes);
  return 1;
}
/* ========================================================================= */