tactime.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Increase or decrease the sample (frame) times in regional or blood/plasma",
  "TAC files.",
  "Samples with negative times are not saved by default.",
  "Time must be given in the same units that are used in the datafile.",
  " ",
  "Usage: @P [Options] tacfile [-]time outputfile",
  " ",
  "Time can be given as positive or negative value directly in the command",
  "line, or in an ASCII file which contains a line with the time change value",
  "in the following format: 'time_difference := time'.",
  "Alternatively, time specified in that file with 'Ta' or 'start_time' is",
  "subtracted from sample times.",
  " ",
  "Options:",
  " -decay",
  "     Physical decay correction is changed with change in sample times;",
  "     without this option, radioactivity values are not changed.",
  "     Note that this option will provide correct result only if time unit",
  "     setting in datafile is correct.",
  " -i=<isotope>",
  "     If tacfile does not contain the isotope, then decay correction cannot",
  "     applied unless isotope is given with this option.",
  "     Accepted isotope codes are for example F-18, C-11, and O-15.",
  "     Isotope code can also be specified in input file in format",
  "     '# isotope := Isotope'.",
  " -keepnegat",
  "     Samples with negative sample times are not removed from output.",
  " -keeptimes",
  "     While correction for physical decay is changed with option -decay,",
  "     sample times will not be changed.",
  " -nogap",
  "     Possible gap between time zero and first sample is filled.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: tacframe, tacdecay, injdifft, tacunit, fitdelay, simdisp, taccut",
  " ",
  "Keywords: TAC, modelling, simulation, late scan, input, time delay",
  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, fj, ri, n, ret;
  int change_decay=0;
  int keep_negat=0;
  int keep_times=0;
  int fill_gap=0;
  int time_sign=1;
  char *cptr, dftfile1[FILENAME_MAX], dftfile2[FILENAME_MAX], tmp[64];
  char isotope[64];
  double time_diff=0.0;
  double halflife=0.0, lambda, dcf, f;
  int time_diff_unit=TUNIT_UNKNOWN;
  DFT dft;
  IFT ift;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dftfile1[0]=dftfile2[0]=isotope[0]=(char)0; 
  iftInit(&ift); 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;
    cptr=argv[ai]+1;
    if(strcasecmp(cptr, "DECAY")==0) {
      change_decay=1; continue;
    } else if(strncasecmp(cptr, "D=", 2)==0) { // deprecated option
      change_decay=1;
      cptr+=2; if(!*cptr) continue;
      halflife=hlFromIsotope(cptr);
      if(halflife<=0.0) {
        fprintf(stderr, "Error: invalid isotope '%s'.\n", cptr);
        return(1);
      }
      strcpy(isotope, hlCorrectIsotopeCode(cptr));
      continue;
    } else if(strncasecmp(cptr, "I=", 2)==0) {
      //change_decay=1;
      cptr+=2; if(!*cptr) continue;
      halflife=hlFromIsotope(cptr);
      if(halflife<=0.0) {
        fprintf(stderr, "Error: invalid isotope '%s'.\n", cptr);
        return(1);
      }
      strcpy(isotope, hlCorrectIsotopeCode(cptr));
      continue;
    } else if(strncasecmp(cptr, "KEEPNEGATIVE", 5)==0) {
      keep_negat=1; continue;
    } else if(strncasecmp(cptr, "KEEPTIMES", 5)==0) {
      keep_times=1; continue;
    } else if(strncasecmp(cptr, "NOGAP", 5)==0) {
      fill_gap=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 */
  /* get data file name */
  if(ai<argc) {
    if(access(argv[ai], 0) == -1) {
      fprintf(stderr, "Error: file '%s' does not exist.\n", argv[ai]);
      return(1);
    }
    strcpy(dftfile1, argv[ai]); ai++;
  }
  /* get time difference */
  if(ai<argc) {
    /* maybe time is given in a file? */
    ret=iftRead(&ift, argv[ai], 1);
    if(ret==0) { /* yes it seems so */
      strcpy(tmp, "scan_start_time_difference");
      n=iftGet(&ift, tmp); time_sign=+1;
      if(n==-1) {
        strcpy(tmp, "time_difference"); n=iftGet(&ift, tmp); time_sign=+1;
      }
      if(n==-1) {
        strcpy(tmp, "Ta"); n=iftGet(&ift, tmp); time_sign=-1;
      }
      if(n==-1) {
        strcpy(tmp, "start_time"); n=iftGet(&ift, tmp); time_sign=-1;
      }
      if(n<0) {
        fprintf(stderr, 
          "Error: keyword for time difference not found in '%s'.\n", argv[ai]);
        iftEmpty(&ift); return(2);
      }
      ret=atof_with_check(ift.item[n].value, &time_diff);
      if(ret!=0) {
        fprintf(stderr, "Error: invalid format for time in '%s'.\n", argv[ai]);
        iftEmpty(&ift); return(2);
      }
      time_diff*=time_sign;
      /* Try to get its units too */
      if(strstr(ift.item[n].value, "[s]")!=NULL) time_diff_unit=TUNIT_SEC;
      else if(strstr(ift.item[n].value, "[min]")!=NULL) time_diff_unit=TUNIT_MIN;
    } else { /* probably not, try to read it directly as a value */
      ret=atof_with_check(argv[ai], &time_diff);
      if(ret!=0) {
        fprintf(stderr, "Error: invalid argument for time: '%s'.\n", argv[ai]);
        return(1);
      }
    }
    if(verbose>1) printf("time_diff := %g\n", time_diff);
    if(time_diff==0.0) {
      fprintf(stderr, "Warning: time change is zero.\n");
    }
    ai++;
  }
  /* get output file name */
  if(ai<argc) {
    strcpy(dftfile2, argv[ai]);
    ai++;
  }
  /* check that there are no extra arguments */
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!dftfile1[0]) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n",
      argv[0]);
    return(1);
  }
  if(!dftfile2[0]) {
    fprintf(stderr, "Error: undefined output file name.\n");
    return(1);
  }
 
  /* Check options */
  if(keep_times==1 && change_decay==0) {
    fprintf(stderr, 
    "Error: change in neither sample times or decay correction were requested.\n");
    return(1);
  }
  if(keep_times==1) keep_negat=1;
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dftfile1 := %s\n", dftfile1);
    printf("dftfile2 := %s\n", dftfile2);
    printf("time_diff := %g\n", time_diff);
    printf("change_decay := %d\n", change_decay);
    printf("halflife := %g\n", halflife);
    if(isotope[0]) printf("isotope := %s\n", isotope);
    printf("keep_times := %d\n", keep_times);
    printf("keep_negat := %d\n", keep_negat);
    printf("fill_gap := %d\n", fill_gap);
  }
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading %s\n", dftfile1);
  if(dftRead(dftfile1, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dftfile1, dfterrmsg);
    return(3);
  }
  if(verbose>2) {
    printf("voiNr := %d\n", dft.voiNr);
    printf("frameNr := %d\n", dft.frameNr);
  }
 
  /* Get halflife from isotope setting in data file */
  if(change_decay && strlen(dft.isotope)>0) {
    f=hlFromIsotope(dft.isotope);
    if(f<=0.0) {
      fprintf(stderr, "Error: invalid isotope '%s' in '%s'.\n",
        dft.isotope, dftfile1);
      dftEmpty(&dft); return(3);
    }
    /* check that command-line and file isotope settings are matching */
    if(halflife>0.0) {
      if(halflife!=f) {
        fprintf(stderr, "Error: invalid isotope in '%s' or in command line.\n",
          dftfile1);
        dftEmpty(&dft); return(1);
      }
    } else {
      halflife=f;
    }
  }
 
  /* Check that we now have the halflife, if decay needs to be changed */
  if(change_decay && halflife<=0.0) {
    fprintf(stderr, "Error: isotope not specified.\n");
    dftEmpty(&dft); return(1);
  }
 
  /* Set isotope code in the file, if necessary */
  if(isotope[0] && strlen(dft.isotope)<1) {
    strcpy(dft.isotope, isotope);
    //strcat(dft.comments, "# isotope := "); strcat(dft.comments, isotope);
  }
  if(!isotope[0]) strcpy(isotope, dft.isotope);
  if(verbose>1) printf("isotope := %s\n", isotope);
 
  /* If necessary, change the halflife time unit based on time unit in data */
  if(halflife>0.0 && dft.timeunit==TUNIT_SEC) {
    halflife*=60.0;
  }
  if(verbose>1) printf("halflife := %g\n", halflife);
 
 
  /* Convert time diff units to the same as in DFT file */
  if(time_diff_unit==TUNIT_MIN && dft.timeunit==TUNIT_SEC) {
    time_diff*=60.0; time_diff_unit=TUNIT_SEC;
  } else if(time_diff_unit==TUNIT_SEC && dft.timeunit==TUNIT_MIN) {
    time_diff/=60.0; time_diff_unit=TUNIT_MIN;
  }
  if(verbose>1) printf("final_time_diff := %g\n", time_diff);
 
  /*
   *  Change the sample times
   */
  if(keep_times==0) {
    for(fi=0; fi<dft.frameNr; fi++) {
      dft.x[fi]+=time_diff; dft.x1[fi]+=time_diff; dft.x2[fi]+=time_diff;
    }
  }
  /* Remove samples with negative times */
  if(keep_negat==0) {
    for(fi=0; fi<dft.frameNr; fi++) if(dft.x[fi]>=0.0) break;
    n=fi; if(n>0) for(fj=fi; fj<dft.frameNr; fj++) {
      for(ri=0; ri<dft.voiNr; ri++) dft.voi[ri].y[fj-n]=dft.voi[ri].y[fj];
      dft.x[fj-n]=dft.x[fj]; dft.x1[fj-n]=dft.x1[fj]; dft.x2[fj-n]=dft.x2[fj];
      dft.w[fj-n]=dft.w[fj];
    }
    dft.frameNr-=n;
    if(dft.frameNr<=0) {
      fprintf(stderr, "Error: No positive times left.\n");
      dftEmpty(&dft); return(4);
    }
  }
  if(verbose>10) dftPrint(&dft);
 
 
  /*
   *  Change the correction for physical decay, if required
   */
  if(change_decay) {
    lambda=hl2lambda(halflife); if(time_diff<0.0) lambda=-lambda;
    dcf=hlLambda2factor(lambda, fabs(time_diff), -1.0);
    if(verbose>0)
      printf("correction factor for physical decay: %g\n", dcf);
    for(fi=0; fi<dft.frameNr; fi++) for(ri=0; ri<dft.voiNr; ri++)
      dft.voi[ri].y[fi]*=dcf;
  }
 
 
  /*
   *  If sample times were not changed, and datafile contains injection time,
   *  then change the injection time accordingly
   */
  if(keep_times!=0 && strlen(dft.injectionTime)>18) {
    time_t t_injection;
    struct tm tm_injection;
    /* Convert time string into time_t */
    ret=get_datetime(dft.injectionTime, &tm_injection, verbose-2);
    if(ret==0) {
      t_injection=timegm(&tm_injection);
      if(t_injection<0) ret=-100;
    }
    if(ret!=0) {
      fprintf(stderr, "Error: invalid injection time.\n");
      dftEmpty(&dft); return(6);
    }
    /* Add time_diff in seconds */
    if(dft.timeunit==TUNIT_SEC) tmAdd(-time_diff, &tm_injection);
    else if(dft.timeunit==TUNIT_MIN) tmAdd(-60.0*time_diff, &tm_injection);
    snprintf(dft.injectionTime, 20, "%.4d-%.2d-%.2d %.2d:%.2d:%.2d",
      tm_injection.tm_year+1900, tm_injection.tm_mon+1, tm_injection.tm_mday,
      tm_injection.tm_hour, tm_injection.tm_min, tm_injection.tm_sec);
    if(verbose>0) printf("updated_injectiontime := %s\n", dft.injectionTime);
  }
 
  /*
   *  Make sure that there is no initial gap
   */
  if(fill_gap!=0) {
    if(verbose>1) printf("checking for initial gap\n");
    ret=dftFillInitialGap(&dft); if(ret!=0) {
      fprintf(stderr, "Error: invalid injection time.\n");
      dftEmpty(&dft); return(9);
    } 
  }
 
 
  /*
   *  Save the modified data
   */
  if(verbose>1) printf("saving modified data in %s\n", dftfile2);
  /* Set comments */
  strcpy(dft.comments, ""); dftSetComments(&dft);
  /* Write the file */
  if(dftWrite(&dft, dftfile2)) {
    fprintf(stderr, "Error in writing '%s': %s\n", dftfile2, dfterrmsg);
    dftEmpty(&dft); return(11);
  }
  /* Quit and tell what was done */
  dftEmpty(&dft);
  if(verbose>=0) {
    if(keep_times==0) fprintf(stdout, "Sample times changed by %g", time_diff);
    else fprintf(stdout, "Injection time changed by %g", -time_diff);
    if(change_decay) fprintf(stdout, " and values with half-life %g", halflife);
    fprintf(stdout, " in %s\n", dftfile2);
  }
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/