imglhbdv.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 image of blood and distribution volume (Vb and Vt)",
  "from dynamic PET image in ECAT, NIfTI, or Analyze format applying",
  "one-tissue compartmental model with arterial blood input.",
  "The compartmental models are transformed to general linear least squares",
  "functions, which are solved using Lawson-Hanson non-negative least squares",
  "(NNLS) algorithm (1). Vt is estimated directly without division (2, 3, 4).",
  " ",
  "Dynamic PET image and plasma time-activity curve (PTAC) must be corrected",
  "for decay to the tracer injection time.",
  " ",
  "Usage: @P [Options] btacfile imgfile vbfile vtfile",
  " ",
  "Options:",
  " -thr=<threshold%>",
  "     Pixels with AUC less than (threshold/100 x BTAC AUC) are set to zero",
  "     default is 0%",
  " -end=<Fit end time (min)>",
  "     Use data from 0 to end time; by default, model is fitted to all frames.",
  " -max=<Max value>",
  "     Upper limit for Vt values; by default 1000.",
  " -withVb",
  "     Vt includes Vb; by default not.",
  " -irrz",
  "     Set Vt to zero in pixels which seem to have irreversible uptake.",
  " -bc=<filename>",
  "     Dynamic image with Vb correction is saved.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "The unit of voxel values in the Vt image is (ml blood)/(ml tissue).",
  " ",
  "Example:",
  "  @P ua3818ab.kbq ua3818dy1.v ua3818vb.v ua3818vt.v",
  " ",
  "References:",
  "1. Lawson CL & Hanson RJ. Solving least squares problems.",
  "   Prentice-Hall, 1974, ISBN 0-89871-356-0.",
  "2. Zhou Y, Brasic J, Endres CJ, Kuwabara H, Kimes A, Contoreggi C, Maini A,",
  "   Ernst M, Wong DF. Binding potential image based statistical mapping for",
  "   detection of dopamine release by [11C]raclopride dynamic PET.",
  "   NeuroImage 2002;16(3):S91.",
  "3. Zhou Y, Brasic JR, Ye W, Dogan AS, Hilton J, Singer HS, Wong DF.",
  "   Quantification of cerebral serotonin binding in normal controls and",
  "   subjects with Tourette's syndrome using [11C]MDL 100,907 and",
  "   (+)[11C]McN 5652 dynamic PET with parametric imaging approach.",
  "   NeuroImage 2004;22(Suppl 2):T98.",
  "4. Hagelberg N, Aalto S, Kajander J, Oikonen V, Hinkka S, Någren K,",
  "   Hietala J, Scheinin H. Alfentanil increases cortical dopamine D2/D3",
  "   receptor binding in healthy subjects. Pain 2004;109:86-93.",
  " ",
  "See also: imglhdv, imglhk3, lhsol, imgdv, imgratio, img2tif",
  " ",
  "Keywords: image, modelling, distribution volume, vascular fraction, 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;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char    *cptr;
  char     inpfile[FILENAME_MAX], petfile[FILENAME_MAX];
  char     vtfile[FILENAME_MAX], vbfile[FILENAME_MAX];
  char     bcfile[FILENAME_MAX];
  float    calcThreshold=0.0;
  double   upperLimit=1000.0;
  double   fittime=-1.0;
  int      weight=0;
  int      irrevZero=0;
  int      withoutVb=1; // 0= Vt includes Vb, 1= Vt does not include Vb
  int      ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=petfile[0]=vtfile[0]=vbfile[0]=bcfile[0]=(char)0;
  /* Get 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==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) {
      fittime=atof_dpi(cptr+4); if(fittime>0.0) continue;
    } else if(strncasecmp(cptr, "WITHVB=", 5)==0) {
      withoutVb=0; continue;
    } else if(strncasecmp(cptr, "WITHOUTVB=", 5)==0) {
      withoutVb=1; continue;
    } else if(*cptr=='w' || *cptr=='W') {
      weight=1; continue;
    } else if(strncasecmp(cptr, "MAX=", 4)==0) {
      upperLimit=atof_dpi(cptr+4); if(upperLimit>0.0) continue;
    } else if(strcasecmp(cptr, "IRRZ")==0) {
      irrevZero=1; continue;
    } else if(strncasecmp(cptr, "BC=", 3)==0) {
      strlcpy(bcfile, cptr+3, FILENAME_MAX); 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(vbfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(vtfile, 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(!vtfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(strcasecmp(petfile, vtfile)==0 || strcasecmp(petfile, vbfile)==0 ||
     strcasecmp(vtfile, vbfile)==0 || strcasecmp(inpfile, vbfile)==0)
  {
    fprintf(stderr, "Error: check the output filenames.\n");
    return(1);
  }
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("inpfile := %s\n", inpfile);
    printf("petfile := %s\n", petfile);
    printf("vbfile := %s\n", vbfile);
    printf("vtfile := %s\n", vtfile);
    if(bcfile[0]) printf("bcfile := %s\n", bcfile);
    printf("calcThreshold :=%g\n", calcThreshold);
    if(upperLimit>0.0) printf("upperLimit := %g\n", upperLimit);
    printf("weight := %d\n", weight);
    printf("withoutVb := %d\n", withoutVb);
    if(fittime>0.0) printf("required_fittime := %g min\n", fittime);
    printf("irrevZero := %d\n", irrevZero);
    fflush(stdout);
  }
  if(verbose>8) IMG_TEST=verbose-8; else IMG_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;
  {
    char tmp[512];
    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 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>0) 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");
 
 
  /*
   *  Allocate result images, for Vb and Vt
   */
  if(verbose>1) fprintf(stdout, "allocating memory for parametric image data\n");
  IMG vbimg; imgInit(&vbimg);
  IMG vtimg; imgInit(&vtimg);
  ret=imgAllocateWithHeader(&vbimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&vtimg, 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(&vbimg); imgEmpty(&vtimg); return(4);
  }
  /* and set the more or less necessary things */
  vbimg.start[0]=vtimg.start[0]=0.0; 
  vbimg.end[0]=vtimg.end[0]=60.*fittime;
  vbimg.unit=vtimg.unit=CUNIT_ML_PER_ML; /* mL/mL */
 
 
  /*
   *  Allocate memory required by NNLS
   */
  if(verbose>1) fprintf(stdout, "allocating memory for NNLS\n");
  int      n, m, 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;
  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(&vbimg); imgEmpty(&vtimg);
    return(4);
  }
  for(n=0, dptr=nnls_mat; 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(weight==1 && img.isWeight==0) {
    for(m=0; m<nnls_m; m++) img.weight[m]=img.end[m]-img.start[m];
    img.isWeight=1;
  }
  /* Compute weights for NNLS */
  if(img.isWeight) {
    for(m=0; m<nnls_m; m++) tac.w[m]=img.weight[m];
    for(m=0; m<nnls_m; m++) {
      if(tac.w[m]<=1.0e-20) tac.w[m]=0.0; else tac.w[m]=sqrt(tac.w[m]);
    }
  }
 
 
  /*
   *  Compute pixel-by-pixel
   */
  /* Determine the threshold */
  double threshold=calcThreshold*tac.voi[0].y2[dataNr-1];
  if(verbose>1) printf("threshold_AUC := %g\n", threshold);
  if(verbose>0) {fprintf(stdout, "computing NNLS pixel-by-pixel\n"); fflush(stdout);}
  nnls_n=3;
  for(int pi=0; pi<img.dimz; pi++) {
    if(verbose>1) printf("computing plane %d\n", img.planeNumber[pi]);
    else if(verbose>0 && img.dimz>2) fprintf(stdout, ".");
    if(verbose>0) fflush(stdout);
    for(int yi=0; yi<img.dimy; yi++) {
      for(int xi=0; xi<img.dimx; xi++) {
        /* Set pixel results to zero */
        vbimg.m[pi][yi][xi][0]=0.0; 
        vtimg.m[pi][yi][xi][0]=0.0;
        /* Copy and integrate pixel curve */
        for(int fi=0; fi<tac.frameNr; fi++)
          tac.voi[1].y[fi]=img.m[pi][yi][xi][fi];
        ret=petintegral(tac.x1, tac.x2, tac.voi[1].y, tac.frameNr, tac.voi[1].y2, tac.voi[1].y3);
        if(ret) continue;
        /* if AUC at the end is less than threshold value, then do nothing */
        if(tac.voi[1].y2[dataNr-1]<threshold) continue;
 
        /* Fill the NNLS data matrix for estimating Vb */
        for(m=0; m<nnls_m; m++) {
          nnls_a[0][m]=tac.voi[0].y[m];
          nnls_a[1][m]=tac.voi[0].y2[m];
          nnls_a[2][m]=-tac.voi[1].y2[m];
          nnls_b[m]=tac.voi[1].y[m];
        }
        if(img.isWeight) {
          for(m=0; m<nnls_m; m++) {
            nnls_b[m]*=tac.w[m];
            for(n=0; n<nnls_n; n++) nnls_a[n][m]*=tac.w[m];
          }
        }
        /* NNLS */
        ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm, nnls_wp, nnls_zz, nnls_index);
        if(ret>1) { /* no solution is possible */
          for(n=0; n<nnls_n; n++) nnls_x[n]=0.0;
          nnls_rnorm=0.0;
        } else if(ret==1) { /* max iteration count exceeded */ }
        // if(verbose>10) printf("  %g %g %g\n", nnls_x[0], nnls_x[1], nnls_x[2]);
        /* Get Vb as the 1st linear coefficient */
        vbimg.m[pi][yi][xi][0]=nnls_x[0];
#if(1)
        /* Subtract vascular contribution from tissue data */
        for(int i=0; i<tac.frameNr; i++) tac.voi[1].y[i]-=nnls_x[0]*tac.voi[0].y[i];
        for(int i=0; i<tac.frameNr; i++) tac.voi[1].y2[i]-=nnls_x[0]*tac.voi[0].y2[i];
        /* Fill the NNLS data matrix for estimating Vt */
        for(m=0; m<nnls_m; m++) {
          nnls_a[0][m]=tac.voi[0].y2[m];
          nnls_a[1][m]=-tac.voi[1].y[m];
          nnls_b[m]=tac.voi[1].y2[m];
        }
        if(img.isWeight) {
          for(m=0; m<nnls_m; m++) {
            nnls_b[m]*=tac.w[m];
            for(n=0; n<2; n++) nnls_a[n][m]*=tac.w[m];
          }
        }
        /* NNLS */
        ret=nnls(nnls_a, nnls_m, 2, nnls_b, nnls_x, &nnls_rnorm, nnls_wp, nnls_zz, nnls_index);
        if(ret>1) { /* no solution is possible */
          for(n=0; n<2; n++) nnls_x[n]=0.0;
          nnls_rnorm=0.0;
        } else if(ret==1) { /* max iteration count exceeded */ }
        // if(verbose>10) printf("  %g %g\n", nnls_x[0], nnls_x[1]);
        /* Get Vt as the 1st linear coefficient */
        vtimg.m[pi][yi][xi][0]=nnls_x[0];
        /* If the second linear component was not needed in NNLS, that means that uptake is irreversible */
        if(nnls_x[1]<1.0E-20 && nnls_x[0]>1.0E-03) {
          if(irrevZero) vtimg.m[pi][yi][xi][0]=0.0; else vtimg.m[pi][yi][xi][0]=upperLimit;
        }
        /* Add Vb if requested */
        if(withoutVb==0) {
          vtimg.m[pi][yi][xi][0]+=vbimg.m[pi][yi][xi][0];
          if(vtimg.m[pi][yi][xi][0]<0.0) vtimg.m[pi][yi][xi][0]=0.0;
        }
 
 
#else
        /* If K1+Vb*k2 (2nd coefficient) is not > 0, then Vt is 0 */
        if(!(nnls_x[1]>1.0E-10)) {
          vtimg.m[pi][yi][xi][0]=0.0;
          continue;
        }
        /* If both K1+Vb*k2 and k2 (last coefficient) are small, then Vt is 0 */
        if((!(nnls_x[1]>1.0E-04) && !(nnls_x[2]>1.0E-04)) || !((nnls_x[1]+nnls_x[2])>1.0E-03)) {
          vtimg.m[pi][yi][xi][0]=0.0;
          continue;
        }
        /* If k2 (last coefficient) is not > 0, then Vt is infinite or zero */
        if(!(nnls_x[2]>1.0E-06)) {
          /* if K1+Vb+k2 is small, then Vt is zero, otherwise 'infinite' */
          if(nnls_x[1]>1.0E-03) vtimg.m[pi][yi][xi][0]=upperLimit;
          else vtimg.m[pi][yi][xi][0]=0.0;
          continue;
        }
        /* Fill the NNLS data matrix for estimating Vt */
        for(m=0; m<nnls_m; m++) {
          nnls_a[0][m]=tac.voi[0].y2[m];
          nnls_a[1][m]=tac.voi[0].y[m];
          nnls_a[2][m]=-tac.voi[1].y[m];
          nnls_b[m]=tac.voi[1].y2[m];
        }
        if(img.isWeight) {
          for(m=0; m<nnls_m; m++) {
            nnls_b[m]*=tac.w[m];
            for(n=0; n<nnls_n; n++) nnls_a[n][m]*=tac.w[m];
          }
        }
        /* NNLS */
        ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm, nnls_wp, nnls_zz, nnls_index);
        if(ret>1) { /* no solution is possible */
          for(n=0; n<nnls_n; n++) nnls_x[n]=0.0;
          nnls_rnorm=0.0;
        } else if(ret==1) { /* max iteration count exceeded */ }
        // if(verbose>10) printf("  %g %g %g\n", nnls_x[0], nnls_x[1], nnls_x[2]);
        /* Get Vt+Vb as the 1st linear coefficient */
        vtimg.m[pi][yi][xi][0]=nnls_x[0];
        /* Subtract Vb */
        if(withoutVb) {
          vtimg.m[pi][yi][xi][0]-=vbimg.m[pi][yi][xi][0];
          if(vtimg.m[pi][yi][xi][0]<0.0) vtimg.m[pi][yi][xi][0]=0.0;
        }
#endif
 
        /* Remove values that are higher than the user-defined limit */
        if(vtimg.m[pi][yi][xi][0]>upperLimit) vtimg.m[pi][yi][xi][0]=upperLimit;
 
 
      } /* next column */
    } /* next row */
  } /* next plane */
  if(verbose>0) {fprintf(stdout, " done.\n"); fflush(stdout);}
  /* No need for NNLS data anymore */
  free(nnls_mat);
 
 
  /*
   *  Save Vb corrected dynamic image, if requested
   */
  if(bcfile[0]) {
    IMG bcimg; imgInit(&bcimg);
    ret=imgAllocateWithHeader(&bcimg, img.dimz, img.dimy, img.dimx, img.dimt, &img);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot allocate memory.\n");
      imgEmpty(&img); dftEmpty(&tac); imgEmpty(&vtimg); imgEmpty(&vbimg);
      return(8);
    }
    for(int pi=0; pi<img.dimz; pi++) {
      for(int yi=0; yi<img.dimy; yi++) {
        for(int xi=0; xi<img.dimx; xi++) {
          for(int ti=0; ti<img.dimt; ti++) {
            double v=img.m[pi][yi][xi][ti]-vbimg.m[pi][yi][xi][ti]*tac.voi[0].y[ti];
            if(v<0.0) v=0.0;
            bcimg.m[pi][yi][xi][ti]=v;
          }
        }
      }
    }
    if(imgWrite(bcfile, &bcimg)) {
      fprintf(stderr, "Error: %s\n", bcimg.statmsg);
      imgEmpty(&img); dftEmpty(&tac); imgEmpty(&vtimg); imgEmpty(&vbimg); imgEmpty(&bcimg);
      return(13);
    }
    if(verbose>0) {fprintf(stdout, "Vb corrected image %s saved.\n", bcfile); fflush(stdout);}
    imgEmpty(&bcimg);
  }
 
 
  /* No need for dynamic image or input data anymore */
  imgEmpty(&img);
  dftEmpty(&tac);
 
 
  /*
   *  Save parametric images
   */
  if(imgWrite(vbfile, &vbimg)) {
    fprintf(stderr, "Error: %s\n", vbimg.statmsg);
    imgEmpty(&vtimg); imgEmpty(&vbimg);
    return(11);
  }
  if(verbose>0) {fprintf(stdout, "Vb image %s saved.\n", vbfile); fflush(stdout);}
 
  if(imgWrite(vtfile, &vtimg)) {
    fprintf(stderr, "Error: %s\n", vtimg.statmsg);
    imgEmpty(&vtimg); imgEmpty(&vbimg);
    return(12);
  }
  if(verbose>0) {fprintf(stdout, "Vt image %s saved.\n", vtfile); fflush(stdout);}
 
  imgEmpty(&vtimg); imgEmpty(&vbimg);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/