dftavg.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculate a volume weighted average TAC of specified regions in a TAC file.",
  " ",
  "Usage: @P [options] tacfile [tacid(s)]",
  " ",
  "Options:",
  " -rm",
  "     Remove the original TACs after averaging.",
  " -hem[isphere]",
  "     Average of hemispheres (dx and sin as the 2nd name field).",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1: Average TAC of all TACs named as 'cer':",
  "     @P b123.dat cer",
  " ",
  "Example 2: Average TACs over planes of all regions separately;",
  "           remove the original TACs with option -rm.",
  "     @P -rm b123.dat",
  " ",
  "Example 3: Average of TAC numbers 4, 6 and 7:",
  "     @P b123.dat 4 6 7",
  " ",
  "Example 4: Average of hemispheres (dx and sin) of all regions separately;",
  "           remove the original TACs with option -rm.",
  "     @P -hemisphere -rm b123.dat",
  " ",
  "See also: taclist, tacdel, dftrmdpl, tacadd, avgttac, taccalc, tacmean",
  " ",
  "Keywords: TAC, modelling, tool, ROI, average",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int    ai, help=0, version=0, verbose=1;
  int    fi, ri, rj;
  int    n, voi, ret, rmOrig=0, hemisphere=0, tacnamenr=0, operation=0;
  DFT    dft;
  char  *cptr, dfile[FILENAME_MAX];
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dfile[0]=(char)0;
  dftInit(&dft);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "HEMISPHERE", 3)==0) {
      hemisphere=1; continue;
    } else if(strcasecmp(cptr, "RM")==0 || strcasecmp(cptr, "DEL")==0) {
      rmOrig=1; continue;
    }
    fprintf(stderr, "Error: invalid option %s\n", argv[ai]); 
    return(1);
  } else break;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  for(; ai<argc; ai++) {
    /* The first 'non-option' argument is the datafile */
    if(!dfile[0]) {strcpy(dfile, argv[ai]); continue;}
    /* The following ones are the TAC names */
    if(tacnamenr==0) tacnamenr=ai;
  }
 
  /* Is something missing? */
  if(!dfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* decide what to do */
  if(hemisphere>0) {
    operation=1; /* separate avg over hemispheres for each region */
  } else if(tacnamenr==0) {
    operation=0; /* separate avg over planes for each region */
  } else {
    operation=2; /* avg of specified TACs */
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", dfile);
    printf("rmOrig := %d\n", rmOrig);
    printf("hemisphere := %d\n", hemisphere);
    printf("operation := %d\n", operation);
    if(tacnamenr>0) {
      printf("TACs :=");
      for(ai=tacnamenr; ai<argc; ai++)
        if(*argv[ai]!='-') printf(" '%s'", argv[ai]);
      printf("\n");
    }
  }
 
 
  /*
   *  Read TAC file
   */
  if(verbose>1) printf("reading %s\n", dfile);
  if(dftRead(dfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dfile, dfterrmsg);
    return(2);
  }
  if(verbose>1) {
    printf("  tacNr := %d\n", dft.voiNr);
    printf("  sampleNr := %d\n", dft.frameNr);
  }
 
  /*
   *  Find the TACs which are to be averaged with each other
   */
  /* first 'deselect' all */
  for(ri=0; ri<dft.voiNr; ri++)
    dft.voi[ri].sw=dft.voi[ri].sw2=dft.voi[ri].sw3=-1;
  n=0;
  if(operation==0) {
    /* For each region, find the corresponding TACs with the same voiname */
    /* and hemisphere */
    for(ri=0; ri<dft.voiNr; ri++) { 
      /* not if already selected */
      if(dft.voi[ri].sw2>=0) continue;
      /* not if this is not a single plane */
      if(strncasecmp(dft.voi[ri].place, "Pl", 2)) continue;
      /* this is the first and it gets its own tac nr */
      dft.voi[ri].sw2=ri; dft.voi[ri].sw3=1; n++;
      /* then search the following tacs */
      for(rj=ri+1; rj<dft.voiNr; rj++) { 
        /* not if already selected */
        if(dft.voi[rj].sw2>=0) continue;
        /* not if this is not a single plane */
        if(strncasecmp(dft.voi[rj].place, "Pl", 2)) continue;
        /* if voiname and hemisphere do match, then set the tac nr */
        if(!strcasecmp(dft.voi[ri].voiname, dft.voi[rj].voiname) &&
           !strcasecmp(dft.voi[ri].hemisphere, dft.voi[rj].hemisphere)) {
          dft.voi[rj].sw2=ri;
          dft.voi[ri].sw3++;
        }
      }
    }
    if(n<=0) {
      fprintf(stderr, "Warning: no planes were found.\n");
      dftEmpty(&dft); return(0);
    }
  } else if(operation==1) {
    /* Find the dx&sin of separate regions */
    /* Select the specified regions, or all if none was specified */
    if(tacnamenr>0) {
      for(ai=tacnamenr; ai<argc; ai++) if(*argv[ai]!='-') {
        /* check if number was given; if not then find the names */
        voi=atoi(argv[ai])-1;
        if(voi>=0 && voi<dft.voiNr) dft.voi[voi].sw=1;
        else dftSelect(&dft, argv[ai]);
      }
    } else {
      for(ri=0; ri<dft.voiNr; ri++) dft.voi[ri].sw=1;
    }
    /* Set sw=1 if dx, sw=2 if sin, and deselect if something else */
    for(ri=0; ri<dft.voiNr; ri++) if(dft.voi[ri].sw) {
      /* check for dx and sin */
      if(!strcasecmp(dft.voi[ri].hemisphere, "dx")) dft.voi[ri].sw=1;
      else if(!strcasecmp(dft.voi[ri].hemisphere, "sin")) dft.voi[ri].sw=2;
      else if(!strcasecmp(dft.voi[ri].hemisphere, "right")) dft.voi[ri].sw=1;
      else if(!strcasecmp(dft.voi[ri].hemisphere, "left")) dft.voi[ri].sw=2;
      else if(!strcasecmp(dft.voi[ri].hemisphere, "r")) dft.voi[ri].sw=1;
      else if(!strcasecmp(dft.voi[ri].hemisphere, "l")) dft.voi[ri].sw=2;
      else if(!strcasecmp(dft.voi[ri].place, "dx")) dft.voi[ri].sw=1;
      else if(!strcasecmp(dft.voi[ri].place, "sin")) dft.voi[ri].sw=2;
      else if(!strcasecmp(dft.voi[ri].place, "right")) dft.voi[ri].sw=1;
      else if(!strcasecmp(dft.voi[ri].place, "left")) dft.voi[ri].sw=2;
      else if(!strcasecmp(dft.voi[ri].place, "r")) dft.voi[ri].sw=1;
      else if(!strcasecmp(dft.voi[ri].place, "l")) dft.voi[ri].sw=2;
      else dft.voi[ri].sw=0;
    }
    /* Find the corresponding sin for each dx */
    for(ri=n=0; ri<dft.voiNr; ri++) if(dft.voi[ri].sw==1) { /* loop for dx's */
      if(dft.voi[ri].sw2>=0) continue;
      for(rj=0; rj<dft.voiNr; rj++) if(dft.voi[rj].sw==2) { /* loop for sin's */
        if(dft.voi[rj].sw2>=0) continue;
        /* check for region name and plane */
        if(strcasecmp(dft.voi[ri].voiname, dft.voi[rj].voiname)) continue;
        if(strcasecmp(dft.voi[ri].place, dft.voi[rj].place)
           && strcasecmp(dft.voi[ri].hemisphere, dft.voi[rj].hemisphere)) 
          continue;
        /* ok, this is a match */
        if(ri<rj) dft.voi[ri].sw2=dft.voi[rj].sw2=ri;
        else dft.voi[ri].sw2=dft.voi[rj].sw2=rj;
        dft.voi[ri].sw3=2; n++;
      }
    }
    if(n<1) {
      fprintf(stderr, "Error: no hemispheres to average were found.\n");
      dftEmpty(&dft); return(2);
    }
  } else {
    /* Find all TACs that match the specified names / numbers */
    for(ai=tacnamenr; ai<argc; ai++) if(*argv[ai]!='-') {
      /* check if number was given; if not then find the names */
      voi=atoi(argv[ai])-1;
      if(voi>=0 && voi<dft.voiNr) {
        dft.voi[voi].sw2=1; n=1;
      } else {
        n=dftSelect(&dft, argv[ai]);
        for(ri=0; ri<dft.voiNr; ri++) if(dft.voi[ri].sw) dft.voi[ri].sw2=1;
      }
    }
    for(ri=0, rj=-1; ri<dft.voiNr; ri++) if(dft.voi[ri].sw2==1) {
      if(rj<0) {rj=ri; dft.voi[ri].sw3=0;}
      dft.voi[ri].sw2=rj; dft.voi[rj].sw3++;
    }
    n=dft.voi[rj].sw3;
    if(n<=0) {
      fprintf(stderr, "Error: region '%s' not found.\n", argv[ai]);
      dftEmpty(&dft); return(2);
    }
  }
  if(verbose>1) {
    printf("%d average TACs will be produced.\n", n);
    for(ri=0; ri<dft.voiNr; ri++) if(dft.voi[ri].sw3>=0) {
      printf("Calculate avg from %d TACs:\n", dft.voi[ri].sw3);
      printf("  %s %s %s %f\n", dft.voi[ri].voiname, dft.voi[ri].hemisphere,
        dft.voi[ri].place, dft.voi[ri].size);
      for(rj=ri+1; rj<dft.voiNr; rj++) if(dft.voi[ri].sw2==dft.voi[rj].sw2)
        printf("  %s %s %s %f\n", dft.voi[rj].voiname, dft.voi[rj].hemisphere,
          dft.voi[rj].place, dft.voi[rj].size);
    }
  }
 
  /*
   *  Calculate average TACs
   */
  /* Allocate memory for new TACs */
  if(dftAddmem(&dft, n)) {
    fprintf(stderr, "Error in reallocating memory.\n");
    dftEmpty(&dft); return(3);
  }
  /* One average TAC at a time */
  /* find the first tac to be averaged */
  for(ri=0; ri<dft.voiNr; ri++) if(dft.voi[ri].sw2==ri) {
    /* tell the user what we are doing */
    if(verbose>0) {
      printf("  averaging TACs:\n");
      printf("    %-6.6s %-6.6s %-6.6s   %g\n", dft.voi[ri].voiname,
        dft.voi[ri].hemisphere, dft.voi[ri].place, dft.voi[ri].size);
    }
    /* start the summing with the first tac */
    if(dft.voi[ri].size>0) {
      dft.voi[dft.voiNr].size=dft.voi[ri].size;
      for(fi=0; fi<dft.frameNr; fi++)
        dft.voi[dft.voiNr].y[fi]=dft.voi[ri].size*dft.voi[ri].y[fi];
    } else {
      dft.voi[dft.voiNr].size=1.0;
      for(fi=0; fi<dft.frameNr; fi++)
        dft.voi[dft.voiNr].y[fi]=dft.voi[ri].y[fi];
    }
    /* copy the names to start with something */
    strcpy(dft.voi[dft.voiNr].voiname, dft.voi[ri].voiname);
    strcpy(dft.voi[dft.voiNr].hemisphere, dft.voi[ri].hemisphere);
    if(dft.voi[ri].sw3>1) strcpy(dft.voi[dft.voiNr].place, dft.voi[ri].place);
    else strcpy(dft.voi[dft.voiNr].place, "All");
    /* Sum all tacs with the same number */
    for(rj=ri+1; rj<dft.voiNr; rj++) if(dft.voi[rj].sw2==ri) {
      if(verbose>0)
        printf("    %-6.6s %-6.6s %-6.6s   %g\n", dft.voi[rj].voiname,
               dft.voi[rj].hemisphere, dft.voi[rj].place, dft.voi[rj].size);
      if(dft.voi[rj].size>0) {
        dft.voi[dft.voiNr].size+=dft.voi[rj].size;
        for(fi=0; fi<dft.frameNr; fi++)
          dft.voi[dft.voiNr].y[fi]+=dft.voi[rj].size*dft.voi[rj].y[fi];
      } else {
        dft.voi[dft.voiNr].size+=1.0;
        for(fi=0; fi<dft.frameNr; fi++)
          dft.voi[dft.voiNr].y[fi]+=dft.voi[rj].y[fi];
      }
      /* let the names be if they are the same */
      if(strcasecmp(dft.voi[dft.voiNr].voiname, dft.voi[rj].voiname))
        strcpy(dft.voi[dft.voiNr].voiname, "");
      if(strcasecmp(dft.voi[dft.voiNr].hemisphere, dft.voi[rj].hemisphere))
        strcpy(dft.voi[dft.voiNr].hemisphere, "");
      if(strcasecmp(dft.voi[dft.voiNr].place, dft.voi[rj].place))
        strcpy(dft.voi[dft.voiNr].place, "");
    }
    /* divide with sum size / number */
    for(fi=0; fi<dft.frameNr; fi++)
      dft.voi[dft.voiNr].y[fi]/=dft.voi[dft.voiNr].size;
    if(dft.voi[ri].size<=0) dft.voi[dft.voiNr].size=0.0;
    /* set names for the avg tac (if common names were not found) */
    if(!strcmp(dft.voi[dft.voiNr].voiname, ""))
      strcpy(dft.voi[dft.voiNr].voiname, "Mean");
    if(!strcmp(dft.voi[dft.voiNr].place, ""))
      strcpy(dft.voi[dft.voiNr].place, "All");
    /* make long name */
    rnameCatenate(dft.voi[dft.voiNr].name, MAX_REGIONNAME_LEN, 
      dft.voi[dft.voiNr].voiname, dft.voi[dft.voiNr].hemisphere, dft.voi[dft.voiNr].place, '-');
    /* set switches */
    dft.voi[dft.voiNr].sw=0; dft.voi[dft.voiNr].sw2=-1;
    /* tell user what we did */
    if(verbose>0) 
        printf("    ->    %-6.6s %-6.6s %-6.6s   %g\n",
          dft.voi[dft.voiNr].voiname, dft.voi[dft.voiNr].hemisphere, 
          dft.voi[dft.voiNr].place, dft.voi[dft.voiNr].size);
    /* next region */
    dft.voiNr++;
  }
 
  /*
   *  Remove those original tacs which were averaged
   */
  if(rmOrig) {
    ri=dft.voiNr-1; n=dft.voiNr;
    while(ri>=0) {
      if(dft.voi[ri].sw2>=0) {
        if((ret=dftDelete(&dft, ri))) {
          fprintf(stderr, "Error (%d) in deleting TAC.\n", ret);
          dftEmpty(&dft); return(5);
        }
      }
      ri--;
    }
    n-=dft.voiNr;
    if(n>0 && verbose>0)
      printf("%d TAC(s) are deleted; backup is in %s.bak\n", n, dfile);
  }
  if(verbose>9) dftPrint(&dft);
 
  /*
   *  Save the TACs
   */
  if(verbose>1) printf("writing %s\n", dfile);
  if(dftWrite(&dft, dfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", dfile, dfterrmsg);
    dftEmpty(&dft); return(11);
  }
  if(verbose>0) printf("%s written.\n", dfile);
 
 
  /*
   *  Free memory
   */
  dftEmpty(&dft);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/