imgmyocc.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 4
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Estimate tissue and vascular components in dynamic myocardial PET imaging",
  "using 1TCM and LV cavity BTAC.",
  " ",
  "Usage: @P [Options] bfile imgfile tcimgfile vcimgfile",
  " ",
  "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.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imgpeak, imgaumc",
  " ",
  "Keywords: image, modelling, 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     bfile[FILENAME_MAX], petfile[FILENAME_MAX], tcfile[FILENAME_MAX];
  char     vcfile[FILENAME_MAX];
  float    calcThreshold=0.01;
  double   fittime=-1.0;
  int      ret;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  bfile[0]=petfile[0]=tcfile[0]=vcfile[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;
    }
    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(bfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(petfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(tcfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(vcfile, 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(!vcfile[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("bfile := %s\n", bfile);
    printf("petfile := %s\n", petfile);
    printf("tcfile := %s\n", tcfile);
    printf("vcfile := %s\n", vcfile);
    printf("calcThreshold :=%g\n", calcThreshold);
    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, bfile, 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<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 RV cavity and tissue TACs, too */
  if(verbose>1) fprintf(stdout, "allocating working memory for pixel TACs\n");
  ret=dftAddmem(&tac, 2);
  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, "LV-cavity");
  strcpy(tac.voi[1].name, "RV-cavity");
  strcpy(tac.voi[2].name, "tissue");
 
 
  /* Estimate RV BTAC by moving LV BTAC to left */
  {
    double tmov[tac.frameNr];
    for(int i=0; i<tac.frameNr; i++) tmov[i]=tac.x[i]-7.0/60;
    ret=interpolate4pet(tmov, tac.voi[0].y, tac.frameNr, tac.x1, tac.x2, tac.voi[1].y, NULL, NULL, tac.frameNr);
  }
  if(ret) {
    fprintf(stderr, "Error: cannot interpolate data.\n");
    imgEmpty(&img); dftEmpty(&tac); return(2);
  }
 
 
  /* 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 tcimg; imgInit(&tcimg);
  IMG vcimg; imgInit(&vcimg);
  ret=imgAllocateWithHeader(&tcimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&vcimg, 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(&tcimg); imgEmpty(&vcimg);
    return(4);
  }
  /* set 'frame time' for parametric images */
  vcimg.start[0]=tcimg.start[0]=0.0;
  vcimg.end[0]=tcimg.end[0]=60.*fittime;
  /* set units in parametric images */
  tcimg.unit=CUNIT_UNITLESS;
  vcimg.unit=CUNIT_UNITLESS;
 
 
  /* Fitting */
 
  int fittedNr=0, fittedokNr=0, thresholdNr=0;
 
  /*
   *  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(&tcimg); imgEmpty(&vcimg);
    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, *clv, *clvi, *crv;
  clv=tac.voi[0].y; clvi=tac.voi[0].y2;
  crv=tac.voi[1].y;
  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++) {
        /* Set pixel results to zero */
        tcimg.m[pi][yi][xi][0]=0.0; 
        vcimg.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 
         */
 
        /* Fill the NNLS data matrix */
        for(int m=0; m<nnls_m; m++) {
          nnls_a[0][m]=clvi[m];
          nnls_a[1][m]=clv[m];
          nnls_a[2][m]=crv[m];
          nnls_a[3][m]=-cti[m];
          nnls_b[m]=ct[m];
        }
        if(/*verbose>5 &&*/ pi==img.dimz/2 && yi==img.dimy/2 && xi==img.dimx/2) {
          printf("Matrix A                                 Array B\n");
          for(m=0; m<nnls_m; m++) {
            printf("%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_b[m]);
          }
        }
        /* NNLS */
        int nnls_index[nnls_n];
        double nnls_x[nnls_n], nnls_wp[nnls_n];
        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++; fittedokNr++;
 
        /* Put results to output images */
        vcimg.m[pi][yi][xi][0]=nnls_x[1]+0.5*nnls_x[2];
        if(vcimg.m[pi][yi][xi][0]>2.0) vcimg.m[pi][yi][xi][0]=2.0;
 
        tcimg.m[pi][yi][xi][0]=nnls_x[0]-nnls_x[3]*(nnls_x[1]+0.9*nnls_x[2]);;
        double b=nnls_x[1]+0.9*nnls_x[2];
        if(b>0.0 && b<0.67) tcimg.m[pi][yi][xi][0]/=(1.0-b);
        if(tcimg.m[pi][yi][xi][0]>1.0) tcimg.m[pi][yi][xi][0]=1.0;
        if(tcimg.m[pi][yi][xi][0]<0.0) tcimg.m[pi][yi][xi][0]=0.0;
 
      } /* 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);
  }
 
 
  /* No need for dynamic image any more */
  imgEmpty(&img); dftEmpty(&tac);
 
 
  /*
   *  Save parametric images
   */
 
  ret=imgWrite(tcfile, &tcimg);
  if(!ret) ret=imgWrite(vcfile, &vcimg);
  imgEmpty(&tcimg); imgEmpty(&vcimg);
  if(ret) {
    fprintf(stderr, "Error: cannot write parametric image.\n");
    return(11);
  }
  if(verbose>0) fprintf(stdout, "Parametric image(s) saved.\n");
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/