logan.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"
#include "tpcpar.h"
#include "tpcli.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "NOT YET FUNCTIONAL",
  " ",
  "Calculate the distribution volume (Vt) or distribution volume ratio (DVR)",
  "using multiple-time graphical analysis (MTGA) for reversible PET",
  "ligands (Logan plot) (1,2,3) from regional PET time-activity curves (TACs).",
  " ",
  "Usage: @P [options] tacfile input starttime endtime resultfile",
  " ",
  "Options:",
  " -k2=<reference region k2>",
  "     With reference region input, the population average of reference",
  "     region k2 (or k2/(1+k5/k6) in 3-compartment model) is set.",
  " -BPnd",
  "     With reference input, BPnd (=DVR-1) is reported instead of DVR",
  " -BPnd=<reference region name>",
  "     With plasma input, BPnd can be calculated as DVroi/DVref-1",
  " -DVR=<reference region name>",
  "     With plasma input, DVR can be calculated as DVroi/DVref",
  " -BPp=<reference region name>",
  "     With plasma input, BPp can be calculated as DVroi-DVref",
  " -sd=<y|n>",
  "     Standard deviations are saved (y, default) or not saved (n) in results.",
  " -mid[=<y|n>]",
  "     Mid frame times are used (y) or not used (n, default) even if frame",
  "     start and end times are available.",
  " -svg=<Filename>",
  "     Plots are written in specified file in Scalable Vector Graphics (SVG) 1.1",
  "     format; specification in https://www.w3.org/TR/SVG/",
  " -plotdata=<Filename>",
  "     Data for plots is written in specified file in XHTML table format for",
  "     easy importing in Excel or OpenOffice spreadsheet, where the data can",
  "     be viewed; if filename extension is .dft, data is written in DFT format",
  " -C | -M | -P | -R",
  "     Options for selecting the least-squares line fit method:",
  "     -C for Traditional regression model",
  "     -M for Median of two-point slopes and intercepts (Cornish-Bowden)",
  "     -P for Perpendicular regression model (4)",
  "     -R for Iterative method (York 1966, Lybanon 1984, Reed 1992); default",
  "     If tissue file contains weights, the iterative method (-R) is weighted.",
  "     With other fitting methods the weights are not used.",
  "     With options -C and -R program can automatically find the linear plot",
  "     range, if fit start time is set to zero.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1:",
  "Tissue curves are in ut2345.dft and plasma curve in ut2345ap.dat;",
  "fitted data range is from 10 to 60 min; standard deviations are not needed:",
  "     @P -nosd ut2345.dft ut2345ap.dat 10 60 ut2345.res",
  " ",
  "Example 2:",
  "Tissue curves in ut1234.dft, including reference region 'cer'; reference",
  "tissue k2 is assumed to equal 0.163; plot is saved for viewing in",
  "ut1234pplot.svg:",
  "     @P -svg=ut2345lplot.svg ut2345.dft cer -k2=0.163 20 60 ut2345.res",
  " ",
  "References:",
  "1. Logan J, Fowler JS, Volkow ND, Wolf AP, Dewey SL, Schlyer DJ,",
  "   MacGregor RR, Hitzemann R, Bendriem B, Gatley SJ, Christman DR.",
  "   Graphical analysis of reversible radioligand binding from time-activity",
  "   measurements applied to [N-11C-methyl]-(-)-cocaine PET studies in human",
  "   subjects. J Cereb Blood Flow Metab 1990; 10: 740-747.",
  "2. Logan J, Fowler JS, Volkow ND, Wang GJ, Ding YS, Alexoff DL.",
  "   Distribution volume ratios without blood sampling from graphical",
  "   analysis of PET data. J Cereb Blood Flow Metab. 1996; 16: 834-840.",
  "3. Logan J. Graphical analysis of PET data applied to reversible and",
  "   irreversible tracers. Nucl Med Biol 2000; 27:661-670.",
  "4. Varga J & Szabo Z. Modified regression model for the Logan plot.",
  "   J Cereb Blood Flow Metab 2002; 22:240-244.",
  " ",
  "See also: imgdv, fitk4, lhsoldv, patlak, rescoll, res2html",
  " ",
  "Keywords: TAC, modelling, Vt, DVR",
  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], inpfile[FILENAME_MAX], resfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  double k2; // reference tissue k2; NaN if not applied
  int save_stat=1; // 0=do not save, 1=save SD and/or other fit statistics
  double tstart=nan(""), tstop=nan(""); // fit start and end times
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=inpfile[0]=resfile[0]=svgfile[0]=(char)0;
  k2=tstart=tstop=nan("");
  /* 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, "SVG=", 4)==0) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); if(strlen(svgfile)>0) continue;
    } else if(strncasecmp(cptr, "K2=", 3)==0) {
      k2=atofVerified(cptr+3); if(!isnan(k2)) continue;
    } else if(strncasecmp(cptr, "SD=", 3)==0) {
      save_stat=tpcYesNo(cptr+3); if(save_stat>=0) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'\n", argv[ai]);
    return(1);
  } else break; // later arguments 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);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) strlcpy(tacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(inpfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(atofCheck(argv[ai], &tstart)) {
      fprintf(stderr, "Error: invalid start time '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    if(atofCheck(argv[ai], &tstop)) {
      fprintf(stderr, "Error: invalid stop time '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) strlcpy(resfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!resfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* Is something wrong? */
  if(!isnan(tstart) && !isnan(tstop) && 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("inpfile := %s\n", inpfile);
    printf("resfile := %s\n", resfile);
    printf("svgfile := %s\n", svgfile);
    if(!isnan(k2)) printf("k2 := %g\n", k2);
    printf("tstart := %g\n", tstart);
    printf("tstop := %g\n", tstop);
    printf("save_stat := %d\n", save_stat);
  }
 
 
  /*
   *  Read tissue and input data
   */
  if(verbose>0) printf("reading tissue and input data\n");
  TAC tac, input; tacInit(&tac); tacInit(&input);
  int fitSampleNr;
  double fitdur;
  if(tstop>0.01) fitdur=tstop; else fitdur=1.0E+10;
 
  if(tacReadModelingData(tacfile, inpfile, NULL, NULL, &fitdur, 0,
                          &fitSampleNr, &tac, &input, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); tacFree(&input); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("tac.sampleNr := %d\n", tac.sampleNr);
    printf("input.sampleNr := %d\n", input.sampleNr);
    printf("fitSampleNr := %d\n", fitSampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    printf("fitdur := %g s\n", fitdur);
  }
 
 
 
  tacFree(&tac); tacFree(&input);
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
/*****************************************************************************/
 
/*****************************************************************************/