dftsuv.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 standardized uptake values (SUV, DUR, DAR) or PID/L from PET",
  "time-activity curves (TAC), injected radioligand dose and subject weight.",
  "SUV is calculated as a mean value between specified sample times, but",
  "optionally TACs can be saved in SUV units.",
  "Information on SUV in:",
  "https://www.turkupetcentre.net/petanalysis/model_suv.html",
  " ",
  "Usage: @P [Options] tacfile starttime endtime dose weight [resultfile]",
  " ",
  "Options:",
  " -curve=<Filename for SUV curve>",
  "     Save regional SUVs from whole measurement time range.",
  " -density=<Tissue density (g/ml)>",
  "     Calculate results per tissue mass instead of tissue volume.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "SUV calculation start and stop time must be entered in minutes;",
  "set both to zero to use the whole time range from TAC data.",
  " ",
  "Injected dose must be given in units MBq at the time of injection.",
  " ",
  "Subject weight must be given in kg or liters.",
  "Instead of SUV, the percentage of injected dose per tissue volume ",
  "(PID/L, %i.d./L) is calculated, if the subject weight is set to 0 (or below).",
  " ",
  "TAC file must be correct for physical decay to the injection time.",
  "If the units of radioactivity concentrations and sample times are not",
  "specified inside the TAC file, the units are assumed to be in kBq/mL and min;",
  "units can be specified by adding line(s) to the end of the TAC file,",
  "for example:",
  "# unit := Bq/cc",
  "# time_unit := sec",
  " ",
  "Results are saved in result file format (.res) by default, but if",
  "file name extension is set to .dft, .csv, .dat, or .html, results are saved",
  "in those formats.",
  " ",
  "Example 1:",
  "Calculate SUV40-60 and save SUV TAC with command",
  "     @P -curve=uia15suv.dat uia15.dat 40 60 330 77 uia15suv40-60.res",
  "Example 2:",
  "Calculate SUV from available time range with command",
  "     @P uia15.dat 0 0 330 77 uia15suv.res",
  " ",
  "See also: imgsuv, regfur, tactime, taccalc, tacunit, tac2suv, rescoll",
  " ",
  "Keywords: TAC, SUV, DUR, DAR, PID, dose, modelling",
  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     ri, fi, ret;
  char    dfile[FILENAME_MAX], rfile[FILENAME_MAX], sfile[FILENAME_MAX];
  char   *cptr, tmp[128];
  DFT     dft, suv;
  RES     res;
  double  time1=-1.0, time2=-1.0, weight=nan(""), dose=-1.0, density=-1.0, f;
  double  t1, t2;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dfile[0]=rfile[0]=sfile[0]=(char)0;
  dftInit(&dft); dftInit(&suv); 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, "CURVE=", 6)==0) {
      cptr+=6; if(strlen(cptr)>0) {strlcpy(sfile, cptr, FILENAME_MAX); continue;}
    } else if(strncasecmp(cptr, "C=", 2)==0) { // for compatibility, remove later
      cptr+=2; if(strlen(cptr)>0) {strlcpy(sfile, cptr, FILENAME_MAX); continue;}
    } else if(strncasecmp(cptr, "DENSITY=", 8)==0) {
      cptr+=8; density=atof_dpi(cptr); if(density>0.0 && density<3.0) continue;
    } else if(strncasecmp(cptr, "D=", 2)==0) { // for compatibility, remove later
      cptr+=2; density=atof_dpi(cptr); if(density>0.0 && density<3.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 */
  if(ai<argc) {strlcpy(dfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    ret=atof_with_check(argv[ai], &time1);
    if(ret!=0 || !(time1>=0.0)) {
      fprintf(stderr, "Error: invalid start time '%s'.\n", argv[ai]);
      return(1);
    }
    ai++;
  }
  if(ai<argc) {
    ret=atof_with_check(argv[ai], &time2);
    if(ret!=0 || !(time2>=time1)) {
      fprintf(stderr, "Error: invalid end time '%s'.\n", argv[ai]);
      return(1);
    }
    ai++;
  }
  if(ai<argc) {
    ret=atof_with_check(argv[ai], &dose);
    if(ret!=0 || !(dose>0.0)) {
      fprintf(stderr, "Error: invalid dose '%s'.\n", argv[ai]);
      return(1);
    }
    ai++;
  }
  if(ai<argc) {
    double v;
    ret=atof_with_check(argv[ai], &v);
    if(ret!=0 || !(v>=0.0)) {
      fprintf(stderr, "Error: invalid weight '%s'.\n", argv[ai]);
      return(1);
    }
    if(v>1.0E-100) weight=v;
    ai++;
  }
  if(ai<argc) {strlcpy(rfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
#if(0)
  /* Process other arguments, starting from the first non-option */
  for(; ai<argc; ai++) {
    if(!dfile[0]) {
      strlcpy(dfile, argv[ai], FILENAME_MAX); continue;
    } else if(time1<0) {
      ret=atof_with_check(argv[ai], &time1);
      if(ret==0 && time1>=0.0) continue;
    } else if(time2<0) {
      ret=atof_with_check(argv[ai], &time2);
      if(ret==0 && time2>=time1) continue;
    } else if(dose<0) {
      ret=atof_with_check(argv[ai], &dose);
      if(ret==0 && dose>0.0) continue;
    } else if(isnan(weight)) {
      ret=atof_with_check(argv[ai], &weight);
      if(ret==0) continue;
    } else if(!rfile[0]) {
      strlcpy(rfile, argv[ai], FILENAME_MAX); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]); return(1);
  }
 
  /* Is something missing? */
  if(!rfile[0]) {
    fprintf(stderr, "Error: missing result file name.\n");
    return(1);
  }
  if(weight<1.0E-100) weight=nan("");
#endif
 
  /* Is something missing? */
  if(!(dose>0.0)) {
    fprintf(stderr, "Error: missing dose.\n");
    return(1);
  }
  if(!rfile[0] && !sfile[0]) {
    fprintf(stderr, "Error: missing result file name.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dfile := %s\n", dfile);
    printf("rfile := %s\n", rfile);
    printf("sfile := %s\n", sfile);
    printf("Time range := %g-%g min\n", time1, time2);
    printf("dose := %g MBq\n", dose);
    if(!isnan(weight)) printf("weight := %g kg\n", weight);
    printf("density := %g g/ml\n", density);
  }
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading TAC data in %s\n", dfile);
  if(dftRead(dfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dfile, dfterrmsg);
    dftEmpty(&dft); return(2);
  }
  if(verbose>1) {
    printf("tacNr := %d\n", dft.voiNr);
    printf("sampleNr := %d\n", dft.frameNr);
  }
  /* Check if file contains NAs (missing values) */
  if(dft_nr_of_NA(&dft)>0) {
    fprintf(stderr, "Error: missing values in %s\n", dfile);
    dftEmpty(&dft); return(2);
  }
  /* Sort data by increasing sample time */
  dftSortByFrame(&dft);
 
  /* Make sure that there is no overlap in frame times */
  if(dft.timetype==DFT_TIME_STARTEND) {
    if(verbose>2) fprintf(stdout, "checking frame overlap in %s\n", dfile);
    ret=dftDeleteFrameOverlap(&dft);
    if(ret) {
      fprintf(stderr, "Error: %s has overlapping frame times.\n", dfile);
      dftEmpty(&dft); return(2);
    }
  }
 
  /* Check the radioactivity units */
  if(petCunitId(dft.unit)==CUNIT_UNKNOWN) {
    /* assuming kBq/mL */
  } else if(petCunitId(dft.unit)==CUNIT_KBQ_PER_ML) {
    /* ok as such */
  } else {
    /* Try to convert to kBq/mL */
    ret=dftUnitConversion(&dft, CUNIT_KBQ_PER_ML);
    if(ret!=0) {
      fprintf(stderr, "Error: activity unit %s is not supported.\n", dft.unit);
      dftEmpty(&dft); return(3);
    }
    if(verbose>0) printf("Data converted to %s\n", dft.unit);
  }
    
  /* Check the time units; convert start and end times if necessary */
  if(time2>0.0) {
    if(dft.timeunit==TUNIT_SEC) {
      time1*=60.; time2*=60.;
      if(verbose>0) fprintf(stderr, "Note: start and end times converted to seconds.\n");
    } else if(dft.timeunit==TUNIT_UNKNOWN) {
      fprintf(stderr, "Warning: unknown data time units.\n");
      fprintf(stderr, "Warning: assuming that data and time range are in same units.\n");
    }
  }
 
  /* Get the TAC start and end times */
  if(dft.timetype==DFT_TIME_STARTEND) {
    t1=dft.x1[0]; t2=dft.x2[dft.frameNr-1];
  } else {
    t1=dft.x[0]; t2=dft.x[dft.frameNr-1];
  }
  if(verbose>1) {
    printf("tac_starttime := %g\n", t1);
    printf("tac_endtime := %g\n", t2);
  }
  if(time2<1.0E-08) {
    /* If user did not specify start and end times, then get those from data */
    time1=t1; time2=t2;
  } else {
    /* otherwise check */
    double allow_dt;
    allow_dt=0.19*(t2-t1);
    if(time2>t2+allow_dt) {
      fprintf(stderr, "Error: time range %g-%g is outside measurement range %g-%g.\n",
              time1, time2, t1, t2);
      dftEmpty(&dft); return(3);
    }
  }
 
 
  /*
   *  Calculate SUV or PID curve
   */
  f=dose;
  if(!isnan(weight)) {
    if(verbose>1) fprintf(stdout, "calculating SUV TAC\n");
    f/=weight; /* Act / (Dose/Weight) */
  } else {
    if(verbose>1) fprintf(stdout, "calculating PID TAC\n");
    /* Percent of injected dose / mL */
    f*=1000.0/100.0;   /* 100% * Act / (Dose*1000) */
    /* per L instead of mL */
    f/=1000.0;
  }
  if(density>0.0) f*=density; /* Act / Density */
  if(verbose>2) printf("conversion factor := %g\n", f);
  /* Calculate curve */
  for(ri=0; ri<dft.voiNr; ri++)
    for(fi=0; fi<dft.frameNr; fi++) 
      dft.voi[ri].y[fi]/=f;
  /* Set calibration units */
  if(weight>0) {
    strcpy(dft.unit, "g/");
    if(density>0) strcat(dft.unit, "g"); else strcat(dft.unit, "mL");
  } else {
    strcpy(dft.unit, "%i.d./");
    if(density>0) strcat(dft.unit, "kg"); else strcat(dft.unit, "L");
  }
 
  /* Remove possible path from TAC file name, to be saved in results */
  filenameRmPath(dfile);
 
 
  /*
   *  Save the curves, if requested
   */
  if(sfile[0]) {
    if(verbose>1) printf("saving converted TAC(s) in %s\n", sfile);
    tpcProgramName(argv[0], 1, 0, tmp, 128);
    if(isnan(weight)) f=0.0; else f=weight;
    sprintf(dft.comments, "# %s: %s %g %g %g %g\n", tmp, dfile, time1, time2, dose, f);
    if(dftWrite(&dft, sfile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", sfile, dfterrmsg);
      dftEmpty(&dft);
      return(11);
    }
    if(!isnan(weight)) strcpy(tmp, "SUV"); else strcpy(tmp, "%i.d.");
    if(verbose>0) fprintf(stdout, "%s curves written in %s\n", tmp, sfile);
  }
 
  /*
   *  Stop, if SUV mean was not requested
   */
  if(!rfile[0]) {dftEmpty(&dft); return(0);}
 
 
  /*
   *  Calculate the TAC average over user-specified time range
   */
  if(verbose>1) printf("calculating average\n");
  ret=dftTimeIntegral(&dft, time1, time2, &suv, 1, tmp, verbose-3);
  if(ret!=0) {
    fprintf(stderr, "Error: %s.\n", tmp);
    if(verbose>2) 
      printf("dftTimeIntegral(suvtac, %g, %g, suvavg, 1, tmp) := %d\n", time1, time2, ret);
    dftEmpty(&dft); dftEmpty(&suv); 
    return(3);
  }
  if(verbose>1) fprintf(stdout, "%s.\n", tmp);
 
  /* Original TAC is not needed any more */
  dftEmpty(&dft);
 
 
  /*
   *  Save the SUV or PID mean
   */
  if(verbose>1) printf("saving average\n");
  /* If filename extension is .dft, .csv, or .dat, then save in DFT, CSV-UK, or 
     simple format, respectively, otherwise as RES */ 
  cptr=strrchr(rfile, '.'); 
  if(cptr!=NULL && strcasecmp(cptr, ".DFT")==0) {
    suv._type=DFT_FORMAT_STANDARD;
    dftSetComments(&suv);
    if(dftWrite(&suv, rfile)) {
      fprintf(stderr, "Error in writing %s: %s\n", rfile, dfterrmsg);
      dftEmpty(&suv); return(21);
    }
  } else if(cptr!=NULL && strcasecmp(cptr, ".CSV")==0) {
    suv._type=DFT_FORMAT_CSV_UK;
    suv.timetype=DFT_TIME_MIDDLE;
    if(dftWrite(&suv, rfile)) {
      fprintf(stderr, "Error in writing %s: %s\n", rfile, dfterrmsg);
      dftEmpty(&suv); return(22);
    }
  } else if(cptr!=NULL && strcasecmp(cptr, ".DAT")==0) {
    suv._type=DFT_FORMAT_PLAIN;
    suv.timetype=DFT_TIME_MIDDLE;
    if(dftWrite(&suv, rfile)) {
      fprintf(stderr, "Error in writing %s: %s\n", rfile, dfterrmsg);
      dftEmpty(&suv); return(23);
    }
  } else {
    /* Copy DFT data into RES */
    ret=dftToResult(&suv, &res, tmp);
    if(ret!=0) {
      fprintf(stderr, "Error: %s.\n", tmp);
      dftEmpty(&suv); resEmpty(&res); return(24);
    }
    /* Set more header information */
    tpcProgramName(argv[0], 1, 1, res.program, 256);
    strcpy(res.datafile, dfile);
    sprintf(res.reffile, "Dose %g MBq", dose);
    if(!isnan(weight)) {
      sprintf(tmp, " ; %g kg", weight);
      strcat(res.reffile, tmp);
    }
    if(density>0.0) res.density=density;
    /* Set the string for parameter names */
    if(isnan(weight) && density<=0.0) {
      strcpy(res.parname[0], "PID/L"); strcpy(res.parunit[0], "");
    } else if(isnan(weight) && density>0.0) {
      strcpy(res.parname[0], "PID/kg"); strcpy(res.parunit[0], "");
    } else if(density<=0.0) {
      strcpy(res.parname[0], "SUV"); strcpy(res.parunit[0], "mL/g");
    } else {
      strcpy(res.parname[0], "SUV"); strcpy(res.parunit[0], "mL/mL");
    }
    res.isweight=0;
    /* Save RES */
    ret=resWrite(&res, rfile, verbose-3); 
    if(ret) {
      fprintf(stderr, "  Error (%d) in writing file %s\n", ret, rfile);
      dftEmpty(&suv); resEmpty(&res); return(25);
    }
    resEmpty(&res);
  }
  if(verbose>0) {
    if(!isnan(weight)) strcpy(tmp, "SUV"); else strcpy(tmp, "%i.d.");
    if(verbose>0) fprintf(stdout, "%s written in %s\n", tmp, rfile);
  }
 
  dftEmpty(&suv);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/