imgflowd.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 "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodel.h"
#include "libtpcmodext.h"
/*****************************************************************************/
#define NNLS_N 3
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Estimation of rate constants K1 (perfusion), k2, and Va from dynamic",
  "radiowater PET image in ECAT 6.3, ECAT 7.x, NIfTI-1, or Analyze 7.5 file",
  "format using the two-compartment model equations",
  "  dCt(t)/dt = K1*Ca(t) - k2*Ct(t)",
  "  Croi(t) = Va*Ca(t) + Ct(t)",
  "are linearized (1) into equation",
  "  Croi(t) = Va*Ca(t) + (K1+Va*k2)*Integral[Ca(t)] - k2*Integral[Croi(t)]", 
  "from which the model parameters are solved using Lawson-Hanson non-negative",
  "least squares (NNLS) method (2).",
  " ",
  "Arterial blood curve (BTAC) and dynamic PET image must be corrected for physical decay.",
  "BTAC sample times must be in seconds.",
  " ",
  "Since time delay varies inside PET image, this program requires a delay map",
  "which contains the time delay for each image pixel.",
  " ",
  "Usage: @P [Options] btacfile imgfile delayfile flowfile",
  " ",
  "Options:",
  " -end=<Fit end time (sec)>",
  "     Use data from 0 to end time; by default, all of it.",
  " -k2=<filename>",
  "     Parametric k2 image is saved; in some situations perfusion calculation",
  "     from k2 can be more accurate than the default assumption of f=K1.",
  "     Perfusion can be calculated from k2 using equation f=k2*pH2O, where",
  "     pH2O is the physiological partition coefficient of water in tissue.",
  " -Va=<filename> | -Va=0",
  "     Parametric Va image is saved, or set Va to 0 if image is pre-corrected",
  "     for arterial blood volume; by default Va is fitted.",
  " -max=<Max value>",
  "     Upper limit for blood flow (K1) values.",
  " -noneg",
  "     Pixels with negative K1 estimates are set to zero.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "The units of pixel values in the blood flow (K1) image is",
  "(mL blood)/((mL tissue) * min), in k2 image 1/min, and",
  "in Va image (mL blood/mL tissue).",
  " ",
  "References:",
  "1. Blomqvist G. On the construction of functional maps in positron",
  "   emission tomography. J Cereb Blood Flow Metab. 1984;4:629-632.",
  "2. Lawson CL & Hanson RJ. Solving least squares problems.",
  "   Prentice-Hall, 1974.",
  " ",
  "See also: imgdelay, imgflowm, imgbfh2o, imgcbv",
  " ",
  "Keywords: image, modelling, perfusion, blood flow, radiowater, NNLS",
  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;
/*****************************************************************************/
 
/*****************************************************************************/
/*
 *  main()
 */
int main(int argc, char *argv[])
{
  int     ai, help=0, version=0, verbose=1;
  int     weight=0, dataNr=0;
  char    inpfile[FILENAME_MAX], petfile[FILENAME_MAX], delayfile[FILENAME_MAX];
  char    k2file[FILENAME_MAX], vafile[FILENAME_MAX], k1file[FILENAME_MAX];
  double  upperLimit=nan("");
  double  lowerLimit=nan("");
  double  fittime=nan("");
  int     fitVa=1;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=petfile[0]=k1file[0]=k2file[0]=vafile[0]=delayfile[0]=(char)0;
  /* Get 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==NULL) continue;
    if(strcasecmp(cptr, "W")==0) {
      weight=1; continue;
    } else if(strncasecmp(cptr, "k2=", 3)==0) {
      strlcpy(k2file, cptr+3, FILENAME_MAX); if(strlen(k2file)) continue;
    } else if(strcasecmp(cptr, "VA=0")==0) {
      fitVa=0; continue;
    } else if(strncasecmp(cptr, "VA=", 3)==0 || strncasecmp(cptr, "VB=", 3)==0) {
      strlcpy(vafile, cptr+3, FILENAME_MAX); if(strlen(vafile)) continue;
    } else if(strncasecmp(cptr, "noneg", 2)==0) {
      lowerLimit=0.0; continue;
    } else if(strncasecmp(cptr, "MAX=", 4)==0) {
      upperLimit=atof_dpi(cptr+4); if(upperLimit>0.0) continue;
    } else if(strncasecmp(cptr, "END=", 4)==0) {
      fittime=atof_dpi(cptr+4)/60.; if(fittime>0.0) continue;
    }
    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) strlcpy(delayfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(k1file, 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(!k1file[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(!(fittime>0.0)) fittime=1.0E+020;
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("inpfile := %s\n", inpfile);
    printf("petfile := %s\n", petfile);
    printf("delayfile := %s\n", delayfile);
    printf("k1file := %s\n", k1file);
    if(k2file[0]) printf("k2file := %s\n", k2file);
    printf("fitVa := %d\n", fitVa);
    if(vafile[0]) printf("vafile := %s\n", vafile);
    if(!isnan(upperLimit)) printf("upperLimit := %f\n", upperLimit);
    if(!isnan(lowerLimit)) printf("lowerLimit := %f\n", lowerLimit);
    if(isfinite(fittime)) printf("fittime := %g [sec]\n", 60.0*fittime);
    printf("weight := %d\n", weight);
    fflush(stdout);
  }
  if(verbose>10) IMG_TEST=verbose-10; else IMG_TEST=0;
 
 
 
  /*
   *  Read PET image and input TAC
   */
  if(verbose>1) printf("reading data files\n");
  DFT inp; dftInit(&inp);
  DFT tac; dftInit(&tac);
  IMG img; imgInit(&img);
  {
    char tmp[FILENAME_MAX+1];
    int ret=imgReadModelingData(
      petfile, NULL, inpfile, NULL, NULL, &fittime, &dataNr, &img,
      &inp, &tac, 1, stdout, verbose-2, tmp);
    if(ret!=0) {
      fprintf(stderr, "Error: %s.\n", tmp);
      if(verbose>1) printf("  ret := %d\n", ret);
      return(2);
    }
    if(imgNaNs(&img, 1)>0)
      if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  }
  /* Let time units be in sec because delay map is in sec */
  if(verbose>1) {
    printf("fittimeFinal := %g s\n", 60.0*fittime);
    printf("dataNr := %d\n", dataNr);
  }
  /* Check that image is dynamic */
  if(dataNr<3) {
    fprintf(stderr, "Error: too few time frames for fitting.\n");
    if(verbose>0) imgInfo(&img);
    imgEmpty(&img); dftEmpty(&tac); dftEmpty(&inp); return(2);
  }
 
  /*
   *  Read delay map
   */
  if(verbose>1) printf("reading delay map\n");
  IMG dtmap; imgInit(&dtmap);
  {
    int ret=imgRead(delayfile, &dtmap);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot read delay map %s.\n", delayfile);
      if(verbose>1) printf("  ret := %d\n", ret);
      imgEmpty(&img); dftEmpty(&tac); dftEmpty(&inp); return(2);
    }
    if(imgMatchMatrixSize(&img, &dtmap)) { 
      fprintf(stderr, "Error: bad delay map size.\n");
      imgEmpty(&img); dftEmpty(&tac); dftEmpty(&inp); imgEmpty(&dtmap); return(2);
    }
  }
 
 
  /* Allocate memory for tissue TAC and integral */
  {
    int ret=dftAddmem(&tac, 1);
    if(ret!=0) {
      fprintf(stderr, "Error (%d) in allocating memory.\n", ret);
      imgEmpty(&img); dftEmpty(&tac); dftEmpty(&inp); imgEmpty(&dtmap);
      return(3);
    }
    strcpy(tac.voi[0].voiname, "input");
    strcpy(tac.voi[1].voiname, "tissue");
  }
 
 
 
  /*
   *  Allocate memory for result images and fill the header info
   */
  IMG k1img, k2img, vaimg;
  imgInit(&k1img); imgInit(&k2img); imgInit(&vaimg);
  {
    int ret=0;
    if((ret=imgAllocateWithHeader(&k1img, img.dimz, img.dimy, img.dimx, 1, &img))) {
      fprintf(stderr, "Error: cannot allocate memory.\n");
    } else {
      for(int zi=0; zi<img.dimz; zi++)
        for(int yi=0; yi<img.dimy; yi++)
          for(int xi=0; xi<img.dimx; xi++)
            k1img.m[zi][yi][xi][0]=0.0;
      k1img.unit=IMGUNIT_ML_PER_ML_PER_MIN;
      k1img.decayCorrection=IMG_DC_NONCORRECTED; k1img.isWeight=0;
      //k1img.start[0]=img.start[0]; k1img.end[0]=img.end[dataNr-1];
      k1img.start[0]=k1img.end[0]=0.0;
    }
    /* optional maps */
    if(ret==0 && k2file[0]) {
      if((ret=imgAllocateWithHeader(&k2img, img.dimz, img.dimy, img.dimx, 1, &k1img))) {
        fprintf(stderr, "Error: cannot allocate memory.\n");
      } else {
        for(int zi=0; zi<img.dimz; zi++)
          for(int yi=0; yi<img.dimy; yi++)
            for(int xi=0; xi<img.dimx; xi++)
              k2img.m[zi][yi][xi][0]=0.0;
        k2img.unit=IMGUNIT_PER_MIN;
      }
    }
    if(ret==0 && vafile[0]) {
      if((ret=imgAllocateWithHeader(&vaimg, img.dimz, img.dimy, img.dimx, 1, &k1img))) {
        fprintf(stderr, "Error: cannot allocate memory.\n");
      } else {
        for(int zi=0; zi<img.dimz; zi++)
          for(int yi=0; yi<img.dimy; yi++)
            for(int xi=0; xi<img.dimx; xi++)
              vaimg.m[zi][yi][xi][0]=0.0;
        vaimg.unit=IMGUNIT_ML_PER_ML;
      }
    }
    /* quit if error */
    if(ret) {
      imgEmpty(&img); dftEmpty(&tac); dftEmpty(&inp); imgEmpty(&dtmap);
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vaimg);
      return(4);
    }
  }
 
 
  /*
   *  Allocate memory required by NNLS
   */
  int     nnls_n, nnls_m, nnls_index[NNLS_N];
  double *nnls_a[NNLS_N], *nnls_b, *nnls_zz, nnls_x[NNLS_N], *nnls_mat,
          nnls_wp[NNLS_N], *dptr, nnls_rnorm;
  if(verbose>1) printf("allocating memory for NNLS\n");
  nnls_n=NNLS_N; nnls_m=dataNr;
  nnls_mat=(double*)malloc(((nnls_n+2)*nnls_m)*sizeof(double));
  if(nnls_mat==NULL) {
    fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
    imgEmpty(&img); dftEmpty(&tac); dftEmpty(&inp); imgEmpty(&dtmap);
    imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vaimg);
    return(5);
  }
  dptr=nnls_mat;
  for(int n=0; n<nnls_n; n++) {nnls_a[n]=dptr; dptr+=nnls_m;}
  nnls_b=dptr; dptr+=nnls_m; nnls_zz=dptr;
 
  /* Copy weights if available */
  /* or set them to frame lengths */
  if(verbose>2) printf("working with NNLS weights\n");
  if(weight==1 && img.isWeight==0) {
    for(int m=0; m<nnls_m; m++) img.weight[m]=img.end[m]-img.start[m];
    img.isWeight=1;
  }
  /* Compute NNLS weights */
  if(img.isWeight) {
    for(int m=0; m<nnls_m; m++) {
      tac.w[m]=img.weight[m]; if(tac.w[m]<=1.0e-20) tac.w[m]=0.0;
    }
  }
 
 
  /*
   *  Compute pixel-by-pixel
   */
  long long okNr=0;
  double *cp, *cpi, *ct, *cti;
  cp=tac.voi[0].y; cpi=tac.voi[0].y2;
  ct=tac.voi[1].y; cti=tac.voi[1].y2;
  if(verbose>0) fprintf(stdout, "computing pixel-by-pixel\n");
  for(int zi=0; zi<img.dimz; zi++) {
    if(verbose>2) printf("computing plane %d\n", img.planeNumber[zi]);
    else if(img.dimz>1 && verbose>0) {fprintf(stdout, "."); fflush(stdout);}
    for(int yi=0; yi<img.dimy; yi++) {
      for(int xi=0; xi<img.dimx; xi++) {
 
        /* Initiate pixel output values */
        k1img.m[zi][yi][xi][0]=0.0;
        if(k2file[0]) k2img.m[zi][yi][xi][0]=0.0;
        if(vafile[0]) vaimg.m[zi][yi][xi][0]=0.0;
 
        /* Check the delay map */
        if(!isfinite(dtmap.m[zi][yi][xi][0])) continue;
 
        /* Copy and integrate pixel curve */
        for(int m=0; m<nnls_m; m++) ct[m]=img.m[zi][yi][xi][m];
        petintegral(tac.x1, tac.x2, ct, nnls_m, cti, NULL);
        /* if AUC at the end is not above zero, then forget this pixel */
        if(!(cti[nnls_m-1]>0.0)) continue;
 
        /* Calculate delay-corrected BTAC and its AUC */
        double dx[inp.frameNr];
        for(int i=0; i<inp.frameNr; i++) dx[i]=inp.x[i]+dtmap.m[zi][yi][xi][0];
        int ret=interpolate4pet(dx, inp.voi[0].y, inp.frameNr, tac.x1, tac.x2, cp, cpi, NULL, dataNr);
        if(ret) {
          if(verbose>3)
            fprintf(stderr, "Error: cannot interpolate/integrate input for pixel %d,%d,%d\n", zi, yi, xi);
          continue;
        }
 
        /*
         *  Estimate K1, k2 and Va
         */
        nnls_m=dataNr; nnls_n=NNLS_N; if(fitVa==0) nnls_n--;
        /* Fill NNLS A matrix: */
        /* function #1: tissue integral x -1 */
        for(int m=0; m<nnls_m; m++) nnls_a[0][m]=-cti[m];
        /* function #2: integral of input */
        for(int m=0; m<nnls_m; m++) nnls_a[1][m]=cpi[m];
        /* function #3: input curve */
        if(nnls_n>2) for(int m=0; m<nnls_m; m++) nnls_a[2][m]=cp[m];
        /* Fill NNLS B array: tissue */
        for(int m=0; m<nnls_m; m++) nnls_b[m]=ct[m];
        /* Apply data weights */
        if(img.isWeight) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, tac.w);
        /* NNLS */
        if(nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm,
                nnls_wp, nnls_zz, nnls_index)>1) continue; /* no solution is possible */
        /* Set pixel values */
        double Va=0.0;
        if(nnls_n>2) Va=nnls_x[2];
        k1img.m[zi][yi][xi][0]=nnls_x[1];
        if(nnls_n>2) k1img.m[zi][yi][xi][0]-=Va*nnls_x[0];
        if(k2file[0]) k2img.m[zi][yi][xi][0]=nnls_x[0];
        if(vafile[0]) vaimg.m[zi][yi][xi][0]=Va;
        /* Convert from per sec to per min */
        k1img.m[zi][yi][xi][0]*=60.0;
        if(k2file[0]) k2img.m[zi][yi][xi][0]*=60.0;
 
        /* Enforce K1 limits, if requested */
        if(!isnan(upperLimit) && k1img.m[zi][yi][xi][0]>upperLimit)
          k1img.m[zi][yi][xi][0]=upperLimit;
        if(!isnan(lowerLimit) && k1img.m[zi][yi][xi][0]<lowerLimit)
          k1img.m[zi][yi][xi][0]=lowerLimit;
 
        /* Va limit */
        if(vafile[0] && vaimg.m[zi][yi][xi][0]>1.0) vaimg.m[zi][yi][xi][0]=1.0;
 
        okNr++;
      } /* next column */
    } /* next row */
  } /* next plane */
  if(verbose>0) fprintf(stdout, "done.\n");
 
  /* No need for dynamic image, NNLS matrix, delay map, or curves any more */
  free(nnls_mat); imgEmpty(&img); imgEmpty(&dtmap); dftEmpty(&tac); dftEmpty(&inp);
 
  /*
   *  Check that not all pixels failed
   */
  if(verbose>1) printf("calculation successful for %lld pixels.\n", okNr);
  if(okNr==0) {
    fprintf(stderr, "Error: failed calculation.\n");
    imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vaimg);
    return(9);
  }
 
 
  /*
   *  Save parametric image(s)
   */
  if(verbose>0) printf("writing parametric images\n");
  {
    /* K1 */
    int ret=imgWrite(k1file, &k1img);
    if(ret) {
      fprintf(stderr, "Error: %s\n", k1img.statmsg);
      if(verbose>1) printf("ret := %d\n", ret);
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vaimg);
      return(13);
    }
    if(verbose>0) fprintf(stdout, "Flow image %s saved.\n", k1file);
    /* k2 */
    if(k2file[0]) {
      ret=imgWrite(k2file, &k2img);
      if(ret) {
        fprintf(stderr, "Error: %s\n", k2img.statmsg);
        imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vaimg);
        return(14);
      }
      if(verbose>0) fprintf(stdout, "k2 image %s saved.\n", k2file);
    }
    /* Va */
    if(vafile[0]) {
      ret=imgWrite(vafile, &vaimg);
      if(ret) {
        fprintf(stderr, "Error: %s\n", vaimg.statmsg);
        imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vaimg);
        return(16);
      }
      if(verbose>0) fprintf(stdout, "Va image %s saved.\n", vafile);
    }
  }
 
  if(verbose>2) printf("before quitting, free memory\n");
  imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vaimg);
  fflush(stdout);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/