tac2suv.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 "tpctac.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Converts TAC radioactivity concentrations into standardized uptake values",
  "(SUV, DUR, DAR) or PID/L, using administered radioligand dose and subject",
  "weight. Injected dose must be given in units MBq at the time of injection.",
  "Subject weight must be given in kg or litres.",
  "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.",
  " ",
  "Usage: @P [Options] tacfile dose weight [outputfile]",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "TAC file must be correct for physical decay to the injection time.",
  "If the units of radioactivity concentrations are not specified inside",
  "the TAC file, the units are assumed to be in kBq/mL;",
  "units can be specified by adding line(s) to the end of the TAC file,",
  "for example:",
  "# unit := Bq/cc",
  "# time_unit := sec",
  " ",
  "If name for output file is not given, then TAC file is overwritten.",
  " ",
  "Example:",
  "     @P uia15.dat 330 77 uia15suv.dat",
  " ",
  "See also: tacunit, tactime, taccalc, dftinteg, dftsuv, imgsuv",
  " ",
  "Keywords: TAC, SUV, 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;
  char    tacfile[FILENAME_MAX], outfile[FILENAME_MAX];
  double  weight=nan(""), dose=nan("");
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=outfile[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;
    fprintf(stderr, "Error: invalid option '%s'\n", argv[ai]);
    return(1);
  } else break; 
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
  
  /* 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(tacfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    if(atofCheck(argv[ai], &dose) || !(dose>0.0)) {
      fprintf(stderr, "Error: invalid dose '%s'.\n", argv[ai]);
      return(1);
    }
    ai++;
  }
  if(ai<argc) {
    double v;
    if(atofCheck(argv[ai], &v)) {
      fprintf(stderr, "Error: invalid weight '%s'.\n", argv[ai]);
      return(1);
    }
    if(v>1.0E-100) weight=v;
    ai++;
  } else {fprintf(stderr, "Error: missing weight.\n"); return(1);}
  if(ai<argc) {strlcpy(outfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!(dose>0.0)) {
    fprintf(stderr, "Error: missing dose.\n");
    return(1);
  }
  if(!outfile[0]) strcpy(outfile, tacfile);
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("outfile := %s\n", outfile);
    printf("dose := %g MBq\n", dose);
    if(!isnan(weight)) printf("weight := %g kg\n", weight);
    fflush(stdout);
  }
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading TAC data in %s\n", tacfile);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  if(verbose>1) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("unit := %s\n", unitName(tac.cunit));
    printf("timeunit := %s\n", unitName(tac.tunit));
    fflush(stdout);
  }
 
  /* Check the radioactivity units */
  if(tac.cunit==UNIT_UNKNOWN) {
    tac.cunit=UNIT_KBQ_PER_ML;
    printf("Warning: assuming unit %s\n", unitName(tac.cunit));
  }
  if(tac.cunit!=UNIT_KBQ_PER_ML && tac.cunit!=UNIT_KBQ_PER_G) {
    int ret=1, oldunit=tac.cunit;
    if(unitDividerHasVolume(tac.cunit)) ret=tacYUnitConvert(&tac, UNIT_KBQ_PER_ML, &status);
    else if(unitDividerHasMass(tac.cunit)) ret=tacYUnitConvert(&tac, UNIT_KBQ_PER_G, &status);
    if(ret) {
      fprintf(stderr, "Error: cannot convert units %s\n", unitName(tac.tunit));
      return(3);
    }
    if(verbose>1) {
      printf("%s: units converted from '%s' to '%s'\n",
             tacfile, unitName(oldunit), unitName(tac.cunit));
      fflush(stdout);
    }
  }
 
 
  /*
   *  Calculate SUV or PID curve
   */
  double f=dose;
  if(!isnan(weight)) {
    if(verbose>1) printf("calculating SUV TAC\n");
    f/=weight; /* Act / (Dose/Weight) */
  } else {
    if(verbose>1) printf("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(verbose>2) {printf("conversion factor := %g\n", f); fflush(stdout);}
  /* Calculate curve */
  for(int i=0; i<tac.tacNr; i++) for(int j=0; j<tac.sampleNr; j++) tac.c[i].y[j]/=f;
  /* Set calibration units */
  if(unitDividerHasVolume(tac.cunit)) {
    if(!isnan(weight)) tac.cunit=UNIT_G_PER_ML; else tac.cunit=UNIT_PID_PER_L;
  } else if(unitDividerHasMass(tac.cunit)) {
    if(!isnan(weight)) tac.cunit=UNIT_UNITLESS; else tac.cunit=UNIT_PID_PER_KG;
  }
 
 
  /* 
   *  Save modified data
   */
  {
    if(verbose>1) printf("writing %s\n", outfile);
    FILE *fp; fp=fopen(outfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&tac); return(11);
    }
    int ret=tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(12);
    }
    tacFree(&tac);
  }
  if(verbose>0) {
    char tmp[32];
    if(!isnan(weight)) strcpy(tmp, "SUV"); else strcpy(tmp, "%i.d.");
    printf("%s written in %s\n", tmp, outfile);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/