/*plotavg.c - translate a GnuroScan averaged EEG file into PostScript code.
  Copyright (c) 1996 Matthew Belmonte

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
  as published by the Free Software Foundation; either version 2
  of the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

  If you find this program useful, please send mail to Matthew Belmonte.
  <mkb4@Cornell.edu>.  If you base a publication on data processed by this
  program, please notify Matthew Belmonte and include the following citation
  in your publication:

	Matthew Belmonte, `A Software System for Analysis of
	Steady-State Evoked Potentials', Association for Computing
	Machinery SIGBIO Newsletter 17:1:9-14 (April 1997).
*/

#include <stdio.h>
#include <unistd.h>
#include <math.h>
#ifdef NOFLOATS
#  define fceil(x) ((float)(ceil((double)x)))
#  define ffloor(x) ((float)(floor((double)x)))
#endif
#include <errno.h>
#include "intel.h"
#include "gnuro.h"
#include "fastshft.h" /*include this iff you want to plot FASTSHFT data*/

#define CH_MAX 16

int fontsize, num_samples;
float lmargin, rmargin, tmargin, bmargin, page_height, page_width, plot_height, plot_leading, title_height, tick_length;
float xmin, xmax, image_width;

#ifdef ST_MSECS_PER_TICK
  long baseline_position;
  int baseline_length;
# define MAX_BASELINE_LEN (1+(int)(((long)(ST_LONGEST_ISI))*AV_MAX_SRATE/1000L))

# include <math.h>
# include "fmedian.h"

  void bootstrap(data, n, bootstrap_distrib, iter)
  float *data;
  int n;
  float *bootstrap_distrib;
  int iter;
    {
    register double sum;
    register int i;
    register btree *boot;
    boot = create_btree(iter);
    while(iter--)
      {
      sum = 0.0;
      for(i = 0; i != n; i++)
	sum += data[random()%n];
      insert_newest(sum/n, boot);
      }
    traverse_btree(boot, bootstrap_distrib);
    destroy_btree(boot);
    }

  /*Return the number of items in distrib[0..k-1], which is sorted in
    nondecreasing order, that are strictly less than 'x'.*/
  int rank(x, distrib, k)
  float x, *distrib;
  int k;
    {
    register int lb, ub;
    lb = 0;
    ub = k;
    while(lb < ub)
    /*inv: distrib[0..lb-1] < x <= distrib[ub..k-1]*/
      {
      k = (lb+ub)/2;
      if(distrib[k] < x)
	lb = k+1;
      else
	ub = k;
      }
    return(lb);
    }
#endif

void draw_axes(channel, label, ymin, ymax)
int channel;
char *label;
float ymin, ymax;
  {
  float x_axis_elevation, y_axis_elevation, hundreds_of_msecs;
  x_axis_elevation = page_height-tmargin-title_height-channel*plot_height-(ymax<=0.0? 0.0: ymin>=0.0? 1.0: ymax/(ymax-ymin))*(plot_height-plot_leading)-plot_leading/2.0;
  y_axis_elevation = (xmax<=0.0)? page_width-rmargin: (xmin>=0.0)? lmargin: lmargin-(xmin/(xmax-xmin))*image_width;
  printf("newpath\n2 setlinewidth\n0.75 setgray\n");
  printf("%f %f moveto\n", lmargin, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0);
  printf("(%s) show\n", label);
  printf("%f %f moveto\n", lmargin, x_axis_elevation);
  printf("%f 0 rlineto\n", image_width);
  printf("%f %f moveto\n", y_axis_elevation, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0);
  printf("0 %f rlineto\n", plot_leading-plot_height);
  printf("(%.0fuV) show\n", ymin);
  printf("%f %f moveto\n", y_axis_elevation, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0);
  printf("(%.0fuV) show\n", ymax);
  for(hundreds_of_msecs = fceil(10.0*xmin); hundreds_of_msecs != ffloor(10.0*xmax); hundreds_of_msecs += 1.0)
    printf("%f %f moveto\n0 %f rlineto\n", y_axis_elevation+(hundreds_of_msecs/(10.0*(xmax-xmin)))*image_width, x_axis_elevation-tick_length/2.0, tick_length);
  printf("stroke\n0 setgray\n");
  }

void plotpoints(fp, n, scale)
FILE *fp;
int n;
float scale;
  {
  register int sample, num_valid_samples_in_baseline;
  register double cur_y, prev_y;
#ifdef ST_MSECS_PER_TICK
# define BOOTSTRAP_SIZE 1000
  register float cur_p, prev_p, grey_level, prev_grey_level;
  register int significance_flag;
  float baseline[MAX_BASELINE_LEN], baseline_bootstrap[BOOTSTRAP_SIZE];
  fseek(fp, baseline_position*sizeof(float), SEEK_CUR);
  num_valid_samples_in_baseline = 0;
  for(sample = 0; sample != baseline_length; sample++)
    {
    baseline[num_valid_samples_in_baseline] = read_Intel_float(fp);
    if(baseline[num_valid_samples_in_baseline] != 0.0)
      num_valid_samples_in_baseline++;
    }
  fseek(fp, -(baseline_position+baseline_length)*sizeof(float), SEEK_CUR);
  bootstrap(baseline, num_valid_samples_in_baseline, baseline_bootstrap, BOOTSTRAP_SIZE);
  significance_flag = 1;
#endif
  printf("1 setlinewidth\n");
  cur_y = read_Intel_float(fp);
#ifdef ST_MSECS_PER_TICK
  cur_p = rank(cur_y, baseline_bootstrap, BOOTSTRAP_SIZE)/(float)BOOTSTRAP_SIZE;
  grey_level = 0.0;
#endif
  printf("0 %f rmoveto\n", scale*cur_y);
  for(sample = 1; sample != n; sample++)
    {
    prev_y = cur_y;
    cur_y = read_Intel_float(fp);
#ifdef ST_MSECS_PER_TICK
    prev_p = cur_p;
    cur_p = rank(cur_y, baseline_bootstrap, BOOTSTRAP_SIZE)/(float)BOOTSTRAP_SIZE;
    prev_grey_level = grey_level;
    if(significance_flag && !(((prev_p <= 0.025) && (cur_p <= 0.025)) || ((prev_p >= 0.975) && (cur_p >= 0.975))))
      {
      significance_flag = 0;
      grey_level = 0.5;
      }
    else if((prev_p <= 0.025) && (cur_p <= 0.025))
      {
      significance_flag = 1;
      grey_level = 10.0*(prev_p+cur_p);
      }
    else if((prev_p >= 0.975) && (cur_p >= 0.975))
      {
      significance_flag = 1;
      grey_level = 10.0*(2.0-prev_p-cur_p);
      }
    if(grey_level != prev_grey_level)
      printf("currentpoint\nstroke\nmoveto\n%f setgray\n", grey_level);
#endif
    printf("%f %f rlineto\n", image_width/(n-1), (cur_y-prev_y)*scale);
    }
  printf("stroke\n0 setgray\n");
  }

void plotvariance(fp, vfp, n, scale)
FILE *fp, *vfp;
int n;
float scale;
  {
  register int sample;
  register double cur_y, prev_y, cur_vy, prev_vy;
  fseek(fp, (long)(-sizeof(float)*n), SEEK_CUR);
  printf("1 setlinewidth\n0.75 setgray\n");
  cur_y = read_Intel_float(fp);
  cur_vy = read_Intel_float(vfp);
  printf("0 %f rmoveto\n", scale*(cur_y+cur_vy));
  for(sample = 1; sample != n; sample++)
    {
    prev_y = cur_y;
    prev_vy = cur_vy;
    cur_y = read_Intel_float(fp);
    cur_vy = read_Intel_float(vfp);
    printf("%f ", image_width/(n-1));
    printf("%f rlineto\n", (cur_y+cur_vy-prev_y-prev_vy)*scale);
    }
  fseek(fp, (long)(-sizeof(float)*n), SEEK_CUR);
  fseek(vfp, (long)(-sizeof(float)*n), SEEK_CUR);
  printf("%f %f rmoveto\n", -image_width, -cur_y-cur_vy);
  cur_y = read_Intel_float(fp);
  cur_vy = read_Intel_float(vfp);
  printf("0 %f rmoveto\n", scale*(cur_y-cur_vy));
  for(sample = 1; sample != n; sample++)
    {
    prev_y = cur_y;
    prev_vy = cur_vy;
    cur_y = read_Intel_float(fp);
    cur_vy = read_Intel_float(vfp);
    printf("%f ", image_width/(n-1));
    printf("%f rlineto\n", (cur_y-cur_vy-prev_y+prev_vy)*scale);
    }
  printf("stroke\n");
  printf("0 setgray\n");
  }

void find_range(fp, n, minp, maxp, calib)
FILE *fp;
int n;
float *minp, *maxp, calib;
  {
  register int sample;
  register double cur_y;
  for(sample = 0; sample != n; sample++)
    {
    cur_y = read_Intel_float(fp);
    if(cur_y*calib < *minp)
      *minp = cur_y*calib;
    else if(cur_y*calib > *maxp)
      *maxp = cur_y*calib;
    }
  }

void main(argc, argv)
int argc;
char **argv;
  {
  register int channel, file, i;
  int num_channels, variance_flag;
  double min_range, round_step;
  FILE *avg, *avg_tail;
  char erp[packed_sizeof_SETUP], ch_hdr[CH_MAX][packed_sizeof_ELECTLOC];
  float calibration[CH_MAX], yscale[CH_MAX], ymin[CH_MAX], ymax[CH_MAX];

  fprintf(stderr, "Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  fontsize = 12;
  title_height = fontsize;
  lmargin = 72*0.5;
  rmargin = 72*0.5;
  tmargin = 72*0.5;
  bmargin = 72*0.5;
  page_height = 72*11.0;
  page_width = 72*8.5;
  plot_leading = 12.0;
  tick_length = 5.0;
  image_width = page_width-lmargin-rmargin;

  if(argc < 5)
    {
    fprintf(stderr, "usage: %s [-v] <label> <min. range> <rounding step> <file.avg>...<file.avg>\n", *argv);
    exit(1);
    }
  variance_flag = (argv[1][0] == '-') && (argv[1][1] == 'v');
  if(variance_flag)
    {
    argv++;
    argc--;
    }
  sscanf(argv[2], "%lf", &min_range);
  sscanf(argv[3], "%lf", &round_step);
  if((min_range <= 0.0) || (round_step <= 0.0))
    {
    fprintf(stderr, "crazy parameter values (%f, %f)\n", min_range, round_step);
    exit(1);
    }
  argv += 4;
  argc -= 4;
  if((avg = fopen(*argv,
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(*argv);
    exit(errno);
    }
/*read # of channels, # of samples, variance flag, and real time bounds*/
  fread(erp, packed_sizeof_SETUP, 1, avg);
  num_samples = S_pnts(erp);
  num_channels = S_nchannels(erp);
  if(num_channels > CH_MAX)
    {
    fprintf(stderr, "Recompile %s with CH_MAX=%d\n", *argv, num_channels);
    fclose(avg);
    exit(1);
    }
  plot_height = (page_height-tmargin-bmargin)/num_channels;
  variance_flag = variance_flag && S_variance(erp);
  xmin = S_xmin(erp);
  xmax = S_xmax(erp);

#ifdef ST_MSECS_PER_TICK
/******A ONE-TICK BASELINE IS TOO SHORT TO BE PRACTICAL.
  baseline_length = (ST_MSECS_PER_TICK*S_rate(erp)+500)/1000;
  baseline_position = ((int)((xmin/(xmin-xmax))*num_samples-0.5))-baseline_length;
  if(baseline_position < 0L)
    {
    baseline_length += baseline_position;
    baseline_position = 0L;
    }
************/
  baseline_length = (int)((xmin/(xmin-xmax))*num_samples-0.5);
  baseline_position = 0L;
  if(baseline_length > MAX_BASELINE_LEN)
    {
    fclose(avg);
    fprintf(stderr, "recompile with MAX_BASELINE_LENGTH=%d\n", baseline_length);
    exit(1);
    }
#endif

/*read channel labels and calibrations*/
  fseek(avg, (long)packed_sizeof_SETUP, SEEK_SET);
  for(channel = 0; channel != num_channels; channel++)
    {
    fread(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, avg);
    calibration[channel] = EL_calib(ch_hdr[channel]);
    }
/*find ranges of all channels*/
  fseek(avg, (long)(packed_sizeof_SETUP+num_channels*packed_sizeof_ELECTLOC), SEEK_SET);
  for(channel = 0; channel != num_channels; channel++)
    {
    fseek(avg, 5L, SEEK_CUR);
    ymin[channel] = 1e+38;
    ymax[channel] = -1e+38;
    find_range(avg, num_samples, ymin+channel, ymax+channel, calibration[channel]);
    }
  fclose(avg);
  for(file = 1; file < argc; file++)
    {
    if((avg = fopen(argv[file],
#ifdef MSDOS
      "rb"
#else
      "r"
#endif
      )) == NULL)
      {
      perror(argv[file]);
      exit(errno);
      }
    fread(erp, packed_sizeof_SETUP, 1, avg);
    if((num_samples != S_pnts(erp)) || (num_channels != S_nchannels(erp)) || (xmin != S_xmin(erp)) || (xmax != S_xmax(erp)))
      {
      fclose(avg);
      fprintf(stderr, "%s (%d points, %d channels, time bounds %f..%f) doesn't match %s (%d points, %d channels, time bounds %f..%f)\n", argv[file], S_pnts(erp), S_nchannels(erp), S_xmin(erp), S_xmax(erp), *argv, num_samples, num_channels, xmin, xmax);
      exit(1);
      }
    variance_flag = variance_flag && S_variance(erp);
    for(channel = 0; channel != num_channels; channel++)
      {
      int c;
      i = 0;
      c = getc(avg);
      while((i != EL_lab_length-1) && (c != '\0') && (c == ch_hdr[channel][EL_lab_offset+i]))
	{
	i++;
	c = getc(avg);
	}
      if((i != EL_lab_length-1) && (c != ch_hdr[channel][EL_lab_offset+i]))
	{
	fclose(avg);
	fprintf(stderr, "Channel label %d in file %s doesn't match channel label %d (%s) in file %s\n", channel, argv[file], channel, EL_lab(ch_hdr[channel]), *argv);
	exit(1);
	}
      fread(ch_hdr[channel]+i+1, packed_sizeof_ELECTLOC-i-1, 1, avg);
      calibration[channel] = EL_calib(ch_hdr[channel]);
      }
  /*find ranges of all channels*/
    fseek(avg, (long)(packed_sizeof_SETUP+num_channels*packed_sizeof_ELECTLOC), SEEK_SET);
    for(channel = 0; channel != num_channels; channel++)
      {
      fseek(avg, 5L, SEEK_CUR);
      find_range(avg, num_samples, ymin+channel, ymax+channel, calibration[channel]);
      }
    fclose(avg);
    }
/*round ranges to nearest 10.0uV, and set scaling factors*/
  for(channel = 0; channel != num_channels; channel++)
    if(ymin[channel] != ymax[channel])
      {
      ymin[channel] = (float)(floor(ymin[channel]/round_step)*round_step);
      ymax[channel] = (float)(ceil(ymax[channel]/round_step)*round_step);
      if(ymax[channel]-ymin[channel] < min_range)
	ymax[channel] = ymin[channel]+min_range;
      yscale[channel] = (plot_height-plot_leading)/(ymax[channel]-ymin[channel]);
      }
    else
      yscale[channel] = 0.0;
/*set font and print title banner*/
  printf("%%!PS\n/Helvetica findfont\n%d scalefont\nsetfont\n", fontsize);
  printf("%f %f moveto\n", lmargin, page_height-tmargin);
  printf("(%s) show\n", argv[-3]);
  for(channel = 0; channel != num_channels; channel++)
    draw_axes(channel, EL_lab(ch_hdr[channel]), ymin[channel], ymax[channel]);
  for(file = 0; file != argc; file++)
    {
  /*if standard deviations have been appended, use a duplicate file pointer*/
    avg = fopen(argv[file],
#ifdef MSDOS
      "rb"
#else
      "r"
#endif
      );
    if(variance_flag)
      {
      avg_tail = fopen(argv[file],
#ifdef MSDOS
	"rb"
#else
	"r"
#endif
	);
      fseek(avg_tail, (long)(packed_sizeof_SETUP+num_channels*(packed_sizeof_ELECTLOC+5L+num_samples*sizeof(float))), SEEK_SET);
      }
  /*read calibrations for current file*/
    for(channel = 0; channel != num_channels; channel++)
      {
      fseek(avg, (long)(packed_sizeof_SETUP+channel*packed_sizeof_ELECTLOC+EL_calib_offset), SEEK_SET);
      calibration[channel] = read_Intel_float(avg);
      }
  /*set a new line type for this file*/
    if(file > 0)
      printf("[%d 1] 0 setdash\n", argc-file);
  /*plot the curves*/
    for(channel = 0; channel != num_channels; channel++)
      {
      fseek(avg, 5L, SEEK_CUR);
      printf("newpath\n%f %f moveto\n", lmargin, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0-ymax[channel]*yscale[channel]);
      plotpoints(avg, num_samples, yscale[channel]*calibration[channel]);
      if(variance_flag)
	{
	printf("newpath\n%f %f moveto\n", lmargin, page_height-tmargin-title_height-channel*plot_height-plot_leading/2.0-ymax[channel]*yscale[channel]);
	plotvariance(avg, avg_tail, num_samples, yscale[channel]*calibration[channel]);
	}
      }
    if(variance_flag)
      fclose(avg_tail);
    fclose(avg);
    }
  printf("showpage\n");
  exit(0);
  }
