imglhk1.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 "libtpccurveio.h"
#include "libtpcmodext.h"
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
/*****************************************************************************/
#define NNLS_N 3
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Computation of parametric images of K1 and Vb, and optionally k2, from",
  "dynamic PET image in ECAT, NIfTI, or Analyze format applying",
  "(ir)reversible on1-tissue compartmental model with arterial plasma input.",
  "The compartmental models are transformed to general linear least squares",
  "functions (1), which are solved using Lawson-Hanson non-negative least",
  "squares (NNLS) algorithm (2).",
  " ",
  "Dynamic PET image and plasma time-activity curve (PTAC) must be corrected",
  "for decay to the tracer injection time.",
  " ",
  "Usage: @P [Options] ptacfile imgfile k1file vbfile [k2file]",
  " ",
  "Options:",
  " -thr=<threshold%>",
  "     Pixels with AUC less than (threshold/100 x PTAC AUC) are set to zero;",
  "     default is 1%.",
  " -end=<Fit end time (min)>",
  "     Use data from 0 to end time; by default, model is fitted to all frames.",
//" -NNLS or -BVLS",
//"     LLSQ method.",
//" -w1, -wf, -wfa",
//"     By default, all weights are set to 1.0 (no weighting, option -w1); option -wf",
//"     sets weights based on frame lengths, and option -wfa based on both frame lengths",
//"     and mean activity during each frame.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "The unidirectional back-transport rate k2 is considered in the model setting",
  "only if file name for k2 image is given.",
  " ",
  "Note that this model can correctly estimate Vb only if",
  "1) plasma does not contain any labelled metabolites, and",
  "2) plasma and blood curves are similar in shape.",
  "Note also that Cpet is modelled as Cpet=Vb*Cb + Ct, thus K1 may need to",
  "be corrected by factor 1/(1-Vb).",
  " ",
  "The units of pixel values in the parametric images are ml/(min*ml) for K1",
  "1/min for k2, and ml/ml for Vb.",
  " ",
  "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, ISBN 0-89871-356-0.",
  " ",
  "See also: imglhk3, imgcbv, imgki, fitk3",
  " ",
  "Keywords: image, modelling, irreversible uptake, 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;
/*****************************************************************************/
 
/*****************************************************************************/
enum {METHOD_UNKNOWN, METHOD_NNLS, METHOD_BVLS};
static char *method_str[] = {"unknown", "NNLS", "BVLS", 0};
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char     inpfile[FILENAME_MAX], petfile[FILENAME_MAX], k1file[FILENAME_MAX];
  char     k2file[FILENAME_MAX], vbfile[FILENAME_MAX];
  float    calcThreshold=0.01;
  double   fittime=-1.0;
  int      weights=2; // 1=frame lengths and activity, 2=frame lengths, 2=no weighting
  int      method=METHOD_NNLS;
  int      ret;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=petfile[0]=k1file[0]=k2file[0]=vbfile[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(strncasecmp(cptr, "THR=", 4)==0) {
      double v; ret=atof_with_check(cptr+4, &v);
      if(!ret && v>=0.0 && v<=200.0) {calcThreshold=(float)(0.01*v); continue;}
    } else if(strncasecmp(cptr, "END=", 4)==0) {
      ret=atof_with_check(cptr+4, &fittime); if(!ret && fittime>0.0) continue;
    } else if(strcasecmp(cptr, "WFA")==0) {
      weights=0; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "W1")==0) {
      weights=2; continue;
    } else if(strcasecmp(cptr, "NNLS")==0) {
      method=METHOD_NNLS; continue;
    } else if(strcasecmp(cptr, "BVLS")==0) {
      method=METHOD_BVLS; 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(k1file, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(vbfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(k2file, 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(!vbfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n"); return(1);}
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("inpfile := %s\n", inpfile);
    printf("petfile := %s\n", petfile);
    printf("k1file := %s\n", k1file);
    printf("vbfile := %s\n", vbfile);
    if(k2file[0]) printf("k2file := %s\n", k2file);
    printf("calcThreshold :=%g\n", calcThreshold);
    printf("weights := %d\n", weights);
    printf("method := %s\n", method_str[method]);
    if(fittime>0.0) printf("required_fittime := %g min\n", fittime);
  }
  if(verbose>8) IMG_TEST=verbose-8; else IMG_TEST=0;
  if(verbose>20) ECAT63_TEST=ECAT7_TEST=verbose-20; else ECAT63_TEST=ECAT7_TEST=0;
 
 
  /*
   *  Read PET image and input TAC
   */
  if(verbose>0) printf("reading data files\n");
  DFT tac; dftInit(&tac); 
  IMG img; imgInit(&img);
  int dataNr=0;
  char errmsg[512];
  ret=imgReadModelingData(
    petfile, NULL, inpfile, NULL, NULL, &fittime, &dataNr, &img,
    NULL, &tac, 1, stdout, verbose-2, errmsg);
  if(ret!=0) {
    fprintf(stderr, "Error: %s.\n", errmsg);
    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 fi=0; fi<tac.frameNr; fi++) tac.voi[0].y2[fi]/=60.0;
    for(int fi=0; fi<tac.frameNr; fi++) tac.voi[0].y3[fi]/=3600.0;
  }
  ret=dftTimeunitConversion(&tac, TUNIT_MIN);
  if(verbose>1) {
    printf("fittimeFinal := %g min\n", fittime);
    printf("dataNr := %d\n", dataNr);
  }
  /* Check that image is dynamic */
  if(dataNr<4) {
    fprintf(stderr, "Error: too few time frames for fitting.\n");
    if(verbose>1) imgInfo(&img);
    imgEmpty(&img); dftEmpty(&tac); return(2);
  }
 
 
  /* Add place for tissue TACs too */
  if(verbose>1) fprintf(stdout, "allocating working memory for pixel TACs\n");
  ret=dftAddmem(&tac, 1);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    if(verbose>0) printf("ret := %d\n", ret);
    imgEmpty(&img); dftEmpty(&tac); return(3);
  }
  strcpy(tac.voi[0].name, "input");
  strcpy(tac.voi[1].name, "tissue");
 
  /* Set weights as requested */
  if(imgSetWeights(&img, weights, verbose-5)!=0) {
    fprintf(stderr, "Error: cannot calculate weights.\n");
    imgEmpty(&img); dftEmpty(&tac); return(3);
  }
  for(int i=0; i<dataNr; i++) tac.w[i]=img.weight[i];
 
  /* Determine the threshold */
  double threshold=calcThreshold*tac.voi[0].y2[dataNr-1];
  if(verbose>1) printf("threshold_AUC := %g\n", threshold);
 
 
  /*
   *  Allocate result images (allocate all, even if user did not want to save those)
   */
  if(verbose>1) fprintf(stdout, "allocating memory for parametric image data\n");
  IMG k1img; imgInit(&k1img);
  IMG k2img; imgInit(&k2img);
  IMG vbimg; imgInit(&vbimg);
  ret=imgAllocateWithHeader(&k1img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&k2img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&vbimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory for result image.\n");
    imgEmpty(&img); dftEmpty(&tac); imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vbimg);
    return(4);
  }
  /* set 'frame time' for parametric images */
  k1img.start[0]=k2img.start[0]=vbimg.start[0]=0.0;
  k1img.end[0]=k2img.end[0]=vbimg.end[0]=60.*fittime;
  /* set units in parametric images */
  k1img.unit=CUNIT_ML_PER_ML_PER_MIN;
  vbimg.unit=CUNIT_ML_PER_ML;
  k2img.unit=CUNIT_PER_MIN;
  /* and set the more or less necessary things */
  k1img.decayCorrection=k2img.decayCorrection=vbimg.decayCorrection=IMG_DC_NONCORRECTED;
  k1img.isWeight=k2img.isWeight=vbimg.isWeight=0;
 
 
  /* Fitting */
 
  int fittedNr=0, fittedokNr=0, thresholdNr=0;
 
  if(method==METHOD_NNLS) {
 
    /* Compute weights for NNLS */
    for(int m=0; m<dataNr; m++) {
      if(tac.w[m]<=1.0e-20) tac.w[m]=0.0; else tac.w[m]=sqrt(tac.w[m]);
    }
 
    /*
     *  Allocate memory required by NNLS
     */
    if(verbose>1) fprintf(stdout, "allocating memory for NNLS\n");
    int nnls_n=NNLS_N; if(!k2file[0]) nnls_n--;
    int nnls_m=dataNr;
    int      n, m, nnls_index[nnls_n];
    double  **nnls_a, nnls_b[nnls_m], nnls_zz[nnls_m], nnls_x[nnls_n], *nnls_mat,
             nnls_wp[nnls_n], nnls_rnorm;
    nnls_mat=(double*)malloc((nnls_n*nnls_m)*sizeof(double));
    nnls_a=(double**)malloc(nnls_n*sizeof(double*));
    if(nnls_mat==NULL || nnls_a==NULL) {
      fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
      imgEmpty(&img); dftEmpty(&tac); imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vbimg);
      return(4);
    }
    for(n=0; n<nnls_n; n++) nnls_a[n]=nnls_mat+n*nnls_m;
 
 
    /*
     *  Compute pixel-by-pixel
     */
    if(verbose>0) fprintf(stdout, "computing NNLS pixel-by-pixel\n");
    double *ct, *cti, *cp, *cpi;
    cp=tac.voi[0].y; cpi=tac.voi[0].y2;
    ct=tac.voi[1].y; cti=tac.voi[1].y2;
    for(int pi=0; pi<img.dimz; pi++) {
      if(verbose>0) {fprintf(stdout, "."); fflush(stdout);}
      for(int yi=0; yi<img.dimy; yi++) {
        for(int xi=0; xi<img.dimx; xi++) {
          double K1=0.0, k2=0.0, Vb=0.0;
          /* Set pixel results to zero */
          k1img.m[pi][yi][xi][0]=0.0; 
          k2img.m[pi][yi][xi][0]=0.0; 
          vbimg.m[pi][yi][xi][0]=0.0; 
          /* Copy and integrate pixel curve */
          for(int fi=0; fi<tac.frameNr; fi++) ct[fi]=img.m[pi][yi][xi][fi];
          ret=petintegral(tac.x1, tac.x2, ct, tac.frameNr, cti, NULL);
          if(ret) continue;
          /* if AUC at the end is less than threshold value, then do nothing more */
          if(cti[dataNr-1]<threshold) {thresholdNr++; continue;}
 
          /*
           *  Estimate parameters K1, Vb, and possibly k2
           */
 
          /* Fill the NNLS data matrix */
          for(int m=0; m<nnls_m; m++) {
            nnls_a[0][m]=cpi[m];
            nnls_a[1][m]=cp[m];
            if(k2file[0]) nnls_a[2][m]=-cti[m];
            nnls_b[m]=ct[m];
          }
          if(img.isWeight)
            for(int m=0; m<nnls_m; m++) {
              nnls_b[m]*=tac.w[m];
              for(int n=0; n<nnls_n; n++) nnls_a[n][m]*=tac.w[m];
            }
          if(verbose>5 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
            printf("Matrix A                                 Array B\n");
            for(m=0; m<nnls_m; m++) {
              printf("%12.3f %12.3f %12.3f     %12.3f\n",
                nnls_a[0][m], nnls_a[1][m], nnls_a[2][m], nnls_b[m]);
            }
          }
          /* NNLS */
          ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm, nnls_wp, nnls_zz, nnls_index);
          if(ret!=0) { /* no solution is possible */
            if(verbose>3) printf("no solution possible (%d)\n", ret);
            if(verbose>4) printf("nnls_n=%d nnls_m=%d\n", nnls_n, nnls_m);
            for(int n=0; n<nnls_n; n++) nnls_x[n]=0.0;
            nnls_rnorm=0.0;
            continue;
          }
          fittedNr++;
 
          Vb=nnls_x[1];
          if(k2file[0]) k2=nnls_x[2];
          K1=nnls_x[0]-Vb*k2; if(K1>=0.0) fittedokNr++; else K1=k2=0.0;
          /* Put results to output images */
          k1img.m[pi][yi][xi][0]=K1;
          vbimg.m[pi][yi][xi][0]=Vb;
          k2img.m[pi][yi][xi][0]=k2;
        } /* next column */
      } /* next row */
    } /* next plane */
    if(verbose>0) {fprintf(stdout, " done.\n"); fflush(stdout);}
    free(nnls_mat); free(nnls_a);
    /* Show statistics on how we succeeded */
    n=(int)img.dimx*(int)img.dimy*(int)img.dimz;
    if(verbose>0) {
      fprintf(stdout, "%d out of %d pixels were fitted; %d pixels ok.\n", fittedokNr, n, fittedNr);
      fprintf(stdout, "%d pixels were thresholded.\n", thresholdNr);
    }
 
  } else {
 
    fprintf(stderr, "Error: selected method not available.");
      imgEmpty(&img); dftEmpty(&tac); imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vbimg);
    return(1);
 
  }
 
  /* No need for dynamic image any more */
  imgEmpty(&img); dftEmpty(&tac);
 
 
  /*
   *  Save parametric images
   */
 
  ret=imgWrite(k1file, &k1img);
  if(!ret) ret=imgWrite(vbfile, &vbimg);
  if(!ret && k2file[0]) ret=imgWrite(k2file, &k2img);
  imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vbimg);
  if(ret) {
    fprintf(stderr, "Error: cannot write parametric image.\n");
    return(11);
  }
  if(verbose>0) fprintf(stdout, "Parametric image(s) saved.\n");
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/