fitdt.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"
#include "tpctacmod.h"
#include "tpcbfm.h"
#include "tpcstatist.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Estimate the time delay (dT) between PET tissue and blood curves.",
  "Initial part of tissue curve is fitted with a range of input curves moved",
  "in time (1, 2, 3), using linearised one or two tissue compartment model",
  "with blood volume term. The time that provides the best fit",
  "(smallest sum-of-squares) is selected. Tissue heterogeneity, dispersion of",
  "the blood curve, metabolites, and variable plasma-to-blood ratio may",
  "only minimally impact the time delay estimate, if only the initial part",
  "of the data is fitted with relatively complex model.",
  " ",
  "Usage: @P [options] btacfile ttacfile parfile [dcttacfile]",
  " ",
  "Options:",
  " -min=<Time (sec)> and -max=<Time (sec)>",
  "     The range of time delays to be tested; by default -10 - +50 s.",
  " -end=<Fit end time (sec)>",
  "     Use data from 0 to end time; by default, 300 s; at least 120 s.",
  "     In case of short PET scan, end time may need to be reduced so that",
  "     BTAC with negative delay extends to end time.",
  " -step=<dT step size (sec)>",
  "     The step length for moving the input curve; by default 1 s.",
//" -model=<R1TCM | I2TCM | R2TCM | DMS>", // -DMS currently does not work reliably.
  " -model=<R1TCM | I2TCM | R2TCM>",
  "     Reversible one-tissue (R1TCM, default), irreversible 2-tissue model",
  "     (I2TCM), or reversible 2-tissue model (R2TCM) can be selected.", 
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Delay times for each tissue curve are reported in parameter file.",
  "File may also contain a median time delay as 'time_difference := time',",
  "depending on the file format.",
  "Positive delay time means that tissue curve is delayed as compared to",
  "the input curve, and vice versa. Thus, input curve needs to be moved",
  "by the delay time to match the tissue curve.",
  " ",
  "Optionally, delay corrected TTACs can be written in specified file.",
  " ",
  "References:",
  " 1. Iida H et al. Evaluation of regional differences of tracer appearance", 
  "    time in cerebral tissues using [15O]water and dynamic positron emission",
  "    tomography. J Cereb Blood Flow Metab. 1988; 8:285-288.",
  " 2. Meyer. Simultaneous correction for tracer arrival delay and dispersion",
  "    in CBF measurements by the H215O autoradiographic method and dynamic PET.",
  "    J Nucl Med 1989; 30:1069-1078.",
  " 3. van den Hoff et al. Accurate local blood flow measurements with",
  "    dynamic PET: fast determination of input function delay and dispersion",
  "    by multilinear minimization. J Nucl Med 1993; 34:1770-1777.",
  " ",
  "See also: fitdelay, imgdelay, tacmove, tactime, simdisp, tacframe",
  " ",
  "Keywords: TAC, modelling, time delay",
  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     ttacfile[FILENAME_MAX], btacfile[FILENAME_MAX],
           parfile[FILENAME_MAX], cttacfile[FILENAME_MAX];
  double   endtime=300.0; // fit end time in seconds
  double   endtimemin=90.0; // shortest fit end time allowed
  double   dtrange[2]={-10,+50}; // dT range
  double   dtstep=1.0; // dT step size
  int      mode=0; // Delay estimation mode
  int      model=0; // 0=R1TCM, 1=I2TCM, 2=R2TCM, 3=DMS
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ttacfile[0]=btacfile[0]=parfile[0]=cttacfile[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(strncasecmp(cptr, "END=", 4)==0) {
      int ret=atofCheck(cptr+4, &endtime); if(!ret && endtime>=endtimemin) continue;
    } else if(strncasecmp(cptr, "MIN=", 4)==0) {
      if(atofCheck(cptr+4, &dtrange[0])==0) continue;
    } else if(strncasecmp(cptr, "MAX=", 4)==0) {
      if(atofCheck(cptr+4, &dtrange[1])==0) continue;
    } else if(strncasecmp(cptr, "STEP=", 5)==0) {
      int ret=atofCheck(cptr+5, &dtstep); if(!ret && dtstep>0.0) continue;
    } else if(strncasecmp(cptr, "MODE=", 5)==0) {
      if(atoiCheck(cptr+5, &mode)==0) continue;
    } else if(strncasecmp(cptr, "MODEL=", 6)==0) {
      if(strcasecmp(cptr+6, "R1TCM")==0) {model=0; continue;}
      if(strcasecmp(cptr+6, "I2TCM")==0) {model=1; continue;}
      if(strcasecmp(cptr+6, "R2TCM")==0) {model=2; continue;}
      if(strcasecmp(cptr+6, "DMS")==0) {model=3; 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(btacfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(ttacfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(parfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(cttacfile, argv[ai], FILENAME_MAX); ai++;}
  /* check that there are no extra arguments */
  if(ai<argc) {fprintf(stderr, "Error: extra command-line argument.\n"); return(1);}
  /* Is something missing? */
  if(!parfile[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ttacfile := %s\n", ttacfile);
    printf("btacfile := %s\n", btacfile);
    printf("parfile := %s\n", parfile);
    printf("cttacfile := %s\n", cttacfile);
    printf("endtime := %g\n", endtime);
    printf("dtrange := %g - %g s\n", dtrange[0], dtrange[1]);
    printf("step := %g\n", dtstep);
    printf("mode := %d\n", mode);
    printf("model := %d\n", model);
    fflush(stdout);
  }
 
  /* Check that TAC files do exist */
  if(!fileExist(ttacfile)) {
    fprintf(stderr, "Error: file '%s' does not exist.\n", ttacfile);
    return(1);
  }
  if(!fileExist(btacfile)) {
    fprintf(stderr, "Error: file '%s' does not exist.\n", btacfile);
    return(1);
  }
 
  /* Check options */
  if(dtrange[0]>dtrange[1]) {double f=dtrange[0]; dtrange[0]=dtrange[1]; dtrange[1]=f;}
  int moveNr=1+(dtrange[1]-dtrange[0])/dtstep;
  if(verbose>1) {
    printf("step_number := %d\n", moveNr);
    fflush(stdout);
  }
  if(moveNr<5 || moveNr>2000) {
    fprintf(stderr, "Error: invalid dT range or step size.\n");
    return(1);
  } 
  if((moveNr<10 || moveNr>500) && verbose>0)
    fprintf(stderr, "Warning: non-optimal dT range or step size.\n");
 
 
 
  /*
   *  Read the data
   */
  if(verbose>1) printf("reading TACs\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC ttac, btac; tacInit(&ttac); tacInit(&btac);
 
  if(tacRead(&ttac, ttacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  if(verbose>2) {
    printf("ttac.fileformat := %s\n", tacFormattxt(ttac.format));
    printf("ttacNr := %d\n", ttac.tacNr);
    printf("ttac.sampleNr := %d\n", ttac.sampleNr);
    fflush(stdout);
  }
 
  if(tacRead(&btac, btacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  if(verbose>2) {
    printf("btac.fileformat := %s\n", tacFormattxt(btac.format));
    printf("btacNr := %d\n", btac.tacNr);
    printf("btac.sampleNr := %d\n", btac.sampleNr);
    fflush(stdout);
  }
  if(btac.tacNr>1) {
    if(verbose>0) fprintf(stderr, "Warning: BTAC file contains more than one TAC.\n");
    btac.tacNr=1;
  }
 
  if(ttac.sampleNr<7 || btac.sampleNr<7) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  /* Check NaNs */
  if(tacNaNs(&ttac)>0 || tacNaNs(&btac)>0) {
    fprintf(stderr, "Error: data contains missing values.\n");
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  /* Sort data by sample time */
  tacSortByTime(&ttac, &status);
  tacSortByTime(&btac, &status);
  /* Fix any gaps and overlap in data if possible and if not then quit with error */
  if(tacSetXContiguous(&ttac)!=TPCERROR_OK || tacSetXContiguous(&btac)!=TPCERROR_OK) {
    fprintf(stderr, "Error: invalid curve data.\n");
    tacFree(&ttac); tacFree(&btac); return(2);
  }
 
  /* Convert sample times into seconds */
  if(tacXUnitConvert(&ttac, UNIT_SEC, &status)!=TPCERROR_OK ||
     tacXUnitConvert(&btac, UNIT_SEC, &status)!=TPCERROR_OK) 
  {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&ttac); return(2);
  }
  /* Convert BTAC concentrations into TTAC units */
  if(tacYUnitConvert(&btac, ttac.cunit, &status)!=TPCERROR_OK) {
    if(verbose>0) fprintf(stderr, "Warning: check and set the data units.\n");
  }
 
 
  /*
   *  Prepare the room for results
   */
  if(verbose>2) printf("allocate memory for PAR\n");
  PAR par; parInit(&par);
  if(parAllocate(&par, 1, ttac.tacNr)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot allocate memory for results.\n");
    tacFree(&ttac); tacFree(&btac); return(3);
  }
  par.tacNr=ttac.tacNr;
  par.parNr=1;
  /* Copy TAC names */
  for(int i=0; i<par.tacNr; i++) strcpy(par.r[i].name, ttac.c[i].name);
  /* Set the parameter name and unit */
  strcpy(par.n[0].name, "dT");
  par.n[0].unit=UNIT_SEC;
  /* Set titles & file names */
  {
    char buf[256];
    time_t t=time(NULL);
    tpcProgramName(argv[0], 1, 1, buf, 256);
    iftPut(&par.h, "program", buf, 0, NULL);
    iftPut(&par.h, "bloodfile", btacfile, 0, NULL);
    iftPut(&par.h, "datafile", ttacfile, 0, NULL);
    iftPut(&par.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
    for(int i=0; i<par.tacNr; i++) {
      par.r[i].start=0.0;
      par.r[i].end=endtime;
    }
  }
 
 
  /*
   *  In case of function fits (instead of compartmental models), fit first the delay of BTAC
   */
  double btacDelaytime=0.0;
  if(model>=3) {
    if(verbose>1) {printf("BTAC delay fitting\n"); fflush(stdout);}
    int ret=0, fitSampleNr=0;
    for(fitSampleNr=0; fitSampleNr<btac.sampleNr; fitSampleNr++) if(btac.x[fitSampleNr]>endtime) break;
    double *x1, *x2;
    if(btac.isframe) {x1=btac.x1; x2=btac.x2;} else {x1=btac.x; x2=NULL;}
    //printf("fitSampleNr=%d lastx=%g endtime=%g\n", fitSampleNr, btac.x[fitSampleNr-1], endtime);
    ret=spectralDMSurge(x1, x2, btac.c[0].y, NULL, fitSampleNr, 1.0E-10, 1.0E+00, 1000, 
            dtrange[0], dtrange[1], dtstep, NULL, NULL, &btacDelaytime, NULL, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); parFree(&par);
      return(4);
    }
    if(verbose>3) printf("BTAC delay := %g\n", btacDelaytime);
  }
 
 
  /* 
   *  Estimate the delay time
   */
  if(verbose>1) {printf("delay fitting\n"); fflush(stdout);}
  if(model<3) { // CM fits
    if(ttac.tacNr==1) { // just one tissue TAC, no need for temp data or median delay
      double delaytime=nan("");
      int ret=tacDelayCMFit(&btac, &ttac, 0, dtrange[0], dtrange[1], endtime, dtstep, &delaytime, 
                          mode, model, NULL, &status);
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&ttac); tacFree(&btac); parFree(&par);
        return(4);
      }
      if(verbose>3) printf("delay := %g\n", delaytime);
      par.r[0].p[0]=delaytime;
    } else { // more than one tissue TAC
      DELAYCMFITDATA ddata;
      initDelayCMFitData(&ddata);
      int ret=0;
      for(int ci=0; ci<ttac.tacNr; ci++) {
        if(verbose>2) {printf("%s\n", ttac.c[ci].name); fflush(stdout);}
        double delaytime=nan("");
        if(ci==0) 
          ret=tacDelayCMFit(&btac, &ttac, ci, dtrange[0], dtrange[1], endtime, dtstep, &delaytime, 
                          mode, model, &ddata, &status);
        else
          ret=tacDelayCMFit(NULL, &ttac, ci, 0.0, 0.0, 0.0, 0.0, &delaytime, mode, model, &ddata, &status);
        if(ret!=TPCERROR_OK) {
          fprintf(stderr, "Error: %s\n", errorMsg(status.error));
          freeDelayCMFitData(&ddata);
          tacFree(&ttac); tacFree(&btac); parFree(&par);
          return(4);
        }
        if(verbose>3) printf("delay := %g\n", delaytime);
        par.r[ci].p[0]=delaytime;
      }
      freeDelayCMFitData(&ddata);
    }
  } else { // Function fits
    int fitSampleNr=0;
    for(fitSampleNr=0; fitSampleNr<ttac.sampleNr; fitSampleNr++) if(ttac.x[fitSampleNr]>endtime) break;
    double *x1, *x2;
    if(ttac.isframe) {x1=ttac.x1; x2=ttac.x2;} else {x1=ttac.x; x2=NULL;}
    //printf("fitSampleNr=%d lastx=%g endtime=%g\n", fitSampleNr, ttac.x[fitSampleNr-1], endtime);
    for(int ci=0; ci<ttac.tacNr; ci++) {
      if(verbose>2) {printf("%s\n", ttac.c[ci].name); fflush(stdout);}
      double delaytime=nan("");
      int ret=spectralDMSurge(x1, x2, ttac.c[ci].y, NULL, fitSampleNr, 1.0E-10, 1.0E+00, 1000, 
                      dtrange[0], dtrange[1], dtstep, NULL, NULL, &delaytime, NULL, &status);
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&ttac); tacFree(&btac); parFree(&par);
        return(4);
      }
      if(verbose>3) printf("delay := %g (%g-%g)\n", delaytime-btacDelaytime, delaytime, btacDelaytime);
      delaytime-=btacDelaytime;
      par.r[ci].p[0]=delaytime;
    }
  }
 
  /* Set median delay in results */
  {
    double delaytime=nan("");
    double t[par.tacNr];
    int n=0;
    for(int i=0; i<par.tacNr; i++) if(isfinite(par.r[i].p[0])) t[n++]=par.r[i].p[0];
    if(n==0) {
      fprintf(stderr, "Error: delay estimation failed.\n");
      tacFree(&ttac); tacFree(&btac); parFree(&par);
      return(5);
    } else {
      delaytime=statMedian(t, n);
    }
    iftPutDouble(&par.h, "time_difference", delaytime, 1, NULL);
  }
 
 
  /* Print and save the parameters */
  if(verbose>0) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 0, NULL);
  if(parfile[0]) {
    par.format=parFormatFromExtension(parfile);
    if(verbose>2) printf("parameter file format := %s\n", parFormattxt(par.format));
    if(par.format==PAR_FORMAT_UNKNOWN) par.format=PAR_FORMAT_TSV_UK;
    /* Save file */
    if(verbose>1) printf("  saving %s\n", parfile);
    FILE *fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&ttac); tacFree(&btac); parFree(&par); return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); parFree(&par); return(12);
    }
    if(verbose>1) printf("parameters saved in %s\n", parfile);
  }
 
 
  /*
   *  Optionally, delay correction for the TTACs
   */
  if(cttacfile[0]) {
    if(verbose>1) {printf("delay correcting TTACs\n"); fflush(stdout);}
    int origframe=ttac.isframe; ttac.isframe=0; // simple interpolation
    for(int i=0; i<ttac.tacNr; i++) {
      if(tacDelay(&ttac, -par.r[i].p[0], i, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&ttac); tacFree(&btac); parFree(&par); return(21);
      }
    }
    ttac.isframe=origframe;
    if(verbose>1) {printf("writing delay correcting TTACs\n"); fflush(stdout);}
    FILE *fp=fopen(cttacfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&ttac); tacFree(&btac); parFree(&par); return(22);
    }
    int ret=tacWrite(&ttac, fp, TAC_FORMAT_UNKNOWN, 0, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); parFree(&par); return(23);
    }
    if(verbose>1) printf("delay corrected TTACs saved in %s\n", cttacfile);
  }
 
  tacFree(&ttac); tacFree(&btac); parFree(&par);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/