arlkup.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 "tpcift.h"
#include "tpctac.h"
#include "tpcli.h"
#include "tpccm.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculation of perfusion look-up table from [O-15]H2O blood curve",
  "for usage with PET in vivo autoradiographic (ARG) method.",
  " ",
  "Usage: @P [options] BTAC p fMax start dur lkupfile",
  " ",
  "Options:",
  " -static=<y|N>",
  "     The look-up table integral can be corrected for physical decay as",
  "     single frame when static PET scan was performed (y), or dynamically",
  "     (N, the default).",
  " -nr=<value>",
  "     Set the size of look-up table; 5000 by default.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Blood TAC (BTAC) must be calibrated and corrected for decay. Sample times",
  "must be in sec unless correctly specified inside the BTAC file.",
  "Concentrations must be in the same units as the units of PET tissue data.",
  "Partition coefficient of water (p) must be given in units mL/mL,",
  "and maximal perfusion (fMax) in units mL/(min*100mL).",
  "Integration start time and duration is given in sec.",
  " ",
  "Look-up table will contain 2 columns: TTAC integral [sec*kBq/mL] and",
  "blood flow [mL/(min*100mL)].",
  " ",
  "Example:",
  "     @P -nr=2000 s1456_blo.fit 0.8 50 0 120 s1456.lkup",
  " ",
  "References:",
  "1. Raichle ME. Quantitative in vivo autoradiography with positron emission",
  "   tomography. Brain Res Rev. 1979;1:47-68.",
  "2. Herscovitch P, Markham J, Raichle ME. Brain blood flow measured with",
  "   intravenous H215O. I. Theory and error analysis.",
  "   J Nucl Med. 1983;24:782-789.",
  "3. Raichle ME, Martin WRW, Herscovitch P, Mintun MA, Markham J. Brain blood",
  "   flow measured with intravenous H215O. II. Implementation and validation.",
  "   J Nucl Med. 1983;24:790-798.",
  "4. Ruotsalainen U, Raitakari M, Nuutila P, Oikonen V, Sipila H, Teras M,",
  "   Knuuti J, Bloomfield PM, Iida H. Quantitative blood flow measurement of",
  "   skeletal muscle using oxygen-15-water and PET. ",
  "   J Nucl Med. 1997; 38:314-319.",
  " ",
  "See also: fitdelay, imginteg, imglkup, taclkup, tacunit, imgbfbp, imgflow",
  " ",
  "Keywords: perfusion, blood flow, radiowater, autoradiography, ARG, look-up table",
  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    btacfile[FILENAME_MAX], lkupfile[FILENAME_MAX];
  char   *cptr;
  double  pWater, maxFlow, start, dur, end;
  int     tableSize=5000;
  int     singleFrame=0;
  int     ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  btacfile[0]=lkupfile[0]=(char)0;
  /* 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) continue;
    if(strncasecmp(cptr, "NR=", 3)==0) {
      if(atoiCheck(cptr+3, &tableSize)==0 && tableSize>5) continue;
    } else if(strncasecmp(cptr, "STATIC=", 7)==0) {
      cptr+=7;
      if(strncasecmp(cptr, "YES", 1)==0) {singleFrame=1; continue;}
      else if(strncasecmp(cptr, "NO", 1)==0) {singleFrame=0; 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-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(btacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(atofCheck(argv[ai], &pWater) || pWater<=0.0 || pWater>1.25) {
      fprintf(stderr, "Error: invalid partition coefficient '%s'.\n", argv[ai]);
      return(1);
    }
    ai++;
  }
  if(ai<argc) {
    if(atofCheck(argv[ai], &maxFlow) || maxFlow<=0.0) {
      fprintf(stderr, "Error: invalid maximum blood flow '%s'.\n", argv[ai]);
      return(1);
    }
    ai++;
  }
  if(ai<argc) {
    if(atofCheck(argv[ai], &start)) {
      fprintf(stderr, "Error: invalid integration start time '%s'.\n", argv[ai]);
      return(1);
    }
    ai++;
  }
  if(ai<argc) {
    if(atofCheck(argv[ai], &dur) || dur<=0.0) {
      fprintf(stderr, "Error: invalid integration duration '%s'.\n", argv[ai]);
      return(1);
    }
    ai++;
  }
  if(ai<argc) strlcpy(lkupfile, 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(!lkupfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  end=start+dur;
  if(end<=0.0) {
    fprintf(stderr, "Error: invalid integration times.\n");
    return(1);
  }
  if(tableSize<50) {
    fprintf(stderr, 
      "Warning: look-up table may be too small to provide accurate results.\n"); 
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("btacfile := %s\n", btacfile);
    printf("pWater := %g\n", pWater);
    printf("maxFlow := %g\n", maxFlow);
    printf("start := %g\n", start);
    printf("dur := %g\n", dur);
    printf("end := %g\n", end);
    printf("tableSize := %d\n", tableSize);
    printf("lkupfile := %s\n", lkupfile);
    printf("singleFrame := %d\n", singleFrame);
  }
 
 
  /*
   *  Read BTAC
   */
  if(verbose>1) printf("reading %s\n", btacfile);
  TAC btac; tacInit(&btac);
  ret=tacRead(&btac, btacfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&btac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(btac.format));
    printf("tacNr := %d\n", btac.tacNr);
    printf("sampleNr := %d\n", btac.sampleNr);
    printf("xunit := %s\n", unitName(btac.tunit));
    printf("yunit := %s\n", unitName(btac.cunit));
  }
  if(btac.tacNr>1) {
    fprintf(stderr, "Warning: only first TAC in input file is used.\n");
    btac.tacNr=1;
  }
 
  /* Sort by sample times */
  tacSortByTime(&btac, NULL);
 
  /* Check for missing values */
  ret=tacNaNs(&btac);
  if(ret>0) {
    if(verbose>1) printf("missing blood concentrations.\n");
    /* Try to fix missing concentrations */
    ret=tacFixNaNs(&btac);
    if(ret!=0) {
      fprintf(stderr, "Error: missing concentrations in %s.\n", btacfile);
      tacFree(&btac); return(2);
    }
  }
 
 
  /* Check and convert time units */
  double xmin, xmax;
  ret=tacXRange(&btac, &xmin, &xmax);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid time range in btac file.\n");
    tacFree(&btac); return(3);
  }
  if(verbose>2) {
    printf("btac.xmin := %g\n", xmin);
    printf("btac.xmax := %g\n", xmax);
  }
  if(btac.tunit==UNIT_UNKNOWN) {
    if(verbose>0) printf("Note: time units had to be guessed.\n");
    if(xmin<1.0 && xmax<30.0) btac.tunit=UNIT_MIN;
    else btac.tunit=UNIT_SEC;
  }  
  ret=tacXUnitConvert(&btac, UNIT_SEC, &status);
  if(ret!=TPCERROR_OK && verbose>0) {
    fprintf(stderr, "Error: invalid time units.\n");
    tacFree(&btac); return(3);
  }
  /* Check that BTAC extends to the end of requested integration period */
  if(xmax<0.99*end) {
    fprintf(stderr, "Error: blood data only upto %g s.\n", xmax);
    tacFree(&btac); return(3);
  }
  /* Check if integration end time seems to be too short */
  if(end<=xmin) {
    fprintf(stderr, "Error: invalid integration time.\n");
    tacFree(&btac); return(3);
  }
  if(end<(0.25*xmin+0.75*xmax)) {
    fprintf(stderr, "Warning: integration time is short compared to BTAC.\n");
  }
  if(btac.sampleNr<5) {
    fprintf(stderr, "Error: only %d samples in %s\n", btac.sampleNr, btacfile);
    tacFree(&btac); return(2);
  } else if(btac.sampleNr<50 && btac.isframe==0) {
    // no warning if image-derived input containing frame times
    fprintf(stderr, "Warning: small blood sample nr in %s\n", btacfile);
  }
 
 
  /* Check and convert concentration units */
  double ymin, ymax;
  ret=tacYRange(&btac, 0, &ymin, &ymax, NULL, NULL, NULL, NULL);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid concentration range in btac file.\n");
    tacFree(&btac); return(4);
  }
  if(verbose>2) {
    printf("btac.ymin := %g\n", ymin);
    printf("btac.ymax := %g\n", ymax);
  }
  if(btac.cunit==UNIT_UNKNOWN) {
    fprintf(stderr, "Warning: concentration units are not known.\n");
  }
 
  /* Add extra column to BTAC data struct to store simulated tissue TAC */
  ret=tacAllocateMore(&btac, 1);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    tacFree(&btac); return(5);
  }
 
 
  /*
   *  Make the look-up table
   */
 
  /* Allocate memory for the lookup table */
  if(verbose>1) printf("allocating memory for table\n");
  TAC lkup; tacInit(&lkup);
  ret=tacAllocate(&lkup, tableSize, 1);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(ret));
    tacFree(&btac); tacFree(&lkup); return(6);
  }
  lkup.sampleNr=tableSize; lkup.tacNr=1;
 
  /* Set look-up table units if possible */
  lkup.cunit=UNIT_ML_PER_DL_MIN;
  if(btac.cunit==UNIT_KBQ_PER_ML) lkup.tunit=UNIT_SEC_KBQ_PER_ML;
  else if(btac.cunit==UNIT_BQ_PER_ML) lkup.tunit=UNIT_SEC_BQ_PER_ML;
  else lkup.tunit=UNIT_UNKNOWN;
 
  /* Calculate the K1 (blood flow) in the look-up table */
  if(verbose>1) printf("setting the perfusion in look-up table\n");
  {
    double K1step=maxFlow/(double)(tableSize-1);
    if(verbose>2) printf("K1step := %E\n", K1step);
    lkup.c[0].y[0]=0.0;
    for(int i=1; i<tableSize; i++) lkup.c[0].y[i]=lkup.c[0].y[i-1]+K1step;
    if(verbose>2) printf("K1 range: %g - %g\n", 
                         lkup.c[0].y[0], lkup.c[0].y[lkup.sampleNr-1]);
  }
 
  /* Calculate the table */
  if(verbose>1) printf("calculating the tissue integrals in look-up table\n");
  {
    double dc=1.0;
    if(singleFrame!=0) {
      /* If PET data is collected as a single frame (static scan), then 
         calculate decay correction factor during the frame */
      dc=decayCorrectionFactorFromIsotope(ISOTOPE_O_15, start/60.0, dur/60.0);
    }
 
    double K1, k2, intx[2], inty[2];
    intx[0]=start; intx[1]=end;
    lkup.x[0]=0.0; 
    ret=0;
    for(int i=1; i<tableSize && !ret; i++) {
      K1=lkup.c[0].y[i]/6000.0; // mL/(min*100mL) -> mL/(sec*mL)
      k2=K1/pWater;
      /* Simulate the tissue curve with these K1 and k2 values */
      ret=simC1(btac.x, btac.c[0].y, btac.sampleNr, K1, k2, btac.c[1].y);
      if(ret) break;
      /* Remove decay correction from tissue curve, if static PET scan */
      if(singleFrame!=0) {
        for(int j=0; j<btac.sampleNr; j++)
          btac.c[1].y[j]/=decayCorrectionFactorFromIsotope(
                                            ISOTOPE_O_15, btac.x[j]/60.0, 0.0);
      }
      /* Integrate the simulated tissue curve */
      ret=liInterpolate(btac.x, btac.c[1].y, btac.sampleNr,
                        intx, NULL, inty, NULL, 2, 4, 1, 0);
      if(ret) break;
      /* Set the table integral values */
      lkup.x[i]=inty[1]-inty[0];
      /* Correct for physical decay again, as single frame, if required */
      if(singleFrame!=0) lkup.x[i]*=dc;
    } // next look-up table row
    if(ret) {
      fprintf(stderr, "Error: cannot compute simulated tissue AUC.\n");
      tacFree(&btac); tacFree(&lkup); return(7);
    }
  }
 
  /* BTAC is not needed anymore */
  tacFree(&btac);  
 
 
  /*
   *  Write look-up table
   */
  if(verbose>1) printf("writing %s\n", lkupfile);
  FILE *fp; fp=fopen(lkupfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", lkupfile);
    tacFree(&lkup); return(11);
  }
  ret=tacWrite(&lkup, fp, TAC_FORMAT_PMOD, 0, &status);
  fclose(fp); tacFree(&lkup);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    return(12);
  }
  if(verbose>0) printf("Look-up table saved in %s.\n", lkupfile);
 
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/