imgflowm.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) must be corrected for decay and time delay,",
  "with sample times in seconds. Dynamic PET image must be corrected for decay.",
  " ",
  "Since time delay varies inside PET image, this program requires a mask image",
  "which specifies to which pixels the model is applied. If output file(s) do",
  "not exist, the output files are created with zero pixel values where mask",
  "image has value 0. If output file(s) exist, the pre-existing pixel values",
  "outside the mask are not changed.",
  " ",
  "Usage: @P [Options] btacfile imgfile maskfile 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).",
  " ",
  "Example 1. Calculation of perfusion and arterial blood volume image,",
  "           stopping fit at 180 s:",
  "  @P -Va=p234va.v -end=180 p234d1.bld p234.nii p234_mask1.nii p234flow.nii",
  " ",
  "Example 2. Dynamic image is pre-corrected for vascular activity,",
  "           and full data length is used:",
  "  @P -Va=0 p345dc2.bld p345_vacorr.nii p345flow.nii",
  " ",
  "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.",
  "3. Ohta S, Meyer E, Thompson CJ, Gjedde A. Oxygen consumption of the",
  "   living human brain measured after a single inhalation of positron",
  "   emitting oxygen. J Cereb Blood Flow Metab. 1992;12:179-192.",
  " ",
  "See also: imgflow, imgdelay, fitdelay, imgcbv, fit_h2o, imgcalc",
  " ",
  "Keywords: image, modelling, perfusion, blood flow, radiowater, mask",
  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], maskfile[FILENAME_MAX];
  char    k2file[FILENAME_MAX], vafile[FILENAME_MAX], k1file[FILENAME_MAX];
  double  upperLimit=nan("");
  double  lowerLimit=nan("");
  double  fittime=nan("");
  int     fitVa=1;
  double *cp, *cpi, *ct, *cti, Va;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=petfile[0]=k1file[0]=k2file[0]=vafile[0]=maskfile[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(maskfile, 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("maskfile := %s\n", maskfile);
    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 tac; dftInit(&tac);
  IMG img; imgInit(&img);
  {
    char tmp[FILENAME_MAX+1];
    int ret=imgReadModelingData(
      petfile, NULL, inpfile, NULL, NULL, &fittime, &dataNr, &img,
      NULL, &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");
  }
  /* Set time unit to min, also for integrals in y2[] */
  if(tac.timeunit==TUNIT_SEC) for(int i=0; i<tac.frameNr; i++) tac.voi[0].y2[i]/=60.0;
  dftTimeunitConversion(&tac, TUNIT_MIN);
  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); return(2);
  }
 
  /*
   *  Read mask image
   */
  if(verbose>1) printf("reading mask image\n");
  IMG mask; imgInit(&mask);
  {
    int ret=imgRead(maskfile, &mask);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot read mask %s.\n", maskfile);
      if(verbose>1) printf("  ret := %d\n", ret);
      imgEmpty(&img); dftEmpty(&tac); return(2);
    }
    if(imgMatchMatrixSize(&img, &mask)) { 
      fprintf(stderr, "Error: bad mask for image.\n");
      imgEmpty(&img); dftEmpty(&tac); imgEmpty(&mask); 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); imgEmpty(&mask);
      return(3);
    }
    strcpy(tac.voi[0].voiname, "input");
    strcpy(tac.voi[1].voiname, "tissue");
  }
 
 
 
  /*
   *  Read pre-existing result images, or if not existing,
   *  allocate result images and fill the header info
   */
  IMG k1img, k2img, vaimg;
  imgInit(&k1img); imgInit(&k2img); imgInit(&vaimg);
  {
    int ret=0;
    if(imgRead(k1file, &k1img)==0) {
      if((ret=imgMatchMatrixSize(&img, &k1img)))
        fprintf(stderr, "Error: wrong matrix size in existing K1 map.\n");
    } else {
      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];
      }
    }
    /* optional maps */
    if(ret==0 && k2file[0]) {
      if(imgRead(k2file, &k2img)==0) {
        if((ret=imgMatchMatrixSize(&img, &k2img)))
          fprintf(stderr, "Error: wrong matrix size in existing k2 map.\n");
      } else {
        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(imgRead(vafile, &vaimg)==0) {
        if((ret=imgMatchMatrixSize(&img, &vaimg)))
          fprintf(stderr, "Error: wrong matrix size in existing Va map.\n");
      } else {
        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); imgEmpty(&mask);
      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); imgEmpty(&mask);
    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
   */
  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++) {
 
        /* Check the mask */
        if(!(fabs(mask.m[zi][yi][xi][0])>1.0E-12)) continue;
 
        /* 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;
 
        /*
         *  Estimate K1, k2 and Va
         */
        nnls_m=dataNr; nnls_n=NNLS_N; if(fitVa==0) nnls_n--;
        /* 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 <= zero, then forget this pixel */
        if(cti[nnls_m-1]<=0.0) continue;
        /* 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 */
        if(nnls_n>2) Va=nnls_x[2]; else Va=0.0;
        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;
 
        /* 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;
 
      } /* next column */
    } /* next row */
  } /* next plane */
  if(verbose>0) fprintf(stdout, "done.\n");
 
  /* No need for dynamic image, NNLS matrix or curves any more */
  free(nnls_mat); imgEmpty(&img); dftEmpty(&tac);
 
 
  /*
   *  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); imgEmpty(&mask);
      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); imgEmpty(&mask);
        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); imgEmpty(&mask);
        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); imgEmpty(&mask);
  fflush(stdout);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/