tocr.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 "tpcift.h"
#include "tpcisotope.h"
#include "tpctac.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Convert the head-curve (count-rate) data file from ECAT HR+ (*.r) or",
  "HRRT (*.hc) to TAC file format that is suitable for correction of ",
  "time delay between blood/plasma and tissue TACs.",
  "Some usual errors in count-rate data is automatically corrected.",
  "Output data is corrected for physical decay.",
  " ",
  "Usage: @P [Options] headcurve isotope outputfile",
  " ",
  "Options:",
  " -min | -sec",
  "     Sample times are written in minutes or seconds;",
  "     by default in sec for tracers with short half-life, otherwise in min.",
  " -copy",
  "     If file is not in HR+ or HRRT head-curve format, it is still saved",
  "     with the new name.",
  " -f=<format-id>",
  " -format=<format-id>",
  "     Accepted format-id's for the output count-rate file:",
  "     CSV-INT   - CSV format with semicolons and decimal commas.",
  "     CSV-UK    - CSV format with commas and decimal points.",
  "     TSV-INT   - TSV format with tabs and decimal commas.",
  "     TSV-UK    - TSV format with tabs and decimal points.",
  "     PMOD      - PMOD tac and bld format.",
  "     DFT       - TPC TAC format (default).",
  "     SIMPLE    - txt file with tabs and decimal points.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P -format=simple is129.r C-11 is129.cr",
  " ",
  "See also: fitdelay, imghead, tacmean, dftscale, tac2svg, tacdecay",
  " ",
  "Keywords: ECAT HR+, HRRT, input, time delay, count-rate, head-curve",
  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     ret;
  int     time_unit=UNIT_UNKNOWN;
  int     outputformat=TAC_FORMAT_DFT;
  int     copy_cr=0;    // 0=no, 1=yes
  char   *cptr, petfile[FILENAME_MAX], crfile[FILENAME_MAX];
  int     isotope=ISOTOPE_UNKNOWN;
  TAC     tac;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacInit(&tac);
  petfile[0]=crfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "MIN", 1)==0) {
      time_unit=UNIT_MIN; continue;
    } else if(strncasecmp(cptr, "SEC", 1)==0) {
      time_unit=UNIT_SEC; continue;
    } else if(strncasecmp(cptr, "COPY", 1)==0) {
      copy_cr=1; continue;
    } else if(strncasecmp(cptr, "F=", 2)==0) {
      cptr+=2;
      outputformat=tacFormatIdentify(cptr);
      if(outputformat!=TAC_FORMAT_UNKNOWN) continue;
      /* for compatibility with previous versions */
      if(strcasecmp(cptr, "CR")==0 || strcasecmp(cptr, "HEAD")==0) {
        outputformat=TAC_FORMAT_SIMPLE; continue;}
    } else if(strncasecmp(cptr, "FORMAT=", 7)==0) {
      cptr+=7;
      outputformat=tacFormatIdentify(cptr);
      if(outputformat!=TAC_FORMAT_UNKNOWN) continue;
      /* for compatibility with previous versions */
      if(strcasecmp(cptr, "CR")==0 || strcasecmp(cptr, "HEAD")==0) {
        outputformat=TAC_FORMAT_SIMPLE; continue;}
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
  
  /* 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);}
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
 
  /* Process other arguments, starting from the first non-option */
  for(; ai<argc; ai++) {
    if(!petfile[0]) {
      strcpy(petfile, argv[ai]); continue;
    } else if(isotope==ISOTOPE_UNKNOWN) {
      isotope=isotopeIdentify(argv[ai]);
      if(isotope!=ISOTOPE_UNKNOWN) continue;
      fprintf(stderr, "Error: invalid isotope '%s'\n", argv[ai]); 
      return(1);
    } else if(!crfile[0]) {
      strcpy(crfile, argv[ai]); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!petfile[0] || isotope==ISOTOPE_UNKNOWN) {
    tpcPrintUsage(argv[0], info, stdout); return(1);}
  /* Make output filename if necessary */
  if(!crfile[0]) {
    strcpy(crfile, petfile);
    cptr=filenameGetExtension(crfile); if(cptr!=NULL) *cptr=(char)0;
    if(outputformat==TAC_FORMAT_SIMPLE) {
      strcat(crfile, ".cr"); // for compatibility with old scripts
    } else {
      strcat(crfile, tacDefaultExtension(outputformat));
    }
  }
  /* If time unit was not set with option, then set it based on half-life */
  if(time_unit==UNIT_UNKNOWN) {
    if(isotopeHalflife(isotope)<3.0) time_unit=UNIT_SEC; 
    else time_unit=UNIT_MIN;
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("crfile := %s\n", crfile);
    printf("isotope := %s\n", isotopeName(isotope));
    printf("time_unit := %s\n", unitName(time_unit));
    printf("copy_cr := %d\n", copy_cr);
    printf("outputformat := %s\n", tacFormattxt(outputformat));
  }
 
  /*
   *  Read the given datafile
   */
  if(verbose>1) printf("reading %s\n", petfile);
  ret=tacRead(&tac, petfile, &status);
  if(ret) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), petfile);
    tacFree(&tac);
    return(2);
  }
  if(verbose>1) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("columnNr := %d\n", tac.tacNr);
    printf("frameNr := %d\n", tac.sampleNr);
    printf("time_unit := %s\n", unitName(tac.tunit));
    printf("unit := %s\n", unitName(tac.cunit));
  }
 
  /*
   *  Process the data
   */
  if(tac.format==TAC_FORMAT_HRRT_HC) {
 
    if(verbose>0) printf("identified as HRRT format.\n");
    int i, n;
    /* Calculate 'frame' start and end times */
    for(i=0; i<tac.sampleNr; i++) {
      if(i==0) tac.x1[i]=0.0; else tac.x1[i]=tac.x2[i-1];
      tac.x2[i]=tac.x[i];
    }
    /* Reset 'frame' mid time */
    for(i=0; i<tac.sampleNr; i++) {
      tac.x[i]=0.5*(tac.x1[i]+tac.x2[i]);
    }
    /* kcps = prompts-randoms */
    for(i=0; i<tac.sampleNr; i++) {
      tac.c[0].y[i]=tac.c[2].y[i]-tac.c[1].y[i];
    }
    /* Remove negative values */
    i=0; n=0;
    while(i<tac.sampleNr) {
      if(tac.c[0].y[i]<0.0) {tacDeleteSample(&tac, i); n++;} else i++;
    }
    if(n>0 && verbose>1) printf("%d negative point(s) deleted.\n", n);
 
    /* Set tacNr etc */
    tac.tacNr=1;
    strcpy(tac.c[0].name, "CR");
 
  } else if(tac.format==TAC_FORMAT_HRPLUS_HC) {
 
    if(verbose>0) printf("identified as HR+ format.\n");
    int i, n;
    double a;
    /* Calculate 'frame' start and end times */
    for(i=0; i<tac.sampleNr; i++) {
      if(i==0) tac.x1[i]=0.0; else tac.x1[i]=tac.x2[i-1];
      tac.x2[i]=tac.x[i];
    }
    /* Reset 'frame' mid time */
    for(i=0; i<tac.sampleNr; i++) {
      tac.x[i]=0.5*(tac.x1[i]+tac.x2[i]);
    }
    /* Calculate count-rates as p_rate-d_rate */
    /* Consider valid data to end when at least 5 consecutive points are zeroes */
    for(i=n=0; i<tac.sampleNr && n<5; i++) {
      tac.c[0].y[i]=tac.c[2].y[i]-tac.c[3].y[i];
      if(tac.c[0].y[i]>0.0) n=0; else n++;
    }
    if(n==5) {
      if(verbose>1) printf("leaving %d/%d points.\n", i-3, tac.sampleNr);
      tac.sampleNr=i-3;
    }
    /* Remove negative values, if also dtime<=0 */
    i=0; n=0;
    while(i<tac.sampleNr) {
      if(tac.c[0].y[i]<=0.0 && tac.c[4].y[i]<=0.0) {
        tacDeleteSample(&tac, i); n++;
      } else i++;
    }
    if(n>0 && verbose>1) printf("%d negative point(s) deleted.\n", n);
    /* Remove values that are over 1.5 times higher than the average
       of previous and next value, if also dtime<=0 */
    i=1; n=0;
    while(i<tac.sampleNr-1) {
      if(tac.c[4].y[i]>0.0) {i++; continue;}
      a=0.75*(tac.c[0].y[i-1]+tac.c[0].y[i+1]); // 1.5*0.5 = 0.75
      if(tac.c[0].y[i]>a) {tacDeleteSample(&tac, i); n++;} else i++;
    }
    if(n>0 && verbose>1) printf("%d outlier point(s) deleted.\n", n);
    /* If still zeroes in the end, then remove those */
    i=tac.sampleNr-1; while(tac.c[0].y[i]<=0.0) i--;
    if(i<tac.sampleNr-1) {
      if(verbose>1) printf("%d trailing zeroes deleted.\n", tac.sampleNr-i+1);
      tac.sampleNr=i+1;
    }
    /* Set tacNr etc */
    tac.tacNr=1;
    strcpy(tac.c[0].name, "CR");
 
  } else {
    if(copy_cr==0) {
      fprintf(stderr, "Error: not identified as HRRT or HR+ countrate file\n");
      tacFree(&tac); return(3);
    }
    if(verbose>0) printf("not identified as HRRT or HR+ countrate file.\n");
  }
 
 
  /*
   *  Check which samples are smaller than 50% of samples at both side of it
   *  AND smaller than the (average - difference) of these samples
   */
  if(tac.format==TAC_FORMAT_HRRT_HC || tac.format==TAC_FORMAT_HRPLUS_HC) {
    int i, n=0;
    double a, d;
    for(i=1; i<tac.sampleNr-1; i++) {
      d=2.0*tac.c[0].y[i];
      if(d>tac.c[0].y[i-1] || d>tac.c[0].y[i+1]) {
        tac.c[1].y[i]=tac.c[0].y[i]; continue;}
      a=0.5*(tac.c[0].y[i-1]+tac.c[0].y[i+1]);
      d=fabs(tac.c[0].y[i-1]-tac.c[0].y[i+1]);
      if(tac.c[0].y[i]>(a-d)) {tac.c[1].y[i]=tac.c[0].y[i]; continue;}
      tac.c[1].y[i]=nan(""); n++;
    }
    /* replace those with linear interpolation */
    if(n>0) {
      for(i=1; i<tac.sampleNr-1; i++) tac.c[0].y[i]=tac.c[1].y[i];
      tacFixNaNs(&tac);
      if(verbose>1) printf("%d outlier point(s) deleted.\n", n);
    }
  }
 
 
  /* Check whether we have any data left */
  if(tac.sampleNr<1) {
    fprintf(stderr, "Error: invalid data in %s\n", petfile);
    tacFree(&tac); return(3);
  }
 
 
  /* 
   *  Correct for decay 
   */
  if(tac.format==TAC_FORMAT_HRRT_HC || tac.format==TAC_FORMAT_HRPLUS_HC) {
    if(verbose>1) printf("correcting data for physical decay\n");
    /* Convert times to min for decay correction */
    ret=tacXUnitConvert(&tac, UNIT_MIN, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
      tacFree(&tac);
      return(7);
    }
    /* Decay correction */
    for(int i=0; i<tac.sampleNr; i++)
      tac.c[0].y[i]*=decayCorrectionFactorFromIsotope(isotope, 
                         tac.x1[i], tac.x2[i]-tac.x1[i]);
  }
 
  /* Convert times to the final units, sec or min */
  if(tac.tunit==UNIT_UNKNOWN) {
    fprintf(stderr, "Warning: time unit unknown.\n");
  } else {
    ret=tacXUnitConvert(&tac, time_unit, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
      tacFree(&tac);
      return(7);
    }
  }
 
  /* Write the last sample time in file */
  if(tac.format==TAC_FORMAT_HRRT_HC || tac.format==TAC_FORMAT_HRPLUS_HC) {
    char tmp[64];
    sprintf(tmp, "%.1f [%s]", tac.x[tac.sampleNr-1], unitName(tac.tunit));
    iftPut(&tac.h, "cr_duration", tmp, (char)1, NULL);
    if(verbose>0) fprintf(stdout, "cr_duration := %s\n", tmp);
  }
 
  /* Add isotope code */
  tacSetHeaderIsotope(&tac.h, isotopeName(isotope));
  /* and that data is corrected for decay */
  tacSetHeaderDecayCorrection(&tac.h, DECAY_CORRECTED);
 
 
  /* 
   *  Save the data
   */
  if(verbose>1) printf("writing %s\n", crfile);
  FILE *fp; fp=fopen(crfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", crfile);
    tacFree(&tac); return(11);
  }
  ret=tacWrite(&tac, fp, outputformat, 1, &status);
  fclose(fp);
  tacFree(&tac);
  if(ret==TPCERROR_UNSUPPORTED) {
    fprintf(stderr, "Error: writing format %s is not supported.\n",
            tacFormattxt(outputformat));
    return(6);
  }
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    return(7);
  }
  if(verbose>0) printf("CR data written in %s\n", crfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/