imgdv.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculates the distribution volume (Vt, DV) or distribution volume ratio (DVR)",
  "from dynamic PET image in ECAT, NIfTI-1, or Analyze 7.5 format,",
  "applying the multiple time graphical analysis (MTGA) approach for reversible",
  "tracer uptake (Logan plot) (1, 2, 3).",
  "Simple non-iterative perpendicular line fitting algorithm (4) is applied",
  "in determination of the slope of the Logan plot.",
  " ",
  "Usage: @P [Options] input image starttime resultimage",
  " ",
  "Input can be either a TAC file containing arterial plasma PTAC or reference",
  "region TTAC. Data must be in the same concentration units as the data in",
  "image, and times in min; however, if units are specified inside the input",
  "and image file, the program can automatically convert units.",
  "Start time of the linear fit must be given in minutes; set to zero to",
  "let the program automatically determine the start time (not recommended).",
  " ",
  "Options:",
  " -k2=<k2 of reference region>",
  "     With reference region input it may be necessary to specify also the",
  "     population average for reference region k2 (2).",
  " -thr=<threshold%>",
  "     Pixels with AUC less than (threshold/100 x input) are set to zero.",
  "     Default is 0 %.",
  " -max=<Max value>",
  "     Upper limit for Vt or DVR values; by default max is set pixel-wise",
  "     to 10 times the AUC ratio.",
  " -min=<Min value>",
  "     Lower limit for Vt or DVR values; 0 by default.",
  " -filter",
  "     Remove parametric pixel values that are over 4x higher than",
  "     their closest neighbours.",
  " -end=<Fit end time (min)>",
  "     By default line is fitted to the end of data. Use this option to enter",
  "     the fit end time.",
  " -v=<filename>",
  "     Y axis intercepts times -1 are written as an image in specified file.",
  " -n=<filename>",
  "     Numbers of selected plot data points are written as an image.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P ua3818ap.kbq ua3818dy1.v 20 ua3818dv.v",
  " ",
  "The unit of pixel values in the parametric Vt image is",
  "(mL plasma)/(mL tissue). Data in DVR image are unitless.",
  " ",
  "In theory, DVR equals BPnd+1, and therefore a binding potential image can",
  "be computed from DVR image by subtraction of unity, e.g. using program",
  "imgcalc. However, this may lead to a large number of pixels with negative",
  "values that may hamper the further statistical or visual analysis.",
  " ",
  "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: imgbfbp, imglhdv, img2tif, logan, imgki, imgcalc",
  " ",
  "Keywords: image, MTGA, Logan plot, modelling, Vt, DVR, distribution volume",
  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;
  linefit_range fit_range;
  int           fi, ret;
  int           param_filt=0;
  double        k2=-1.0;
  double        startTime=-1.0, lowerLimit=0.0, upperLimit=-1.0, fittime=-1.0;
  float         calcThreshold=0.0;
  char          inpfile[FILENAME_MAX], petfile[FILENAME_MAX];
  char          outfile[FILENAME_MAX], icfile[FILENAME_MAX];
  char          nrfile[FILENAME_MAX];
  char          tmp[1024], *cptr;
  IMG           img, out, icimg, nrimg;
  DFT           input;
  int           startSample=0, endSample=0, fitdimt=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=petfile[0]=outfile[0]=icfile[0]=nrfile[0]=(char)0;
  dftInit(&input);
  imgInit(&img); imgInit(&out); imgInit(&icimg); imgInit(&nrimg);
  /* Get options */
  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, "V=", 2)==0) {
      strcpy(icfile, cptr+2); if(strlen(icfile)>1) continue;
    } else if(strncasecmp(cptr, "N=", 2)==0) {
      strcpy(nrfile, cptr+2); if(strlen(nrfile)>1) continue;
    } else if(strncasecmp(cptr, "K2=", 3)==0) {
      k2=atof_dpi(cptr+3); if(k2>0.0) continue;
    } else if(strncasecmp(cptr, "MIN=", 4)==0) {
      if(atof_with_check(cptr+4, &lowerLimit)==0) continue;
    } else if(strncasecmp(cptr, "MAX=", 4)==0) {
      upperLimit=atof_dpi(cptr+4); if(upperLimit>0.0) continue;
    } else if(strncasecmp(cptr, "E=", 2)==0) {
      fittime=60.0*atof_dpi(cptr+2); if(fittime>0.0) continue;
    } else if(strncasecmp(cptr, "END=", 4)==0) {
      fittime=60.0*atof_dpi(cptr+4); if(fittime>0.0) continue;
    } else if(strncasecmp(cptr, "FILTER", 4)==0) {
      param_filt=1; continue;
    } else if(strncasecmp(cptr, "THR=", 4)==0) {
      double v; ret=atof_with_check(cptr+4, &v);
      if(ret==0) {calcThreshold=0.01*v; continue;} // negative is ok
    }
    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(inpfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(petfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {startTime=atof_dpi(argv[ai++]); if(startTime<0.0) startTime=0.0;}
  if(ai<argc) strlcpy(outfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    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(upperLimit>0.0 && upperLimit<=lowerLimit) {
    fprintf(stderr, "Error: invalid -min or -max setting.\n");
    return(1);
  }
 
  /* 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("startTime := %g min\n", startTime);
    if(k2>0.0) printf("k2 := %g\n", k2);
    printf("lowerLimit := %g\n", lowerLimit);
    if(upperLimit>0.0) printf("upperLimit := %g\n", upperLimit);
    printf("icfile := %s\n", icfile);
    printf("nrfile := %s\n", nrfile);
    printf("fittime := %g\n", fittime);
    printf("param_filt := %d\n", param_filt);
    printf("calcThreshold := %g%%\n", 100.*calcThreshold);
  }
  if(verbose>9) IMG_TEST=verbose-9; else IMG_TEST=0;
 
 
  /*
   *  Read PET image and input TAC
   */
  fittime/=60.; 
  ret=imgReadModelingData(
    petfile, NULL, inpfile, NULL, NULL, &fittime, &fitdimt, &img,
    &input, NULL, 1, NULL, verbose-2, tmp);
  if(ret!=0) {
    fprintf(stderr, "Error in reading data: %s.\n", tmp);
    imgEmpty(&img); dftEmpty(&input);
    return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  fittime*=60.;
  endSample=fitdimt-1;
  if(verbose>2) {
    printf("fittimeFinal := %g [min]\n", fittime/60.0);
    printf("fitdimt := %d\n", fitdimt);
    printf("endSampleIndex := %d\n", endSample);
    printf("endSample := %d\n", endSample+1);
  }
  /* Find out the first sample to use in the line fit */
  for(fi=startSample=0; fi<fitdimt; fi++)
    if((img.mid[fi]/60.0)>startTime) break; else startSample++;
  if(verbose>2) {
    printf("startSampleIndex := %d\n", startSample);
    printf("startSample := %d\n", startSample+1);
    printf("startTimeFinal := %g [min]\n", img.start[startSample]/60.0);
  }
  if((fitdimt-startSample)<2) {
    fprintf(stderr, "Error: too few time frames to fit.\n");
    imgEmpty(&img); dftEmpty(&input);
    return(3);
  }
 
#if(0)
  /*
   *  Thresholding
   */   
  if(verbose>1) fprintf(stdout, "thresholding\n");
  ret=imgThresholding(&img, calcThreshold, &n);
  if(ret!=0) {
    fprintf(stderr, "Error in thresholding the dynamic image: %s\n", img.statmsg);
    imgEmpty(&img); dftEmpty(&input);
    return(4);
  }
  if(verbose>1) {
    fprintf(stdout, "threshold_cutoff_nr := %d / %d\n", n, img.dimx*img.dimy*img.dimz);
    fflush(stdout);
  }
#endif
 
  /*
   *  Calculate Logan plot
   */
  /* Fixed or free time range */
  if(startTime>0.0) fit_range=PRESET; else fit_range=EXCLUDE_BEGIN;
  if(verbose>1) printf("fit_range := %d\n", (int)fit_range);
  /* Set optional output IMG struct pointers to NULL, if not needed */
  IMG *picimg=NULL, *pnrimg=NULL;
  if(icfile[0]) picimg=&icimg;
  if(nrfile[0]) pnrimg=&nrimg;
  /* MTGA */
  if(verbose>0) fprintf(stdout, "computing MTGA pixel-by-pixel\n");
  clock_t fitStart, fitFinish;
  fitStart=clock();
  ret=img_logan(&input, &img, startSample, endSample, fit_range, calcThreshold, k2, 
                &out, picimg, pnrimg, tmp, verbose-5);
  if(ret!=0) {
    fprintf(stderr, "Error (%d): %s\n", ret, tmp);
    imgEmpty(&img); dftEmpty(&input); return(5);
  }
  fitFinish=clock();
  if(verbose>0) {
    if(fitFinish>fitStart) fprintf(stdout, "done in %g seconds.\n",
      (double)(fitFinish - fitStart) / (double)CLOCKS_PER_SEC );
    else fprintf(stdout, "done.\n");
  }
  /* No need for the dynamic image or input data anymore */
  imgEmpty(&img); dftEmpty(&input);
 
 
 
  /*
   *  Filter the parametric images if required
   */
  /* Remove values that are higher than the user-defined limit */
  if(upperLimit>0.0) {
    if(verbose>0) fprintf(stdout, "filtering out pixel values exceeding %g\n", upperLimit);
    imgCutoff(&out, (float)upperLimit, 0);
  }
  /* Remove values that are lower than the user-defined or default limit */
  if(verbose>0)
    fprintf(stdout, "filtering out pixel values lower than %g\n", lowerLimit);
  imgCutoff(&out, (float)lowerLimit, 1);
  /*
   *  Filter out pixels that are over 4x higher than their
   *  closest neighbours
   */
  if(param_filt>0) {
    if(verbose>0) printf("filtering extreme values from parametric images\n");
    imgOutlierFilter(&out, 4.0);
    if(icfile[0]) imgOutlierFilter(&icimg, 4.0);
  }
 
 
  /*
   *  Save parametric Vt image
   */
  ret=backupExistingFile(outfile, NULL, tmp); if(ret!=0) {
    fprintf(stderr, "Error: %s\n", tmp); 
    imgEmpty(&out); imgEmpty(&icimg); imgEmpty(&nrimg); return(11);
  }
  if(verbose>1) printf("%s\n", tmp);
  ret=imgWrite(outfile, &out); 
  if(ret) {
    fprintf(stderr, "Error: %s\n", out.statmsg); 
    imgEmpty(&out); imgEmpty(&icimg); imgEmpty(&nrimg); return(12);
  }
  if(verbose>0) {
    fprintf(stdout, "Vt or DVR image %s saved.\n", outfile);
  }
 
  /*
   *  Save parametric Ic image
   */
  if(icfile[0]) {
    ret=backupExistingFile(icfile, NULL, tmp); if(ret!=0) {
      fprintf(stderr, "Error: %s\n", tmp); 
      imgEmpty(&out); imgEmpty(&icimg); imgEmpty(&nrimg); return(13);
    }
    if(verbose>1) printf("%s\n", tmp);
    ret=imgWrite(icfile, &icimg); if(ret) {
      fprintf(stderr, "Error: %s\n", icimg.statmsg); 
      imgEmpty(&out); imgEmpty(&icimg); imgEmpty(&nrimg); return(14);
    }
    if(verbose>0) fprintf(stdout, "Intercept image %s saved.\n", icfile);
  }
 
  /*
   *  Save parametric nr image
   */
  if(nrfile[0]) {
    ret=backupExistingFile(nrfile, NULL, tmp); if(ret!=0) {
      fprintf(stderr, "Error: %s\n", tmp); 
      imgEmpty(&out); imgEmpty(&icimg); imgEmpty(&nrimg); return(15);
    }
    if(verbose>1) printf("%s\n", tmp);
    ret=imgWrite(nrfile, &nrimg); if(ret) {
      fprintf(stderr, "Error: %s\n", nrimg.statmsg); 
      imgEmpty(&out); imgEmpty(&icimg); imgEmpty(&nrimg); return(16);
    }
    if(verbose>0) fprintf(stdout, "Nr image %s saved.\n", nrfile);
  }
 
  /* Free memory */
  imgEmpty(&out); imgEmpty(&icimg); imgEmpty(&nrimg);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/