imgbfk2.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 MAX_N 2
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Estimation of rate constants K1, k2, and Vb from dynamic PET image in",
  "ECAT 6.3, ECAT 7.x, NIfTI-1, or Analyze 7.5 file format using",
  "one-tissue compartment model, solved using the basis function approach.",
  "The radioactivity concentration of each voxel is assumed to consist of",
  "tissue (Ct) and blood (Cb):",
  "  Cpet(t) = Ct(t) + Vb*Cb",
  "and differential equation for Ct is:",
  "  dCt(t)/dt = K1*Cp(t) - k2*Ct(t)",
  ", where Cp(t) is arterial plasma concentration curve.",
  "Arterial plasma and blood data must be corrected for decay and delay.",
  "Dynamic PET image must be corrected for decay.",
  "Fit time must be given in seconds.",
  " ",
  "Usage: @P [Options] ptacfile btacfile imgfile fittime K1file [k2file]",
  " ",
  "Options:",
  " -Vt=<filename>",
  "     Parametric K1/k2 (Vt) image is saved.",
  " -Vb=<filename>",
  "     Parametric Vb image is saved.",
  " -k2min=<Min k2> and -k2max=<Max k2>",
  "     Enter the minimum and maximum k2 in units 1/min, applying to decay",
  "     corrected data. Defaults are k2min=0 and k2max=10, respectively.",
  " -K1min=<Min K1> and -K1max=<Max K1>",
  "     Enter the minimum and maximum K1 value; defaults are",
  "     0 and 10 mL/(mL*min), respectively.",
  " -Vtmin=<Min Vt> and -Vtmax=<Max Vt>",
  "     Enter the minimum and maximum value for apparent Vt (K1/k2);",
  "     defaults are 0 and 10 mL/mL, respectively.",
  " -Vbmin=<Min Vb> and -Vbmax=<Max Vb>",
  "     Enter the minimum and maximum value for Vb (volume fraction);",
  "     defaults are 0 and 1 mL/mL, respectively.",
  " -nr=<value>",
  "     Set number of basis functions; default is 500, minimum 100.",
  " -wss=<filename>",
  "     Weighted sum-of-squares are written in specified image file.",
  " -thr=<threshold%>",
  "     Pixels with AUC less than (threshold/100 x input AUC) are set to zero;",
  "     default is 0%",
  " -bf=<filename>",
  "     Basis function curves are written in specified file.",
  " -err=<filename>",
  "     Excluded voxels and voxels with any of the result parameters at either",
  "     of its limits is written in the specified image file with value 1,",
  "     otherwise with value 0.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1. Calculation of K1 and k2 images using 0-20 min of",
  "the dynamic image:",
  "  @P P223344apc.bld P223344ab.bld P223344.nii 20 P223344k1.nii P223344k2.nii",
  " ",
  "By default, the units of voxel values in the K1 image are",
  "(mL plasma)/((mL tissue) * min), in k2 image 1/min, in Vb image",
  "mL blood/mL tissue), and in Vt image (mL plasma)/(mL tissue),",
  " ",
  "See also: imgbfh2o, fitk2, imgunit, fitdelay, imgcbv",
  " ",
  "Keywords: image, modelling, compartmental model, basis function method",
  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 bf1TCM(
  DFT *input,
  DFT *tissue,
  DFT *bf,
  int bfNr,
  double k2min,
  double k2max,
  char *status,
  int verbose
) {
  if(verbose>0)
    printf("\n%s(*inp, *tis, *bf, %d, %g, %g, status, %d)\n", __func__, bfNr, k2min, k2max, verbose);
 
  /* Check the parameters */
  if(input==NULL || tissue==NULL || bf==NULL) {
    if(status!=NULL) strcpy(status, "program error");
    return(1);
  }
  if(input->frameNr<3 || input->voiNr<1) {
    if(status!=NULL) strcpy(status, "no input data");
    return(2);
  }
  if(tissue->frameNr<1) {
    if(status!=NULL) strcpy(status, "no pet data");
    return(3);
  }
  if(input->timeunit!=tissue->timeunit) {
    if(status!=NULL) strcpy(status, "invalid time units");
    return(4);
  }
  if(bfNr<2) {
    if(status!=NULL) strcpy(status, "invalid nr of basis functions");
    return(5);
  }
  if(k2min<1.0E-10) k2min=1.0E-10; // range calculation does not work otherwise
  if(k2min>=k2max || k2min<0.0) {
    if(status!=NULL) strcpy(status, "invalid k2 range");
    return(6);
  }
  if(verbose>1) {
    printf("input timerange: %g - %g\n", input->x[0], input->x[input->frameNr-1]);
    printf("tissue timerange: %g - %g\n", tissue->x[0], tissue->x[tissue->frameNr-1]);
  }
 
  int bi, fi, ret;
  double a, b, c;
 
  /* Allocate memory for basis functions */
  if(verbose>1) printf("allocating memory for basis functions\n");
  ret=dftSetmem(bf, tissue->frameNr, bfNr);
  if(ret) {
    if(status!=NULL) strcpy(status, "out of memory");
    return(10);
  }
 
  /* Copy and set information fields */
  bf->voiNr=bfNr; bf->frameNr=tissue->frameNr;
  bf->_type=tissue->_type;
  dftCopymainhdr2(tissue, bf, 1);
  for(bi=0; bi<bf->voiNr; bi++) {
    snprintf(bf->voi[bi].voiname, 6, "B%5.5d", bi+1);
    strcpy(bf->voi[bi].hemisphere, ".");
    strcpy(bf->voi[bi].place, ".");
    strcpy(bf->voi[bi].name, bf->voi[bi].voiname);
  }
  for(fi=0; fi<bf->frameNr; fi++) {
    bf->x[fi]=tissue->x[fi];
    bf->x1[fi]=tissue->x1[fi];
    bf->x2[fi]=tissue->x2[fi];
  }
  
  /* Compute the range of k2 values to size fields */
  if(verbose>1) printf("computing k2 values\n");
  a=log10(k2min); b=log10(k2max); c=(b-a)/(double)(bfNr-1);
  if(verbose>20) printf("a=%g b=%g, c=%g\n", a, b, c);
  for(bi=0; bi<bf->voiNr; bi++) {
    bf->voi[bi].size=pow(10.0, (double)bi*c+a);
  }
  if(verbose>2) {
    printf("final BF k2 range: %g - %g\n", bf->voi[0].size, bf->voi[bf->voiNr-1].size);
  }
  
  /* Allocate memory for simulated TAC */
  double *sim=(double*)malloc(input->frameNr*sizeof(double));
  if(sim==NULL) {
    if(status!=NULL) strcpy(status, "out of memory");
    dftEmpty(bf); return(11);
  }
  
  /* Calculate the basis functions at input time points */
  if(verbose>1) printf("computing basis functions at input sample times\n");
  for(bi=0; bi<bf->voiNr; bi++) {
    a=bf->voi[bi].size;
    ret=simC1(input->x, input->voi[0].y, input->frameNr, 1.0, a, sim);
    if(ret) {
      if(status!=NULL) strcpy(status, "simulation problem");
      free(sim); dftEmpty(bf);
      return(20);
    }
    if(verbose>100) {
      printf("\nk2 := %g\n", a);
      printf("simulated TAC:\n");
      for(fi=0; fi<input->frameNr; fi++)
        printf("  %12.6f  %12.3f\n", input->x[fi], sim[fi]);
    }
    /* interpolate to PET time frames */
    if(tissue->timetype==DFT_TIME_STARTEND)
      ret=interpolate4pet(input->x, sim, input->frameNr, tissue->x1, tissue->x2,
                          bf->voi[bi].y, NULL, NULL, bf->frameNr);
    else
      ret=interpolate(input->x, sim, input->frameNr, tissue->x,
                      bf->voi[bi].y, NULL, NULL, bf->frameNr);
    if(ret) {
      if(status!=NULL) strcpy(status, "simulation problem");
      free(sim); dftEmpty(bf);
      return(20);
    }
 
  } // next basis function
 
  free(sim);
  if(verbose>1) printf("%s() done.\n\n", __func__);
  if(status!=NULL) strcpy(status, "ok");
  return(0);
}
 
 
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  int      fitdimt;
  char     apfile[FILENAME_MAX], abfile[FILENAME_MAX], petfile[FILENAME_MAX];
  char     k1file[FILENAME_MAX], k2file[FILENAME_MAX], vbfile[FILENAME_MAX];
  char     vtfile[FILENAME_MAX], wssfile[FILENAME_MAX], errfile[FILENAME_MAX];
  char     bfsfile[FILENAME_MAX];
  int      bfNr=500, *bf_opt_nr;
  float    threshold, calcThreshold=0.0;
  double   fittime=0.0;
  double   k1min=0.0, k1max=10.0; // mL/(min*mL)
  double   k2min=0.0, k2max=10.0; // 1/min
  double   vtmin=0.0, vtmax=100.0; // mL/mL
  double   vbmin=0.0, vbmax=1.0; // mL/mL
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  apfile[0]=abfile[0]=petfile[0]=k1file[0]=k2file[0]=vbfile[0]=vtfile[0]=(char)0;
  wssfile[0]=errfile[0]=bfsfile[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, "VB=", 3)==0 || strncasecmp(cptr, "VA=", 3)==0) {
      strlcpy(vbfile, cptr+3, FILENAME_MAX); if(strlen(vbfile)) continue;
    } else if(strncasecmp(cptr, "VT=", 3)==0 || strncasecmp(cptr, "VD=", 3)==0) {
      strlcpy(vtfile, cptr+3, FILENAME_MAX); if(strlen(vtfile)) continue;
    } else if(strncasecmp(cptr, "K1min=", 6)==0 && strlen(cptr)>6) {
      if(atof_with_check(cptr+6, &k1min)==0 && k1min>=0.0) continue;
    } else if(strncasecmp(cptr, "K1max=", 6)==0 && strlen(cptr)>6) {
      if(atof_with_check(cptr+6, &k1max)==0 && k1max>=0.0) continue;
    } else if(strncasecmp(cptr, "k2min=", 6)==0) {
      if(atof_with_check(cptr+6, &k2min)==0 && k2min>=0.0) continue;
    } else if(strncasecmp(cptr, "k2max=", 6)==0) {
      if(atof_with_check(cptr+6, &k2max)==0 && k2max>=0.0) continue;
    } else if(strncasecmp(cptr, "Vtmin=", 6)==0) {
      if(atof_with_check(cptr+6, &vtmin)==0 && vtmin>=0.0) continue;
    } else if(strncasecmp(cptr, "Vtmax=", 6)==0) {
      if(atof_with_check(cptr+6, &vtmax)==0 && vtmax>=0.0) continue;
    } else if(strncasecmp(cptr, "Vbmin=", 6)==0) {
      if(atof_with_check(cptr+6, &vbmin)==0 && vbmin>=0.0) continue;
    } else if(strncasecmp(cptr, "Vbmax=", 6)==0) {
      if(atof_with_check(cptr+6, &vbmax)==0 && vbmax>=0.0) continue;
    } else if(strncasecmp(cptr, "NR=", 3)==0) {
      bfNr=atoi(cptr+3); if(bfNr>5E+04) bfNr=5E+04;
      if(bfNr>=100) continue;
    } else if(strncasecmp(cptr, "BF=", 3)==0) {
      strlcpy(bfsfile, cptr+3, FILENAME_MAX); if(strlen(bfsfile)>0) continue;
    } else if(strncasecmp(cptr, "WSS=", 4)==0) {
      strlcpy(wssfile, cptr+4, FILENAME_MAX); if(strlen(wssfile)>0) continue;
    } else if(strncasecmp(cptr, "ERR=", 4)==0) {
      strlcpy(errfile, cptr+4, FILENAME_MAX); if(strlen(errfile)>0) continue;
    } else if(strncasecmp(cptr, "THR=", 4)==0 && strlen(cptr)>4) {
      cptr+=4; if(isdigit(*cptr) || *cptr=='+' || *cptr=='-') {
        calcThreshold=0.01*atof_dpi(cptr);
        if(calcThreshold>=0.0 && calcThreshold<=2.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(apfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(abfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(petfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(atof_with_check(argv[ai], &fittime)) {
      fprintf(stderr, "Error: invalid fit time '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) strlcpy(k1file, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(k2file, 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<1.0E-04) fittime=1.0E+020;
  /* Check if Vb is constrained to zero; then abfile is not needed */
  if(vbmin<1.0E-05 && vbmax<1.0E-05) {vbmin=vbmax=0.0; abfile[0]=(char)0;}
  /* Check if Vb is constrained to a certain value; it would be stupid to save vbfile */
  if(vbmin==vbmax && vbfile[0]) {
    fprintf(stderr, "Warning: Vb image will not be created.\n");
    vbfile[0]=(char)0;
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("apfile := %s\n", apfile);
    if(abfile[0]) printf("abfile := %s\n", abfile);
    printf("k1file := %s\n", k1file);
    if(k2file[0]) printf("k2file := %s\n", k2file);
    if(vbfile[0]) printf("vbfile := %s\n", vbfile);
    if(vtfile[0]) printf("vtfile := %s\n", vtfile);
    if(wssfile[0]) printf("wssfile := %s\n", wssfile);
    if(errfile[0]) printf("errfile := %s\n", errfile);
    if(bfsfile[0]) printf("bfsfile := %s\n", bfsfile);
    printf("bfNr := %d\n", bfNr);
    printf("calcThreshold := %g\n", calcThreshold);
    printf("requested_fittime := %g [min]\n", fittime);
    printf("k1min := %g\n", k1min);
    printf("k1max := %g\n", k1max);
    printf("k2min := %g\n", k2min);
    printf("k2max := %g\n", k2max);
    printf("vbmin := %g\n", vbmin);
    printf("vbmax := %g\n", vbmax);
    printf("vtmin := %g\n", vtmin);
    printf("vtmax := %g\n", vtmax);
  }
  if(verbose>8) {IMG_TEST=verbose-8; SIF_TEST=verbose-8;} else IMG_TEST=SIF_TEST=0;
 
  /* Check user-defined parameter ranges and calculate range of k2 */
  if(k1min>=k1max) {
    fprintf(stderr, "Error: invalid range for K1 (%g - %g).\n", k1min, k1max);
    return(1);
  }
  if(k2min>=k2max) {
    fprintf(stderr, "Error: invalid range for k2 (%g - %g).\n", k2min, k2max);
    return(1);
  }
  if(vbmin>vbmax || vbmax>1.0) {
    fprintf(stderr, "Error: invalid range for Vb (%g - %g).\n", vbmin, vbmax);
    return(1);
  }
  if(vtmin>=vtmax) {
    fprintf(stderr, "Error: invalid range for Vt (%g - %g).\n", vtmin, vtmax);
    return(1);
  }
  if(k1min/vtmax>k2min) k2min=k1min/vtmax;
  if(vtmin>0.0 && k1max/vtmin<k2max) k2max=k1max/vtmin;
  if(k2max<=k2min || k2min<0.) {
    fprintf(stderr, "Error: invalid range for K1 and Vt.\n");
    return(1);
  }
 
 
  /*
   *  Read PET image and input TAC
   */
  if(verbose>1) printf("reading data files\n");
  DFT input, tac; dftInit(&input); dftInit(&tac);
  IMG img; imgInit(&img);
  {
    char tmp[256];
    int ret=imgReadModelingData(
      petfile, NULL, apfile, abfile, NULL, &fittime, &fitdimt, &img,
      &input, &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");
  }
  //printf("last x2 = %g\n", tac.x2[tac.frameNr-1]);
  /* Set time unit to min */
  dftTimeunitConversion(&input, TUNIT_MIN);
  dftTimeunitConversion(&tac, TUNIT_MIN);
  if(verbose>1) {
    printf("fittimeFinal := %g min\n", fittime);
    printf("fitdimt := %d\n", fitdimt);
    fflush(stdout);
  }
  /* Check that the image is dynamic */
  if(fitdimt<3) {
    fprintf(stderr, "Error: too few time frames for fitting.\n");
    if(verbose>2) imgInfo(&img);
    imgEmpty(&img); dftEmpty(&input); dftEmpty(&tac); return(2);
  }
  /* Allocate memory for tissue data and integrals */
  if(dftAddmem(&tac, 1)!=0) {
    fprintf(stderr, "Error: cannot allocate memory.\n"); fflush(stderr);
    imgEmpty(&img); dftEmpty(&input); dftEmpty(&tac); return(3);
  }
  /* And set indices for input and tissue TACs */
  int iAPC=0, iAB=0, iT=1; if(tac.voiNr==2) {iAB++; iT++;}
  strcpy(tac.voi[iAPC].voiname, "input");
  strcpy(tac.voi[iT].voiname, "tissue");
  if(iAB>0) strcpy(tac.voi[iAB].voiname, "blood");
  if(verbose>50) dftPrint(&tac);
 
  /* Determine the threshold */
  if(verbose>2)
    printf("input_AUC[%g] := %g\n", tac.x2[tac.frameNr-1], tac.voi[iAPC].y2[tac.frameNr-1]/60.0);
  threshold=calcThreshold*tac.voi[iAPC].y2[tac.frameNr-1]/60.0;
  if(verbose>2) {printf("threshold_AUC = %g\n", threshold); fflush(stdout);}
 
 
  /*
   *  Determine the weights for the fit
   */
  if(verbose>1) printf("setting weights based on frame lengths\n");
  if(imgSetWeights(&img, 1, verbose-5)!=0) {
    fprintf(stderr, "Warning: cannot calculate weights.\n"); fflush(stderr);
    /* set weights to 1 */
    for(int i=0; i<img.dimt; i++) img.weight[i]=1.0;
    img.isWeight=1;
  }
 
 
  /*
   *  Calculate the basis functions
   */
  if(verbose>1) {fprintf(stdout, "calculating basis functions\n"); fflush(stdout);}
  DFT bf; dftInit(&bf);
  {
    char tmp[256];
    int ai=tac.frameNr; tac.frameNr=fitdimt;
    int ret=bf1TCM(&input, &tac, &bf, bfNr, k2min, k2max, tmp, verbose-2);
    tac.frameNr=ai;
    if(ret) {
      fprintf(stderr, "Error: cannot calculate basis functions.\n");
      if(verbose>1) fprintf(stderr, "bf1TCM(): %s.\n", tmp);
      fflush(stderr);
      imgEmpty(&img); dftEmpty(&input); dftEmpty(&tac); return(6);
    }
  }
  /* Note that basis functions may need to be saved later (in bfsfile), after it is known 
     in how many image pixels each basis function was found to give the best fit */
 
 
  /*
   *  Allocate result images and fill the header info
   */
  if(verbose>1) {printf("allocating memory for parametric images\n"); fflush(stdout);}
  IMG k1img; imgInit(&k1img);
  IMG k2img; imgInit(&k2img);
  IMG vbimg; imgInit(&vbimg);
  IMG vtimg; imgInit(&vtimg);
  IMG wssimg; imgInit(&wssimg);
  IMG errimg; imgInit(&errimg);
  if(imgAllocateWithHeader(&k1img, img.dimz, img.dimy, img.dimx, 1, &img)!=0) {
    fprintf(stderr, "Error: out of memory.\n");
    imgEmpty(&img); dftEmpty(&tac); dftEmpty(&bf);
    return(8);
  }
  k1img.start[0]=0.0; k1img.end[0]=fittime*60.0;
  k1img.decayCorrection=IMG_DC_NONCORRECTED;
  k1img.isWeight=0;
  k1img.unit=IMGUNIT_ML_PER_ML_PER_MIN;
  {
    int ret=0;
    if(k2file[0]) {
      ret=imgAllocateWithHeader(&k2img, img.dimz, img.dimy, img.dimx, 1, &k1img);
      if(ret==0) k2img.unit=IMGUNIT_PER_MIN;
    }
    if(ret==0 && vtfile[0]) {
      ret=imgAllocateWithHeader(&vtimg, img.dimz, img.dimy, img.dimx, 1, &k1img);
      if(ret==0) vtimg.unit=IMGUNIT_UNITLESS;
    }
    if(ret==0 && vbfile[0]) {
      ret=imgAllocateWithHeader(&vbimg, img.dimz, img.dimy, img.dimx, 1, &k1img);
      if(ret==0) vbimg.unit=IMGUNIT_UNITLESS;
    }
    if(ret==0 && wssfile[0]) {
      ret=imgAllocateWithHeader(&wssimg, img.dimz, img.dimy, img.dimx, 1, &k1img);
      if(ret==0) wssimg.unit=IMGUNIT_UNITLESS;
    }
    if(ret==0 && errfile[0]) {
      ret=imgAllocateWithHeader(&errimg, img.dimz, img.dimy, img.dimx, 1, &k1img);
      if(ret==0) errimg.unit=IMGUNIT_UNITLESS;
    }
    if(ret) {
      fprintf(stderr, "Error: out of memory.\n"); fflush(stderr);
      imgEmpty(&img); dftEmpty(&input); dftEmpty(&tac); dftEmpty(&bf);
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vbimg); imgEmpty(&vtimg);
      imgEmpty(&wssimg); imgEmpty(&errimg);
      return(8);
    }
  }
 
 
  /*
   *  Allocate memory for QR
   */
  int /*m, n,*/ M, N;
  double **mem, **A, *B, X[MAX_N], *tau, *residual, RNORM, *chain;
  double *qrweight, **wws, *ws, *wwschain;
 
  if(verbose>1) {fprintf(stdout, "allocating memory for QR\n"); fflush(stdout);}
  M=tac.frameNr; N=2; if(vbmin==vbmax) N--;
  chain=(double*)malloc((M+1)*N*bf.voiNr * sizeof(double));
  mem=(double**)malloc(bf.voiNr * sizeof(double*));
  A=(double**)malloc(M * sizeof(double*));
  B=(double*)malloc(M*sizeof(double));
  residual=(double*)malloc(M*sizeof(double));
  qrweight=(double*)malloc(M*sizeof(double));
  wwschain=(double*)malloc((M*N+2*M)*sizeof(double));
  wws=(double**)malloc(M * sizeof(double*));
  if(chain==NULL || B==NULL || A==NULL || residual==NULL || qrweight==NULL || wwschain==NULL || wws==NULL)
  {
    fprintf(stderr, "Error: out of memory.\n"); fflush(stderr);
    imgEmpty(&img); dftEmpty(&input); dftEmpty(&tac); dftEmpty(&bf);
    imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vbimg); imgEmpty(&vtimg);
    imgEmpty(&wssimg); imgEmpty(&errimg);
    return(8);
  }
  for(int bi=0; bi<bf.voiNr; bi++) mem[bi]=chain+bi*(M+1)*N;
  for(int m=0; m<M; m++) wws[m]=wwschain+m*N;
  ws=wwschain+M*N;
 
  /* Pre-compute QR weights for faster execution */
  for(int m=0; m<M; m++) {
    if(img.weight[m]<=1.0e-20) qrweight[m]=0.0;
    else qrweight[m]=sqrt(img.weight[m]);
  }
 
  /* Make A matrix, and QR decomposition for it, for all pixels
     beforehand for faster execution */
  if(verbose>1) {fprintf(stdout, "calculating QR decomposition\n"); fflush(stdout);}
  for(int bi=0; bi<bf.voiNr; bi++) {
 
    /* Define memory site for coefficient matrix and vector tau */
    for(int m=0; m<M; m++) A[m]=mem[bi]+m*N;
    tau=mem[bi]+M*N;
 
    /* Initiate matrix  (A = mem[bi]) */
    for(int m=0; m<M; m++) {
      A[m][0]=bf.voi[bi].y[m];
      if(N>1) A[m][1]=tac.voi[iAB].y[m]; // blood TAC for Vb estimation
    }
 
    /* Apply data weights */
    for(int m=0; m<M; m++) for(int n=0; n<N; n++) A[m][n]*=qrweight[m];
 
    /* Compute QR decomposition of the coefficient matrix */
    int ret=qr_decomp(A, M, N, tau, wws, ws);
 
    if(ret>0) { /* Decomposition failed */
      free(chain); free(B); free(residual); 
      free(A); free(wwschain); free(wws); free(qrweight); free(mem);
      imgEmpty(&img); dftEmpty(&input); dftEmpty(&tac); dftEmpty(&bf);
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vbimg); imgEmpty(&vtimg);
      imgEmpty(&wssimg); imgEmpty(&errimg);
      return(9);
    } 
  } /* next BF */
 
 
  /*
   *  Compute pixel-by-pixel
   */
  if(verbose>0) {fprintf(stdout, "computing QR pixel-by-pixel\n"); fflush(stdout);}
  int nosolution_nr=0, thresholded_nr=0;
  double *ct=tac.voi[iT].y;
  double *cti=tac.voi[iT].y2;
  /* Allocate memory for BF counters on how often each BF is found to provide the optimal fit */
  bf_opt_nr=(int*)malloc(bfNr*sizeof(int));
  for(int bi=0; bi<bf.voiNr; bi++) bf_opt_nr[bi]=0.0;
  /* pixel-by-pixel */
  for(int zi=0; zi<img.dimz; zi++) {
    if(img.dimz>1 && verbose>0) {
      if(verbose>6) printf("computing plane %d\n", img.planeNumber[zi]); else fprintf(stdout, ".");
      fflush(stdout);
    }
    for(int yi=0; yi<img.dimy; yi++) {
      if(verbose>6 && yi==4*img.dimy/10) printf("  computing row %d\n", yi+1);
      for(int xi=0; xi<img.dimx; xi++) {
        //printf("%d,%d,%d\n", zi, yi, xi); fflush(stdout);
        if(verbose>6 && yi==4*img.dimy/10 && xi==4*img.dimx/10)
          printf("    computing column %d\n", xi+1);
        /* To start with, set result voxel values to zero */
        k1img.m[zi][yi][xi][0]=0.0;
        if(k2file[0]) k2img.m[zi][yi][xi][0]=0.0;
        if(vbfile[0]) vbimg.m[zi][yi][xi][0]=0.0;
        if(vtfile[0]) vtimg.m[zi][yi][xi][0]=0.0;
        if(wssfile[0]) wssimg.m[zi][yi][xi][0]=0.0; 
        if(errfile[0]) errimg.m[zi][yi][xi][0]=0.0;
 
        /* if end AUC is less than threshold value, then let results be zero and quit this pixel */
        /* Calculate pixel integral */
        for(int m=0; m<M; m++) {ct[m]=img.m[zi][yi][xi][m];}
        petintegral(tac.x1, tac.x2, ct, tac.frameNr, cti, NULL);
        //printf("last x2 = %g\n", tac.x2[tac.frameNr-1]); fflush(stdout);
        //printf("     Pixel (%d,%d,%d), int= %f, threshold= %g\n", zi, yi, xi, cti[M-1], threshold); fflush(stdout);
        if(verbose>6 && yi==4*img.dimy/10 && xi==4*img.dimx/10) {
          printf("     Pixel (%d,%d,%d), int= %f, threshold= %g\n", zi, yi, xi, cti[M-1], threshold); 
          if(verbose>7) {
            for(int m=0; m<M; m++)
              printf("     %02d:\t%g\t%g\n", m, img.m[zi][yi][xi][m], tac.voi[iAPC].y[m]);
          }
        }
        if(!(cti[M-1]>=threshold)) {
          if(errfile[0]) errimg.m[zi][yi][xi][0]=1.0;
          thresholded_nr++;
          continue;
        }
 
        /* If Vb is constrained to a value higher than 0, then subtract BTAC from tissue */
        if(vbmin==vbmax && vbmin>0.0) {
          for(int m=0; m<M; m++) {ct[m]=img.m[zi][yi][xi][m]-vbmin*tac.voi[iAB].y[m];}
          /* re-calculate the integral */
          petintegral(tac.x1, tac.x2, ct, tac.frameNr, cti, NULL);
          /* if end AUC is less than threshold value, then set values to 0 */
          if(!(cti[M-1]>0.0)) {
            if(errfile[0]) errimg.m[zi][yi][xi][0]=1.0;
            thresholded_nr++;
            continue;
          }
        }
 
        /* Go through all basis functions */
        int bi_min=-1; double rnorm_min=1.0E80;
        double p1=0.0, p2=0.0, p3=0.0;
        for(int bi=0; bi<bf.voiNr; bi++) {
        //printf("  bi=%d\n", bi); fflush(stdout);
 
          /* Define memory site for present coefficient matrix and vector tau */
          for(int m=0; m<M; m++) {A[m]=mem[bi]+ m*N;}
          tau=mem[bi]+M*N;
 
          /* Get data vector, and apply data weights */
          for(int m=0; m<M; m++) {
            B[m]=img.m[zi][yi][xi][m]*qrweight[m];
          }
 
          /* Compute solution */
          int ret=qr_solve(A, M, N, tau, B, X, residual, &RNORM, wws, ws);
          if(ret>0) { /* no solution is possible */
            for(int n=0; n<N; n++) X[n]=0.0; 
            RNORM=1.0E80;
          }
          //printf("  %g,%g %g\n", X[0], X[1], RNORM); fflush(stdout);
          /* Check if this was the best fit so far */
          if(!(RNORM<rnorm_min)) continue; // not
          /* Ok, fit is best so far, but accept it only if parameters stay inside limits */
          /* Vb */
          if(N>1 && X[1]>=1.0) { /* If Vb>1 then set Vb=1 and K1=0, and accept; will result to k2=0 */
            X[1]=1.0; X[0]=0.0;
          } else if(N>1 && (X[1]<vbmin || X[1]>vbmax)) continue;
          /* K1 */ if(X[0]<0.0 || X[0]>k1max) continue;
          /* K1/k2 */ if(!(X[0]/bf.voi[bi].size <= vtmax)) continue;
          /* Fine, we accept it for now */
          rnorm_min=RNORM; bi_min=bi;
          /* K1 */ p1=X[0];
          /* Vb */ if(N>1) p2=X[1]; else p2=vbmin;
          /* k2 */ p3=bf.voi[bi_min].size;
                   if(p1<k1min+1.0E-05) p3=0.0; // if K1==0 then k2=0
        } /* next basis function */
        //printf("  %g,%g,%g %g\n", p1, p2, p3, rnorm_min); fflush(stdout);
 
 
        /* Check that we got an acceptable result */
        if(rnorm_min>=1.0E60) { // if not...
          nosolution_nr++;
          if(errfile[0]) errimg.m[zi][yi][xi][0]=2.0;
          continue; // ... then next pixel
        }
 
        bf_opt_nr[bi_min]+=1;
        //printf("      Pixel (%d,%d,%d), K1=%g Vb=%g k2=%g\n", zi, yi, xi, p1, p2, p3); fflush(stdout);
 
        if(verbose>6 && yi==4*img.dimy/10 && xi==4*img.dimx/10) {
          printf("      Pixel (%d,%d,%d), K1=%g Vb=%g k2=%g\n", zi, yi, xi, p1, p2, p3);
          if(verbose>10) dftPrint(&tac);
        }
 
        /* Calculate WSS */
        double wss=0.0;
        for(int m=0; m<M; m++) {
          double f=p1*bf.voi[bi_min].y[m];
          if(N>1) f+=p2*tac.voi[iAB].y[m];
          f-=img.m[zi][yi][xi][m];
          wss+=img.weight[m]*f*f;
        }
        /* Put results to output images */
        k1img.m[zi][yi][xi][0]=p1;
        if(vbfile[0]) vbimg.m[zi][yi][xi][0]=p2;
        if(k2file[0]) k2img.m[zi][yi][xi][0]=p3;
        if(vtfile[0]) {double f=p1/p3; if(f>=vtmin && f<=vtmax) vtimg.m[zi][yi][xi][0]=f;}
        if(wssfile[0]) wssimg.m[zi][yi][xi][0]=wss;
        if(errfile[0]) {
          if(bi_min==0 || bi_min==bf.voiNr-1) errimg.m[zi][yi][xi][0]=2.0;
          else errimg.m[zi][yi][xi][0]=0.0;
        }
 
      } /* next column */
    } /* next row */
  } /* next plane */
  if(verbose>0) {fprintf(stdout, "done.\n"); fflush(stdout);}
 
  if(verbose>1 || thresholded_nr>0) {
    double f;
    f=(double)thresholded_nr/((double)(k1img.dimx*k1img.dimy*k1img.dimz));
    f*=100.; if(f<3.0) printf("%g%%", f); else printf("%.0f%%", f);
    printf(" of pixels were not fitted due to threshold.\n");
    if(verbose>2) printf("thresholded %d pixels\n", thresholded_nr);
    fflush(stdout);
  }
  if(verbose>1 || nosolution_nr>0) {
    fprintf(stdout, "no QR solution for %d pixels.\n", nosolution_nr); fflush(stdout);}
 
  /* No need for dynamic image or input tac anymore */
  imgEmpty(&img); dftEmpty(&input); dftEmpty(&tac);
  /* Free memory of QR */
  free(chain);  free(B); free(residual); free(A); free(wwschain); 
  free(wws); free(qrweight); free(mem);
 
  /*
   *  Save basis functions if required;
   *  this is done not before, so that also the number of optimal fits
   *  achieved with each BF can be saved as the "size".
   */
  if(bfsfile[0]) {
    for(int bi=0; bi<bf.voiNr; bi++) sprintf(bf.voi[bi].place, "%d", bf_opt_nr[bi]);
    if(dftWrite(&bf, bfsfile)) {
      fprintf(stderr, "Error in writing %s: %s\n", bfsfile, dfterrmsg); fflush(stderr);
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vbimg); imgEmpty(&vtimg);
      imgEmpty(&wssimg); imgEmpty(&errimg);
      dftEmpty(&bf); free(bf_opt_nr); return(11);
    }
    if(verbose>0) {fprintf(stdout, "basis functions written in %s\n", bfsfile); fflush(stdout);}
  }
 
  /* No need for basis functions any more */
  dftEmpty(&bf); free(bf_opt_nr);
 
 
  /*
   *  Save parametric image(s)
   */
  if(verbose>1) {printf("Saving parametric images\n"); fflush(stdout);}
  {
    int ret=imgWrite(k1file, &k1img);
    if(ret) {fprintf(stderr, "Error: %s\n", k1img.statmsg); fflush(stderr);}
    else if(verbose>0) {fprintf(stdout, "K1 image %s saved.\n", k1file); fflush(stdout);}
 
    if(ret==0 && vbfile[0]) {
      ret=imgWrite(vbfile, &vbimg);
      if(ret) {fprintf(stderr, "Error: %s\n", vbimg.statmsg); fflush(stderr);}
      else if(verbose>0) {fprintf(stdout, "Vb image %s saved.\n", vbfile); fflush(stdout);}
    }
    if(ret==0 && vtfile[0]) {
      ret=imgWrite(vtfile, &vtimg);
      if(ret) {fprintf(stderr, "Error: %s\n", vtimg.statmsg); fflush(stderr);}
      else if(verbose>0) {fprintf(stdout, "Vt image %s saved.\n", vtfile); fflush(stdout);}
    }
    if(ret==0 && k2file[0]) {
      ret=imgWrite(k2file, &k2img);
      if(ret) {fprintf(stderr, "Error: %s\n", k2img.statmsg); fflush(stderr);}
      else if(verbose>0) {fprintf(stdout, "k2 image %s saved.\n", k2file); fflush(stdout);}
    }
    if(ret==0 && wssfile[0]) {
      ret=imgWrite(wssfile, &wssimg);
      if(ret) {fprintf(stderr, "Error: %s\n", wssimg.statmsg); fflush(stderr);}
      else if(verbose>0) {fprintf(stdout, "WSS image %s saved.\n", wssfile); fflush(stdout);}
    }
    if(ret==0 && errfile[0]) {
      ret=imgWrite(errfile, &errimg);
      if(ret) {fprintf(stderr, "Error: %s\n", errimg.statmsg); fflush(stderr);}
      else if(verbose>0) {fprintf(stdout, "Error image %s saved.\n", errfile); fflush(stdout);}
    }
    if(ret) {
      imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vbimg); imgEmpty(&vtimg);
      imgEmpty(&wssimg); imgEmpty(&errimg);
      return(ret);
    }
  }
  imgEmpty(&k1img); imgEmpty(&k2img); imgEmpty(&vbimg); imgEmpty(&vtimg);
  imgEmpty(&wssimg); imgEmpty(&errimg);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/