annotate saitools/saisummary.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.
*/
/*---------------------------------------------------------------------------*/
/* 
* saisummary
*
* builds a summary .sai 
* Output is one channel.
*
* M. Akeroyd.   May 1993.  Version 2.0
*               Revised Winter 1994.
* revised Autumn 1993.
* variuos fiddles added Summer 1994.
* (and some new options)
* and some more Spring 1995.
* and more Summer 1995.
*/


#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 writedata_output (FILE *outputfp, int samples_to_write);
FILE *open_file (char filefn[], FILE *dir_default, int streamtype);
void copytooutput_single(int no_peaks, int frame, int channel, int trigger_peak);

int findpeaksreverse();
int findbigpeaks(int total_peaks);
void integratebigpeaks(int total_bigpeaks);

void changeheader(int no_frames, int new_pwidth, int new_nwidth, int framewidth_samples_output);
void changeheader_B(int new_mincf, int new_maxcf, int new_channels);

void write_matrix_hori_time(int framewidth, FILE *outputfigurefp);



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


short inputdata[MAX_DATA];
short outputdata[MAX_DATA];
short outputfiguredata[MAX_DATA];
short outputgrounddata[MAX_DATA];
float summary[MAX_DATA];
float framesummary[MAX_DATA];

short interval[MAX_CHANNELS][MAX_DATA];   /* required for matrix.c */
short summation[MAX_CHANNELS];            /* required for matrix.c */ 

int clip[MAX_DATA];

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


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

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



/* .sai parameters */
/*-------------------------------------------------------------------------*/

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 summationflag = OFF;
int frameflag = OFF;                 /* -f */
int asciiflag = OFF;
int greyscaleflag = OFF;
int divideflag = OFF;

int oppositearchflag = OFF;          /* -oparch : see saigraph.c for info */

int nstopflag = ON;
int nwarnflag = OFF;
int n32760flag = OFF;
int n0flag = OFF;
int headeroutputflag = ON;
int singleframeflag = OFF;
int frameaverageflag = ON;
int specificchannelflag = OFF;
int cutchannels[MAX_CHANNELS];
int nchannels=0;

int framestart;   
int frameend;
int cut_framestart;
int cut_frameend;


/* ..............           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;
  int total_bigpeaks = 0;
  int trueframeheight = 0;
  int header_bytes = 0;

  double temp;
 
  char inputfn[MAX_STRING_LENGTH], 
       outputbasefn[MAX_STRING_LENGTH];

  FILE *inputfp, *outputfigurefp;

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

  for (n=0; n<MAX_CHANNELS; n++)
    cutchannels[n] = OFF;

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

/*    for (n=0; n<MAX_CHANNELS; n++)
      if (cutchannels[n] == ON)
	fprintf(stderr, "%d\n", n);
*/
  if (specificchannelflag == OFF)
    for (n=0; n<MAX_CHANNELS; n++)
      cutchannels[n] = ON;
  
  if (strcmp(outputbasefn, "") == 0) 
    strcpy(outputfigurefn, "");
  else {
    strcpy(outputfigurefn, outputbasefn);
  /*  strcat(outputfigurefn, OUTPUT_EXT);*/
  }

  /* 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);

  /* do some error-checking on the header, and set the fread 
   * pointer to the right place */
  checkheader();

  if (frameflag == ON) {
    cut_frameend = frameend;
    cut_framestart = framestart; }
  else {
    cut_frameend = no_frames;
    cut_framestart = 1; }

  if (frameflag == ON)
    changeheader((cut_frameend - cut_framestart + 1), pwidth, nwidth, framewidth_samples);

  if (singleframeflag == ON)
    changeheader(1, pwidth, nwidth, framewidth_samples);

  /* reset the variables */
  /* arguments: mincf  maxcf  channels */
  changeheader_B(0, 0, 1);  
  
  /* write header */
  if (singleframeflag == OFF) {
    if ((asciiflag == OFF) && (headeroutputflag == ON))
      writeheader(outputfigurefp);
    else ;}
  else {
    if ((asciiflag == OFF) && (headeroutputflag == ON))
      writeheader(outputfigurefp);
    else if (asciiflag == ON){
      /* begining of ascii header. The rest depends upon the output format. */
      fprintf(outputfigurefp, "# saisummary output. %s\n", inputfn);
      fprintf(outputfigurefp, "# \n");
      fprintf(outputfigurefp, "# samplerate=%i  pwidth=%d  nwidth=%d  frameheight=%i\n", samplerate, pwidth, nwidth, frameheight);
      fprintf(outputfigurefp, "# \n");
    }}


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

  /* Main loop */


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

  for(sample=0; sample<framewidth_samples; sample++) 
    framesummary[sample] = 0.0;

  for (frame=1; frame<=no_frames; frame++) {
    if  (verboseflag==ON) {
      fprintf(stderr, " %i ", frame);
      fflush(stderr);}

    for(sample=0; sample<framewidth_samples; sample++) {
      interval[1][sample] = outputfiguredata[sample] = 0;
      summary[sample] = 0.0;
      clip[sample] = OFF;}

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

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

      /* clear arrays */
      for(sample=0; sample<framewidth_samples; sample++)
        inputdata[sample] = 0;

      /* load data */
      fread (inputdata, 2, (size_t) framewidth_samples, inputfp);

      /* This next is a simple input check: */
      /* revised for various warnings */
      for(sample=0; sample<framewidth_samples; sample++) {
/*	  fprintf(stderr, "data: %i  sample: %i channel: %i \n", inputdata[sample], sample, channel);*/
        if (inputdata[sample] < 0 )  {
	  if (nstopflag == ON) {
	    if (verboseflag == ON) fprintf(stderr, "\n");
	    fprintf(stderr, "%s: negative data: ", progname);
	    fprintf(stderr, "%i  sample: %i channel: %i \n", inputdata[sample], sample, channel);
	    fprintf(stderr, "bye.\n");
	    exit(-1);}
	  if (nwarnflag == ON){
	    if (verboseflag == ON) fprintf(stderr, "\n");
	    fprintf(stderr, "warning: negative data: %i  (sample: %i channel: %i)\n", inputdata[sample], sample, channel);}
	  if (n0flag == ON)
	    inputdata[sample] = 0;
	  else if (n32760flag == ON)
	    inputdata[sample] = 32760;}}

      /* add data to the summary autoco: outputfigurefp */
      if ((frame >= cut_framestart) && (frame <= cut_frameend)){
	if (cutchannels[channel] == ON){
	  nchannels ++;
	  /*fprintf(stderr, "c%d ", channel);*/
	  for(sample=0; sample<framewidth_samples; sample++)
	    summary[sample] += (double) inputdata[sample];}
	/*fprintf(stderr, "%d hello\n", channel);*/
      }
    } /* channel */

/*    fprintf(stderr, "c%d ", nchannels);*/
    /* divide everything by no. of channels */
    if (divideflag == ON) {
      for(sample=0; sample<framewidth_samples; sample++){
	temp = (double) summary[sample] / nchannels;
	if ( ( temp > 32760.0)||( temp < 0)) {
	  if (clip[sample] == OFF) {
	    fprintf(stderr, "\nclipping: time interval %i\n", sample);
	    clip[sample] = ON; }
	  temp = 32761.0; }
	interval[1][sample] = outputfiguredata[sample] = (short) temp; }}
    else {
      for(sample=0; sample<framewidth_samples; sample++){
	temp = (double) summary[sample];
	if (( temp > 32760.0) || (temp < 0.0)){
	  if (clip[sample] == OFF) {
	    fprintf(stderr, "\nclipping: time interval %i\n", sample);
	    clip[sample] = ON; }
	  temp = 32761.0; }
	interval[1][sample] = outputfiguredata[sample] = (short) temp; }}

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

    /* if required, sum the values */
    if (summationflag == ON) {
       /* clear ... */
       summation[1] = 0;
       clip[0] = OFF;
       /* set ... */
       for (sample=0; sample < framewidth_samples; sample++) {
         summation[1] += interval[1][sample];
	 /* clip check */
	 if (( summation[1] > 32760 )||(summation[1] < 0 )) {
	   if (clip[0] == OFF) {
	     fprintf(stderr, "\nclipping: summation\n");
	     clip[0] = ON; }
	   summation[1] = 32761;}}

    }

    /* add to frame average */
    if ((frame >= cut_framestart) && (frame <= cut_frameend)){
      for (sample=0; sample < framewidth_samples; sample++) 
	framesummary[sample] += outputfiguredata[sample];
    }

    /* write output: outputfiguredata[] for the .sai
     *               interval[] & summation[] for the ascii
     * Note hack for frameheight
     * BUT: not if only 1 frame is to be output 
     */

    if ((frame >= cut_framestart) && (frame <= cut_frameend)){    
      if (singleframeflag == OFF) {
	if (asciiflag == OFF)
	  writedata_output(outputfigurefp, framewidth_samples);
	else {
	  trueframeheight=frameheight;
	  frameheight=1;
	  write_matrix_hori_time(framewidth_samples, outputfigurefp);
	  frameheight=trueframeheight; }}
      else ;
    }

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

  }  /* frame */

    if (singleframeflag == ON) {
      if (verboseflag == ON) {
	fprintf(stderr, " averaging ");
	fflush(stderr);}
      clip[0] = OFF;
      for (sample=0; sample < framewidth_samples; sample++) {
	if (frameaverageflag == ON)
	  temp = framesummary[sample]/(cut_frameend - cut_framestart + 1);
	else
	  temp = framesummary[sample];
	/* clip check */
	if (temp > 32760) {
	    if (clip[0] == OFF) {
	      fprintf(stderr, "\nclipping: frameaverage at sample %i\n", sample);
	      clip[0] = ON; }
	    temp = 32760;}
	outputfiguredata[sample]= (short)temp;
	interval[1][sample] = outputfiguredata[sample];}
      if (asciiflag == OFF)
	writedata_output(outputfigurefp, framewidth_samples);
      else {
	trueframeheight=frameheight;
	frameheight=1;
	write_matrix_hori_time(framewidth_samples, outputfigurefp);
	frameheight=trueframeheight; }}
    else ;


  /* Tidy up and exit */
  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 outputfn[])
{
  int x=1, helpflag = OFF;
  int group;

  while (x < argc){
    if      (!strcmp(argv[x], "-o")) {strcpy(outputfn, argv[x+1]); x+=2;}
    else if (!strcmp(argv[x], "-output")) {strcpy(outputfn, 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], "-s")) {summationflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-g")) {greyscaleflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-gray")) {greyscaleflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-grey")) {greyscaleflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-ar")) {asciiflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-d")) {divideflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-channels")) {specificchannelflag = ON; 
					     x+=1;
					     while (x<argc){
					       if ((int) strspn(argv[x], "-") == 0)
						 cutchannels[atoi(argv[x])]=ON;
					       else break;
					       x++;}}
    else if (!strcmp(argv[x], "-oneframe")) {singleframeflag = ON; frameaverageflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-sumframe")) {singleframeflag = ON; frameaverageflag = OFF; x+=1;}
    else if (!strcmp(argv[x], "-oparch")) {oppositearchflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-nstop")) {nstopflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-header")) {headeroutputflag = ON; x+=1;}
    else if (!strcmp(argv[x], "-noheader")) {headeroutputflag = OFF; x+=1;}
    else if (!strcmp(argv[x], "-nwarn")) {nwarnflag = ON; nstopflag = OFF; x+=1;}
    else if (!strcmp(argv[x], "-nignore")) {nwarnflag = OFF; nstopflag = OFF; x+=1;}
    else if (!strcmp(argv[x], "-n0")) {n0flag = ON; n32760flag = OFF; nstopflag = OFF; x+=1;}
    else if (!strcmp(argv[x], "-n32760")) {n0flag = OFF; n32760flag = ON; nstopflag = OFF; x+=1;}
    else if (!strcmp(argv[x], "-frames")) {frameflag=ON; framestart=atoi(argv[x+1]); frameend=atoi(argv[x+2]); x+=3;}

    else {fprintf(stderr, "%s: unknown option %s\n", progname, argv[x]);
	  exit(-1);}
  }
  
  if (helpflag == ON)
    {
      fprintf(stderr, "\n");
      fprintf(stderr, "                         saisummary\n");
      fprintf(stderr, " -----------------------------------------------------------\n");
      fprintf(stderr, " Works out a summary image (cf a summary autocorrelogram).\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, " -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 cpu\n");
      fprintf(stderr, "\n");
      fprintf(stderr, " -frames <start>  <end>   (inclusive) (default=all)\n");
      fprintf(stderr, " -channels <> <> .....   specify channel numbers to be included (1++)\n");
      fprintf(stderr, " -d               divide by no. of channels (Warning: precision may be lost).\n");
      fprintf(stderr, " -ar              ASCII output: rows=time intervals (requires -oneframe)\n");
      fprintf(stderr, " -s               include summations (ascii only)\n");
      fprintf(stderr, "\n");
      fprintf(stderr, " -header          include AIM header in output file (default)\n");
      fprintf(stderr, " -noheader        don't  ...  ...  ... \n");
      fprintf(stderr, " -v               verbose output (to stderr)\n");
      fprintf(stderr, " -g               change .sai format to 'grayscale'\n");
      fprintf(stderr, " -oneframe        average frames together; one (1) output frame\n");
      fprintf(stderr, " -sumframe        sum frames together; one (1) output frame\n");
      fprintf(stderr, "\n");
      fprintf(stderr, " -nstop           if negative inputs found, stop processing (default)\n");
      fprintf(stderr, " -nwarn           if negative inputs found, continue (but warn user) \n");
      fprintf(stderr, " -nignore         if negative inputs found, continue \n");
      fprintf(stderr, " -n0              if negative inputs found, replace with 0 \n");
      fprintf(stderr, " -n32760          if negative inputs found, replace with 32760 \n");
      fprintf(stderr, "\n\n");
      exit(-1);}
}  

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