imgmtrap.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"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Computation of parametric images of K1 and Vb from dynamic PET image in",
  "ECAT, NIfTI, or Analyze format, applying NNLS.",
  "The multilinear model is either",
  "  Cpet(t) = Vb*Cb(t) + ((1-Vb)*K1+Vb*k2)*Integral[Cb(t)]",
  "            - k2*Integral[Cpet(t)]",
  "or if plasma data is available",
  "  Cpet(t) = Vb*Cb(t) + (1-Vb)*K1*Integral[(1-HCT)*Cp(t)]",
  " ",
  "The program needs the blood curve and dynamic PET image, and optionally",
  "plasma curve, multiplied by (1-HCT).",
  "The name for the resulting parametric K1 image must be given.", 
  " ",
  "Usage: @P [Options] btacfile imgfile k1file [vbfile]",
  " ",
  "Options:",
  " -k2=<value> | -k2=median",
  "     Parameter k2 is fixed to a known value, or to a median k2 from",
  "     initial fitting of unconstrained model.",
  " -ptac=<filename>",
  "     Use plasma curve, multiplied by (1-HCT), as input function;",
  "     the -k2 options are not available with this.",
  " -thr=<threshold%>",
  "     Pixels with last frame concentration less than (threshold/100 x BTAC)",
  "     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.",
  " -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.",
  " -p2=<filename>",
  "     Image with pixel values of the second fitted parameter, (1-Vb)*K1.",
  " -p3=<filename>",
  "     Image with pixel values of the third fitted parameter, HCT*kc;",
  "     only available with the first model setting (with unknown PTAC).",
  " -cbv=<filename>",
  "     Image with pixel values corrected for estimated blood volume,",
  "     Cpet(t) - (1-Vb)*Cb(t).",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "The units of pixel values in the parametric images are ml/(min*ml) for K1",
  "and ml/ml for Vb.",
  " ",
  "See also: fitmtrap, b2ptrap, img2dft",
  " ",
  "Keywords: image, modelling, irreversible uptake, perfusion",
  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     ptacfile[FILENAME_MAX], btacfile[FILENAME_MAX], petfile[FILENAME_MAX];
  char     k1file[FILENAME_MAX], vbfile[FILENAME_MAX], cbvfile[FILENAME_MAX];
  char     p2file[FILENAME_MAX], p3file[FILENAME_MAX];
  float    calcThreshold=0.01;
  double   fittime=-1.0;
  double   fixedk2=nan("");
  int      fixk2=0; // Fix k2: 0=no; 1=fix to value given by user; 2=fix to median from fit
  int      weights=2; // 1=frame lengths and activity, 2=frame lengths, 2=no weighting
  int      ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ptacfile[0]=btacfile[0]=petfile[0]=k1file[0]=vbfile[0]=p2file[0]=p3file[0]=cbvfile[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(strncasecmp(cptr, "K2=MEDIAN", 4)==0) {
      fixk2=2; fixedk2=nan(""); continue;
    } else if(strncasecmp(cptr, "K2=", 3)==0) {
      ret=atof_with_check(cptr+3, &fixedk2);
      if(!ret && fixedk2>=0.0) {fixk2=1; continue;}
    } else if(strncasecmp(cptr, "PTAC=", 5)==0) {
      strlcpy(ptacfile, cptr+5, FILENAME_MAX); if(strlen(ptacfile)>0) continue;
    } else if(strncasecmp(cptr, "P2=", 3)==0) {
      strlcpy(p2file, cptr+3, FILENAME_MAX); if(strlen(p2file)>0) continue;
    } else if(strncasecmp(cptr, "P3=", 3)==0) {
      strlcpy(p3file, cptr+3, FILENAME_MAX); if(strlen(p3file)>0) continue;
    } else if(strncasecmp(cptr, "CBV=", 4)==0) {
      strlcpy(cbvfile, cptr+4, FILENAME_MAX); if(strlen(cbvfile)>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(btacfile, 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) {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);}
 
  /* Check that the options do not clash */
  if(ptacfile[0] && fixk2>0) {
    fprintf(stderr, "Error: -k2 option not available with PTAC.\n");
    return(1);
  }
  if(p3file[0]) {
    if(ptacfile[0]) {
      fprintf(stderr, "Warning: p3 not available with PTAC.\n");
      p3file[0]=(char)0;
    }
    if(fixk2==1) {
      fprintf(stderr, "Warning: p3 not available with k2.\n");
      p3file[0]=(char)0;
    }
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("btacfile := %s\n", btacfile);
    if(ptacfile[0]) printf("ptacfile := %s\n", ptacfile);
    printf("petfile := %s\n", petfile);
    printf("k1file := %s\n", k1file);
    if(vbfile[0]) printf("vbfile := %s\n", vbfile);
    if(p2file[0]) printf("p2file := %s\n", p2file);
    if(p3file[0]) printf("p3file := %s\n", p3file);
    if(cbvfile[0]) printf("cbvfile := %s\n", cbvfile);
    printf("calcThreshold :=%g\n", calcThreshold);
    printf("weights := %d\n", weights);
    printf("fixk2 := %d\n", fixk2);
    if(!isnan(fixedk2)) printf("fixedk2 := %g\n", fixedk2);
    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 TACs
   */
  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, btacfile, ptacfile, 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");
  strcpy(tac.voi[0].name, "blood");
  if(tac.voiNr>1) strcpy(tac.voi[1].name, "plasma");
  /* Set time unit to min, also for integrals in y2[] */
  if(tac.timeunit==TUNIT_SEC) {
    for(int ti=0; ti<tac.voiNr; ti++) {
      for(int fi=0; fi<tac.frameNr; fi++) tac.voi[ti].y2[fi]/=60.0;
      for(int fi=0; fi<tac.frameNr; fi++) tac.voi[ti].y3[fi]/=3600.0;
    }
  }
  ret=dftTimeunitConversion(&tac, TUNIT_MIN);
  if(verbose>1) {
    printf("fittimeFinal := %g min\n", fittime);
    printf("dataNr := %d\n", dataNr);
    printf("image_matrix_size := %dx%dx%d\n", img.dimx, img.dimy, img.dimz);
    printf("frames := %d\n", img.dimt);
  }
  /* 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[tac.voiNr].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].y[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 vbimg; imgInit(&vbimg);
  IMG p2img; imgInit(&p2img);
  IMG p3img; imgInit(&p3img);
  ret=imgAllocateWithHeader(&k1img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&vbimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&p2img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&p3img, 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(&vbimg);
    imgEmpty(&p2img); imgEmpty(&p3img);
    return(4);
  }
  /* set 'frame time' for parametric images */
  k1img.start[0]=vbimg.start[0]=p2img.start[0]=p3img.start[0]=0.0;
  k1img.end[0]=vbimg.end[0]=p2img.end[0]=p3img.end[0]=60.*fittime;
  /* set units in parametric images */
  k1img.unit=CUNIT_ML_PER_ML_PER_MIN;
  vbimg.unit=CUNIT_ML_PER_ML;
  p2img.unit=p3img.unit=CUNIT_PER_MIN;
  /* and set the more or less necessary things */
  k1img.decayCorrection=vbimg.decayCorrection=p2img.decayCorrection=p3img.decayCorrection=
    IMG_DC_NONCORRECTED;
  k1img.isWeight=vbimg.isWeight=p2img.isWeight=p3img.isWeight=0;
 
 
  /* Fitting */
 
  int fittedNr=0, fittedokNr=0, thresholdNr=0;
 
  /* 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]);
  }
 
  /*
   *  Compute full 3-parameter model, if needed
   */
  if(!ptacfile[0] && fixk2!=1) {
    if(verbose>0) fprintf(stdout, "fitting 3-parameter model\n");
    int nnls_n=3;
 
    /* Allocate memory required by NNLS */
    if(verbose>1) fprintf(stdout, "allocating memory for NNLS\n");
    int nnls_m=dataNr;
    int      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(&vbimg);
      imgEmpty(&p2img); imgEmpty(&p3img);
      return(4);
    }
    for(int n=0; n<nnls_n; n++) nnls_a[n]=nnls_mat+n*nnls_m;
 
    double *ct, *cti, *cb, *cbi;
    cb=tac.voi[0].y; cbi=tac.voi[0].y2; 
    ct=tac.voi[tac.voiNr].y; cti=tac.voi[tac.voiNr].y2;
    /* Median of p3 from reasonable pixel results will be computed */
    double *p3list=(double*)malloc(sizeof(double)*img.dimz*img.dimy*img.dimx);
    if(p3list==NULL) {
      fprintf(stderr, "Error: out of memory.\n");
      imgEmpty(&img); dftEmpty(&tac); imgEmpty(&k1img); imgEmpty(&vbimg);
      imgEmpty(&p2img); imgEmpty(&p3img);
      return(4);
    }
    int p3nr=0;
    /* pixel-by-pixel */
    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, Vb=0.0, p2=0.0, p3=0.0;
          /* Set pixel results to zero */
          k1img.m[pi][yi][xi][0]=0.0; 
          vbimg.m[pi][yi][xi][0]=0.0; 
          p2img.m[pi][yi][xi][0]=0.0; 
          p3img.m[pi][yi][xi][0]=0.0; 
          /* if value at the end is less than threshold value, then do nothing more */
          if(img.m[pi][yi][xi][dataNr-1]<threshold) {thresholdNr++; continue;}
          /* Copy pixel curve */
          for(int fi=0; fi<tac.frameNr; fi++) ct[fi]=img.m[pi][yi][xi][fi];
          /* Integrate it */
          ret=petintegral(tac.x1, tac.x2, ct, tac.frameNr, cti, NULL);
          if(ret) continue;
 
          /* Fill the NNLS data matrix */
          for(int m=0; m<nnls_m; m++) {
            nnls_a[0][m]=cb[m];
            nnls_a[1][m]=cbi[m];
            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];
            }
          /* 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[0];
          p2=nnls_x[1];
          p3=nnls_x[2];
          if(Vb>0.98) {
            K1=0.0; Vb=1.0;
          } else {
            K1=(p2-Vb*p3)/(1.0-Vb);
            if(K1>0.1 && Vb<0.6 && p3>0.1 && p3<1.0) p3list[p3nr++]=p3;
            if(K1>0.0 && K1<8.0) fittedokNr++;
          }
          if(K1<0.0) K1=0.0; else if(K1>8.0) K1=8.0;
          if(p2>5.0) p2=5.0;
          /* Put results to output images */
          vbimg.m[pi][yi][xi][0]=Vb;
          p2img.m[pi][yi][xi][0]=p2;
          p3img.m[pi][yi][xi][0]=p3;
          k1img.m[pi][yi][xi][0]=K1;
 
          /* If cbvfile was requested, and we do not run 2-parameter model */
          if(cbvfile[0] && fixk2==0) {
            for(int ti=0; ti<img.dimt; ti++) {
              img.m[pi][yi][xi][ti]-=nnls_x[0]*cb[ti];
              if(img.m[pi][yi][xi][ti]<0.0) img.m[pi][yi][xi][ti]=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 */
    {
      int 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", fittedNr, n, fittedokNr);
        fprintf(stdout, "%d pixels were thresholded.\n", thresholdNr);
      }
    }
    /* Calculate p3 median */
    if(p3nr>0) fixedk2=dmedian(p3list, p3nr);
    free(p3list);
    if(!isnan(fixedk2)) {
      if(verbose>0) printf("Median k2 := %g\n", fixedk2);
      if(verbose>1) printf("Median calculated from %d pixels\n", p3nr);
    } else if(fixk2==2) {
      fprintf(stderr, "Error: reasonable k2 median could not be calculated.\n");
      imgEmpty(&k1img); imgEmpty(&vbimg); imgEmpty(&p2img); imgEmpty(&p3img);
      imgEmpty(&img); dftEmpty(&tac);
      return(6);
    }
    /* Save p3 image, if requested */
    if(p3file[0]) {
      if(imgWrite(p3file, &p3img)) {
        fprintf(stderr, "Error: cannot write p3 image.\n");
        imgEmpty(&k1img); imgEmpty(&vbimg); imgEmpty(&p2img); imgEmpty(&p3img);
        imgEmpty(&img); dftEmpty(&tac);
        return(11);
      }
      if(verbose>0) fprintf(stdout, "P3 image saved.\n");
    }
 
    /* If 2-parameter fitting is not needed, then save results and quit */
    if(fixk2==0) {
      /* Save cbv image if requested */
      if(cbvfile[0]) {
        if(imgWrite(cbvfile, &img)) fprintf(stderr, "Error: cannot write cbv image.\n");
      }
      imgEmpty(&img); dftEmpty(&tac);
      ret=imgWrite(k1file, &k1img);
      if(!ret && vbfile[0]) ret=imgWrite(vbfile, &vbimg);
      if(!ret && p2file[0]) ret=imgWrite(p2file, &p2img);
      imgEmpty(&k1img); imgEmpty(&vbimg); imgEmpty(&p2img); imgEmpty(&p3img);
      if(ret) {
        fprintf(stderr, "Error: cannot write parametric image.\n");
        return(12);
      }
      if(verbose>0) fprintf(stdout, "Parametric image(s) saved.\n");
      return(0);
    }
    imgEmpty(&p3img);
  }
 
 
 
  /*
   *  If still here, compute 2-parameter model
   */
  if(verbose>0) fprintf(stdout, "fitting 2-parameter model\n");
  int nnls_n=2;
 
  /* Allocate memory required by NNLS */
  if(verbose>1) fprintf(stdout, "allocating memory for NNLS\n");
  int nnls_m=dataNr;
  int      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(&vbimg);
    imgEmpty(&p2img); imgEmpty(&p3img);
    return(4);
  }
  for(int n=0; n<nnls_n; n++) nnls_a[n]=nnls_mat+n*nnls_m;
 
  double *ct, *cti, *cb, *cbi, *cpi=NULL;
  cb=tac.voi[0].y; cbi=tac.voi[0].y2; 
  if(ptacfile[0]) cpi=tac.voi[1].y2;
  ct=tac.voi[tac.voiNr].y; cti=tac.voi[tac.voiNr].y2;
  /* pixel-by-pixel */
  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, Vb=0.0, p2=0.0;
        /* Set pixel results to zero */
        k1img.m[pi][yi][xi][0]=0.0; 
        vbimg.m[pi][yi][xi][0]=0.0; 
        p2img.m[pi][yi][xi][0]=0.0; 
        /* if value at the end is less than threshold value, then do nothing more */
        if(img.m[pi][yi][xi][dataNr-1]<threshold) {thresholdNr++; continue;}
        /* Copy pixel curve */
        for(int fi=0; fi<tac.frameNr; fi++) ct[fi]=img.m[pi][yi][xi][fi];
        /* Integrate it */
        ret=petintegral(tac.x1, tac.x2, ct, tac.frameNr, cti, NULL);
        if(ret) continue;
 
        /*
         *  Estimate parameters K1 and Vb
         */
 
        /* Fill the NNLS data matrix */
        if(ptacfile[0]) {
          for(int m=0; m<nnls_m; m++) {
            nnls_a[0][m]=cb[m];
            nnls_a[1][m]=cpi[m];
            nnls_b[m]=ct[m];
          }
        } else {
          for(int m=0; m<nnls_m; m++) {
            nnls_a[0][m]=cb[m]+fixedk2*cbi[m];
            nnls_a[1][m]=cbi[m];
            nnls_b[m]=ct[m]+fixedk2*cti[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];
          }
        /* 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[0];
        p2=nnls_x[1];
        if(Vb>0.98) {
          K1=p2; Vb=1.0; //K1=0.0; Vb=1.0;
        } else {
          K1=p2/(1.0-Vb);
          if(K1>0.0 && K1<8.0) fittedokNr++;
        }
        if(K1<0.0) K1=0.0; else if(K1>8.0) K1=8.0;
        if(p2>5.0) p2=5.0;
        /* Put results to output images */
        vbimg.m[pi][yi][xi][0]=Vb;
        p2img.m[pi][yi][xi][0]=p2;
        k1img.m[pi][yi][xi][0]=K1;
 
        /* If cbvfile was requested, then calculate it */
        if(cbvfile[0]) {
          for(int ti=0; ti<img.dimt; ti++) {
            img.m[pi][yi][xi][ti]-=nnls_x[0]*cb[ti];
            if(img.m[pi][yi][xi][ti]<0.0) img.m[pi][yi][xi][ti]=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 */
  {
    int 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", fittedNr, n, fittedokNr);
      fprintf(stdout, "%d pixels were thresholded.\n", thresholdNr);
    }
  }
 
  /*
   *  Save cbv image if requested
   */
  if(cbvfile[0]) {
    if(imgWrite(cbvfile, &img)) fprintf(stderr, "Error: cannot write cbv image.\n");
  }
 
 
  /* No need for dynamic data any more */
  imgEmpty(&img); dftEmpty(&tac);
 
 
  /*
   *  Save parametric images
   */
 
  ret=imgWrite(k1file, &k1img);
  if(!ret && vbfile[0]) ret=imgWrite(vbfile, &vbimg);
  if(!ret && p2file[0]) ret=imgWrite(p2file, &p2img);
  imgEmpty(&k1img); imgEmpty(&vbimg); imgEmpty(&p2img);
  if(ret) {
    fprintf(stderr, "Error: cannot write parametric image.\n");
    return(13);
  }
  if(verbose>0) fprintf(stdout, "Parametric image(s) saved.\n");
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/