imglhk3.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[] = {
  "Computation of parametric image of k3 from dynamic PET image in ECAT, NIfTI,",
  "or Analyze format applying irreversible two-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 k3file",
  " ",
  "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.",
  " -dv=<filename>",
  "     Parametric K1/(k2+k3) image is saved.",
  " -K1=<filename>",
  "     Parametric K1 image is saved.",
  " -k2=<filename>",
  "     Parametric k2 image is saved.",
  " -Ki=<filename>",
  "     Parametric Ki image is saved.",
  " -Vb=<filename>",
  "     Vb is fitted, and parametric Vb image is saved; without this,",
  "     Vb is assumed to be zero (or precorrected).",
  " -k2k3=<filename>",
  "     Parametric k2+k3 image is saved.",
//" -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 vascular volume is considered in the model setting only if file name for",
  "Vb image is given. Otherwise the contribution of vasculature to the total",
  "radioactivity concentration is assumed to be negligible, or pre-corrected.",
  "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.",
  "Vascular volume can be pre-corrected with imgcbv.", 
  " ",
  "The units of pixel values in the parametric images are 1/min for k3,",
  "ml/(min*ml) for K1 and Ki, and ml/ml for DV and Vb.",
  " ",
  "Example 1b: Vb is assumed negligible:",
  "     @P ua2917ap.kbq ua2917dy1.v ua2917k3.v",
  "Example 1b: Vb is fitted as one of the model parameters:",
  "     @P -Vb=ua2917vb.v ua2917ap.kbq ua2917dy1.v ua2917k3.v",
  "Example 1c: Vb is pre-corrected:",
  "     imgcbv ua2917dy1.v ua2917ab.kbq 0.04 ua2917dy1_cbv.v",
  "     @P ua2917ap.kbq ua2917dy1_cbv.v ua2917k3.v",
  "Example 2: K1, DV, and k3 images are saved:",
  "     @P -K1=ua2917k1.v -DV=ua2917dv.v ua2917ap.kbq ua2917dy1.v ua2917k3.v",
  " ",
  "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: imgcbv, imgki, fitk3",
  " ",
  "Keywords: image, modelling, irreversible uptake, Ki, 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], k3file[FILENAME_MAX];
  char     k1file[FILENAME_MAX], k2file[FILENAME_MAX], vbfile[FILENAME_MAX];
  char     dvfile[FILENAME_MAX], kifile[FILENAME_MAX], k2k3file[FILENAME_MAX];
  float    calcThreshold=0.01;
  int      fitVb=0; // 0=not fitted, 1=fitted
  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]=k3file[0]=kifile[0]=(char)0;
  vbfile[0]=k1file[0]=k2file[0]=k2k3file[0]=dvfile[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, "K1=", 3)==0) {
      strlcpy(k1file, cptr+3, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "K2=", 3)==0) {
      strlcpy(k2file, cptr+3, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "K2K3=", 5)==0) {
      strlcpy(k2k3file, cptr+5, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "KI=", 3)==0) {
      strlcpy(kifile, cptr+3, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "VB=", 3)==0) {
      strlcpy(vbfile, cptr+3, FILENAME_MAX); fitVb=1; continue;
    } else if(strncasecmp(cptr, "DV=", 3)==0) {
      strlcpy(dvfile, cptr+3, FILENAME_MAX); continue;
    } else 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(k3file, 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(!k3file[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n"); return(1);}
  if(strcasecmp(petfile, k3file)==0 || strcasecmp(petfile, k1file)==0 ||
     strcasecmp(petfile, vbfile)==0 || strcasecmp(inpfile, dvfile)==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("k3file := %s\n", k3file);
    if(vbfile[0]) printf("vbfile := %s\n", vbfile);
    if(k1file[0]) printf("k1file := %s\n", k1file);
    if(k2file[0]) printf("k2file := %s\n", k2file);
    if(k2k3file[0]) printf("k2k3file := %s\n", k2k3file);
    if(kifile[0]) printf("kifile := %s\n", kifile);
    if(dvfile[0]) printf("dvfile := %s\n", dvfile);
    printf("calcThreshold :=%g\n", calcThreshold);
    printf("weights := %d\n", weights);
    printf("fitVb := %d\n", fitVb);
    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 k3img; imgInit(&k3img);
  IMG k1img; imgInit(&k1img);
  IMG k2img; imgInit(&k2img);
  IMG k2k3img; imgInit(&k2k3img);
  IMG kiimg; imgInit(&kiimg);
  IMG dvimg; imgInit(&dvimg);
  IMG vbimg; imgInit(&vbimg);
  ret=imgAllocateWithHeader(&k3img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) 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(&k2k3img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&kiimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&dvimg, 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(&k3img); imgEmpty(&kiimg);
    imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vbimg); imgEmpty(&k2k3img);
    return(4);
  }
  /* set 'frame time' for parametric images */
  k3img.start[0]=k1img.start[0]=k2img.start[0]=kiimg.start[0]=vbimg.start[0]=
    dvimg.start[0]=k2k3img.start[0]=0.0;
  k3img.end[0]=k1img.end[0]=k2img.end[0]=kiimg.end[0]=vbimg.end[0]=
    dvimg.end[0]=k2k3img.end[0]=60.*fittime;
  /* set units in parametric images */
  k3img.unit=k2img.unit=k2k3img.unit=CUNIT_PER_MIN;
  dvimg.unit=vbimg.unit=CUNIT_ML_PER_ML;
  k1img.unit=kiimg.unit=CUNIT_ML_PER_ML_PER_MIN;
  /* and set the more or less necessary things */
  k3img.decayCorrection=k1img.decayCorrection=k2img.decayCorrection=kiimg.decayCorrection=
    vbimg.decayCorrection=dvimg.decayCorrection=k2k3img.decayCorrection=IMG_DC_NONCORRECTED;
  k3img.isWeight=k1img.isWeight=k2img.isWeight=kiimg.isWeight=
    vbimg.isWeight=dvimg.isWeight=k2k3img.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(fitVb==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(&k3img); imgEmpty(&kiimg);
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vbimg); imgEmpty(&k2k3img);
      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 K1, k3, Vb, Ki, K1k2k3, k2k3;
    double *ct, *cti, *cp, *cpi, *cpii;
    cp=tac.voi[0].y; cpi=tac.voi[0].y2; cpii=tac.voi[0].y3;
    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++) {
          nnls_n=NNLS_N; if(fitVb==0) nnls_n--;
          /* Set pixel results to zero */
          k3img.m[pi][yi][xi][0]=0.0; 
          k1img.m[pi][yi][xi][0]=0.0; 
          k2img.m[pi][yi][xi][0]=0.0; 
          k2k3img.m[pi][yi][xi][0]=0.0; 
          kiimg.m[pi][yi][xi][0]=0.0; 
          dvimg.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 and possibly Vb
           */
 
          /* Fill the NNLS data matrix */
          for(int m=0; m<nnls_m; m++) {
            nnls_a[0][m]=-cti[m];
            nnls_a[1][m]=cpii[m];
            nnls_a[2][m]=cpi[m];
            if(fitVb!=0) nnls_a[3][m]=cp[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;
          }
 
          if(fitVb!=0) Vb=nnls_x[3]; else Vb=0.0;
          if(Vb>=0.99) {
            fittedNr++; 
            vbimg.m[pi][yi][xi][0]=Vb;
            continue;  // leave other parameters to zero
          }
          k2k3=nnls_x[0]; // k2+k3
          K1=nnls_x[2]-Vb*k2k3;
          if(Vb<0.99) k3=nnls_x[1]/K1; else k3=0.0;
          if(Vb<0.99) K1k2k3=K1/k2k3; else K1k2k3=0.0;
          if(Vb<0.99) Ki=nnls_x[1]/k2k3; else Ki=0.0;
          if(Vb<0.9) { // If Vb was fitted, then scale to tissue volume without blood
            K1/=(1.0-Vb); 
            Ki/=(1.0-Vb); 
            K1k2k3/=(1.0-Vb);
          }
          fittedNr++;
 
          /* Subtract vascular contribution from tissue data */
          if(fitVb!=0 && Vb>0.0 && Vb<1.0) {
            for(int fi=0; fi<tac.frameNr; fi++) ct[fi]-=Vb*cp[fi];
            for(int fi=0; fi<tac.frameNr; fi++) ct[fi]/=(1.0-Vb);
            ret=petintegral(tac.x1, tac.x2, ct, tac.frameNr, cti, NULL);
            if(ret) continue;
          }
 
          /*
           *  Estimate Ki and Vd=K1/(k2+k3)
           */
          if(fitVb!=0) {
            nnls_n--; /* Vb is not fitted, already corrected */
            for(int n=0; n<nnls_n; n++) nnls_a[n]=nnls_mat+n*nnls_m;
          }
          /* Fill the NNLS data matrix */
          for(m=0; m<nnls_m; m++) {
            nnls_a[0][m]=-ct[m];
            nnls_a[1][m]=cpii[m];
            nnls_a[2][m]=cpi[m];
            nnls_b[m]=cti[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];
          }
          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\n");
            for(n=0; n<nnls_n; n++) nnls_x[n]=0.0;
            nnls_rnorm=0.0;
            continue;
          }
          Ki=nnls_x[1]; K1k2k3=nnls_x[2];
 
          /*
           *  Estimate k3
           */
          k3=0.0;
          /* Fill the NNLS data matrix */
          for(m=0; m<nnls_m; m++) {
            nnls_a[0][m]=-cpii[m];
            nnls_a[1][m]=cti[m];
            nnls_a[2][m]=ct[m];
            nnls_b[m]=cpi[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];
          }
          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\n");
            for(n=0; n<nnls_n; n++) nnls_x[n]=0.0;
            nnls_rnorm=0.0;
            continue;
          }
          k3=nnls_x[0]; fittedNr++;
 
          /* Put results to output images */
          k3img.m[pi][yi][xi][0]=k3; if(k3>0.0) fittedokNr++;
          k1img.m[pi][yi][xi][0]=K1;
          vbimg.m[pi][yi][xi][0]=Vb;
          dvimg.m[pi][yi][xi][0]=K1k2k3;
          kiimg.m[pi][yi][xi][0]=Ki;
          k2k3img.m[pi][yi][xi][0]=k2k3;
          if(Ki>0.0) k2img.m[pi][yi][xi][0]=k3*(K1/Ki-1.0);
          else if(K1k2k3>0.0) k2img.m[pi][yi][xi][0]=K1/K1k2k3;
        } /* 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(&k3img); imgEmpty(&kiimg);
    imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vbimg); imgEmpty(&k2k3img);
    return(1);
 
  }
 
  /* No need for dynamic image any more */
  imgEmpty(&img); dftEmpty(&tac);
 
 
  /*
   *  Save parametric images
   */
 
  ret=imgWrite(k3file, &k3img);
  if(!ret && k1file[0]) ret=imgWrite(k1file, &k1img);
  if(!ret && k2file[0]) ret=imgWrite(k2file, &k2img);
  if(!ret && k2k3file[0]) ret=imgWrite(k2k3file, &k2k3img);
  if(!ret && kifile[0]) ret=imgWrite(kifile, &kiimg);
  if(!ret && vbfile[0]) ret=imgWrite(vbfile, &vbimg);
  if(!ret && dvfile[0]) ret=imgWrite(dvfile, &dvimg);
  imgEmpty(&k3img); imgEmpty(&kiimg); imgEmpty(&k1img); 
  imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vbimg); imgEmpty(&k2k3img);
  if(ret) {
    fprintf(stderr, "Error: cannot write parametric image.\n");
    return(11);
  }
  if(verbose>0) fprintf(stdout, "Parametric image(s) saved.\n");
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/