tactime.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 "tpcfileutil.h"
#include "tpcift.h"
#include "tpccsv.h"
#include "tpcisotope.h"
#include "tpctac.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.",
  "If 'peak' is given as time, TAC(s) are moved to start at the peak time",
  "or at start of peak frame.", 
  "If 'gap' is given as time, TAC(s) are moved to remove initial gap time.",
  " ",
  "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>",
  "     Isotope must be known for the decay correction.",
  "     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 := C-11'.",
  " -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, tacmove, taccut",
  " ",
  "Keywords: TAC, simulation, late scan, input, time delay, peak",
  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;
  isotope  isot=ISOTOPE_UNKNOWN;
  int      change_decay=0;
  int      keep_negat=0;
  int      keep_times=0;
  int      fill_gap=0;
  char     tacfile1[FILENAME_MAX], tacfile2[FILENAME_MAX];
  int      peak_start=0; // 1=move TAC to start at peak
  int      gap_start=0; // 1=move TAC so that initial gap is removed
  int      time_sign=1;
  double   time_diff=0.0;
  int      time_diff_unit=UNIT_UNKNOWN;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile1[0]=tacfile2[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, "DECAY")==0) {
      change_decay=1; continue;
    } else if(strncasecmp(cptr, "I=", 2)==0) {
      cptr+=2;
      if((isot=isotopeIdentify(cptr))==ISOTOPE_UNKNOWN) {
        fprintf(stderr, "Error: invalid isotope '%s'.\n", cptr);
        return(1);
      }
      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; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-3;
  
  /* 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);}
 
  /* The first argument (non-option) is the file name */
  if(ai<argc) {strlcpy(tacfile1, argv[ai], FILENAME_MAX); ai++;}
  else {fprintf(stderr, "Error: missing file name.\n"); return(1);}
  if(!fileExist(tacfile1)) {
    fprintf(stderr, "Error: file '%s' does not exist.\n", tacfile1);
    return(1);
  }
 
  /* The second argument is the time difference */
  if(ai<argc) {
    if(verbose>3) printf("argument for time difference: '%s'\n", argv[ai]);
    /* maybe user wants to move TAC to start from the peak? */
    if(!strcasecmp("PEAK", argv[ai]) || !strcasecmp("-PEAK", argv[ai])) {
      peak_start=1; keep_times=0;
    } else if(!strcasecmp("GAP", argv[ai])) {
      gap_start=1;
    } else if(fileExist(argv[ai])) { /* maybe time is given in a file? */
      IFT ift; iftInit(&ift);
      FILE *fp=fopen(argv[ai], "r");
      if(fp==NULL) {
        fprintf(stderr, "Error: cannot open file %s\n", argv[ai]);
        return(2);
      }
      int ret=iftRead(&ift, fp, 1, 1, &status); fclose(fp);
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "Error: invalid contents in %s\n", argv[ai]);
        iftFree(&ift); return(2);
      }
      /* Find the specific key string */
      int i=iftFindKey(&ift, "scan_start_time_difference", 0); if(i>=0) time_sign=+1;
      if(i<0) {i=iftFindKey(&ift, "time_difference", 0); if(i>=0) time_sign=+1;}
      if(i<0) {i=iftFindKey(&ift, "Ta", 0); if(i>=0) time_sign=-1;}
      if(i<0) {i=iftFindKey(&ift, "start_time", 0); if(i>=0) time_sign=-1;}
      if(i<0) {
        fprintf(stderr, "Error: keyword for time difference not found in '%s'.\n", argv[ai]);
        iftFree(&ift); return(2);
      }
      /* Read the time and its unit */
      if(iftGetDoubleWithUnit(&ift, i, &time_diff, &time_diff_unit)!=0) {
        fprintf(stderr, "Error: invalid format for time in '%s'.\n", argv[ai]);
        iftFree(&ift); return(2);
      }
      time_diff*=time_sign;
      iftFree(&ift);
    } else { /* then try to read it directly as a value */
      if(atofCheck(argv[ai], &time_diff)!=0) {
        fprintf(stderr, "Error: invalid argument for time: '%s'.\n", argv[ai]);
        return(1);
      }
    }
    if(verbose>7) printf("time_diff := %g\n", time_diff);
    if(time_diff==0.0 && !peak_start) {
      fprintf(stderr, "Warning: time change is zero.\n"); fflush(stderr);
    }
    ai++;
  } else {
    fprintf(stderr, "Error: missing time.\n");
    return(1);
  }
 
  /* Third argument, if it exists, is the output file name */
  if(ai<argc) {strlcpy(tacfile2, argv[ai], FILENAME_MAX); ai++;}
  else {strcpy(tacfile2, tacfile1);}
 
  /* check that there are no extra arguments */
  if(ai<argc) {fprintf(stderr, "Error: extra command-line argument.\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("tacfile1 := %s\n", tacfile1);
    printf("tacfile2 := %s\n", tacfile2);
    printf("keep_times := %d\n", keep_times);
    printf("keep_negat := %d\n", keep_negat);
    printf("fill_gap := %d\n", fill_gap);
    printf("change_decay := %d\n", change_decay);
    if(isot!=ISOTOPE_UNKNOWN) printf("isotope := %s\n", isotopeName(isot));
    printf("peak_start := %d\n", peak_start);
    printf("gap_start := %d\n", gap_start);
    if(time_diff!=0.0) {
      printf("time_diff := %g\n", time_diff);
      if(time_diff_unit!=UNIT_UNKNOWN) printf("time_diff_unit := %s\n", unitName(time_diff_unit));
    }
    fflush(stdout);
  }
 
 
  /* 
   *  Read the file
   */
  if(verbose>1) printf("reading %s\n", tacfile1);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile1, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(3);
  }
  if(verbose>1) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    if(tac.isframe) printf("frames := yes\n"); else printf("frames := no\n");
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    fflush(stdout);
  }
  if(tacSortByTime(&tac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(3);
  }
 
  /* Get and check TAC data range */
  {
    double xmin, xmax;
    if(tacXRange(&tac, &xmin, &xmax)) {
      fprintf(stderr, "Error: invalid sample range.\n");
      tacFree(&tac); return(3);
    }
    if(verbose>1) printf("x_range: %g - %g\n", xmin, xmax);
    if(!(xmax>xmin)) {fprintf(stderr, "Warning: check the sample times.\n"); fflush(stderr);}
    /* If initial gap is to be removed, then set the time difference */
    if(gap_start==1) {
      time_diff=-xmin;
      if(verbose>1) printf("gap_time_diff := %g\n", time_diff);
    }
  }
  int ymax_index;
  {
    double ymin, ymax;
    if(tacYRange(&tac, -1, &ymin, &ymax, NULL, &ymax_index, NULL, NULL)) {
      fprintf(stderr, "Error: invalid concentration range.\n");
      tacFree(&tac); return(3);
    }
    if(verbose>1) printf("y_range: %g - %g\n", ymin, ymax);
    if(!(ymax>ymin)) {fprintf(stderr, "Warning: check the data.\n"); fflush(stderr);}
  }
  /* Set time_diff based on the peak */
  if(peak_start) {
    if(tac.isframe) time_diff=-tac.x1[ymax_index]; else time_diff=-tac.x[ymax_index];
    time_diff_unit=tac.tunit;
    if(verbose>1) {
      printf("peak_index := %d\n", ymax_index);
      printf("time_diff := %g\n", time_diff);
      if(time_diff_unit!=UNIT_UNKNOWN) printf("time_diff_unit := %s\n", unitName(time_diff_unit));
    }
  }
 
 
  /* If datafile contains valid isotope, then check that it is the same as given by user */
  {
    isotope fisot=tacGetIsotope(&tac);
    if(verbose>3) {printf("tac.isotope :=  %s\n", isotopeName(fisot)); fflush(stdout);}
    if(isot==ISOTOPE_UNKNOWN && fisot!=ISOTOPE_UNKNOWN) {
      isot=fisot;
      if(verbose>1) printf("isotope := %s\n", isotopeName(isot));
    } else if(isot!=ISOTOPE_UNKNOWN && fisot==ISOTOPE_UNKNOWN) {
      fisot=isot;
      tacSetIsotope(&tac, fisot);
    } else if(isot!=ISOTOPE_UNKNOWN && fisot!=ISOTOPE_UNKNOWN && isot!=fisot) {
      fprintf(stderr, "Error: different isotope in %s\n", tacfile1);
      tacFree(&tac); return(2);
    }
  }
 
  /* Check that we now have the half life, if decay needs to be changed */
  if(change_decay && isot==ISOTOPE_UNKNOWN) {
    fprintf(stderr, "Error: isotope not specified.\n");
    tacFree(&tac); return(1);
  }
 
  /* Check that we know the time unit, if decay needs to be changed */
  if(change_decay && !unitIsTime(tac.tunit)) {
    fprintf(stderr, "Error: sample time unit unknown.\n");
    tacFree(&tac); return(1);
  }
 
  /* Convert time difference units to the same as in TAC file */
  if(unitIsTime(time_diff_unit) && tac.tunit!=UNIT_UNKNOWN) {
    double f=unitConversionFactor(time_diff_unit, tac.tunit);
    if(isnan(f)) {
      fprintf(stderr, "Error: incompatible time units.\n");
      tacFree(&tac); return(2);
    }
    time_diff*=f; time_diff_unit=tac.tunit;
  }
  if(verbose>1) printf("final_time_diff := %g %s\n", time_diff, unitName(tac.tunit));
 
 
 
  /*
   *  Change the sample times
   */
  if(keep_times==0) {
    for(int i=0; i<tac.sampleNr; i++) {
      tac.x[i]+=time_diff; tac.x1[i]+=time_diff; tac.x2[i]+=time_diff;
    }
  }
  /* Remove samples with negative times */
  if(keep_negat==0) {
    int i=0;
    while(tac.sampleNr>0 && tac.x[i]<0.0) tacDeleteSample(&tac, i); 
    if(tac.sampleNr<=0) {
      fprintf(stderr, "Error: No positive times left.\n");
      tacFree(&tac); return(5);
    }
  }
  /* Make sure that there is no initial time gap, if requested */
  if(fill_gap!=0) {
    if(verbose>1) printf("filling initial time gap\n");
    if(tacAddZeroSample(&tac, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot fill time gap.\n");
      tacFree(&tac); return(5);
    }
  }
 
 
  /*
   *  Change decay correction for the time difference, if requested
   */
  if(change_decay && time_diff!=0.0) {
    if(verbose>3) printf("decay correcting for time %g\n", time_diff);
    /* Calculate decay correction factor */
    double tdiff=fabs(time_diff)*unitConversionFactor(tac.tunit, UNIT_MIN);
    if(verbose>5) printf("  tdiff := %g\n", tdiff);
    double dcf=decayCorrectionFactorFromIsotope(isot, tdiff, 0.0);
    if(time_diff<0.0) dcf=1.0/dcf;
    if(isnan(dcf)) {
      fprintf(stderr, "Error: cannot calculate decay correction factor.\n");
      tacFree(&tac); return(1);
    }
    if(verbose>1) printf("decay_correction_factor := %g\n", dcf);
    /* Correct all TACs */
    for(int j=0; j<tac.tacNr; j++) for(int i=0; i<tac.sampleNr; i++) tac.c[j].y[i]*=dcf;
    /* Opposite correction for weights, if available */
    if(tacIsWeighted(&tac)) {
      for(int i=0; i<tac.sampleNr; i++) tac.w[i]/=dcf; 
      tacWeightNorm(&tac, NULL);
    }
  }
 
 
  /*
   *  If sample times were not changed, and datafile contains injection time,
   *  then change the injection time accordingly
   */
  if(keep_times!=0 && tacGetHeaderInjectiontime(&tac.h, NULL, NULL)==TPCERROR_OK) {
    char buf[20];
    tacGetHeaderInjectiontime(&tac.h, buf, &status);
    if(verbose>2) printf("injection_time: %s\n", buf);
    double tdiff=-time_diff*unitConversionFactor(tac.tunit, UNIT_SEC);
    if(verbose>5) printf("  injection_time_tdiff := %g s\n", tdiff);
    if(strDateTimeAdd((int)round(tdiff), buf) || tacSetHeaderInjectiontime(&tac.h, buf)!=TPCERROR_OK)
    {
      fprintf(stderr, "Error: cannot correct injection time.\n");
      tacFree(&tac); return(1);
    }
    if(verbose>0) printf("updated_injection_time := %s\n", buf);
  }
 
 
 
  /*
   *  Save the modified data
   */
  if(verbose>1) printf("saving modified data in %s\n", tacfile2);
   {
    FILE *fp; fp=fopen(tacfile2, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing\n");
      tacFree(&tac); return(11);
    }
    int ret=tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(12);
    }
  
    /* Tell user what we did */
    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 decay of %s", isotopeName(isot));
      fprintf(stdout, " in %s\n", tacfile2);
      fflush(stdout);
    }
  }
  tacFree(&tac);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/