dftratio.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[] = {
  "Calculate the regional tissue to reference tissue ratio, which is often",
  "used as a robust and model independent index of receptor availability.",
  "By default, ratio of AUCs is calculated as",
  " ",
  "         Regional AUC(t1..t2) ",
  "Ratio = --------------------- ",
  "        Reference AUC(t1..t2) ",
  " ",
  "Data concentration units are cancelled out from the ratio. Therefore",
  "the result is the same if it is calculated from original radioactivity",
  "concentrations or SUV TACs (SUV ratio, SUR).",
  " ",
  "Usage: @P [options] ttacfile reference t1 t2 resultfile",
  " ",
  "Reference can be given either as the name or number of the reference TAC",
  "inside the TTAC file or the name of file containing the reference TAC.",
  "Start and end times (t1 and t2) must be given in minutes, if TAC file(s)",
  "do not contain time units. To use the whole time range of ttacfile, enter",
  "zeroes as t1 and t2.",
  " ",
  "Options:",
  " -bound",
  "     Bound/free-ratio is calculated instead of Total/free;",
  "     reference region TAC is first subtracted from other regional TACs.",
  " -max",
  "     Bound TAC maximum value is searched between start and end times,",
  "     and Bound/free ratio is calculated at this time point; this should",
  "     not be used with noisy data.",
  " -rat=<Filename>",
  "     Ratio curve is saved for plotting.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P ut1234.tac occip 40 60 ut1234ratio.res",
  " ",
  "See also: imgratio, taccalc, dftinteg, interpol, tacunit, dftsuv",
  " ",
  "Keywords: TAC, modelling, AUC, ratio",
  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      *cptr, tmp[512];
  char       tacfile[FILENAME_MAX], resfile[FILENAME_MAX], ratfile[FILENAME_MAX];
  char       refname[FILENAME_MAX], reffile[FILENAME_MAX];
  double     tstart, tstop;
  int        ratmode=0; // 0=normal, 1=bound/ref, 2=bound/ref at bound max
 
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=resfile[0]=ratfile[0]=reffile[0]=refname[0]=(char)0;
  tstart=tstop=nan("");
  /* Get options first, because it affects what arguments are read */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strncasecmp(cptr, "BOUND", 1)==0) {
      ratmode=1; continue;
    } else if(strncasecmp(cptr, "MAX", 1)==0) {
      ratmode=2; continue;
    } else if(strncasecmp(cptr, "RAT=", 4)==0 && strlen(cptr)>4) {
      strlcpy(ratfile, cptr+4, FILENAME_MAX); continue;
    }
    fprintf(stderr, "Error: invalid option '%s'\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
  
  /* 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);}
  if(ai<argc) {strlcpy(reffile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    if(atof_with_check(argv[ai], &tstart)!=0) {
      fprintf(stderr, "Error: invalid t1 '%s'\n", argv[ai]); return(1);
    }
    ai++;
  }
  if(ai<argc) {
    if(atof_with_check(argv[ai], &tstop)!=0) {
      fprintf(stderr, "Error: invalid t2 '%s'\n", argv[ai]); return(1);
    }
    ai++;
  }
  if(ai<argc) {strlcpy(resfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: extra command-line argument.\n");
    return(1);
  }
  /* Check that we got what we need */
  if(!resfile[0]) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  if(tstop<tstart) {
    fprintf(stderr, "Error: illegal time range %g-%g\n", tstart, tstop);
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("reference := %s\n", reffile);
    printf("t1 := %g\n", tstart);
    printf("t2 := %g\n", tstop);
    printf("resfile := %s\n", resfile);
    printf("ratmode := %d\n", ratmode);
    if(ratfile[0]) printf("ratfile := %s\n", ratfile);
  }
 
 
  /*
   *  Read tissue data
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  DFT dft; dftInit(&dft);
  if(dftRead(tacfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", tacfile, dfterrmsg);
    return(2);
  }
  /* Check if file contains NAs (missing values) */
  if(dft_nr_of_NA(&dft)>0) {
    fprintf(stderr, "Error: missing values in %s\n", tacfile);
    dftEmpty(&dft); return(2);
  }
  /* Sort data by increasing sample time */
  dftSortByFrame(&dft);
 
  /* Convert time range to seconds, if data is given in seconds */
  if(dft.timeunit==TUNIT_SEC) {
    tstart*=60.0; tstop*=60.0;
    if(verbose) printf("time range is converted to seconds.\n");
  } else if(dft.timeunit==TUNIT_UNKNOWN) {
    fprintf(stderr, "Warning: unknown data time units.\n");
    fprintf(stderr, "Assuming that data and time range are in same units.\n");
  }
 
  /* If user gave zeroes for the start and end times, then get the times from tissue data */
  if(tstart==0.0 && tstop==0.0) {
    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];}
  }
 
  /* Check the time range */
  int n=0;
  if(dft.timetype==DFT_TIME_STARTEND) {
    if(ratmode==2 && tstop>dft.x2[dft.frameNr-1]) tstop=dft.x2[dft.frameNr-1];
    if(dft.x2[dft.frameNr-1] < 0.5*(tstart+tstop) || dft.x1[0]> 0.5*(tstart+tstop) ) {
      fprintf(stderr, "Error: data does not contain specified time range.\n");
      dftEmpty(&dft); return(2);
    }
    if(dft.x2[dft.frameNr-1] < tstop) {
      fprintf(stderr, "Warning: specified time range exceeds the data.\n");
    }
    for(int fi=0; fi<dft.frameNr; fi++)
      if(dft.x2[fi]>tstart && dft.x1[fi]<tstop) n++;
  } else {
    if(ratmode==2 && tstop>dft.x[dft.frameNr-1]) tstop=dft.x[dft.frameNr-1];
    if(dft.x[dft.frameNr-1] < 0.5*(tstart+tstop) || dft.x[0]> 0.5*(tstart+tstop) ) {
      fprintf(stderr, "Error: data does not contain the specified time range.\n");
      dftEmpty(&dft); return(2);
    }
    if(dft.x[dft.frameNr-1] < tstop) {
      fprintf(stderr, "Warning: specified time range exceeds the data.\n");
    }
    for(int fi=0; fi<dft.frameNr; fi++) if(dft.x[fi]>=tstart && dft.x[fi]<=tstop) n++;
  }
  if(verbose>1) {
    printf("tstart := %g min\n", tstart);
    printf("tstop := %g\n", tstop);
  }
  /* Verifying how many samples are inside the time range */
  if(n<1) {
    fprintf(stderr, "Error: data does not contain the specified time range.\n");
    dftEmpty(&dft); return(2);
  } else if(n<3) {
    fprintf(stderr, "Warning: only %d sample(s) in specified time range.\n", n);
  }
 
 
  /*
   *  Select reference TAC file or reference VOI in tissue file, 
   *  set ref = its voi number
   */
  if(verbose>1) printf("reading reference %s\n", reffile);
  int inputtype;
  double f, g;
  DFT rdft; dftInit(&rdft);
  if(dftReadinput(&rdft, &dft, reffile, &inputtype, &f, &g, 0, tmp, verbose-1)) {
    fprintf(stderr, "Error in reading '%s': %s\n", reffile, tmp);
    dftEmpty(&dft); dftEmpty(&rdft); return(3);
  }
  if(verbose>9) dftPrint(&rdft); 
  if(verbose>2) printf("inputtype := %d\n", inputtype);
  if(inputtype==5) { // Reference region name was given
    if(verbose>1)
      fprintf(stdout, "selected reference region := %s\n", rdft.voi[0].name);
    if(rdft.voiNr>1)
      fprintf(stderr, "Warning: %s selected of %d reference regions.\n",
        rdft.voi[0].name, rdft.voiNr);
    strcpy(refname, rdft.voi[0].name);
    strcpy(reffile, "");
  } else {
    if(verbose>2) printf("Reference tissue calibration unit := %s\n", rdft.unit);
    if(rdft.voiNr>1)
      fprintf(stderr, "Warning: only the first of reference curves is used.\n");
    if(f>tstart || g<tstop)
      fprintf(stderr, "Warning: time range may not be optimal for reference data.\n");
  }
 
 
  /*
   *  Subtract the reference region TAC, if required
   */
  if(ratmode!=0) {
    if(verbose>1) printf("Subtracting the ref TAC\n");
    for(int ri=0; ri<dft.voiNr; ri++)
      for(int fi=0; fi<dft.frameNr; fi++)    
        dft.voi[ri].y[fi]-=rdft.voi[0].y[fi];
    if(verbose>8) dftPrint(&dft);
  }
 
 
  /*
   *  Prepare the room for results
   */
  if(verbose>1) printf("initializing result data\n");
  RES res; resInit(&res); 
  if(res_allocate_with_dft(&res, &dft)!=0) {
    fprintf(stderr, "Error: cannot setup memory for results.\n");
    dftEmpty(&dft); dftEmpty(&rdft); return(6);
  }
  /* Copy titles & filenames */
  tpcProgramName(argv[0], 1, 1, res.program, 256);
  strcpy(res.datafile, tacfile);
  if(refname[0]) snprintf(res.refroi, 63, "%s", refname);
  if(reffile[0]) strcpy(res.reffile, reffile); 
  /* Set data range */
  sprintf(res.datarange, "%g - %g %s", tstart, tstop, petTunit(dft.timeunit));
  /* Set current time to results */
  res.time=time(NULL);
  res.isweight=0; res.Vb=-1.0;
  /* Set parameter number, including also the extra "parameters" */
  /* Set the string for parameter names */
  res.parNr=2;
  if(ratmode==2) {
    int pi;
    pi=0; strcpy(res.parname[pi], "B/F"); strcpy(res.parunit[pi], "");
    pi++; strcpy(res.parname[pi], "tMax");
          strcpy(res.parunit[pi], petTunit(dft.timeunit));
  } else {
    int pi;
    if(ratmode==1) {
      pi=0; strcpy(res.parname[pi], "B/F"); strcpy(res.parunit[pi], "");
    } else {
      pi=0; strcpy(res.parname[pi], "T/F"); strcpy(res.parunit[pi], "");
    }
    pi++; strcpy(res.parname[pi], "AUC"); strcpy(res.parunit[pi], "");
  }
 
 
  /*
   *  Compute the ratios
   */
  if(verbose>1) printf("computing the ratios\n");
  if(ratmode!=2) { /* AUC ratio */
 
 
#if(1) // 2022-05-09 (takes into account the frame lengths)
 
    double aucref;
    DFT auc; dftInit(&auc);
 
    /*
     *  Compute the AUC in reference region
     */
    if(dftTimeIntegral(&rdft, tstart, tstop, &auc, 0, NULL, verbose-3)!=0) {
      fprintf(stderr, "Error in integration of reference region TAC.\n");
      dftEmpty(&dft); dftEmpty(&rdft); resEmpty(&res); dftEmpty(&auc); return(7);
    }
    aucref=auc.voi[0].y[0];
    dftEmpty(&auc);
    if(verbose>2) printf("%s AUC = %g\n", rdft.voi[0].name, aucref);
    if(aucref<=0.0) {
      fprintf(stderr, "Error: invalid reference region AUC; check the data.\n");
      dftEmpty(&dft); dftEmpty(&rdft); resEmpty(&res); return(7);
    }
 
    /*
     *  Compute the AUC and ratio in regions of interest
     */
    if(dftTimeIntegral(&dft, tstart, tstop, &auc, 0, NULL, verbose-3)!=0) {
      fprintf(stderr, "Error in integration of regional TACs.\n");
      dftEmpty(&dft); dftEmpty(&rdft); resEmpty(&res); dftEmpty(&auc); return(8);
    }
    for(int ri=0; ri<dft.voiNr; ri++) {
      if(verbose>2) printf("%s %s %s AUC = %g\n",
        dft.voi[ri].voiname, dft.voi[ri].hemisphere, dft.voi[ri].place, auc.voi[ri].y[0]);
      /* Calculate the ratio */
      res.voi[ri].parameter[0]=auc.voi[ri].y[0]/aucref;
      res.voi[ri].parameter[1]=auc.voi[ri].y[0];
    } /* next region */
    dftEmpty(&auc);
 
#else // before 2022-05-09 (did not work with static data at all)
 
    /*
     *  Compute the AUC in reference region
     */
    double x[2], yi[2], aucref, aucroi;
    x[0]=tstart; x[1]=tstop;
    if(interpolate(dft.x, rdft.voi[0].y, dft.frameNr, x, NULL, yi, NULL, 2)!=0) {
      fprintf(stderr, "Error in integration of reference region TAC.\n");
      dftEmpty(&dft); dftEmpty(&rdft); resEmpty(&res); return(7);
    }
    aucref=yi[1]-yi[0];
    if(verbose>2) printf("%s AUC%g - AUC%g = %g\n", rdft.voi[0].name, yi[1], yi[0], aucref);
    if(aucref<=0.0) {
      fprintf(stderr, "Error: invalid reference region AUC; check the data.\n");
      dftEmpty(&dft); dftEmpty(&rdft); resEmpty(&res); return(7);
    }
 
    /*
     *  Compute the AUC and ratio in regions of interest
     */
    for(int ri=0; ri<dft.voiNr; ri++) {
      /* compute the AUC */
      if(interpolate(dft.x, dft.voi[ri].y, dft.frameNr, x, NULL, yi, NULL, 2)!=0) {
        fprintf(stderr, "Error in integration of region %d.\n", ri+1);
        dftEmpty(&dft); dftEmpty(&rdft); resEmpty(&res); return(8);
      }
      aucroi=yi[1]-yi[0];
      if(verbose>2) printf("%s %s %s AUC%g - AUC%g = %g\n",
        dft.voi[ri].voiname, dft.voi[ri].hemisphere, dft.voi[ri].place, yi[1], yi[0], aucroi);
      /* Calculate the ratio */
      res.voi[ri].parameter[0]=aucroi/aucref;
      res.voi[ri].parameter[1]=aucroi;
    } /* next region */
#endif
    
  } else { /* Ratio at max bound */
  
    double aucroi;
    for(int ri=0; ri<dft.voiNr; ri++) {
      /* Bound TAC is already computed; search its max */
      aucroi=0.0; n=-1;
      for(int fi=0; fi<dft.frameNr; fi++) {
        /* search only inside time range and with sensible ref values */
        if(dft.x[fi]>=tstart && dft.x[fi]<=tstop && rdft.voi[0].y[fi]>0.0) {
          if(dft.voi[ri].y[fi]>aucroi) {
            n=fi; aucroi=dft.voi[ri].y[fi];
          }
        }
      }
      /* Compute the ratio and set the tMax */
      if(n>=0) {
        res.voi[ri].parameter[0]=aucroi/rdft.voi[0].y[n];
        res.voi[ri].parameter[1]=dft.x[n];
      } else {
        res.voi[ri].parameter[0]=0.0;
        res.voi[ri].parameter[1]=0.0;
      }
    } /* next region */
 
  }
 
  /*
   *  Print results on screen
   */
  if(verbose>0) {resPrint(&res); fprintf(stdout, "\n");}
 
 
  /*
   *  Save results
   */
  if(verbose>1) printf("saving results\n");
  if(resWrite(&res, resfile, verbose-3)!=0) {
    fprintf(stderr, "Error in writing '%s': %s\n", resfile, reserrmsg);
    dftEmpty(&dft); dftEmpty(&rdft); resEmpty(&res);
    return(20);
  }
 
 
  /*
   *  Calculate and save the ratio curves, if required
   */
  if(ratfile[0]) {
    if(verbose>1) printf("Calculating ratio curves\n");
    for(int ri=0; ri<dft.voiNr; ri++) {
      for(int fi=0; fi<dft.frameNr; fi++) {
        if(isnan(rdft.voi[0].y[fi]) || isnan(dft.voi[ri].y[fi]) ||
           fabs(rdft.voi[0].y[fi])<1.0E-20) {
          dft.voi[ri].y[fi]=0.0;
        } else {
          dft.voi[ri].y[fi]/=rdft.voi[0].y[fi];
        }
      }
    }
    /* check that ratio is not < 0 or -1 */
    for(int ri=0; ri<dft.voiNr; ri++) for(int fi=0; fi<dft.frameNr; fi++) {
      if(ratmode>0) {if(dft.voi[ri].y[fi]<-1.0) dft.voi[ri].y[fi]=-1.0;}
      else {if(dft.voi[ri].y[fi]<0.0) dft.voi[ri].y[fi]=0.0;}
    }
    /* remove initial high peaks */
    f=g=0.0;
    for(int ri=0; ri<dft.voiNr; ri++) for(int fi=2; fi<dft.frameNr; fi++) {
      if(dft.voi[ri].y[fi]>f) f=dft.voi[ri].y[fi];
      else if(dft.voi[ri].y[fi]<g) g=dft.voi[ri].y[fi];
    }
    for(int ri=0; ri<dft.voiNr; ri++) for(int fi=0; fi<2 && fi<dft.frameNr; fi++) {
      if(dft.voi[ri].y[fi]>f) dft.voi[ri].y[fi]=f;
      else if(dft.voi[ri].y[fi]<g) dft.voi[ri].y[fi]=g;
    }
    /* Correct 'header' information */
    dftUnitToDFT(&dft, CUNIT_UNITLESS);
    dft.comments[0]=(char)0;
    dft.isweight=0;
    /* and write as ratio file */
    if(dftWrite(&dft, ratfile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", ratfile, dfterrmsg);
      dftEmpty(&dft); dftEmpty(&rdft); resEmpty(&res);
      return(22);
    }
  }
 
  dftEmpty(&dft); dftEmpty(&rdft); resEmpty(&res);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/