simframe.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 "tpccsv.h"
#include "tpctac.h"
#include "tpcisotope.h"
#include "tpcli.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "PET TACs are simulated with short sample time intervals.",
  "This program sums up those points within specified time frames to simulate",
  "a measured PET tissue uptake curve.",
  " ",
  "Usage: @P [options] tacfile framefile newfile [isotope]",
  " ",
  "Options:",
  " -sec or -min",
  "     Sample times are known to be in seconds or minutes, but is not",
  "     specified or is wrong in TAC file.", 
  " -mid",
  "     Frame mid time is written in output file instead of start and end times.",
  " -i  Calculates integrals at frame mid times, instead of frame averages;",
  "     not available with [isotope].",
  " -ii Calculates 2nd integrals at frame mid times, instead of frame averages;",
  "     not available with [isotope].",
  " --force",
  "     Program allows extensive extrapolation.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "TAC file can contain more than one TAC.",
  "Time frames data can be in SIF or TAC format. Alternatively, file can consist",
  "of one or two columns of data, containing either 1) frame durations, or",
  "2) frame start times and frame durations; frame time units must be same",
  "as in the datafile. Frames are allowed to overlap.",
  " ",
  "If the isotope is specified, the correction for physical decay is at first",
  "removed, then PET framing is simulated, and after that the framed data",
  "is decay corrected again, based on the frame start time and length as is",
  "the normal procedure when collecting PET image data.",
  "Verify that time units are correct when using this possibility.",
  " ",
  "See also: interpol, tacframe, tac4frpl, tactime, fit2dat, sim_3tcm",
  " ",
  "Keywords: TAC, SIF, simulation, interpolation, time frame",
  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;
  unit       knownTimeunit=UNIT_UNKNOWN;
  isotope    isot=ISOTOPE_UNKNOWN;
  int        mode=0; // 0=mean, 1=integral, 2=2nd integral
  int        forceMode=0; // 0=off, 1=on
  int        midFrame=0; // save with mid frame times (1) or not (0)
  char       datfile[FILENAME_MAX], frafile[FILENAME_MAX], outfile[FILENAME_MAX];
 
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  datfile[0]=frafile[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, "C")==0 || strcasecmp(cptr, "Y")==0) {
      mode=0; continue;
    } else if(strcasecmp(cptr, "I")==0) {
      mode=1; continue;
    } else if(strcasecmp(cptr, "II")==0) {
      mode=2; continue;
    } else if(strncasecmp(cptr, "MINUTES", 3)==0) {
      knownTimeunit=UNIT_MIN; continue;
    } else if(strncasecmp(cptr, "SECONDS", 3)==0) {
      knownTimeunit=UNIT_SEC; continue;
    } else if(strcasecmp(cptr, "MID")==0) {
      midFrame=1; continue;
    } else if(strcasecmp(cptr, "F")==0 || strcasecmp(cptr, "FORCE")==0) {
      forceMode=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-1;
  
  /* 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(datfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(frafile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(outfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if((isot=isotopeIdentify(argv[ai]))==ISOTOPE_UNKNOWN) {
      fprintf(stderr, "Error: invalid isotope argument '%s'.\n", argv[ai]);
      return(1);
    }
    if(mode>0) {
      fprintf(stderr, "Error: decay option cannot be used with integrals.\n");
      return(1);
    }
    ai++;    
  }
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!outfile[0]) {
    fprintf(stderr, "Error: missing file name; use option --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("datfile := %s\n", datfile);
    printf("frafile := %s\n", frafile);
    printf("outfile := %s\n", outfile);
    printf("mode := %d\n", mode);
    printf("forceMode := %d\n", forceMode);
    if(knownTimeunit!=UNIT_UNKNOWN) printf("knownTimeunit := %s\n", unitName(knownTimeunit));
    printf("midFrame := %d\n", midFrame);
    if(isot!=ISOTOPE_UNKNOWN) printf("isotope := %s\n", isotopeName(isot));
    fflush(stdout);
  }
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) {fprintf(stdout, "reading %s\n", datfile); fflush(stdout);}
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, datfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), datfile);
    tacFree(&tac); return(2);
  }
  if(verbose>2) {
    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");
    fflush(stdout);
  }
  tacDeleteMissingSamples(&tac);
  if(tacSortByTime(&tac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: invalid sample times.\n");
    tacFree(&tac); return(2);
  }
  /* Correct the time unit, if given by user */
  if(knownTimeunit!=UNIT_UNKNOWN) {
    tac.tunit=knownTimeunit;
  } else if(!unitIsTime(tac.tunit)) {
    double xmin, xmax;
    if(tacXRange(&tac, &xmin, &xmax)) {
      fprintf(stderr, "Error: invalid sample times.\n");
      tacFree(&tac); return(2);
    }
    if(xmax>30.0 && lambdaFromIsotope(isot)>0.3) {
      tac.tunit=UNIT_SEC;
      fprintf(stderr, "Warning: assuming that sample times are in seconds.\n");
    } else {
      tac.tunit=UNIT_MIN;
      fprintf(stderr, "Warning: assuming that sample times are in minutes.\n");
    }
  }
  if(verbose>2) {
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    fflush(stdout);
  }
 
 
  /* If isotope was given, remove decay correction */
  if(isot!=ISOTOPE_UNKNOWN) {
    if(verbose>1) printf("removing decay correction.\n");
    if(tacDecayCorrection(&tac, isot, 0, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), datfile);
      tacFree(&tac); return(2);
    }
  }
 
 
  /*
   *  Read frame data
   */
  if(verbose>1) printf("reading frame data in %s\n", frafile);
  TAC sif; tacInit(&sif);
  {
    /* Open the file */
    FILE *fp;
    fp=fopen(frafile, "r");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open %s\n", frafile);
      tacFree(&tac); return(3);
    }
    /* Read contents into CSV structure */
    CSV csv; csvInit(&csv);
    if(csvRead(&csv, fp, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), frafile);
      tacFree(&tac); csvFree(&csv); fclose(fp); return(3);
    }
    fclose(fp);
    /* Try to read as SIF */
    if(tacReadSIF(&sif, &csv, NULL, &status)==TPCERROR_OK) {
      if(verbose>2) printf("frame times from SIF\n");
      csvFree(&csv);
    } else {
      if(verbose>2) printf("not SIF, trying to read as TAC file.\n");
      csvFree(&csv);
      if(tacRead(&sif, frafile, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), frafile);
        tacFree(&tac); return(3);
      }
      if(verbose>2) {
        printf("fileformat := %s\n", tacFormattxt(sif.format));
        printf("tacNr := %d\n", sif.tacNr);
        printf("sampleNr := %d\n", sif.sampleNr);
        printf("xunit := %s\n", unitName(sif.tunit));
        printf("isframe := %d\n", sif.isframe);
        if(verbose>3) tacWrite(&sif, stdout, TAC_FORMAT_TSV_UK, 1, NULL);
        fflush(stdout);
      }
 
      if(sif.tacNr>0 && sif.isframe) {
        if(verbose>2) printf("seems to be valid TAC file\n");
      } else if(sif.isframe==0 && (sif.format==TAC_FORMAT_PMOD || sif.format==TAC_FORMAT_DFT)) {
        fprintf(stderr, "Error: no frame start and end times in %s\n", frafile);
        tacFree(&tac); tacFree(&sif); return(3);
      } else if(sif.tacNr==1 && sif.isframe==0 && tacYNaNs(&sif, 0)==sif.sampleNr) {
        if(verbose>2) printf("just one column containing frame lengths\n");
        sif.x1[0]=0.0; sif.x2[0]=sif.x[0]; sif.x[0]=0.5*(sif.x1[0]+sif.x2[0]);
        for(int i=1; i<sif.sampleNr; i++) {
          sif.x1[i]=sif.x2[i-1]; sif.x2[i]=sif.x1[i]+sif.x[i]; sif.x[i]=0.5*(sif.x1[i]+sif.x2[i]);
        }
        sif.isframe=1;
      } else if(sif.tacNr==1 && sif.isframe==0) {
        if(verbose>2) printf("two columns, probably containing frame start times and lengths\n");
        for(int i=0; i<sif.sampleNr; i++) {
          sif.x1[i]=sif.x[i]; sif.x2[i]=sif.x1[i]+sif.c[0].y[i]; sif.x[i]=0.5*(sif.x1[i]+sif.x2[i]);
        }
        sif.isframe=1;
      } else {
        fprintf(stderr, "Error: invalid format in %s\n", frafile);
        tacFree(&tac); tacFree(&sif); return(3);
      }
      if(verbose>3) tacWrite(&sif, stdout, TAC_FORMAT_TSV_UK, 1, NULL);
    }
    /* If time unit is known, convert to same units as in the data to be edited */
    if(sif.tunit!=UNIT_UNKNOWN) {
      if(tacXUnitConvert(&sif, tac.tunit, &status)!=TPCERROR_OK)
        fprintf(stderr, "Warning: %s.\n", errorMsg(status.error));
    }
  }
  /* Warning, if it looks like time units are different in frame frame */
  /* Check for the need to extrapolate */
  {
    double xmin1, xmax1, xmin2, xmax2;
    if(tacXRange(&tac, &xmin1, &xmax1)) {
      fprintf(stderr, "Error: invalid sample times in %s.\n", datfile);
      tacFree(&tac); tacFree(&sif); return(2);
    }
    if(tacXRange(&sif, &xmin2, &xmax2)) {
      fprintf(stderr, "Error: invalid frame times in %s.\n", frafile);
      tacFree(&tac); tacFree(&sif); return(3);
    }
    if(xmax2<0.1*xmax1 || xmax2>10.*xmax1) {
      fprintf(stderr, "Warning: PET frame times may not be in same unit as TAC data.\n");
      fflush(stderr);
    }
    if(verbose>0 && xmax1<xmax2) {
      printf("Note: extrapolation needed from %g to %g\n", xmax1, xmax2); fflush(stdout);}
    if(xmax1<0.95*xmax2 || (isot!=ISOTOPE_UNKNOWN && xmax1<0.99*xmax2)) {
      if(forceMode==0) {
        fprintf(stderr, "Error: required extrapolation is too risky.\n");
        tacFree(&tac); tacFree(&sif); return(3);
      }
    }
    if(forceMode==0 && xmin2<0.95*xmin1) {
      printf("Note: extrapolation needed from %g to %g\n", xmin1, xmin2); fflush(stdout);
      /* If mean value of first sample is <=0 then there's no problem */
      double a=0.0;
      for(int i=0; i<tac.tacNr; i++) a+=tac.c[i].y[0]; 
      a/=(double)tac.tacNr;
      if(a>1.0E-12) {
        fprintf(stderr, "Error: required extrapolation is too risky.\n");
        tacFree(&tac); tacFree(&sif); return(3);
      }
    }
  }
 
 
  /* 
   *  Allocate memory for output data
   */
  TAC tac2; tacInit(&tac2);
  if(tacDuplicate(&tac, &tac2)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot setup new TAC data.\n");
    tacFree(&tac); tacFree(&sif); return(4);
  }
  if(sif.sampleNr>tac.sampleNr && tacAllocateMoreSamples(&tac2, sif.sampleNr-tac.sampleNr)) {
    fprintf(stderr, "Error: cannot setup new TAC data.\n");
    tacFree(&tac); tacFree(&sif); return(4);
  }
  tac2.sampleNr=sif.sampleNr;
  /* Set output frame times */
  tacXCopy(&sif, &tac2, 0, tac2.sampleNr-1);
  if(midFrame) tac2.isframe=0; else tac2.isframe=1; 
  /* Simple file format cannot store frame start and end times; change that, when necessary */
  if(tac2.isframe && tac2.format==TAC_FORMAT_SIMPLE) {
    char *ex=filenameGetExtension(outfile); if(ex!=NULL) ex++;
    tac2.format=tacFormatIdentify(ex);
    if(tac2.format==TAC_FORMAT_UNKNOWN) tac2.format=TAC_FORMAT_PMOD;
  }
  /* Make sure that file format writing is supported */
  if(!tacFormatWriteSupported(tac2.format)) {
    tac2.format=TAC_FORMAT_PMOD;
  }
 
 
  /* 
   *  Interpolate
   */
  if(verbose>1) {printf("interpolating\n"); fflush(stdout);}
  {
    int ret=0;
    for(int i=0; i<tac2.tacNr; i++) {
      double *y=NULL, *yi=NULL, *yii=NULL;
      if(mode==0) y=tac2.c[i].y; else if(mode==1) yi=tac2.c[i].y; else yii=tac2.c[i].y;
      ret=liInterpolateForPET(tac.x, tac.c[i].y, tac.sampleNr, 
                              tac2.x1, tac2.x2, y, yi, yii, tac2.sampleNr, 3, 1, verbose-10);
      if(ret) break;
    }
    if(ret) {
      fprintf(stderr, "Error: cannot interpolate the data.\n");
      if(verbose>1) printf("  ret := %d\n", ret);
      tacFree(&sif); tacFree(&tac); tacFree(&tac2); return(5);
    }
  }
  tacFree(&tac);
 
  /* If necessary, correct again for radioactive decay */
  if(isot!=ISOTOPE_UNKNOWN) {
    //tacWrite(&tac2, stdout, TAC_FORMAT_TSV_UK, 1, NULL);
    if(verbose>1) printf("correcting for decay.\n");
    if(tacDecayCorrection(&tac2, isot, 1, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&sif); tacFree(&tac2); return(6);
    }
  }
 
 
 
 
  /*
   *  Save TACs
   */
  /* Isotope from SIF is not used to do the decay off/on procedure, but it will be written into
     TAC file if it can be found */
  if(isot==ISOTOPE_UNKNOWN) isot=tacGetIsotope(&sif);
  if(isot!=ISOTOPE_UNKNOWN) tacSetIsotope(&tac2, isot);
  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");
      tacFree(&sif); tacFree(&tac2); return(11);
    }
    int ret=tacWrite(&tac2, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&sif); tacFree(&tac2); return(12);
    }
  
    /* Tell user what we did */
    if(verbose>=0) {
      if(mode==0) printf("  %d frames saved.\n", tac2.sampleNr);
      else if(mode==1) printf("  %d integral frames saved.\n", tac2.sampleNr);
      else if(mode==2) printf("  %d 2nd integral frames saved.\n", tac2.sampleNr);
    }
  }
  tacFree(&tac2);
  tacFree(&sif);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/