imgflow.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 "libtpcmisc.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodel.h"
#include "libtpcmodext.h"
/*****************************************************************************/
#define NNLS_N 3
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Estimation of rate constants K1, k2 and Va from dynamic PET image in",
  "ECAT 6.3, ECAT 7.x, NIfTI-1, or Analyze 7.5 file format using",
  "linearized two-compartment model [1].",
  "Lawson-Hanson non-negative least squares (NNLS) method [2] is used to solve",
  "general linear least squares functions.",
  " ",
  "When applied to dynamic [O-15]H2O studies, the resulting K1 image",
  "equals perfusion (blood flow) image. K1 image can be divided by tissue",
  "density (g/mL) (option -density) and multiplied by 100 (option -dL)",
  "to achieve the blood flow image in units (mL blood)/((100 g tissue) * min).",
  " ",
  "When applied to dynamic [O-15]O2 brain studies, the resulting K1 image",
  "can be converted to oxygen consumption image by multiplying it by",
  "arterial oxygen concentration [3] (ml O2 / dL blood) to get the",
  "parametric image in units mL O2 / ((100 ml tissue) * min).",
  "The model assumptions hold only when oxygen consumption is 1-6.7",
  "mL O2/(100g * min) and fit time is set to 300 s or less [3].",
  " ",
  "Arterial blood TAC must be corrected for decay and delay, with sample times",
  "in seconds. Dynamic PET image must be corrected for decay.",
  " ",
  "Usage: @P [Options] btacfile imgfile flowfile",
  " ",
  "Options:",
  " -end=<Fit end time (sec)>",
  "     Use data from 0 to end time; by default, all of it.",
  " -k2=<filename>",
  "     Parametric k2 image is saved; in some situations perfusion calculation",
  "     from k2 can be more accurate than the default assumption of f=K1.",
  "     Perfusion can be calculated from k2 using equation f=k2*pH2O, where",
  "     pH2O is the physiological partition coefficient of water in tissue.",
  " -Va=<filename>",
  "     Parametric Va image is saved.",
  "     Set -Va=0, if Va=0 is assumed; otherwise Va is always fitted.",
//  " -Vd=<filename>", // option to be removed
//  "     Parametric K1/k2 image (apparent pH2O) is saved.",
//  " -wss=<filename>", // option to be removed
//  "     Weighted sum-of-squares are written in specified image file.",
  " -thr=<threshold%>",
  "     Pixels with AUC less than (threshold/100 x max AUC) are set to zero.",
  "     Default is 5%.",
  " -max=<Max value>",
  "     Upper limit for flow values in final units.",
  " -filter",
  "     Remove parametric pixel values that are over 4x higher than",
  "     their closest neighbours.",
  " -noneg",
  "     Pixels where K1 estimates are negative are fitted again with",
  "     constraint Va=0.",
  " -mL or -dL",
  "     Units in flow and Va images will be given per mL or per dL,",
  "     respectively. By default, units are per mL.",
  " -density[=<value>]",
  "     With option -density the flow is calculated per gram or 100g tissue.",
  "     Tissue density can be changed from the default 1.04 g/ml.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "By default, the units of pixel values in the flow (K1) image is",
  "(mL blood)/((mL tissue) * min), in DV image (mL blood)/(mL tissue),",
  "in k2 image 1/min, and in Va image (mL blood/mL tissue),",
  "but the flow and Va units can be changed with above listed options.",
  " ",
  "Example 1. Calculation of perfusion and arterial blood volume image,",
  "           stopping fit at 180 s:",
  "  @P -Va=s2345va.v -end=180 s2345abfit.kbq s2345dy1.v s2345flow.v",
  " ",
  "Example 2. Dynamic image is precorrected for vascular activity,",
  "           and all available data is used:",
  "  @P -Va=0 s2345abfit.kbq s2345dy1_vacorr.v s2345flow.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.",
  "3. Ohta S, Meyer E, Thompson CJ, Gjedde A. Oxygen consumption of the",
  "   living human brain measured after a single inhalation of positron",
  "   emitting oxygen. J Cereb Blood Flow Metab. 1992;12:179-192.",
  " ",
  "See also: imgflowm, imgbfh2o, imgcbv, fit_h2o, b2t_h2o, imglhdv, fitdelay",
  " ",
  "Keywords: image, modelling, perfusion, blood flow, radiowater, 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;
/*****************************************************************************/
 
/*****************************************************************************/
/*
 *  main()
 */
int main(int argc, char *argv[])
{
  int     ai, help=0, version=0, verbose=1;
  int     fi, pi, yi, xi, ret;
  int     weight=0, dataNr=0, param_filt=0;
  char    inpfile[FILENAME_MAX], petfile[FILENAME_MAX], k1file[FILENAME_MAX];
  char    k2file[FILENAME_MAX], vafile[FILENAME_MAX], dvfile[FILENAME_MAX], wssfile[FILENAME_MAX];
  char   *cptr, tmp[FILENAME_MAX+1];
  float   calcThreshold=0.05;
  double  upperLimit=-1.0;
  int     per_gram=0;
  double  density=1.04;
  double  fittime=nan("");
  int     fitVa=1;
  int     k1_nn_constraint=0;
  int     per_dl=0;
  DFT     tac;
  IMG     img, k1img, k2img, vaimg, dvimg, wssimg;
  double *cp, *cpi, *ct, *cti, Va;
  /* nnls */
  int     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;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=petfile[0]=k1file[0]=k2file[0]=vafile[0]=dvfile[0]=wssfile[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(strcasecmp(cptr, "W")==0) {
      weight=1; continue;
    } else if(strncasecmp(cptr, "FILTER", 4)==0) {
      param_filt=1; continue;
    } else if(strncasecmp(cptr, "THR=", 4)==0) {
      double v; ret=atof_with_check(cptr+4, &v);
      if(!ret && v<100.0) {calcThreshold=(float)(0.01*v); continue;}
    } else if(strncasecmp(cptr, "DV=", 3)==0) {
      strlcpy(dvfile, cptr+3, FILENAME_MAX); if(strlen(dvfile)) continue;
    } else if(strncasecmp(cptr, "k2=", 3)==0) {
      strlcpy(k2file, cptr+3, FILENAME_MAX); if(strlen(k2file)) continue;
    } else if(strncasecmp(cptr, "VA=", 3)==0 || strncasecmp(cptr, "VB=", 3)==0) {
      strlcpy(vafile, cptr+3, FILENAME_MAX); if(strlen(vafile)) continue;
    } else if(strncasecmp(cptr, "DV=", 3)==0 || strncasecmp(cptr, "VD=", 3)==0) {
      strlcpy(dvfile, cptr+3, FILENAME_MAX); if(strlen(dvfile)) continue;
    } else if(strncasecmp(cptr, "noneg", 2)==0) {
      k1_nn_constraint=1; continue;
    } else if(strcasecmp(cptr, "DL")==0) {
      per_dl=1; continue;
    } else if(strcasecmp(cptr, "ML")==0) {
      per_dl=0; continue;
    } else if(strcasecmp(cptr, "DENSITY")==0) {
      /* if plain -density, then use default density */
      per_gram=1; continue;
    } else if(strncasecmp(cptr, "DENSITY=", 8)==0) {
      /* or read user-defined density */
      per_gram=1;
      density=atof_dpi(cptr+8); if(density>0.0) continue;
    } else if(strncasecmp(cptr, "MAX=", 4)==0) {
      upperLimit=atof_dpi(cptr+4); if(upperLimit>0.0) continue;
    } else if(strncasecmp(cptr, "END=", 4)==0) {
      fittime=atof_dpi(cptr+4)/60.; if(fittime>0.0) continue;
    } else if(strncasecmp(cptr, "WSS=", 4)==0) {
      strlcpy(wssfile, cptr+4, FILENAME_MAX); if(strlen(wssfile)>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(inpfile, 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) {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);
  }
  if(!(fittime>0.0)) fittime=1.0E+020;
  /* Check if Va=0 is assumed */
  if(strcasecmp(vafile, "NONE")==0 || strcasecmp(vafile, "ZERO")==0 || strcasecmp(vafile, "0")==0) 
    {vafile[0]=(char)0; fitVa=0;}
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("inpfile := %s\n", inpfile);
    printf("petfile := %s\n", petfile);
    printf("k1file := %s\n", k1file);
    if(k2file[0]) printf("k2file := %s\n", k2file);
    printf("fitVa := %d\n", fitVa);
    if(vafile[0]) printf("vafile := %s\n", vafile);
    if(dvfile[0]) printf("dvfile := %s\n", dvfile);
    if(wssfile[0]) printf("wssfile := %s\n", wssfile);
    printf("calcThreshold := %g\n", calcThreshold);
    printf("max := %f\n", upperLimit);
    printf("fittime := %g [sec]\n", 60.0*fittime);
    printf("per_dl := %d\n", per_dl);
    printf("per_gram := %d\n", per_gram);
    printf("density := %g\n", density);
    fflush(stdout);
  }
  if(verbose>10) IMG_TEST=verbose-10; else IMG_TEST=0;
 
 
  /*
   *  Read PET image and input TAC
   */
  if(verbose>1) printf("reading data files\n");
  dftInit(&tac); imgInit(&img);
  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(fi=0; fi<tac.frameNr; fi++) tac.voi[0].y2[fi]/=60.0;
  ret=dftTimeunitConversion(&tac, TUNIT_MIN);
  if(verbose>1) {
    printf("fittimeFinal := %g s\n", 60.*fittime);
    printf("dataNr := %d\n", dataNr);
  }
  /* Check that image is dynamic */
  if(dataNr<3) {
    fprintf(stderr, "Error: too few time frames for fitting.\n");
    if(verbose>0) imgInfo(&img);
    imgEmpty(&img); dftEmpty(&tac); return(2);
  }
 
  /*
   *  Thresholding dynamic image
   */   
  if(verbose>0) fprintf(stdout, "thresholding\n");
  long long tn;
  ret=imgThresholding(&img, calcThreshold, &tn);
  if(ret!=0) {
    fprintf(stderr, "Error in thresholding the dynamic image: %s\n", img.statmsg);
    imgEmpty(&img); dftEmpty(&tac); return(4);
  }
  if(verbose>1)
    fprintf(stdout, "threshold_cutoff_nr := %lld / %lld\n", tn, (long long)img.dimx*img.dimy*img.dimz); 
 
 
  /* Allocate memory for tissue TAC and integral */
  ret=dftAddmem(&tac, 1);
  if(ret!=0) {
    fprintf(stderr, "Error (%d) in allocating memory.\n", ret);
    imgEmpty(&img); dftEmpty(&tac);
    return(4);
  }
  strcpy(tac.voi[0].voiname, "input");
  strcpy(tac.voi[1].voiname, "tissue");
 
 
 
  /*
   *  Allocate result images and fill the header info
   */
  if(verbose>1) printf("allocating memory for parametric images\n");
  imgInit(&k1img); imgInit(&k2img); imgInit(&dvimg); imgInit(&vaimg); imgInit(&wssimg);
  /* K1 image */
  ret=imgAllocateWithHeader(&k1img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret==0) {
    k1img.unit=IMGUNIT_ML_PER_ML_PER_MIN;
    k1img.decayCorrection=IMG_DC_NONCORRECTED; k1img.isWeight=0;
    k1img.start[0]=img.start[0]; k1img.end[0]=img.end[dataNr-1];
  }
  /* Optional parametric images */
  if(ret==0 && k2file[0]) {
    ret=imgAllocateWithHeader(&k2img, img.dimz, img.dimy, img.dimx, 1, &k1img);
    if(ret==0) k2img.unit=IMGUNIT_PER_MIN;
  }
  if(ret==0 && dvfile[0]) {
    ret=imgAllocateWithHeader(&dvimg, img.dimz, img.dimy, img.dimx, 1, &k1img);
    if(ret==0) dvimg.unit=IMGUNIT_ML_PER_ML;
  }
  if(ret==0 && vafile[0]) {
    ret=imgAllocateWithHeader(&vaimg, img.dimz, img.dimy, img.dimx, 1, &k1img);
    if(ret==0) vaimg.unit=IMGUNIT_ML_PER_ML;
  }
  if(ret==0 && wssfile[0]) {
    ret=imgAllocateWithHeader(&wssimg, img.dimz, img.dimy, img.dimx, 1, &k1img);
    if(ret==0) wssimg.unit=IMGUNIT_UNKNOWN;
  }
  if(ret) {
    fprintf(stderr, "Error: cannot setup memory for parametric image.\n");
    imgEmpty(&img); dftEmpty(&tac);
    imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vaimg); imgEmpty(&wssimg);
    return(6);
  }
 
  /*
   *  Allocate memory required by NNLS
   */
  if(verbose>1) printf("allocating memory for NNLS\n");
  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(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vaimg); imgEmpty(&wssimg);
    return(7);
  }
  dptr=nnls_mat;
  for(int n=0; 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(verbose>2) printf("working with NNLS weights\n");
  if(weight==1 && img.isWeight==0) {
    for(int m=0; m<nnls_m; m++) img.weight[m]=img.end[m]-img.start[m];
    img.isWeight=1;
  }
  /* Compute NNLS weights */
  if(img.isWeight) {
    for(int m=0; m<nnls_m; m++) {
      tac.w[m]=img.weight[m]; if(tac.w[m]<=1.0e-20) tac.w[m]=0.0;
    }
  }
 
 
  /*
   *  Compute pixel-by-pixel
   */
  long long nr_of_negative_k1=0;
  cp=tac.voi[0].y; cpi=tac.voi[0].y2;
  ct=tac.voi[1].y; cti=tac.voi[1].y2;
  if(verbose>0) fprintf(stdout, "computing pixel-by-pixel\n");
  for(pi=0; pi<img.dimz; pi++) {
    if(verbose>2) printf("computing plane %d\n", img.planeNumber[pi]);
    else if(img.dimz>1 && verbose>0) {fprintf(stdout, "."); fflush(stdout);}
    for(yi=0; yi<img.dimy; yi++) {
      for(xi=0; xi<img.dimx; xi++) {
 
        if(verbose>6 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
          printf("    computing row %d, column %d\n", yi+1, xi+1);
        }
 
        /* Initiate pixel output values */
        k1img.m[pi][yi][xi][0]=0.0;
        if(k2file[0]) k2img.m[pi][yi][xi][0]=0.0;
        if(vafile[0]) vaimg.m[pi][yi][xi][0]=0.0;
        if(dvfile[0]) dvimg.m[pi][yi][xi][0]=0.0;
 
        /*
         *  Estimate K1, k2 and Va
         */
        nnls_m=dataNr; nnls_n=NNLS_N; if(fitVa==0) nnls_n--;
        /* Copy and integrate pixel curve */
        for(int m=0; m<nnls_m; m++) ct[m]=img.m[pi][yi][xi][m];
        petintegral(tac.x1, tac.x2, ct, nnls_m, cti, NULL);
        /* if AUC at the end is <= zero, then forget this pixel */
        if(cti[nnls_m-1]<=0.0) continue;
        /* Fill NNLS A matrix: */
        /* function #1: tissue integral x -1 */
        for(int m=0; m<nnls_m; m++) nnls_a[0][m]=-cti[m];
        /* function #2: integral of input */
        for(int m=0; m<nnls_m; m++) nnls_a[1][m]=cpi[m];
        /* function #3: input curve */
        if(nnls_n>2) for(int m=0; m<nnls_m; m++) nnls_a[2][m]=cp[m];
        /* Fill NNLS B array: tissue */
        for(int m=0; m<nnls_m; m++) nnls_b[m]=ct[m];
        /* Apply data weights */
        if(img.isWeight) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, tac.w);
        if(verbose>6 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
          printf("Matrix A                                 Array B\n");
          for(int 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>1) continue; /* no solution is possible */
        if(nnls_n>2) Va=nnls_x[2]; else Va=0.0;
        k1img.m[pi][yi][xi][0]=nnls_x[1];
        if(nnls_n>2) k1img.m[pi][yi][xi][0]-=Va*nnls_x[0];
        if(k2file[0]) k2img.m[pi][yi][xi][0]=nnls_x[0];
        if(vafile[0]) vaimg.m[pi][yi][xi][0]=Va;
        if(verbose>6 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
          printf("    Va=%g x[0]=%g x[1]=%g\n", Va, nnls_x[0], nnls_x[1]);
        }
 
        /*
         *  If k1<0, then fit again with assumption Va=0
         */
        if(k1_nn_constraint!=0 && k1img.m[pi][yi][xi][0]<0.0) {
          nr_of_negative_k1++;
          nnls_n=NNLS_N-1;
          /* Fill NNLS A matrix: */
          /* function #1: tissue integral x -1 */
          for(int m=0; m<nnls_m; m++) nnls_a[0][m]=-cti[m];
          /* function #2: integral of input */
          for(int m=0; m<nnls_m; m++) nnls_a[1][m]=cpi[m];
          /* Fill NNLS B array: tissue */
          for(int m=0; m<nnls_m; m++) nnls_b[m]=ct[m];
          /* Apply data weights */
          if(img.isWeight) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, tac.w);
          if(verbose>6 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
            printf("Matrix A                                 Array B\n");
            for(int m=0; m<nnls_m; m++) {
              printf("%12.3f %12.3f                 %12.3f\n",
                nnls_a[0][m], nnls_a[1][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>1) continue; /* no solution is possible */
          Va=0.0;
          k1img.m[pi][yi][xi][0]=nnls_x[1];
          if(k2file[0]) k2img.m[pi][yi][xi][0]=nnls_x[0];
          if(vafile[0]) vaimg.m[pi][yi][xi][0]=Va;
          if(verbose>6 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
            printf("    Va=%g x[0]=%g x[1]=%g\n", Va, nnls_x[0], nnls_x[1]);
          }
        } else if(k1img.m[pi][yi][xi][0]<0.0) {
          nr_of_negative_k1++;
        }
 
        /* Calculate Sum-of-Squares, if user wants to save it */
        if(wssfile[0]) {
          double wss=0.0;
          for(int m=0; m<nnls_m; m++) {
            double f=0.0;
            for(int n=0; n<nnls_n; n++) f+=nnls_x[n]*nnls_a[n][m];
            f-=nnls_b[m]; wss+=f*f;
          }
          wssimg.m[pi][yi][xi][0]=wss;
        }
 
        /*
         *  Estimate K1/k2, if necessary
         */
        if(dvfile[0]) {
          /* If Va was fitted, then subtract Va*Ca from tissue TAC */
          if(nnls_n>2) {
            for(int m=0; m<nnls_m; m++) ct[m]-=Va*cp[m];
            petintegral(tac.x1, tac.x2, ct, nnls_m, cti, NULL);
            //integrate(tac.x, ct, nnls_m, cti);
          }
          nnls_n=NNLS_N-1; /* Vb is already corrected */
          /* Use previous copies of pixel curves */
          /* Fill NNLS A matrix: */
          /* function #1: tissue x -1 */
          for(int m=0; m<nnls_m; m++) nnls_a[0][m]=-ct[m];
          /* function #2: integral of input */
          for(int m=0; m<nnls_m; m++) nnls_a[1][m]=cpi[m];
          /* Fill NNLS B array: tissue integral */
          for(int m=0; m<nnls_m; m++) nnls_b[m]=cti[m];
          /* Apply data weights */
          if(img.isWeight) nnlsWght(nnls_n, nnls_m, nnls_a, nnls_b, tac.w);
          if(verbose>6 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
            printf("Matrix A                                 Array B\n");
            for(int m=0; m<nnls_m; m++) {
              printf("%12.3f %12.3f            %12.3f\n",
                nnls_a[0][m], nnls_a[1][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>1) continue; /* no solution is possible */
          dvimg.m[pi][yi][xi][0]=nnls_x[1];
          if(verbose>6 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
            printf("    x[0]=%g x[1]=%g\n", nnls_x[0], nnls_x[1]);
          }
        }
 
      } /* next column */
    } /* next row */
  } /* next plane */
  if(verbose>0) fprintf(stdout, "done.\n");
 
  /* No need for dynamic image, NNLS matrix or curves any more */
  free(nnls_mat); imgEmpty(&img); dftEmpty(&tac);
 
  /* Tell user how many negative k1 values we got 
     (it is possible when Va is fitted) */
  if(verbose>1 && nr_of_negative_k1>0) {
    fprintf(stdout, "Negative flow in %lld pixels\n", nr_of_negative_k1);
    if(k1_nn_constraint==1) fprintf(stdout, "  corrected by assuming Va=0.\n");
  }
 
 
  /*
   *  Convert units to per dL and/or gram if necessary
   */
  if(per_dl || per_gram) {
    /* K1 */
    float f;
    if(per_dl && per_gram) {
      if(verbose>1) printf("converting to per 100g\n");
      f=100.0/density; k1img.unit=IMGUNIT_ML_PER_DL_PER_MIN;
    } else if(per_dl) {
      if(verbose>1) printf("converting to per 100mL\n");
      f=100.0; k1img.unit=IMGUNIT_ML_PER_DL_PER_MIN;
    } else if(per_gram) {
      if(verbose>1) printf("converting to per g\n");
      f=1.0/density; k1img.unit=IMGUNIT_ML_PER_ML_PER_MIN;
    } else {
      if(verbose>2) printf("no conversion from mL/(min*mL)\n");
      f=1.0;
    }
    ret=imgArithmConst(&k1img, f, '*', 1.0E+06, verbose-4);
  }
  if(per_dl && vafile[0]) {
    /* Va */
    vaimg.unit=IMGUNIT_ML_PER_DL;
    ret=imgArithmConst(&vaimg, 100.0, '*', 1.0E+06, verbose-4);
  }
 
 
  /*
   *  Filter the parametric images if required
   */
  /* Remove flow values that are higher than user-defined limit */
  if(upperLimit>0.0) {
    if(verbose>0) fprintf(stdout, "filtering out flow values exceeding %g\n", upperLimit);
    imgCutoff(&k1img, (float)upperLimit, 0);
  }
  /*
   *  Filter out pixels that are over 4x higher than their
   *  closest neighbours
   */
  if(param_filt>0) {
    if(verbose>0) printf("filtering extreme values from parametric images\n");
    imgOutlierFilter(&k1img, 4.0);
    if(k2file[0]) imgOutlierFilter(&k2img, 4.0);
    if(vafile[0]) imgOutlierFilter(&vaimg, 4.0);
    if(dvfile[0]) imgOutlierFilter(&dvimg, 4.0);
  }
 
 
  /*
   *  Save parametric image(s)
   */
  if(verbose>0) printf("writing parametric images\n");
  /* K1 */
  ret=imgWrite(k1file, &k1img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", k1img.statmsg);
    if(verbose>1) printf("ret := %d\n", ret);
    imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vaimg); imgEmpty(&wssimg);
    return(13);
  }
  if(verbose>0) fprintf(stdout, "Flow image %s saved.\n", k1file);
  /* k2 */
  if(k2file[0]) {
    ret=imgWrite(k2file, &k2img);
    if(ret) {
      fprintf(stderr, "Error: %s\n", k2img.statmsg);
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vaimg); imgEmpty(&wssimg);
      return(14);
    }
    if(verbose>0) fprintf(stdout, "k2 image %s saved.\n", k2file);
  }
  /* DV */
  if(dvfile[0]) {
    ret=imgWrite(dvfile, &dvimg);
    if(ret) {
      fprintf(stderr, "Error: %s\n", dvimg.statmsg);
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vaimg); imgEmpty(&wssimg);
      return(15);
    }
    if(verbose>0) fprintf(stdout, "DV image %s saved.\n", dvfile);
  }
  /* Va */
  if(vafile[0]) {
    ret=imgWrite(vafile, &vaimg);
    if(ret) {
      fprintf(stderr, "Error: %s\n", vaimg.statmsg);
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vaimg); imgEmpty(&wssimg);
      return(16);
    }
    if(verbose>0) fprintf(stdout, "Va image %s saved.\n", vafile);
  }
  /* WSS */
  if(wssfile[0]) {
    ret=imgWrite(wssfile, &wssimg);
    if(ret) {
      fprintf(stderr, "Error: %s\n", wssimg.statmsg);
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vaimg); imgEmpty(&wssimg);
      return(16);
    }
    if(verbose>0) fprintf(stdout, "WSS image %s saved.\n", wssfile);
  }
 
  if(verbose>2) printf("before quitting, free memory\n");
  imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&dvimg); imgEmpty(&vaimg); imgEmpty(&wssimg);
  fflush(stdout);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/