dftinteg.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
#include "libtpcmodel.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculation of an AUC (integral over time) from regional or plasma TACs.",
  " ",
  "Usage: @P [Options] tacfile starttime integrationtime aucfile",
  " ",
  "Enter the start time and integration time in the same units that are used",
  "in the TAC file. Note that 3rd argument is not the integration end time.",
  "Program will automatically set the integration start time and duration",
  "based on the time range in the TAC file, if both are set to zero.", 
  "AUC/mean is by default written in DFT format, or in RES format if filename",
  "has extension .res.",
  " ",
  "Options:",
  " -avg",
  "     Average during specified range is calculated instead of AUC.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1: calculate AUC(32-302) with command",
  "  @P s5998dy1.dft 32 270 s5998int.dft",
  "Example 2: integrate the regional TACs from a static or dynamic study",
  "from the start time of the first frame to the end of the last frame:",
  "  @P a773dy1.dft 0 0 a773int.dft",
  "Example 3: mean between 0 and 10 is calculated and saved in result format:",
  "  @P -avg ec4568.dft 0 10 ec4568_mean0-10.res", 
  " ",
  "See also: interpol, dftratio, dftsuv, regfur, tacunit, imginteg, tac2suv",
  " ",
  "Keywords: TAC, modelling, AUC, autoradiography",
  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;
  char  *cptr, dftfile[FILENAME_MAX], outfile[FILENAME_MAX], tmp[512];
  int    calc_avg=0; // 0=integral, 1=average
  DFT    dft, idft;
  double tstart, tstop, tdur;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dftfile[0]=outfile[0]=(char)0;
  tstart=tdur=tstop=-1.0;
  dftInit(&dft); dftInit(&idft);
  /* 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(strcasecmp(cptr, "AVG")==0) {
      calc_avg=1; continue;
    } else if(strcasecmp(cptr, "MEAN")==0) {
      calc_avg=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 */
  for(; ai<argc; ai++) {
    if(!dftfile[0]) {
      strlcpy(dftfile, argv[ai], FILENAME_MAX); continue;
    } else if(tstart<0) {
      if(atof_with_check(argv[ai], &tstart)==0) continue;
    } else if(tdur<0) {
      if(atof_with_check(argv[ai], &tdur)==0) {tstop=tstart+tdur; continue;}
    } else if(!outfile[0]) {
      strlcpy(outfile, argv[ai], FILENAME_MAX); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!outfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dftfile := %s\n", dftfile);
    printf("outfile := %s\n", outfile);
    printf("calc_avg := %d\n", calc_avg);
    if(tstart>=1.0E-03 || tdur>=1.0E-02) printf("tstart := %g\ntstop := %g\n", tstart, tstop);
  }
 
 
  /*
   *  Read the datafile
   */
  if(dftRead(dftfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dftfile, dfterrmsg);
    return(2);
  }
  /* Check and remove NaNs if possible */
  if(dftNAfill(&dft)!=0) {
    fprintf(stderr, "Error: missing values in %s.\n", dftfile);
    dftEmpty(&dft); return(2);
  }
  /* Make sure that there is no overlap in image frames */
  if(dft.timetype==DFT_TIME_STARTEND) {
    if(verbose>1) fprintf(stdout, "checking frame overlap in %s\n", dftfile);
    ret=dftDeleteFrameOverlap(&dft);
    if(ret) {
      fprintf(stderr, "Error: %s has overlapping frame times.\n", dftfile);
      dftEmpty(&dft); return(2);
    }
  }
  /* If both integration start time and duration were set to zero, then take
     these from the TAC data */
  if(tstart<1.0E-03 && tdur<1.0E-02) {
    if(dft.timetype==DFT_TIME_STARTEND) {
      tstart=dft.x1[0]; tstop=dft.x2[dft.frameNr-1];
    } else {
      tstart=dft.x[0]; tstop=dft.x[dft.frameNr-1];
    }
    tdur=tstop-tstart;
    if(verbose>1) {
      printf("tstart := %g\n", tstart);
      printf("tstop := %g\n", tstop);
      printf("tdur := %g\n", tdur);
    }
  }
  if(verbose>10) dftPrint(&dft);
 
  /*
   *  Calculate the AUC
   */
  ret=dftTimeIntegral(&dft, tstart, tstop, &idft, calc_avg, tmp, verbose-2);
  if(ret!=0) {
    fprintf(stderr, "Error: %s\n", tmp);
    dftEmpty(&dft); return(5);
  }
  dftEmpty(&dft);
 
  /*
   *  Save the AUC file, either as DFT or as RES
   */
  if(verbose>2) {
    if(calc_avg==0) printf("writing AUC in %s\n", outfile);
    else printf("writing mean in %s\n", outfile);
  }
  cptr=strrchr(outfile, '.'); if(cptr!=NULL) cptr++;
  if(strcasecmp(cptr, "RES")==0) { // Save in RES format
    RES res;
    resInit(&res);
    ret=dftToResult(&idft, &res, tmp);
    if(ret!=0) {
      fprintf(stderr, "Error in making results: %s\n", tmp);
      dftEmpty(&idft); resEmpty(&res); return(11);
    }
    if(calc_avg==0) {strcpy(res.parname[0], "AUC");}
    else {strcpy(res.parname[0], "Mean"); strcpy(res.parunit[0], idft.unit);}
    tpcProgramName(argv[0], 1, 1, res.program, 256);
    strcpy(res.datafile, dftfile);
    sprintf(res.datarange, "%g - %g", idft.x1[0], idft.x2[0]);
    res.datanr=0;
    if(verbose>0) resPrint(&res);
    if(resWrite(&res, outfile, verbose-3)) {
      fprintf(stderr, "Error in writing '%s': %s\n", outfile, reserrmsg);
      dftEmpty(&idft); resEmpty(&res); return(11);
    }
    resEmpty(&res);
  } else { // save in DFT format
    dftSetComments(&idft);
    if(dftWrite(&idft, outfile)) {
      fprintf(stderr, "Error in writing %s: %s\n", outfile, dfterrmsg);
      dftEmpty(&idft); dftEmpty(&idft); return(11);
    }
  }
  if(verbose>0) printf("%s written.\n", outfile);
 
  dftEmpty(&idft);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/