parai.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcpar.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculate Asymmetry Index (AI) from PET results inside a parameter file.",
  " ",
  "                        (dx - sin)        ",
  "  Asymmetry index AI = ------------ x 100 ",
  "                       (dx + sin)/2       ",
  " ",
  "The right and left hemisphere is recognized from the following strings:",
  "  dx, dex, dexter, right, rgt, sn, sin, sinist, left, lft",
  " ",
  "Usage: @P [options] file [resultfile]",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: pardiff, paradd, parcoll, parformat, tacren, parget, parrenp",
  " ",
  "Keywords: parameter, tool, ratio",
  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;
  char parfile[FILENAME_MAX], outfile[FILENAME_MAX];
  /* List of identified dx/sin representation in region names */
  static char *dx_string[] = {"dx", "dex", "dexter", "right", "rgt", "d", 0};
  static char *sin_string[] = {"sn", "sin", "sinist", "left", "lft", "s", 0};
 
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  parfile[0]=outfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strcasecmp(cptr, "DRY")==0) {
      continue;
    }
    fprintf(stderr, "Error: invalid option '%s'\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-5;
  
  /* 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 */
  if(ai<argc) strlcpy(parfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(outfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]); return(1);}
 
  /* Is something missing? */
  if(!parfile[0]) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  if(strcasecmp(parfile, outfile)==0) {
    fprintf(stderr, "Error: input file would be overwritten.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("parfile := %s\n", parfile);
    if(outfile[0]) printf("outfile := %s\n", outfile);
    fflush(stdout);
  }
 
 
  /* 
   *  Read the file
   */
  if(verbose>1) printf("reading %s\n", parfile);
  PAR par; parInit(&par);
  if(parRead(&par, parfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    parFree(&par); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", parFormattxt(par.format));
    printf("parNr := %d\n", par.parNr);
    printf("tacNr := %d\n", par.tacNr);
  }
  if(par.tacNr<2) {
    fprintf(stderr, "Error: no data to calculate AI.\n");
    parFree(&par); return(2);
  }
  /* Sort regions by name */
  if(parSortByName(&par, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    parFree(&par); return(2);
  }
 
  /* List the region names */
  if(verbose>10) {
    fflush(stdout);
    for(int i=0; i<par.tacNr; i++) printf("%s\n", par.r[i].name);
    fflush(stdout);
  }
 
 
  /*
   *  Search and count the dx and sin regions
   *  Set sw=1 for dx and sw=2 for sin, and sw=0 for those that are neither
   */
  if(verbose>1) {printf("counting dx and sin regions\n"); fflush(stdout);}
  int ndx=0, nsin=0;
  char test_string[MAX_TACNAME_LEN+1];
 
  for(int ri=0; ri<par.tacNr; ri++) {
    if(verbose>3) printf("--------------\n%s\n", par.r[ri].name);
    par.r[ri].sw=0;
    if(strnlen(par.r[ri].name, 4)<3) continue;
    /* Try to find hemisphere strings */
    int j=0; 
    while(dx_string[j]!=0) {
      /* Try to find in the middle */
      sprintf(test_string, "-%s-", dx_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(strcasestr(par.r[ri].name, test_string)!=NULL) {
        par.r[ri].sw=1; break;
      }
      sprintf(test_string, "_%s_", dx_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(strcasestr(par.r[ri].name, test_string)!=NULL) {
        par.r[ri].sw=1; break;
      }
      sprintf(test_string, " %s ", dx_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(strcasestr(par.r[ri].name, test_string)!=NULL) {
        par.r[ri].sw=1; break;
      }
      /* Try to find in the end */
      unsigned int n=strlen(par.r[ri].name);
      sprintf(test_string, "-%s", dx_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(n>strlen(test_string) && strcasecmp(par.r[ri].name+n-strlen(test_string), test_string)==0) {
        par.r[ri].sw=1; break;
      }
      sprintf(test_string, "_%s", dx_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(n>strlen(test_string) && strcasecmp(par.r[ri].name+n-strlen(test_string), test_string)==0) {
        par.r[ri].sw=1; break;
      }
      sprintf(test_string, " %s", dx_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(n>strlen(test_string) && strcasecmp(par.r[ri].name+n-strlen(test_string), test_string)==0) {
        par.r[ri].sw=1; break;
      }
      j++;
    }
    if(par.r[ri].sw!=0) {
      ndx++;
      if(verbose>3) printf("  found string '%s'\n", test_string);
      (void)strdelstr(par.r[ri].name, test_string);
      continue;
    }
    j=0; 
    while(sin_string[j]!=0) {
      sprintf(test_string, "-%s-", sin_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(strcasestr(par.r[ri].name, test_string)!=NULL) {
        par.r[ri].sw=2; break;
      }
      sprintf(test_string, "_%s_", sin_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(strcasestr(par.r[ri].name, test_string)!=NULL) {
        par.r[ri].sw=2; break;
      }
      sprintf(test_string, " %s ", sin_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(strcasestr(par.r[ri].name, test_string)!=NULL) {
        par.r[ri].sw=2; break;
      }
      /* Try to find in the end */
      unsigned int n=strlen(par.r[ri].name);
      sprintf(test_string, "-%s", sin_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(n>strlen(test_string) && strcasecmp(par.r[ri].name+n-strlen(test_string), test_string)==0) {
        par.r[ri].sw=2; break;
      }
      sprintf(test_string, "_%s", sin_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(n>strlen(test_string) && strcasecmp(par.r[ri].name+n-strlen(test_string), test_string)==0) {
        par.r[ri].sw=2; break;
      }
      sprintf(test_string, " %s", sin_string[j]);
      if(verbose>6) printf("  testing '%s'\n", test_string);
      if(n>strlen(test_string) && strcasecmp(par.r[ri].name+n-strlen(test_string), test_string)==0) {
        par.r[ri].sw=2; break;
      }
      j++;
    }
    if(par.r[ri].sw!=0) {
      nsin++;
      if(verbose>3) printf("  found string '%s'\n", test_string);
      (void)strdelstr(par.r[ri].name, test_string);
      continue;
    }
  }
  if(ndx<1 || nsin<1) {
    fprintf(stderr, "Error: no dx/sin pairs were identified.\n");
    parFree(&par); return(3);
  }
  /* Delete regions where dx/sin was not identified */
  if(verbose>2) {printf("deleting unidentified regions from list\n"); fflush(stdout);}
  {
    int i=0;
    while(i<par.tacNr) {
      if(par.r[i].sw==0) parDeleteTAC(&par, i); else i++;
    }
  }
  /* List the region names */
  if(verbose>3) {
    fflush(stdout);
    for(int i=0; i<par.tacNr; i++) printf("%s\n", par.r[i].name);
    fflush(stdout);
  }
 
  /*
   *  Create a new parameter structure for AIs.
   */
  if(verbose>2) {printf("allocating place for AIs\n"); fflush(stdout);}
  PAR pai; parInit(&pai);
  int pn=ndx; if(nsin<pn) pn=nsin;
  if(parAllocate(&pai, par.parNr, pn)) {
    fprintf(stderr, "Error: cannot allocate memory for AIs.\n");
    parFree(&par); parFree(&pai); return(4);
  }
  /* Set results information */
  pai.tacNr=0; pai.parNr=par.parNr;
  pai.format=par.format;
  for(int i=0; i<pai.parNr; i++) {
    strcpy(pai.n[i].name, par.n[i].name);
    pai.n[i].unit=UNIT_PERCENTAGE;
  }
  {
    char buf[256];
    tpcProgramName(argv[0], 1, 1, buf, 256);
    iftPut(&pai.h, "program", buf, 0, NULL);
    time_t t=time(NULL);
    iftPut(&pai.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
  }
  iftPut(&pai.h, "datafile", parfile, 0, NULL);
 
 
  /*
   *  Compute AIs
   */
  if(verbose>1) {printf("computing AIs\n"); fflush(stdout);}
  for(int ri=0; ri<par.tacNr; ri++) if(par.r[ri].sw==1) {
    /* now search for the matching left hemisphere */
    for(int rj=0; rj<par.tacNr; rj++) if(par.r[rj].sw==2) {
      /* check that the remaining region name is matching */
      if(strcasecmp(par.r[ri].name, par.r[rj].name)!=0) continue;
      /* Set AI region name */
      strcpy(pai.r[pai.tacNr].name, par.r[ri].name);
      /* Calculate AI for all parameters */
      for(int pi=0; pi<par.parNr; pi++) {
        double f=par.r[ri].p[pi]+par.r[rj].p[pi];
        if(f==0.) pai.r[pai.tacNr].p[pi]=0.;
        else pai.r[pai.tacNr].p[pi]=200.*(par.r[ri].p[pi]-par.r[rj].p[pi])/f;
      } /* next parameter */
      pai.tacNr++;
      /* Do not process these regions again */
      par.r[rj].sw=0;
      break;
    }
    /* Do not process this region again */
    par.r[ri].sw=0;
  } /* next region */
  /* Original parameter data no more needed */
  parFree(&par);
  if(pai.tacNr<1) {
    fprintf(stderr, "Error: no data to calculate AI.\n");
    parFree(&par); return(6);
  }
 
 
  /* Print and save */
  if(verbose>0 || !outfile[0]) parWrite(&pai, stdout, PAR_FORMAT_TSV_UK, 0, NULL);
  if(outfile[0]) {
    /* Save file */
    if(verbose>1) printf("  saving %s\n", outfile);
    FILE *fp;
    fp=fopen(outfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      parFree(&pai); return(11);
    }
    int ret=parWrite(&pai, fp, PAR_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      parFree(&pai); return(12);
    }
    if(verbose>0) printf("parameter AIs saved in %s\n", outfile);
  }
 
  parFree(&pai);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/