interpol.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 "tpctac.h"
#include "tpcli.h"
#include "tpcstatist.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Linear (dot-to-dot) interpolation and integration of PET time-activity data.",
  " ",
  "Usage: @P [options] tacfile [outputfile]",
  " ",
  "Options:",
  " -y or -i or -ii",
  "     TACs are interpolated (-y, default), integrated (-i), or the 2nd",
  "     integral (-ii) is calculated.",
  " -X=<x1,x2,...> or -C=<start,stop,step> or -F=<xfile>",
  "     The sample times for interpolated or integral TACs can be specified with",
  "     these options; if none of these is given, the sample times of the",
  "     original datafile are used. Definite sample times can be given after",
  "     -X; regular sample timing can be specified using option -C; with -F the",
  "     sample times can be read from given TAC file.",
  " -A[=<t>]",
  "     Input TACs are interpolated and integrated with sample (0,0) or (t,0)",
  "     as their first measurement points.",
  " -0=<y|F>",
  "     Input TACs are interpolated and integrated with sample (0,y) added",
  "     if TAC starts later than at zero time. With option -0=F the first",
  "     existing sample value is used as y value.",
  " -header=<Yes|no>",
  "     Descriptive title lines (if header information is available) can be",
  "     included (default), or not included.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Input TAC data is considered to be non-framed even if it contains frame",
  "start and end times; if necessary, frame mid times are calculated and used.",
  "The interpolated or integrated values are calculated exactly at the",
  "specified sample times, except when the start and end times of PET frames",
  "are given with -F=<xfile>; then the time frame average values or integrals",
  "at the frame mid time are calculated.",
  " ",
  "If output filename is not specified, the results are written to stdout.",
  " ",
  "If result consists of only one value, output file is not given, and option",
  "-header=no is given, then just the value is printed, without time, for",
  "easier use in scripts. For example in bash:",
  " ",
  "  auc=$(interpol -i -x=50 -header=no plasma.dat)",
  "  printf \"AUC=%s\\n\" $auc ",
  " ",
  "See also: taccalc, dftinteg, simframe, tacframe, ainterp, taccut, extrapol",
  " ",
  "Keywords: TAC, modelling, simulation, tool, AUC, interpolation",
  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          mode=0; /* 0=interpolation; 1=integral; 2=2nd integral */
  char        *cptr, tacfile[FILENAME_MAX], outfile[FILENAME_MAX], 
               xfile[FILENAME_MAX];
  static char *mode_str[] = {"concentration", "integral", "2nd_integral", 0};
  TAC          tac1, tac2;
  double       start0, startY;
  int          xnr; // length of xlist
  double      *xlist=NULL; // remember to free 
  int          header=1; // save header (1) or do not save (0)
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacInit(&tac1); tacInit(&tac2);
  tacfile[0]=outfile[0]=xfile[0]=(char)0;
  xnr=0;
  start0=startY=nan("");
  /* 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(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, "F=", 2)==0) {
      if(xnr>0) {
        fprintf(stderr,"Error: sample times specified with more than one way.\n");
        free(xlist); return(1);
      }
      strcpy(xfile, cptr+2); if(strlen(xfile)>0) continue;
    } else if(strncasecmp(cptr, "C=", 2)==0) {
      if(xnr>0) {
        fprintf(stderr,"Error: sample times specified with more than one way.\n");
        free(xlist); return(1);
      }
      cptr+=2; 
      /* Read x start, end, and stepsize */
      double v[4];
      if(atofList(cptr, ";,_|: \t", v, 4) != 3) {
        fprintf(stderr, "Error: invalid format with option -C.\n"); 
        return(1);
      }
      if(v[0]>=v[1] || v[0]<0.0 || v[2]<=0.0) {
        fprintf(stderr, "Error: invalid values with option -C.\n"); 
        return(1);
      }
      int n;
      n=1+(int)ceil(fabs(v[1]-v[0])/v[2]);
      if(n<3 || n>10000) {
        fprintf(stderr, "Error: invalid step size.\n"); 
        return(1);
      }
      xlist=(double*)malloc(n*sizeof(double));
      if(xlist==NULL) {
        fprintf(stderr, "Error: out of memory.\n"); 
        return(1);
      }
      /* Fill data from start to stop time */
      double f;
      for(xnr=0; xnr<n; xnr++) {
        f=v[0]+(double)xnr*v[2]; if(f>v[1]) break; else xlist[xnr]=f;}
      continue;
    } else if(strncasecmp(cptr, "X=", 2)==0) {
      if(xnr>0) {
        fprintf(stderr,"Error: sample times specified with more than one way.\n");
        free(xlist); return(1);
      }
      cptr+=2;
      /* Read x values */
      xnr=strTokenNr((char*)cptr, ";,_|: \t"); 
      if(xnr<1 || xnr>100) {
        fprintf(stderr, "Error: invalid format with option -X.\n"); 
        return(1);
      }
      xlist=(double*)malloc(xnr*sizeof(double));
      if(xlist==NULL) {
        fprintf(stderr, "Error: out of memory.\n"); 
        xnr=0; return(1);
      }
      xnr=atofList(cptr, ";,_|: \t", xlist, xnr);
      if(xnr<1) {
        fprintf(stderr, "Error: invalid format with option -X.\n"); 
        free(xlist); xnr=0; return(1);
      }
      statSortDouble(xlist, xnr, 0); // make sure that values are ascending
      continue;
    } else if(strcasecmp(cptr, "A")==0) {
      start0=0.0; continue;
    } else if(strncasecmp(cptr, "A=", 2)==0) {
      start0=atofVerified(cptr+2); if(isfinite(start0)) continue;
    } else if(strncasecmp(cptr, "0=", 2)==0) {
      cptr+=2; 
      if(strcasecmp(cptr, "F")==0) {startY=-1.0E+20; continue;}
      startY=atofVerified(cptr); if(isfinite(startY)) continue;
    } else if(strncasecmp(cptr, "HEADER=", 7)==0) {
      cptr+=7;
      if(strncasecmp(cptr, "YES", 1)==0) {
        header=1; continue;
      } else if(strncasecmp(cptr, "NO", 1)==0) {
        header=0; continue;
      }
    } else if(strncasecmp(cptr, "HDR=", 4)==0) {
      cptr+=4;
      if(strncasecmp(cptr, "YES", 1)==0) {
        header=1; continue;
      } else if(strncasecmp(cptr, "NO", 1)==0) {
        header=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-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);}
 
  /* In verbose mode print options */
  if(verbose>1) {
    printf("mode := %s\n", mode_str[mode]);
    if(xfile[0]) printf("xfile := %s\n", xfile);
    if(!isnan(start0)) printf("start0 := %g\n", start0);
    if(!isnan(startY)) printf("startY := %g\n", startY);
    if(xnr>0) {
      printf("xlist :="); 
      for(int i=0; i<xnr; i++) printf(" %g", xlist[i]); 
      printf("\n");
    }
  }
  
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) strlcpy(tacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(outfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!tacfile[0]) {
    fprintf(stderr, "Error: missing file name.\n");
    return(1);
  }
 
  /* In verbose mode print arguments */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    if(outfile[0]) printf("outfile := %s\n", outfile);
  }
 
 
  /*
   *  Read the TAC file
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  ret=tacRead(&tac1, tacfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac1); free(xlist); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac1.format));
    printf("tacNr := %d\n", tac1.tacNr);
    printf("sampleNr := %d\n", tac1.sampleNr);
    printf("xunit := %s\n", unitName(tac1.tunit));
    printf("yunit := %s\n", unitName(tac1.cunit));
  }
  /* Sort data by sample times */
  ret=tacSortByTime(&tac1, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac1); free(xlist); return(2);
  }
  /* Check for missing y values */
  ret=tacNaNs(&tac1);
  if(ret>0) {
    if(verbose>1) printf("missing concentrations.\n");
    /* Try to fix missing concentrations */
    ret=tacFixNaNs(&tac1);
    if(ret!=0) {
      fprintf(stderr, "Error: missing concentrations in %s.\n", tacfile);
      tacFree(&tac1); free(xlist); return(2);
    }
  }
 
  /* 
   *  Add an initial value, if required with option -A
   */
  if(isnormal(start0)) {
    if(verbose>1) printf("adding initial point\n");
    /* Make sure that sample is added, but keep the original first sample, too */
    double x, x1, x2;
    x=tac1.x[0]; x1=tac1.x1[0]; x2=tac1.x2[0];
    tac1.x[0]=tac1.x1[0]=tac1.x2[0]=1.0;
    ret=tacAddZeroSample(&tac1, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      if(verbose>=0) printf("Note: check the data and option -A\n");
      tacFree(&tac1); free(xlist); return(3);
    }
    tac1.x[1]=x; tac1.x1[1]=x1; tac1.x2[1]=x2;
    /* Set the new first sample to required values */
    tac1.x[0]=tac1.x1[0]=tac1.x2[0]=start0;
    for(int j=0; j<tac1.tacNr; j++) tac1.c[j].y[0]=0.0;
    /* Remove samples with earlier time than the new first sample */
    if(tac1.isframe) 
      for(int i=1; i<tac1.sampleNr; i++) tac1.x[i]=0.5*(tac1.x1[i]+tac1.x2[i]);
    ret=TPCERROR_OK;
    while(tac1.sampleNr>1 && tac1.x[1]<start0 && ret==TPCERROR_OK) {
      if(verbose>4) printf("deleting sample x=%g\n", tac1.x[1]);
      ret=tacDeleteSample(&tac1, 1);
    }
    if(ret!=TPCERROR_OK || tac1.sampleNr<2) {
      fprintf(stderr, "Error: invalid data or option -A.\n");
      tacFree(&tac1); free(xlist); return(3);
    }
  }
 
 
  /* 
   *  Add an initial value, if required with option -0 
   */
  if(tac1.isframe) tac1.x[0]=0.5*(tac1.x1[0]+tac1.x2[0]);
  if(isnormal(startY) && tac1.x[0]>0.0) {
    if(verbose>1) printf("adding initial value at zero\n");
    ret=tacAddZeroSample(&tac1, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      if(verbose>=0) printf("Note: check the data and option -0\n");
      tacFree(&tac1); free(xlist); return(4);
    }
    /* Set first sample to zero time */
    tac1.x[0]=tac1.x1[0]=tac1.x2[0]=0.0;
    /* Set the value */
    if(startY<-1.0E+10 && tac1.sampleNr>1)
      for(int j=0; j<tac1.tacNr; j++) tac1.c[j].y[0]=tac1.c[j].y[1];
    else
      for(int j=0; j<tac1.tacNr; j++) tac1.c[j].y[0]=startY;
  }
 
 
  /*
   *  Set up space for output data and set output times
   */
  if(xfile[0]) {
  
    /* If sample times are given in a file, then read it,
       and use it as the basis of output data */
    if(verbose>1) printf("reading %s\n", xfile);
    ret=tacRead(&tac2, xfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac1); tacFree(&tac2); free(xlist); return(5);
    }
    if(verbose>2) {
      printf("fileformat2 := %s\n", tacFormattxt(tac2.format));
      printf("tacNr2 := %d\n", tac2.tacNr);
      printf("sampleNr2 := %d\n", tac2.sampleNr);
      printf("xunit2 := %s\n", unitName(tac2.tunit));
      printf("yunit2 := %s\n", unitName(tac2.cunit));
    }
    /* Convert time units, if necessary/possible */
    ret=tacXUnitConvert(&tac2, tac1.tunit, &status);
    if(ret!=TPCERROR_OK && verbose>0) {
      fprintf(stderr, "Warning: different or unknown time units.\n");
    }
    /* Allocate space for all TACs */
    if(tac2.tacNr<tac1.tacNr) {
      tac2.tacNr=0;
      ret=tacAllocateMore(&tac2, tac1.tacNr);
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(ret));
        tacFree(&tac1); tacFree(&tac2); free(xlist); return(5);
      }
    }
    tac2.tacNr=tac1.tacNr;
    /* Set header */
    tac2.weighting=WEIGHTING_OFF;
    tac2.cunit=tac1.cunit;
    tac2.format=tac1.format;
    iftFree(&tac2.h); iftDuplicate(&tac1.h, &tac2.h);
    for(int j=0; j<tac1.tacNr; j++) tacCopyTacchdr(&tac1.c[j], &tac2.c[j]);
    
  } else if(xnr>0) {
  
    /* x values were given with command-line options -X or -C */
    if(verbose>1) printf("allocating memory for %d samples\n", xnr);
    ret=tacAllocate(&tac2, xnr, tac1.tacNr);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(ret));
      tacFree(&tac1); tacFree(&tac2); free(xlist); return(5);
    }
    tac2.sampleNr=xnr; tac2.tacNr=tac1.tacNr;
    tac2.weighting=WEIGHTING_OFF;
    tac2.cunit=tac1.cunit;
    tac2.tunit=tac1.tunit;
    tac2.format=tac1.format;
    iftFree(&tac2.h); iftDuplicate(&tac1.h, &tac2.h);
    for(int j=0; j<tac1.tacNr; j++) tacCopyTacchdr(&tac1.c[j], &tac2.c[j]);
    tac2.isframe=0;
    for(int j=0; j<xnr; j++) tac2.x[j]=xlist[j];
 
  } else {
  
    /* x values must be retrieved from the input data */
    ret=tacDuplicate(&tac1, &tac2);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(ret));
      tacFree(&tac1); tacFree(&tac2); free(xlist); return(5);
    }
 
  }
  /* list of x values is not needed later */
  free(xlist);
 
  if(verbose>3) {
    printf("_tacNr2 := %d\n", tac2._tacNr);
    printf("_sampleNr2 := %d\n", tac2._sampleNr);
    printf("tacNr2 := %d\n", tac2.tacNr);
    printf("sampleNr2 := %d\n", tac2.sampleNr);
  }
 
 
  /*
   *  Calculate
   */
  ret=0;
  if(mode==0) {
    if(verbose>1) printf("interpolating\n");
    if(tac2.isframe) for(int j=0; j<tac1.tacNr && ret==0; j++) {
      ret=liInterpolateForPET(tac1.x, tac1.c[j].y, tac1.sampleNr,
            tac2.x1, tac2.x2, tac2.c[j].y, NULL, NULL, tac2.sampleNr, 
            3, 1, verbose-2);
    } else for(int j=0; j<tac1.tacNr && ret==0; j++) {
      ret=liInterpolate(tac1.x, tac1.c[j].y, tac1.sampleNr,
            tac2.x, tac2.c[j].y, NULL, NULL, tac2.sampleNr, 
            3, 1, verbose-2);
    }
  } else if(mode==1) {
    if(verbose>1) printf("integrating\n");
    if(tac2.isframe) for(int j=0; j<tac1.tacNr && ret==0; j++) {
      ret=liInterpolateForPET(tac1.x, tac1.c[j].y, tac1.sampleNr,
            tac2.x1, tac2.x2, NULL, tac2.c[j].y, NULL, tac2.sampleNr, 
            3, 1, verbose-2);
    } else for(int j=0; j<tac1.tacNr && ret==0; j++) {
      ret=liInterpolate(tac1.x, tac1.c[j].y, tac1.sampleNr,
            tac2.x, NULL, tac2.c[j].y, NULL, tac2.sampleNr, 
            3, 1, verbose-2);
    }
  } else if(mode==2) {
    if(verbose>1) printf("computing 2nd integrals\n");
    if(tac2.isframe) for(int j=0; j<tac1.tacNr && ret==0; j++) {
      ret=liInterpolateForPET(tac1.x, tac1.c[j].y, tac1.sampleNr,
            tac2.x1, tac2.x2, NULL, NULL, tac2.c[j].y, tac2.sampleNr, 
            3, 1, verbose-2);
    } else for(int j=0; j<tac1.tacNr && ret==0; j++) {
      ret=liInterpolate(tac1.x, tac1.c[j].y, tac1.sampleNr,
            tac2.x, NULL, NULL, tac2.c[j].y, tac2.sampleNr, 
            3, 1, verbose-2);
    }
  }
  if(ret) {
    fprintf(stderr, "Error: cannot interpolate/integrate.\n");
    if(verbose>1) printf("return_code := %d\n", ret);
    tacFree(&tac1); tacFree(&tac2); return(8);
  }
 
 
  /*
   *  Save data 
   */
  if(tac2.isframe && tac2.format==TAC_FORMAT_SIMPLE) {
    /* make sure that frame start and end times are not replaced by middle
       times because of file format */
    tac2.format=TAC_FORMAT_DFT;
  }
  FILE *fp; 
  if(outfile[0]) {
    if(verbose>1) printf("writing %s\n", outfile);
    fp=fopen(outfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", outfile);
      tacFree(&tac1); tacFree(&tac2); return(11);
    }
  } else if(header==0 && tac2.tacNr==1 && tac2.sampleNr==1) {
    /* If user does not want headers, and we have just one value to print to
       stdout, then do not print time(s) but only the value */
    fprintf(stdout, "%g\n", tac2.c[0].y[0]);
    tacFree(&tac1); tacFree(&tac2);
    return(0);
  } else {
    fp=stdout;
  }
  ret=tacWrite(&tac2, fp, TAC_FORMAT_UNKNOWN, header, &status);
  tacFree(&tac1); tacFree(&tac2);
  if(outfile[0]) fclose(fp);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    return(12);
  }
  if(verbose>=0 && outfile[0]) printf("%s saved.\n", outfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/