imglhk2d.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 5
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Computation of parametric images of K1A, VbA, K1B, VbB, and k2, from",
  "dynamic PET image in ECAT, NIfTI, or Analyze format applying",
  "reversible on1-tissue compartmental model with two input functions:",
  "  Cpet(t) = VbA*CbA(t) + VbB*CbB(t) + Ct(t)",
  "  dCt(t) = K1A*CbA(t) + K1B*CbB(t) - k2*Ct(t)",
  " ",
  "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 input TACs must be corrected for decay",
  "to the tracer injection time.",
  " ",
  "Usage: @P [Options] iAfile iBfile imgfile k1Afile k1Bfile vbAfile vbBfile 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 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     inpAfile[FILENAME_MAX], inpBfile[FILENAME_MAX], petfile[FILENAME_MAX];
  char     k1Afile[FILENAME_MAX], vbAfile[FILENAME_MAX];
  char     k1Bfile[FILENAME_MAX], vbBfile[FILENAME_MAX];
  char     k2file[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);}
  inpAfile[0]=inpBfile[0]=petfile[0]=(char)0;
  k1Afile[0]=k1Bfile[0]=k2file[0]=vbAfile[0]=vbBfile[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(inpAfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(inpBfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(petfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(k1Afile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(k1Bfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(vbAfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(vbBfile, 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(!k2file[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("inpAfile := %s\n", inpAfile);
    printf("inpBfile := %s\n", inpBfile);
    printf("petfile := %s\n", petfile);
    printf("k1Afile := %s\n", k1Afile);
    printf("vbAfile := %s\n", vbAfile);
    printf("k1Bfile := %s\n", k1Bfile);
    printf("vbBfile := %s\n", vbBfile);
    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, inpAfile, inpBfile, 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 ri=0; ri<tac.voiNr; ri++) {
      for(int fi=0; fi<tac.frameNr; fi++) tac.voi[ri].y2[fi]/=60.0;
      for(int fi=0; fi<tac.frameNr; fi++) tac.voi[ri].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<6) {
    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, "inputA");
  strcpy(tac.voi[1].name, "inputB");
  strcpy(tac.voi[2].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 k1Aimg; imgInit(&k1Aimg);
  IMG vbAimg; imgInit(&vbAimg);
  IMG k1Bimg; imgInit(&k1Bimg);
  IMG vbBimg; imgInit(&vbBimg);
  IMG k2img; imgInit(&k2img);
  ret=imgAllocateWithHeader(&k1Aimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&vbAimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&k1Bimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&vbBimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&k2img, 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(&k1Aimg); imgEmpty(&vbAimg); imgEmpty(&k1Bimg); imgEmpty(&vbBimg); imgEmpty(&k2img);
    return(4);
  }
  /* set 'frame time' for parametric images */
  k1Aimg.start[0]=k1Bimg.start[0]=vbAimg.start[0]=vbBimg.start[0]=k2img.start[0]=0.0;
  k1Aimg.end[0]=k1Bimg.end[0]=vbAimg.end[0]=vbBimg.end[0]=k2img.end[0]=60.*fittime;
  /* set units in parametric images */
  k1Aimg.unit=k1Bimg.unit=CUNIT_ML_PER_ML_PER_MIN;
  vbAimg.unit=vbBimg.unit=CUNIT_ML_PER_ML;
  k2img.unit=CUNIT_PER_MIN;
 
 
  /* 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;
    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(&k1Aimg); imgEmpty(&vbAimg); imgEmpty(&k1Bimg); imgEmpty(&vbBimg); imgEmpty(&k2img);
      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, *cpA, *cpAi, *cpB, *cpBi;
    cpA=tac.voi[0].y; cpAi=tac.voi[0].y2;
    cpB=tac.voi[1].y; cpBi=tac.voi[1].y2;
    ct=tac.voi[2].y; cti=tac.voi[2].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 K1A=0.0, K1B=0.0, k2=0.0, VbA=0.0, VbB=0.0;
          /* Set pixel results to zero */
          k1Aimg.m[pi][yi][xi][0]=0.0; 
          vbAimg.m[pi][yi][xi][0]=0.0; 
          k1Bimg.m[pi][yi][xi][0]=0.0; 
          vbBimg.m[pi][yi][xi][0]=0.0; 
          k2img.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]=cpAi[m];
            nnls_a[1][m]=cpBi[m];
            nnls_a[2][m]=cpA[m];
            nnls_a[3][m]=cpB[m];
            nnls_a[4][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 %12.3f     %12.3f\n",
                nnls_a[0][m], nnls_a[1][m], nnls_a[2][m], nnls_a[3][m], nnls_a[4][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++;
 
          VbA=nnls_x[2];
          VbB=nnls_x[3];
          k2=nnls_x[4];
          K1A=nnls_x[0]; if(VbA<1.0) K1A-=VbA*k2;
          K1B=nnls_x[1]; if(VbB<1.0) K1B-=-VbB*k2;
          if((K1A+K1B)>=0.0 && (VbA+VbB)>=0.0 && (VbA+VbB)<=1.0) fittedokNr++;
          /* Put results to output images */
          k1Aimg.m[pi][yi][xi][0]=K1A;
          vbAimg.m[pi][yi][xi][0]=VbA;
          k1Bimg.m[pi][yi][xi][0]=K1B;
          vbBimg.m[pi][yi][xi][0]=VbB;
          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(&k1Aimg); imgEmpty(&vbAimg); imgEmpty(&k1Bimg); imgEmpty(&vbBimg); imgEmpty(&k2img);
    return(1);
 
  }
 
  /* No need for dynamic image any more */
  imgEmpty(&img); dftEmpty(&tac);
 
 
  /*
   *  Save parametric images
   */
 
  ret=imgWrite(k1Afile, &k1Aimg);
  if(!ret) ret=imgWrite(vbAfile, &vbAimg);
  if(!ret) ret=imgWrite(k1Bfile, &k1Bimg);
  if(!ret) ret=imgWrite(vbBfile, &vbBimg);
  if(!ret) ret=imgWrite(k2file, &k2img);
  imgEmpty(&k1Aimg); imgEmpty(&vbAimg); imgEmpty(&k1Bimg); imgEmpty(&vbBimg); imgEmpty(&k2img);
  if(ret) {
    fprintf(stderr, "Error: cannot write parametric image.\n");
    return(11);
  }
  if(verbose>0) fprintf(stdout, "Parametric image(s) saved.\n");
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/