extrapol.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcsvg.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Extrapolation of the decreasing tail of PET input curve to given time.",
  "This is accomplished by fitting line to the end-part of the plot of the",
  "natural logarithm of tracer concentration against time.",
  "This approach may be applied to TACs which approach zero; thus this",
  "works best for metabolite-corrected input TACs, and cannot be used with",
  "tracers with very slow clearance, such as radiowater.", 
  " ",
  "Usage: @P [Options] inputfile time outputfile",
  " ",
  "Options:",
  " -e[nd]=<Fit end time>",
  "     Start the search to fit end time; by default, the search for the best",
  "     line fit is started from the last sample.",
  " -minnr=<Minimum nr of samples>",
  "     Set the minimum number of samples used in searching the best fit;",
  "     by default 3.",
  " -maxnr=<Maximum nr of samples>",
  "     Set the maximum number of samples used in searching the best fit;",
  "     by default all samples.",
  " -mintime=<Minimum time>",
  "     Set a minimum time range used in searching the best fit.",
  " -svg=<Filename>",
  "     Measured and extrapolated TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Input TAC file must contain a time column, and one or more concentration",
  "columns. The extrapolation time must be given in same units as are the sample",
  "times in the datafile.",
  " ",
  "See also: paucinf, avgbolus, fit_dexp, fit_feng, fit2dat, tacln",
  " ",
  "Keywords: input, plasma, TAC, modelling, simulation",
  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;
  char inpfile[FILENAME_MAX], outfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  char temp[FILENAME_MAX], *cptr;
  DFT dft, ext;
  double extr_to=-1.0, fittime=-1.0, mintime=-1.0;
  int ri, n, ret, min_nr=3, max_nr=-1;
  FILE *logfp;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=outfile[0]=svgfile[0]=(char)0;
  dftInit(&dft); dftInit(&ext);
 
  /* 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;
    if(strncasecmp(cptr, "E=", 2)==0) {
      fittime=atof_dpi(cptr+2); if(fittime>0.0) continue;
    } else if(strncasecmp(cptr, "END=", 4)==0) {
      fittime=atof_dpi(cptr+4); if(fittime>0.0) continue;
    } else if(strncasecmp(cptr, "MINNR=", 6)==0) {
      min_nr=atoi(cptr+6); if(min_nr>1) continue;
    } else if(strncasecmp(cptr, "MAXNR=", 6)==0) {
      max_nr=atoi(cptr+6); if(max_nr>1) continue;
    } else if(strncasecmp(cptr, "MINTIME=", 8)==0) {
      mintime=atof_dpi(cptr+8); if(mintime>0.0) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0) {
      strcpy(svgfile, cptr+4); if(strlen(svgfile)>0) 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 */
  for(; ai<argc; ai++) {
    if(!inpfile[0]) {
      strcpy(inpfile, argv[ai]); continue;
    } else if(extr_to<0.0) {
      extr_to=atof_dpi(argv[ai]); if(extr_to>0.0) continue;
      fprintf(stderr, "Error: invalid extrapolation time: '%s'.\n", argv[ai]);
      return(1);
    } else if(!outfile[0]) {
      strcpy(outfile, argv[ai]); continue;
    }
    /* we should never get this far */
    fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!outfile[0]) {
    fprintf(stderr, "Error: missing file name for extrapolated data.\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("inpfile := %s\n", inpfile);
    printf("outfile := %s\n", outfile);
    printf("extr_to := %g\n", extr_to);
    printf("svgfile := %s\n", svgfile);
    printf("fittime := %g\n", fittime);
    printf("mintime := %g\n", mintime);
    printf("min_nr := %d\n", min_nr);
    printf("max_nr := %d\n", max_nr);
  }
 
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading %s\n", inpfile);
  if(dftRead(inpfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", inpfile, dfterrmsg);
    return(2);
  }
 
 
  /*
   *  Delete the ROIs that have name 'SD' (produced at least by avgbolus)
   */
  for(ri=0, n=0; ri<dft.voiNr; ri++) {
    dft.voi[ri].sw=0;
    if(strcasecmp(dft.voi[ri].voiname, "SD")==0) {
      dft.voi[ri].sw=1; n++; continue;}
    if(strcasecmp(dft.voi[ri].hemisphere, "SD")==0) {
      dft.voi[ri].sw=1; n++; continue;}
  }
  if(n>0 && n<dft.voiNr) {
    ri=dft.voiNr-1;
    while(ri>=0) {
      if(dft.voi[ri].sw) {
        if(verbose>0) 
          printf("  TAC '%s' is not extrapolated\n", dft.voi[ri].name);
        dftDelete(&dft, ri);
      }
      ri--;
    }
  }
  if(dft.voiNr<1) {
    fprintf(stderr, "Error: file does not contain data for extrapolation.\n");
    dftEmpty(&dft); return(2);
  }
 
 
  /*
   *  Extrapolation
   */
  if(verbose>1) printf("extrapolation\n");
  if(verbose>2) logfp=stdout; else logfp=NULL;
  ret=extrapolate_monoexp(
    &dft, &fittime, &min_nr, max_nr, mintime, extr_to, &ext, logfp, temp
  );
  if(ret!=0) {
    fprintf(stderr, "Error in extrapolation: %s\n", temp);
    dftEmpty(&dft); dftEmpty(&ext);
  }
 
  /*
   *  Writing extrapolated TAC
   */
  if(verbose>1) printf("writing %s\n", outfile);
  /* Set comments */
  dftSetComments(&ext);
  snprintf(temp, FILENAME_MAX-1, "# extrapolated: %s %g\n", inpfile, extr_to);
  strcat(ext.comments, temp);
  /* Write the file */
  if(ext._type<1) ext._type=DFT_FORMAT_STANDARD;
  ret=dftWrite(&ext, outfile);
  if(ret) {
    fprintf(stderr, "Error (%d) in writing '%s': %s\n", ret, outfile, dfterrmsg);
    dftEmpty(&dft); dftEmpty(&ext); return(11);
  }
 
  /*
   *  Save SVG plot of original and extrapolated TAC
   */
  if(svgfile[0]) {
    if(verbose>1) printf("saving SVG plot\n");
    sprintf(temp, "Extrapolation ");
    if(strlen(dft.studynr)>1) strcat(temp, dft.studynr);
    ret=plot_fit_svg(&dft, &ext, temp, svgfile, verbose-2);
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
      dftEmpty(&dft); dftEmpty(&ext);
      return(30+ret);
    }
    if(verbose>=0) printf("Plots written in %s\n", svgfile);
  }
 
  /* Free memory */
  dftEmpty(&dft); dftEmpty(&ext);
 
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/