view saitools/saigraph.c @ 0:5242703e91d3 tip

Initial checkin for AIM92 aimR8.2 (last updated May 1997).

/*
    Copyright (c) Applied Psychology Unit, Medical Research Council. 1993
    ===========================================================================

    Permission to use, copy, modify, and distribute this software without fee 
    is hereby granted for research purposes, provided that this copyright 
    notice appears in all copies and in all supporting documentation, and that 
    the software is not redistributed for any fee (except for a nominal 
    shipping charge). Anyone wanting to incorporate all or part of this 
    software in a commercial product must obtain a license from the Medical 
    Research Council.

    The MRC makes no representations about the suitability of this 
    software for any purpose.  It is provided "as is" without express or 
    implied warranty.
 
    THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING 
    ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL 
    THE A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES 
    OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, 
    WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 
    ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS 
    SOFTWARE.
*/

/*------------------------------------------------------------------------*/

/*  saigraph.c
*  --------------
*  
* Works out the (first-order) temporal intervals between peaks in a .sai file,
* and saves them in a SINGLE-FRAME .sai
* In essence, a histogram.
*
* Uses ALL peaks: if you just want the locallly maxima peaks, then use 
* trbigpeak first.
*
* Output is another .sai file: but only a single frame. The header is
* modified accordingly.
* 
*
*
* M. Akeroyd.   May 1993.   Version 2.0
*
*               Autumn 1993: added grayscaleflag
*               Winter 1994: allowed weighting, and other general revisions.
*/



#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include "tip.h"


/* Function declarations */
/*------------------------------------------------------------*/

void parsecommandline (int argc, char *argv[], char inputfn[], char outputfn[]);

int readheader(FILE *inputfp);
void checkheader();
void writeheader(FILE *outputfp);

void changeheader(int no_frames, int new_pwidth, int new_nwidth, int framewidth_samples_output);

int findpeaksreverse (); 
void find_intervals(int total_peaks, int channel, int max_interval);

void writedata_output (FILE *outputfp, int samples_to_write);
void writedata_output_fg (FILE *figurefp, FILE *groundfp, int samples_to_write);
FILE *open_file (char filefn[], FILE *dir_default, int streamtype);

void copytooutput_fg_single (int total_peaks, int frame, int channel, int trigger_peak);

void write_matrix_hori_time(int framewidth, FILE *outputfigurefp);
void write_matrix_hori_freq(int framewidth, FILE *outputfigurefp);




/* Data arrays: global */
/*-----------------------------------------------------------------*/


short inputdata[MAX_DATA];        /* input .sai: a single channel of a frame */
short outputdata[MAX_DATA];       /* output .sai: a single channel */
short outputfiguredata[MAX_DATA]; /* output figure.sai: a single channel */
short outputgrounddata[MAX_DATA]; /* output ground.sai: a single channel */
short clip[MAX_CHANNELS];

short interval[MAX_CHANNELS][MAX_DATA];      /* histogram */ 
short summation[MAX_CHANNELS];



/* other variables */
/*----------------------------------------------------------------*/

struct Peak peak[MAX_PEAKS];                    /* used in findpeaksreverse():
						 * single channel */
struct NextEvent nextevent[MAX_NEXTEVENTS];     /* used in findpeaksreverse():
						 * single channel */



/* variables read from header */
/*---------------------------------------------------------------*/

char header[MAX_LINES_HEADER][MAX_LINE_LENGTH];
int header_lines;

int no_frames;
int frameheight;             /* number of channels */
int framewidth_samples;      /* pwidth + nwidth * samplerate */
int frameshift_samples;      /* frstep_aid * samplerate */
int pwidth;                  /* in msecs */
int nwidth;                  /* in msecs: NEGATIVE */
int width_win;               /* pixels */
int height_win;              /* pixels */
long samplerate;             /* samples per sec */
int mincf;                   /* Hz */
int maxcf;                   /* Hz */




/* misc */
/*-----------------------------------------------------------*/

char progname[MAX_STRING_LENGTH];
char outputfigurefn[MAX_STRING_LENGTH];


int verboseflag = OFF;               /* -v */
int removeflag = OFF;                /* -r */
double removevalue = 0.1;            
int widthflag = OFF;                 /* -width : if OFF, copy pwdith/nwidth 
				      *      off Input .sai */
int total_width;                     /* -width : total width of new .sai */
int asciiflag = OFF;                 /* -a */
int asciireverseflag = OFF;          /* -ar */
int summationflag = OFF;             /* -s */
int greyscaleflag = OFF;             /* -g */
int weightflag = OFF;             /* -weight */
int oppositearchflag = OFF;          /* -oparch: if ON, then the header-reading
				      * code will load one (1) extra byte. This
				      * is so that Sun .sai/.nap can be read
				      * on a DEC, and equally the reverse.
				      * It is NOT required if reading Sun on
				      * Sun, or DEC on DEC. */



/*...............           Main           ................................*/
/*.........................................................................*/
/*.........................................................................*/



void main (int argc, char *argv[])
{
  int sample;
  int n;

  int frame, channel;
  int trigger_peak = 0;
  int no_peak = 0, total_peaks =0, total_peaks_2;
  int framewidth_samples_input,
      framewidth_samples_output;
  int new_pwidth, 
      new_nwidth;
  int header_bytes = 0;
 
  char inputfn[MAX_STRING_LENGTH], 
       outputbasefn[MAX_STRING_LENGTH];

  FILE *inputfp, *outputfigurefp;

  int clipflag = OFF;

/*-----------------------------------------*/

  strcpy(progname, argv[0]);
  strcpy(inputfn, "");
  strcpy(outputbasefn, "");
  strcpy(outputfigurefn, "");

  /* parse command line */
  parsecommandline(argc, argv, inputfn, outputbasefn);

  if (strcmp(outputbasefn, "") == 0) 
    strcpy(outputfigurefn, "");
  else 
    strcpy(outputfigurefn, outputbasefn);


  /* open files */
  /* default directions are:
   * input = stdin 
   * figure = stdout
   */

  inputfp = open_file(inputfn, stdin, READ);
  outputfigurefp = open_file(outputfigurefn, stdout, WRITE);

  /* read Header */
  header_bytes = readheader(inputfp);
  checkheader();


  /* reset the variables */
  framewidth_samples_input = framewidth_samples;
  if (widthflag == OFF) {
    new_pwidth = + pwidth;
    new_nwidth = -nwidth;
    framewidth_samples_output = framewidth_samples_input;}
  else {
    new_pwidth = + NEW_PWIDTH;
    new_nwidth = (new_pwidth + total_width);
    framewidth_samples_output = (int) ((int) (abs(new_pwidth) + abs (new_nwidth)) * samplerate ) / 1000;}

  if (framewidth_samples_output >= MAX_DATA) {
    fprintf(stderr, "%s: new frame is  too wide at %i: only allowed %i samples\n", progname, framewidth_samples_output, MAX_DATA);
    exit(-1); }


  /* change header */
  changeheader(1, new_pwidth, new_nwidth, framewidth_samples_output);

  /* write header */
  if (asciiflag == OFF) 
    writeheader(outputfigurefp);
  else {
    /* begining of ascii header. The rest depends on the matrix format ... */
    fprintf(outputfigurefp, "# saigraph output.  %s\n", inputfn);
    fprintf(outputfigurefp, "# \n");
    fprintf(outputfigurefp, "# samplerate=%i  pwidth=%d  nwidth=%d  frames=%i ", samplerate, pwidth, nwidth, no_frames);
    fprintf(outputfigurefp, " mincf=%i  maxcf=%i  channels=%i \n", mincf, maxcf,  frameheight);
    fprintf(outputfigurefp, "# \n");
  }

/*------------------------------------------*/

  /* Main loop */

  if  (verboseflag == ON) { 
    fprintf(stderr, " frames "); 
    fflush(stderr);}

  
  /* clear arrays */
  
  for (channel=0; channel<MAX_CHANNELS; channel++) {
    for(sample=0; sample<MAX_DATA; sample++)
      interval[channel][sample] = 0;
    clip[channel]=OFF;
  }


  /* reset this, so everything still works */
  framewidth_samples = framewidth_samples_input;


  for (frame=1; frame<=no_frames; frame++) {

    if  (verboseflag == ON) { 
      fprintf(stderr, " %i ", frame); 
      fflush(stderr);}

    /*--------------------------------------------*/

    /* load the whole frame's worth of data first */

    for (channel=1; channel <=frameheight; channel ++) {

      /* clear arrays */
      for(sample=0; sample<framewidth_samples; sample++) 
	inputdata[sample] = outputdata[sample] = 0;
      for(n=0; n < MAX_PEAKS; n++){
	peak[n].tent = UNSET;
	peak[n].sample=0;}

      /* Nextevents are unnecessary */

      /* load a single channel's worth of data */
      fread (inputdata, 2, (size_t) framewidth_samples, inputfp);

      /* This next is a simple input check: if any numbers < 0, then say so */
      for(sample=0; sample<framewidth_samples; sample++)
	if (inputdata[sample] < 0 )  {
	  fprintf(stderr, "%s: something's gone wrong: the data is negative.\n", progname);
	  exit(-1); }

      /* find peaks. WARNING: in reverse */
      total_peaks = findpeaksreverse();

      /* bin any small peaks ... */
      if ((removeflag == ON) && (total_peaks !=0 ) ) {
        total_peaks_2 = removepeaks(total_peaks);
        for(n=0; n < MAX_PEAKS; n++) {
          peak[n].tent = UNSET;
          peak[n].sample = 0;}
        total_peaks = findpeaksreverse();}

      /* get intervals */
      find_intervals(total_peaks, channel, framewidth_samples_output);

    } /* channel */

  } /* frame */


/*------------------------------------------*/
/*------------------------------------------*/


  fflush(stdout);


  /* If required, sum the values */
  if (summationflag==ON) {

    for (channel = 1; channel <= frameheight; channel ++)
	summation[channel] = 0;

    for (channel = 1; channel <= frameheight; channel ++) {
      clipflag = OFF;
      for (sample = 0; sample < framewidth_samples_output; sample++)
	summation[channel] += interval[channel][sample];

      /* clip check */
      if ((summation[channel] > 32765) || 
	  (summation[channel] < 0)){
	    if (clipflag == OFF) {
	      fprintf(stderr, " \nclipping: channel %i, summations.\n ", channel);
	      clipflag = ON;}
	    summation[channel] = 32765;
	  }
    }
  }

/*-------------------------------------------*/

  /* write output */

  if (asciiflag == OFF) {
    for (channel = 1; channel <= frameheight; channel ++){
      for (sample = 0; sample < framewidth_samples_output; sample++)
	outputfiguredata[sample] = interval[channel][sample];
      writedata_output(outputfigurefp, framewidth_samples_output);
    }
  }  
  else {
    if (asciireverseflag == OFF) 
      write_matrix_hori_freq(framewidth_samples_output, outputfigurefp);
    else
      write_matrix_hori_time(framewidth_samples_output, outputfigurefp);
  }    
  
  /*-------------------------------------------*/

  /* Tidy up and exit */
    
  if (verboseflag == ON) 
    fprintf(stderr, "\n");
  fclose(inputfp); 
  fclose(outputfigurefp); 
 
  if ( ferror(inputfp) != 0) {
    fprintf(stderr, " %s : error closing input file.\n", progname);
    exit(-1);}

  if ( ferror(outputfigurefp) != 0) {
    fprintf(stderr, " %s : error closing figure file.\n", progname);
    exit(-1);}

  exit(0);

} /* Main */






/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/





void parsecommandline(int argc, char *argv[], char inputfn[], char outputbasefn[])
{
  int x=1, helpflag = OFF;

  while (x < argc){
    if      (!strcmp(argv[x], "-o")) {strcpy(outputbasefn, argv[x+1]); x+=2;}
    else if (!strcmp(argv[x], "-output")) {strcpy(outputbasefn, argv[x+1]); x+=2;}
    else if (!strcmp(argv[x], "-i")) {strcpy(inputfn, argv[x+1]); x+=2;}
    else if (!strcmp(argv[x], "-input")) {strcpy(inputfn, argv[x+1]); x+=2;}
    else if (!strcmp(argv[x], "-help")) {helpflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-h")) {helpflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-verbose")) {verboseflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-v")) {verboseflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-weight")) {weightflag = NORMAL_WEIGHTING; x+=1;}
    else if (!strcmp(argv[x], "-weightlog")) {weightflag = LOG_WEIGHTING; x+=1;}
    else if (!strcmp(argv[x], "-width")) {widthflag = ON; total_width=atoi(argv[x+1]); x+=2;}
    else if (!strcmp(argv[x], "-r")) { removeflag = ON; removevalue = atof(argv[x+1]); x+=2;}
    else if (!strcmp(argv[x], "-a")) { asciiflag = ON;x+=1;}
    else if (!strcmp(argv[x], "-ar")) { asciiflag = ON; asciireverseflag=ON; x+=1;}
    else if (!strcmp(argv[x], "-s")) { summationflag=ON; x+=1;}
    else if (!strcmp(argv[x], "-g")) { greyscaleflag=ON; x+=1;}
    else if (!strcmp(argv[x], "-grey")) { greyscaleflag=ON; x+=1;}
    else if (!strcmp(argv[x], "-gray")) { greyscaleflag=ON; x+=1;}
    else if (!strcmp(argv[x], "-oparch")) { oppositearchflag=ON; x+=1;}
    else {fprintf(stderr, "%s: unknown option %s\n", progname, argv[x]);
	  exit(-1);}
  }
  
  if (helpflag == ON)
    {
      fprintf(stderr, "\n");
      fprintf(stderr, "                         %s\n", progname);
      fprintf(stderr, " -----------------------------------------------------------------------------\n");
      fprintf(stderr, " Works out a histogram of peak intervals, from a .sai\n");
      fprintf(stderr, " Default output is a .sai file: if -a or -ar specified,\n");

      fprintf(stderr, " then output is a ASCII matrix.\n");
      fprintf(stderr, "\n");
      fprintf(stderr, " The -oparch option is REQUIRED if reading a Sun .sai files \n");
      fprintf(stderr, " on a DEC cpu or vice-versa.\n");
      fprintf(stderr, " Add .sai to both input and output filenames.\n");

      fprintf(stderr, " -----------------------------------------------------------------------------\n");
      fprintf(stderr, "\n");
      fprintf(stderr, " -i <.sai file>      input file             default = stdin\n");
      fprintf(stderr, " -o <.sai file>      output file            default = stdout\n");
      fprintf(stderr, " -oparch             assume input was generated on an opposing architecture\n");
      fprintf(stderr, "\n");
      fprintf(stderr, " -a                  ASCII output: rows = channels\n");
      fprintf(stderr, " -ar                 ASCII output: rows = time intervals\n");
      fprintf(stderr, " -s                  include summations over time intervals (ASCII options only)\n");
      fprintf(stderr, "\n");
      fprintf(stderr, " -width <int>        width of graph.sai (msecs).\n");
      fprintf(stderr, " -r <fraction>       remove any peaks F(n) < F(n+1) & F(n-1) * fract \n");
      fprintf(stderr, "                                            default = %.3f\n", removevalue);
      fprintf(stderr, " -weight             weight histogram entries by mean height of peaks\n");
      fprintf(stderr, " -weightlog          as per -weight, but use 10log(mean height)\n");
      fprintf(stderr, "\n");
      fprintf(stderr, " -v                  verbose output (to stderr)\n");
      fprintf(stderr, " -g                  change .sai 'view' format to 'grayscale'.\n");
      fprintf(stderr, "\n\n");
      exit(-1);}

}  




/* The End */
/*---------------------------------------------------------------------------*/