regfur.c

Go to the documentation of this file.

/******************************************************************************/
#include "tpcclibConfig.h"
/******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
/******************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculates Fractional Uptake Rate (FUR) or FUR-based Metabolic Rate (MR)",
  "from regional PET TTACs. Information on FUR in:",
  "https://www.turkupetcentre.net/petanalysis/model_fur.html",
  " ",
  "Usage: @P [Options] inputfile ttacfile starttime endtime resultfile",
  " ",
  "FUR calculation start and stop time must be entered in minutes;",
  "set both to zero to use the whole time range from TTAC data.",
  " ",
  "Options:",
  " -Ca=<value>",
  "     Concentration of native substrate in arterial plasma (mM),",
  "     for example plasma glucose in [18F]FDG studies.",
  "     With this option the metabolic rate (umol/(min*100 g)) is calculated.",
  " -LC=<value>",
  "     Lumped Constant in MR calculation; default is 1.0.",
  " -density=<value>",
  "     Tissue density in MR calculation; default is 1.0 g/ml.",
  " -curve=<filename>",
  "     FUR as a function of time is written in specified file; this can be",
  "     used to study the time-dependence of FUR estimates.",
  " -it=<Time (min)>",
  "     Input AUC is normally calculated from 0 to the middle time of FUR",
  "     calculation time; in special cases this option can be used to",
  "     calculate it from 0 to the specified time instead.",
  " -deriv[ative]",
  "     Tentative option for calculating FUR as ratio of tissue derivative and",
  "     plasma. This does not affect the FUR curve made with option -curve.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1. Calculation of FUR from dynamic PET data from 45 to 60 min:",
  "     @P ua2918ap.kbq ua2918dy1.dft 45 60 ua2918fur.res",
  " ",
  "Example 2. Calculation of glucose uptake, when tissue density is 1.04,",
  "plasma glucose concentration is 5.2 mM, lumped constant is 0.52, from",
  "a static (one frame) scan:",
  "     @P -Ca=5.2 -LC=0.52 -d=1.04 a864ap.kbq a864dy1.dft 0 0 a864mrglu.res",
  " ",
  "FUR and MR results are saved in result file format by default, but if",
  "filename extension is set to .dft, .csv, .dat, or .html, results are saved",
  "in those formats.",
  " ",
  "The unit of FUR is (mL plasma)/(min*(mL tissue)) by default, and",
  "umol/(min*100 g) if metabolic rate is calculated.",
  " ",
  "See also: patlak, dftinteg, taccalc, tactime, dftsuv, rescoll, imgfur",
  " ",
  "Keywords: TAC, modelling, FUR, retention index, irreversible uptake",
  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;
/*****************************************************************************/
 
/*****************************************************************************/
#ifndef DEFAULT_LC
#define DEFAULT_LC 1.00
#endif
#ifndef DEFAULT_DENSITY
#define DEFAULT_DENSITY 1.00
#endif
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int        ai, help=0, version=0, verbose=1;
  int        ret;
  char       inpfile[FILENAME_MAX], petfile[FILENAME_MAX];
  char       curfile[FILENAME_MAX], outfile[FILENAME_MAX], tmp[1024], *cptr;
  DFT        input, auc, pet, avg;
  RES        res;
  double     startTime=-1.0, endTime=-1.0, aucTime=-1.0;
  double     LC=-1.0, Ca=-1.0, density=-1.0;
  int        fur_mode=0; // 0=traditional Ct/iCp, 1=derivative dCt/Cp
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=petfile[0]=outfile[0]=curfile[0]=(char)0;
  dftInit(&input); dftInit(&auc); dftInit(&pet); dftInit(&avg); resInit(&res);
  /* Get 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, "CA=", 3)==0) {
      Ca=atof_dpi(cptr+3); if(Ca>0.0) continue;
    } else if(strncasecmp(cptr, "LC=", 3)==0) {
      LC=atof_dpi(cptr+3); if(LC>0.0) continue;
    } else if(strncasecmp(cptr, "D=", 2)==0) {
      density=atof_dpi(cptr+2); if(density>0.0) continue;
    } else if(strncasecmp(cptr, "DENSITY=", 8)==0) {
      density=atof_dpi(cptr+8); if(density>0.0) continue;
    } else if(strncasecmp(cptr, "DERIVATIVE", 5)==0) {
      fur_mode=1; continue;
    } else if(strncasecmp(cptr, "CURVE=", 6)==0 && strlen(cptr)>6) {
      strcpy(curfile, cptr+6); continue;
    } else if(strncasecmp(cptr, "IT=", 3)==0) {
      aucTime=atof_dpi(cptr+3); if(aucTime>0.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(!petfile[0]) {
      strcpy(petfile, argv[ai]); continue;
    } else if(startTime<0.0) {
      startTime=atof_dpi(argv[ai]); if(startTime<0.0) startTime=0.0;
      continue;
    } else if(endTime<0.0) {
      endTime=atof_dpi(argv[ai]); if(endTime<0.0) endTime=0.0;
      continue;
    } else if(!outfile[0]) {
      strcpy(outfile, argv[ai]); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!outfile[0]) {
    fprintf(stderr, "Error: missing result file name.\n");
    return(1);
  }
 
  /* If MR will be calculated, set defaults and give warnings as necessary */
  if(Ca>0.0) {
    if(LC<0.0) {
      LC=DEFAULT_LC;
      fprintf(stderr, "Warning: LC not set, using default %g\n", LC);
    }
    if(density<0.0) {
      density=DEFAULT_DENSITY;
      fprintf(stderr, "Warning: tissue density not set, using default %g\n",
        density);
    } 
  } else { /* Warnings if density or LC set when MR will not be calculated */
    if(LC>0.0) fprintf(stderr, "Warning: LC was set but is not used.\n");
    if(density>0.0) 
      fprintf(stderr, "Warning: tissue density was set but is not used.\n");
  }
  /* Check the time range */
  if(endTime<startTime) {
    fprintf(stderr, "Error: invalid time range.\n");
    return(1);
  }
  if(startTime==endTime && startTime>0.5) {
    startTime-=0.5; endTime+=0.5;
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("inpfile := %s\n", inpfile);
    printf("petfile := %s\n", petfile);
    printf("outfile := %s\n", outfile);
    printf("curfile := %s\n", curfile);
    if(endTime>0.0) {
      printf("startTime := %g min\n", startTime);
      printf("endTime := %g min\n", endTime);
    }
    if(aucTime>0.0) {
      printf("aucTime := %g min\n", aucTime);
    }
    if(Ca>0.0) {
      printf("Ca := %g\n", Ca);
      printf("LC := %g\n", LC);
      printf("density := %g\n", density);
    }
    printf("fur_mode := %d\n", fur_mode);
    fflush(stdout);
  }
 
 
  /*
   *  Read regional data
   */
  if(verbose>1) printf("reading regional data %s\n", petfile);
  if(dftRead(petfile, &pet)) {
    fprintf(stderr, "Error in reading '%s': %s\n", petfile, dfterrmsg);
    return(2);
  }
  if(dft_nr_of_NA(&pet)>0) {  // check if file contains NAs (missing values)
    fprintf(stderr, "Error: missing values in %s\n", petfile);
    dftEmpty(&pet); return(2);
  }
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&pet);
 
  /* Make sure that there is no overlap in frame times */
  if(pet.timetype==DFT_TIME_STARTEND) {
    if(verbose>2) fprintf(stdout, "checking frame overlap in %s\n", petfile);
    ret=dftDeleteFrameOverlap(&pet);
    if(ret) {
      fprintf(stderr, "Error: %s has overlapping frame times.\n", petfile);
      dftEmpty(&pet); return(2);
    }
  }
 
  /* Set time unit to min */
  ret=dftTimeunitConversion(&pet, TUNIT_MIN);
  if(ret) fprintf(stderr, "Warning: check that regional data times are in minutes.\n");
  /* If user did not specify start and end times, then get those from data */
  if(endTime<=1.0E-02) {
    if(pet.timetype==DFT_TIME_STARTEND) {
      startTime=pet.x1[0]; endTime=pet.x2[pet.frameNr-1];
    } else {
      startTime=pet.x[0]; endTime=pet.x[pet.frameNr-1];
    }
    if(verbose>1) {
      printf("startTime := %g min\n", startTime);
      printf("endTime := %g min\n", endTime);
    }
  }
 
 
  /*
   *  Calculate the average (or slope) over specified time range
   */
  if(verbose>1) printf("calculating average\n");
  ret=dftTimeIntegral(&pet, startTime, endTime, &avg, 1, tmp, verbose-3);
  if(ret!=0) {
    fprintf(stderr, "Error: %s.\n", tmp);
    if(verbose>2)
      printf("dftTimeIntegral(pet, %g, %g, avg, 1, tmp) := %d\n", startTime, endTime, ret);
    dftEmpty(&pet); dftEmpty(&avg); return(3);
  }
  if(verbose>1) fprintf(stdout, "%s.\n", tmp);
  if(fur_mode==1) {
    int ri;
    double k, ksd, b, bsd, r, ysd;
    if(verbose>0) {
      fprintf(stdout, "calculating slope\n"); fflush(stdout);}
    for(ri=0; ri<pet.voiNr; ri++) {
      ret=pearson4(pet.x, pet.voi[ri].y, pet.frameNr, startTime, endTime,
                  &k, &ksd, &b, &bsd, &r, &ysd);
      if(ret!=0) break;
      avg.voi[ri].y[0]=k;
    }
    if(ret!=0) {
      fprintf(stderr, "Error: tissue slope calculation not successful.\n");
      if(verbose>1) printf("      ret := %d\n", ret);
      dftEmpty(&pet); dftEmpty(&avg); return(3);
    }
  }
  if(verbose>2) {
    printf("Regional tissue value or derivative\n");
    for(int ri=0; ri<avg.voiNr; ri++)
      printf("%s : %g\n", avg.voi[ri].name, avg.voi[ri].y[0]);
  }
 
  /*
   *  Read input data
   */
  if(verbose>1) printf("Reading input file %s\n", inpfile);
  if(dftRead(inpfile, &input)) {
    fprintf(stderr, "Error in reading '%s': %s\n", inpfile, dfterrmsg);
    dftEmpty(&pet); dftEmpty(&avg); return(4);
  }
  if(input.voiNr>1) {
    fprintf(stderr, "Warning: only first TAC is used as input.\n");
    input.voiNr=1;
  }
  if(dft_nr_of_NA(&input)>0) {  // check if file contains NAs (missing values)
    fprintf(stderr, "Error: missing values in %s\n", inpfile);
    dftEmpty(&pet); dftEmpty(&avg); dftEmpty(&input); return(4);
  }
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&pet);
 
  /* Set time unit to min */
  ret=dftTimeunitConversion(&input, TUNIT_MIN);
  if(ret) fprintf(stderr, "Warning: check that input times are in minutes.\n");
 
  /*
   *  Check the regional and plasma TAC concentration units
   */
  ret=dftUnitConversion(&input, dftUnitId(pet.unit));
  if(ret!=0) {
    fprintf(stderr, "Warning: check the units of input and regional data.\n");
  }
 
 
  /*
   *  Calculate and save FUR curve, if required;
   *  this is not used in calculation of regional FUR
   */
  if(curfile[0]) {
    if(verbose>1) {printf("calculating FUR curve\n"); fflush(stdout);}
    DFT fur;
    int fi, ri;
    dftInit(&fur); ret=dftdup(&pet, &fur);
    if(ret!=0) {
      fprintf(stderr, "Warning: cannot allocate memory for FUR curves\n");
      goto failed;
    }
    fur.frameNr=0;
    for(fi=0; fi<pet.frameNr; fi++) {
      if(pet.x[fi]<startTime || pet.x[fi]>endTime) continue; 
      /* AUC 0-t for each frame */
      if(pet.x[fi]>0.0) {
        ret=dftTimeIntegral(&input, 0.0, pet.x[fi], &auc, 0, tmp, verbose-4);
        if(ret!=0) {fprintf(stderr, "Warning (%d): %s\n", ret, tmp); break;}
        if(auc.voi[0].y[0]<1.0E-006) continue;
      } else continue;
      fur.x1[fur.frameNr]=pet.x1[fi]; fur.x2[fur.frameNr]=pet.x2[fi];
      fur.x[fur.frameNr]=pet.x[fi]; fur.w[fur.frameNr]=pet.w[fi];
      /* Divide each region by AUC */
      for(ri=0; ri<fur.voiNr; ri++)
        fur.voi[ri].y[fur.frameNr]= pet.voi[ri].y[fi] / auc.voi[0].y[0];
      fur.frameNr++;
    }
    if(ret!=0) goto failed;
    /* Write FUR curve */
    dftUnitToDFT(&fur, CUNIT_PER_MIN);
    dftSetComments(&fur);
    if(dftWrite(&fur, curfile))
      fprintf(stderr, "Error in writing %s: %s\n", curfile, dfterrmsg);
    failed:
    dftEmpty(&fur); fflush(stderr);
  }
 
 
  /*
   *  Calculate input integral from 0 to PET middle time point,
   *  or input value at PET middle time point
   */
  if(aucTime<=0.0) aucTime=0.5*(startTime+endTime);
  if(fur_mode==0)
    ret=dftTimeIntegral(&input, 0.0, aucTime, &auc, 0, tmp, verbose-3);
  else
    ret=dftTimeIntegral(&input, startTime, endTime, &auc, 1, tmp, verbose-3);
  if(ret!=0) {
    fprintf(stderr, "Error (%d): %s\n", ret, tmp);
    dftEmpty(&avg); dftEmpty(&input); dftEmpty(&auc); dftEmpty(&pet); 
    return(6);
  }
  if(verbose>1) {
    if(fur_mode==0)
      printf("AUC[%g-%g] := %g\n", auc.x1[0], auc.x2[0], auc.voi[0].y[0]);
    else
      printf("Input[%g-%g] := %g\n", auc.x1[0], auc.x2[0], auc.voi[0].y[0]);
  }
  /* Original input curve is not needed anymore */
  dftEmpty(&input);
 
 
  /*
   *  Divide average regional data by input AUC, or
   *  regional concentration derivative by input
   */
  for(int ri=0; ri<avg.voiNr; ri++) avg.voi[ri].y[0]/=auc.voi[0].y[0];
  dftUnitToDFT(&avg, CUNIT_PER_MIN);
 
  /* Input AUC is not needed anymore */
  dftEmpty(&auc);
 
 
  /*
   *  Calculate metabolic rate, if necessary
   */
  if(Ca>0.0) {
    double MRf;
    MRf=100.*Ca/(density*LC);
    if(verbose>1)
      fprintf(stdout, "converting FUR to metabolic rate with factor %g\n", MRf);
    for(int ri=0; ri<avg.voiNr; ri++) avg.voi[ri].y[0]*=MRf;
    dftUnitToDFT(&avg, CUNIT_UMOL_PER_MIN_PER_100G);
  }
 
 
  /*
   *  Save FUR/MR data
   */
  ret=backupExistingFile(outfile, NULL, tmp); if(ret!=0) {
    fprintf(stderr, "Error: %s\n", tmp); 
    dftEmpty(&avg); dftEmpty(&pet); return(11);
  }
  if(verbose>2) printf("%s\n", tmp);
  /* If filename extension is .dft, .csv, or .dat, then save in DFT, CSV-UK, or 
     simple format, respectively, otherwise as RES */ 
  cptr=strrchr(outfile, '.'); 
  if(cptr!=NULL && strcasecmp(cptr, ".DFT")==0) {
    dftSetComments(&avg);
    if(dftWrite(&avg, outfile)) {
      fprintf(stderr, "Error in writing %s: %s\n", outfile, dfterrmsg);
      dftEmpty(&avg); dftEmpty(&pet); return(12);
    }
  } else if(cptr!=NULL && strcasecmp(cptr, ".CSV")==0) {
    avg._type=DFT_FORMAT_CSV_UK;
    avg.timetype=DFT_TIME_MIDDLE;
    if(dftWrite(&avg, outfile)) {
      fprintf(stderr, "Error in writing %s: %s\n", outfile, dfterrmsg);
      dftEmpty(&avg); dftEmpty(&pet); return(12);
    }
  } else if(cptr!=NULL && strcasecmp(cptr, ".DAT")==0) {
    avg._type=DFT_FORMAT_PLAIN;
    avg.timetype=DFT_TIME_MIDDLE;
    if(dftWrite(&avg, outfile)) {
      fprintf(stderr, "Error in writing %s: %s\n", outfile, dfterrmsg);
      dftEmpty(&avg); dftEmpty(&pet); return(12);
    }
  } else {
    /* Copy DFT data into RES */
    ret=dftToResult(&avg, &res, tmp);
    if(ret!=0) {
      fprintf(stderr, "Error: %s.\n", tmp);
      dftEmpty(&avg); resEmpty(&res); dftEmpty(&pet); return(13);
    }
    /* Set more header information */
    tpcProgramName(argv[0], 1, 1, res.program, 256);
    cptr=strrchr(petfile, '/'); if(cptr==NULL) cptr=strrchr(petfile, '\\');
    if(cptr==NULL) cptr=petfile; else cptr++; strcpy(res.datafile, cptr);
    cptr=strrchr(inpfile, '/'); if(cptr==NULL) cptr=strrchr(inpfile, '\\');
    if(cptr==NULL) cptr=inpfile; else cptr++; strcpy(res.plasmafile, cptr);
    if(Ca>0.0) {
      res.concentration=Ca;
      res.lc=LC;
      res.density=density;
    }
    if(Ca>0.0) strcpy(res.parname[0], "MR"); else strcpy(res.parname[0], "FUR");
    /* Save RES */
    ret=resWrite(&res, outfile, verbose-3); 
    if(ret) {
      fprintf(stderr, "  Error (%d) in writing file %s\n", ret, outfile);
      dftEmpty(&avg); resEmpty(&res); dftEmpty(&pet); return(11);
    }
    resEmpty(&res);
  }
  if(verbose>0) {
    if(Ca<=0.0) fprintf(stdout, "FUR(s) saved in %s.\n", outfile);
    else fprintf(stdout, "MRs saved in %s.\n", outfile);
  }
 
  /* Free memory */
  dftEmpty(&avg);
  dftEmpty(&pet);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/