imgki.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"
/*****************************************************************************/
 
/*****************************************************************************/
#ifndef DEFAULT_LC
#define DEFAULT_LC 1.00
#endif
#ifndef DEFAULT_DENSITY
#define DEFAULT_DENSITY 1.00
#endif
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculates the net influx (uptake) rate constant (Ki) from dynamic PET image",
  "in ECAT, NIfTI, or Analyze format, applying the multiple time graphical",
  "analysis (MTGA) approach for irreversible tracer uptake (Patlak plot).",
  "Simple non-iterative perpendicular line fitting algorithm (4) is applied.",
  " ",
  "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:",
  " -Ca=<value>",
  "     Concentration of native substrate in arterial plasma (mM),",
  "     for example plasma glucose in [18F]FDG studies.",
  "     With this option the metabolic rate (umol/(min*100 g)) is calculated.",
  " -LC=<value>",
  "     Lumped Constant in MR calculation; default is 1.0.",
  " -density=<value>",
  "     Tissue density in MR calculation; default is 1.0 g/ml.",
  " -thr=<threshold%>",
  "     Pixels with AUC less than (threshold/100 x input AUC) are set to zero;",
  "     default is 0%.",
//  "     Pixels with AUC less than (threshold/100 x max AUC) are set to zero.",
//  "     Default is 0 %.",
  " -noneg",
  "     Pixels with negative Ki are set to zero.",
  " -max=<Max value>",
  "     Upper limit for Ki values.",
  " -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 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 1. Calculation of Ki image, starting linear fit from 20 minutes:",
  "     @P ua2918ap.kbq ua2918dy1.v 20 ua2918ki.v",
  " ",
  "Example 2. Calculation of glucose uptake image, when tissue density is 1.02,",
  "           plasma glucose concentration is 5.2 mM, lumped constant is 1.0,",
  "           starting linear fit from 10 minutes after radiotracer injection:",
  "     @P -density=1.02 -Ca=5.2 -LC=1.0 s8642ap.kbq s8642dy1.v 10 s8642ki.v",
  " ",
  "The unit of pixel values in the parametric (Ki) image is",
  "(mL plasma)/(min*(mL tissue)) by default, and umol/(min*100 g) in metabolic",
  "rate image.",
  " ",
  "To calculate Ki images from [18F]FDOPA studies with plasma input:",
  "subtract reference region TAC from dynamic image before using this program.",
  " ",
  "References:",
  "1. Gjedde A. Calculation of cerebral glucose phosphorylation from brain",
  "   uptake of glucose analogs in vivo: a re-examination. Brain Res. 1982;",
  "   257:237-274.",
  "2. Patlak CS, Blasberg RG, Fenstermacher JD. Graphical evaluation of",
  "   blood-to-brain transfer constants from multiple-time uptake data.",
  "   J Cereb Blood Flow Metab 1983;3:1-7.",
  "3. Patlak CS, Blasberg RG. Graphical evaluation of blood-to-brain transfer",
  "   constants from multiple-time uptake data. Generalizations",
  "   J Cereb Blood Flow Metab 1985;5:584-590.",
  "4. Varga J & Szabo Z. Modified regression model for the Logan plot.",
  "   J Cereb Blood Flow Metab 2002; 22:240-244.",
  " ",
  "See also: imgfur, imgbound, imgunit, imglhki, img2tif, patlak ",
  " ",
  "Keywords: image, MTGA, Patlak plot, modelling, Ki, metabolic rate",
  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           fi, ret;
  int           param_filt=0;
  int           noneg=0; // 0=negative values are kept, 1=negative values set to zero
  char          inpfile[FILENAME_MAX], petfile[FILENAME_MAX];
  char          icfile[FILENAME_MAX], nrfile[FILENAME_MAX], outfile[FILENAME_MAX];
  char          tmp[1024], *cptr;
  DFT           input;
  IMG           img, out, icimg, nrimg;
  double        startTime=-1.0, upperLimit=-1.0, fittime=-1.0;
  float         calcThreshold=0.0;
  int           startSample=0, endSample=0, fitdimt=0;
  linefit_range fit_range;
  double        LC=-1.0, Ca=-1.0, density=-1.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]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) 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, "CA=", 3)==0) {
      Ca=atof_dpi(cptr+3); if(Ca>0.0) continue;
    } else if(strncasecmp(cptr, "LC=", 3)==0) {
      LC=atof_dpi(cptr+3); if(LC>0.0) continue;
    } else if(strncasecmp(cptr, "D=", 2)==0) {
      density=atof_dpi(cptr+2); if(density>0.0) continue;
    } else if(strncasecmp(cptr, "DENSITY=", 8)==0) {
      density=atof_dpi(cptr+8); if(density>0.0) continue;
    } else if(strncasecmp(cptr, "NONEGATIVES", 5)==0) {
      noneg=1; 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 MR will be calculated, set defaults and give warnings as necessary */
  if(Ca>0.0) {
    if(LC<0.0) {
      LC=DEFAULT_LC;
      fprintf(stderr, "Warning: LC not set, using default %g\n", LC);
    }
    if(density<0.0) {
      density=DEFAULT_DENSITY;
      fprintf(stderr, "Warning: tissue density not set, using default %g\n", density);
    } 
  } else { /* Warnings if density or LC set when MR will not be calculated */
    if(LC>0.0) fprintf(stderr, "Warning: LC was set but is not used.\n");
    if(density>0.0) fprintf(stderr, "Warning: tissue density was set but is not used.\n");
  }
 
 
  /* 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);
    printf("icfile := %s\n", icfile);
    printf("nrfile := %s\n", nrfile);
    printf("max := %f\n", upperLimit);
    printf("Ca := %g\n", Ca);
    printf("LC := %g\n", LC);
    printf("fittime := %g\n", fittime);
    printf("density := %g\n", density);
    printf("param_filt := %d\n", param_filt);
    printf("calcThreshold := %g%%\n", 100.*calcThreshold);
    printf("noneg := %d\n", noneg);
    fflush(stdout);
  }
  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); 
#endif
 
  /*
   *  Calculate Gjedde-Patlak 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_patlak(&input, &img, startSample, endSample, fit_range, calcThreshold,
                 &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 Ki values that are higher than the user-defined limit */
  if(upperLimit>0.0) {
    if(verbose>0) fprintf(stdout, "filtering out Ki values exceeding %g\n", upperLimit);
    imgCutoff(&out, (float)upperLimit, 0);
  }
  /* Remove negative Ki values, if requested */
  if(noneg) {
    if(verbose>0) fprintf(stdout, "filtering out negative Ki values\n");
    imgCutoff(&out, 0.0, 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);
  }
 
 
  /*
   *  Calculate metabolic rate, if necessary
   */
  if(Ca>0.0) {
    double MRf;
    MRf=100.*Ca/(density*LC);
    if(verbose>1) fprintf(stdout, "converting Ki to metabolic rate with factor %g\n", MRf);
    ret=imgArithmConst(&out, MRf, '*', 1.0E+06, verbose-6);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot calculate metabolic rate.\n"); 
      imgEmpty(&out); imgEmpty(&icimg); imgEmpty(&nrimg);
      return(7);
    }
    out.unit=CUNIT_UMOL_PER_MIN_PER_100G;
  }
 
  /*
   *  Save parametric Ki 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) {
    if(Ca<=0.0) fprintf(stdout, "Ki image %s saved.\n", outfile);
    else fprintf(stdout, "MR 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);
}
/*****************************************************************************/
 
/*****************************************************************************/