imgbfbp.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[] = {
  "Computation of parametric image of binding potential (BPnd) from",
  "dynamic PET images in ECAT, NIfTI, or Analyze format applying simplified",
  "reference tissue model (SRTM) [1]. The model is solved using the basis",
  "function approach [2], similar to the RPM program.",
  "Note that default limits for theta3 are suitable for [C-11]raclopride;",
  "for other tracers the limits need to be set using command-line options."
  " ",
  "Dynamic PET image and reference region TAC must be corrected for decay.",
  " ",
  "Fit is always weighted (unless specific weighting options are used), either",
  "based on the counts in SIF provided by user, or if SIF is not provided,",
  "then weights are estimated based on the mean radioactivity concentration",
  "in the dynamic image.",
  " ",
  "Usage: @P [Options] imgfile rtacfile bpfile [sif]",
  " ",
  "Options:",
  " -R1=<filename>",
  "     Programs computes also an R1 image.",
  " -k2=<filename>",
  "     Programs computes also a k2 image.",
  " -min=<value (1/min)>",
  "     Set minimum value for theta3; it must be >= k2min/(1+BPmax)+lambda.",
  "     Default is 0.06 min-1. Lambda for F-18 is 0.0063 and for C-11 0.034.",
  " -max=<value (1/min)>",
  "     Set maximum value for theta3; it must be <= k2max+lambda.",
  "     Default is 0.60 min-1.",
  " -nr=<value>",
  "     Set number of basis functions; default is 500, minimum 100.",
  " -bf=<filename>",
  "     Basis function curves are written in specified file.",
  " -wss=<filename>",
  "     Weighted sum-of-squares are written in specified image file.",
  " -err=<filename>",
  "     Pixels with their theta3 in its min or max value are written",
  "     in the specified imagefile with values 1 and 2, respectively,",
  "     others with value 0.",
  " -thr=<threshold%>",
  "     Pixels with AUC less than (threshold/100 x ref AUC) are set to zero;",
  "     default is 0%",
  " -DVR",
  "     Instead of BP, program saves the DVR (=BP+1) values.",
  " -noneg",
  "     Pixels with negative BP values are set to zero.", 
  " -end=<Fit end time (min)>",
  "     Use data from 0 to end time; by default, model is fitted to all frames.",
  " -savesif=<filename>",
  "     If weights are created based on image data, then created SIF can be",
  "     saved with this option.",
  " -wf | -w1",
  "     Weights are based only on frame length (-wf), or set to 1.0 (-w1).",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "  @P ua2918dy1.v ua2918cer.dft ua2918bp.v",
  " ",
  "References:",
  "1. Lammertsma AA, Hume SP. Simplified reference tissue model for PET",
  "   receptor studies. NeuroImage 1996;4:153-158.",
  "2. Gunn RN, Lammertsma AA, Hume SP, Cunningham VJ. Parametric imaging of",
  "   ligand-receptor binding in PET using a simplified reference region",
  "   model. NeuroImage 1997;6:279-287.",
  " ",
  "See also: imgunit, eframe, imgweigh, tacweigh, imgdecay, imgsrtm, bfmsrtm",
  " ",
  "Keywords: image, modelling, binding potential, SRTM, reference input",
  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;
  int      fitdimt;
  char     inpfile[FILENAME_MAX], petfile[FILENAME_MAX], bpfile[FILENAME_MAX];
  char     r1file[FILENAME_MAX], k2file[FILENAME_MAX], wssfile[FILENAME_MAX];
  char     errfile[FILENAME_MAX], siffile[FILENAME_MAX], bffile[FILENAME_MAX];
  char     newsiffile[FILENAME_MAX], tmp[FILENAME_MAX+1], *cptr;
  int      bfNr=500, bp_plus_one=0, *bf_opt_nr;
  float    threshold, calcThreshold=0.0;
  double   fittime=-1.0, f;
  double   t3min=0.06, t3max=0.6, lambda;
  double   BP, k2, wss, p1, p2, p3, rnorm_min;
  SIF      sif;
  int      ret, fi, pi, xi, yi, bi, bi_min;
  int      nosolution_nr=0, thresholded_nr=0;
  int      nonnegativity_constraint=0;
  double   WCORR_LIMIT = 100.0;
  clock_t  fitStart, fitFinish;
  int      weight_method=0; // 0=Mazoyer; 1=no weighting (set to 1.0); 
                            // 2=frame duration
 
  DFT tac, bf;
  IMG img, bpimg, r1img, k2img, wssimg, errimg;
  /* qr */
  int m, n, M, N;
  double **mem, **A, *B, X[MAX_N], *tau, *residual, RNORM, *chain;
  double *qrweight, **wws, *ws, *wwschain, *ct, *cti;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=petfile[0]=bpfile[0]=r1file[0]=k2file[0]=wssfile[0]=(char)0;
  errfile[0]=siffile[0]=newsiffile[0]=bffile[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, "DVR")==0 || strcasecmp(cptr, "BP+1")==0) {
      bp_plus_one=1; continue;
    } else if(strncasecmp(cptr, "R1=", 3)==0) {
      strlcpy(r1file, cptr+3, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "k2=", 3)==0) {
      strlcpy(k2file, cptr+3, FILENAME_MAX); 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, "MIN=", 4)==0) {
      t3min=atof_dpi(cptr+4); if(t3min>0.0) continue;
    } else if(strncasecmp(cptr, "MAX=", 4)==0) {
      t3max=atof_dpi(cptr+4); if(t3max>0.0) continue;
    } else if(strncasecmp(cptr, "BF=", 3)==0) {
      strlcpy(bffile, cptr+3, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "WSS=", 4)==0) {
      strlcpy(wssfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "ERR=", 4)==0) {
      strlcpy(errfile, cptr+4, 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) {
      fittime=60.0*atof_dpi(cptr+4); if(fittime>0.0) continue;
    } else if(strncasecmp(cptr, "NONEG", 3)==0) {
      nonnegativity_constraint=1; continue;
    } else if(strncasecmp(cptr, "SAVESIF=", 8)==0) {
      strlcpy(newsiffile, cptr+8, FILENAME_MAX); continue;
    } else if(strcasecmp(cptr, "W1")==0) {
      weight_method=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weight_method=2; 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(petfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(inpfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(bpfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(siffile, 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(!bpfile[0]) {
      fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(t3min>=t3max) {
    fprintf(stderr, "Error: invalid theta3 bounds.\n"); return(1);}
  /* Do not accidentally overwrite sif, in case user forgot to enter previous
     argument */
  if(!siffile[0]) {
    cptr=strrchr(bpfile, '.');
    if(cptr!=NULL && strcasecmp(cptr, ".SIF")==0) {
      fprintf(stderr, "Error: missing command-line argument; use option --help\n");
      return(1);
    }
  }
  /* Give warning if user entered option -savesif, but we're not gonna do it */
  if(weight_method!=0 && newsiffile[0]) {
    fprintf(stderr, "Warning: SIF will not be saved with options -w1 and -wf.\n");
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("inpfile := %s\n", inpfile);
    printf("bpfile := %s\n", bpfile);
    printf("siffile := %s\n", siffile);
    if(r1file[0]) printf("r1file := %s\n", r1file);
    if(k2file[0]) printf("k2file := %s\n", k2file);
    if(wssfile[0]) printf("wssfile := %s\n", wssfile);
    if(errfile[0]) printf("errfile := %s\n", errfile);
    if(newsiffile[0]) printf("newsiffile := %s\n", newsiffile);
    printf("bp_plus_one := %d\n", bp_plus_one);
    printf("t3min := %g\n", t3min);
    printf("t3max := %g\n", t3max);
    printf("bfNr := %d\n", bfNr);
    printf("calcThreshold := %g\n", calcThreshold);
    printf("nonnegativity_constraint := %d\n", nonnegativity_constraint);
    if(fittime>0.0) printf("required_fittime := %g min\n", fittime/60.0);
    printf("weight_method := %d\n", weight_method);
  }
  if(verbose>8) {IMG_TEST=verbose-8; SIF_TEST=verbose-8;} 
  else IMG_TEST=SIF_TEST=0;
 
 
  sifInit(&sif);
 
  /*
   *  Read PET image and reference tissue TAC
   */
  fittime/=60.; imgInit(&img); dftInit(&tac);
  ret=imgReadModelingData(
    petfile, siffile, inpfile, NULL, NULL, &fittime, &fitdimt, &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);
  }
  fittime*=60.;
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  /* Half-life ok? */
  if(img.isotopeHalflife<1.0) {
    fprintf(stderr, "Error: image does not contain isotope halflife.\n");
    imgEmpty(&img); dftEmpty(&tac); return(3);
  }
 
  /*
   *  Remove decay correction from the image and from the reference data
   */
  lambda=hl2lambda(img.isotopeHalflife/60.0);
  if(verbose>1) fprintf(stdout, "lambda := %g min-1\n", lambda);
  /* image */
  if(img.decayCorrection==IMG_DC_UNKNOWN) {
    // Assuming that image is decay corrected
    fprintf(stderr, "Warning: image decay correction is removed.\n");
    img.decayCorrection=IMG_DC_CORRECTED;
  }
  if(img.decayCorrection==IMG_DC_CORRECTED) {
    if(verbose>1) fprintf(stdout, "Removing decay correction from image\n");
    ret=imgDecayCorrection(&img, 0);
    if(ret) {
      fprintf(stderr, "Error: cannot remove decay correction from image.\n");
      imgEmpty(&img); dftEmpty(&tac); return(4);
    }
  } else if(verbose>0) {
    fprintf(stderr, "Note: image is supposed to have no decay correction.\n");
  }
  /* reference tac */
  if(tac.decayCorrected==DFT_DECAY_UNKNOWN) {
    // Assuming that reference TAC is decay corrected
    fprintf(stderr, "Warning: decay correction in reference TAC is removed.\n");
    tac.decayCorrected=DFT_DECAY_CORRECTED;
  }
  if(tac.decayCorrected==DFT_DECAY_CORRECTED) {
    if(verbose>1) 
      fprintf(stdout, "Removing decay correction from reference TAC\n");
    ret=dftDecayCorrection(&tac, img.isotopeHalflife/60.0, 0, 1, 0, 0,
                           tmp, verbose-3);
    if(ret) {
      fprintf(stderr, 
              "Error: cannot remove decay correction from reference TAC.\n");
      fprintf(stderr, "Error: %s.\n", tmp);
      imgEmpty(&img); dftEmpty(&tac); return(4);
    }
  } else if(verbose>0) {
    fprintf(stderr, 
            "Note: reference TAC is supposed to have no decay correction.\n");
  }
  /* Set time unit to min */
  ret=dftTimeunitConversion(&tac, TUNIT_MIN);
  /* Allocate memory for tissue data and integrals */
  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");
  if(verbose>1) {
    printf("fittimeFinal := %g min\n", fittime/60.0);
    printf("fitdimt := %d\n", fitdimt);
  }
 
 
  /* Determine the threshold */
  if(verbose>20) dftPrint(&tac);
  if(verbose>2) printf("ref_AUC[%g] := %g\n",
                  tac.x2[fitdimt-1], tac.voi[0].y2[fitdimt-1]);
  threshold=calcThreshold*tac.voi[0].y2[fitdimt-1];
  if(verbose>1) printf("threshold_AUC = %g\n", threshold);
 
 
  /*
   *  Determine the weights for the fit
   */
  if(verbose>1) printf("setting weights\n");
  /* If there is a SIF file with count data then use it */
  if(weight_method==0 && siffile[0]) {
    /* Read SIF file, this hopefully also contains the weights */
    if(verbose>1) printf("reading %s\n", siffile);
    ret=sifRead(siffile, &sif);
    if(ret) {
      fprintf(stderr, "Error in reading '%s': %s\n", siffile, siferrmsg);
      imgEmpty(&img); dftEmpty(&tac); return(6);
    }
    if(verbose>7) sifPrint(&sif);
    if(sif.frameNr<tac.frameNr) {
      fprintf(stderr, "Error: different frame number in SIF and image.\n");
      imgEmpty(&img); dftEmpty(&tac); sifEmpty(&sif); return(6);
    }
    /* Check that SIF actually contains count data */
    if(sifExistentCounts(&sif)>0) { 
      /* Calculate weights for non-decay corrected data */
      if(verbose>3) printf("compute the SIF weights.\n");
      sifModerateTrues(&sif, WCORR_LIMIT);
      sifWeight(&sif, 0.0); if(verbose>3) sifPrint(&sif);
      /* Copy the weights to tac and img */
      for(fi=0; fi<tac.frameNr; fi++) tac.w[fi]=img.weight[fi]=sif.weights[fi];
      tac.isweight=1; img.isWeight=1;
    } else {
      fprintf(stderr, "Warning: SIF does not contain valid prompts data for weighting.\n");
      strcpy(siffile, ""); // then calculated from image data in the next step
    }
    /* no need for SIF any more */
    sifEmpty(&sif);
  }
  if(weight_method==0 && !siffile[0]) { 
    /* If no SIF, then calculate weights from the image */
    if(verbose>0) fprintf(stdout, "calculating weights for the fit.\n");
    /* Allocate memory for 'sif' data */
    ret=sifSetmem(&sif, fitdimt);
    if(ret) {
      fprintf(stderr, "Error: out of memory.\n");
      imgEmpty(&img); dftEmpty(&tac); return(6);
    }
    if(verbose>3)
      printf("Memory for SIF allocated (sif.frameNr=%d).\n", sif.frameNr);
    /* Set SIF information */
    sif.colNr=4; strcpy(sif.isotope_name, "");
    for(fi=0; fi<fitdimt; fi++) {
      sif.x1[fi]=img.start[fi]; sif.x2[fi]=img.end[fi];}
    sif.scantime=img.scanStart; sif.version=1;
    /* Calculate average curve of all pixels */
    if(verbose>3) printf("calculate image average curve.\n");
    (void)imgAverageTAC(&img, img.weight);
    for(fi=0; fi<fitdimt; fi++) sif.trues[fi]=img.weight[fi];
    sifModerateTrues(&sif, WCORR_LIMIT);
    /* Multiply image data with frame durations and 1E5 to reach count scale */
    if(verbose>3) printf("Convert concentrations into count related scale.\n");
    for(fi=0; fi<sif.frameNr; fi++) {
      f=sif.x2[fi]-sif.x1[fi]; if(f<=1.0e-20) f=1.0;
      sif.trues[fi]=sif.trues[fi]*f*1.0E5;
    }
    /* Fill more of SIF contents */
    for(fi=0; fi<sif.frameNr; fi++) {
      sif.prompts[fi]=sif.trues[fi]; sif.randoms[fi]=0.0;}
    sif.colNr=5;
    /* Compute the weights */
    if(verbose>2) printf("compute the SIF weights.\n");
    sifWeight(&sif, 0.0); sif.colNr++;
    if(verbose>3) sifPrint(&sif);
    /* Save SIF data, if required */
    if(newsiffile[0]) {
      ret=sifWrite(&sif, newsiffile);
      if(ret==0 && verbose>0)
        printf("%s was created based on the applied weights.\n", newsiffile);
    }
    /* Copy the weights to tac and img */
    for(fi=0; fi<tac.frameNr; fi++) tac.w[fi]=img.weight[fi]=sif.weights[fi];
    tac.isweight=1; img.isWeight=1;
    /* no need for SIF any more */
    sifEmpty(&sif);
  }
  if(weight_method==1) {
    if(verbose>0) fprintf(stdout, "setting weights to 1.0 (no weights).\n");
    for(fi=0; fi<tac.frameNr; fi++) tac.w[fi]=img.weight[fi]=1.0;
    tac.isweight=0; img.isWeight=0;
  } 
  if(weight_method==2) {
    if(verbose>0) fprintf(stdout, "setting weights based on frame lengths.\n");
    ret=dftWeightByFreq(&tac);
    if(ret) {
      fprintf(stderr, "Error: cannot weight based on frame durations.\n");
      imgEmpty(&img); dftEmpty(&tac); return(6);
    }
    for(fi=0; fi<tac.frameNr; fi++) tac.w[fi]/=tac.frameNr;
    for(fi=0; fi<tac.frameNr; fi++) img.weight[fi]=tac.w[fi];
    tac.isweight=1; img.isWeight=1;
  }
  if(verbose>1) {
    printf("\nFrame weights:\n");
    for(fi=0; fi<tac.frameNr; fi++) 
      printf("%g %g %g\n", tac.x1[fi], tac.x2[fi], tac.w[fi]);
  }
 
 
  /*
   *  Calculate the basis functions
   */
  if(verbose>1) fprintf(stdout, "calculating basis functions\n");
  dftInit(&bf);
  ret=bf_srtm(tac.x, tac.voi[0].y, tac.frameNr, bfNr, t3min, t3max, &bf);
  if(ret) {
    fprintf(stderr, "Error: cannot calculate basis functions (%d).\n", ret);
    imgEmpty(&img); dftEmpty(&tac); return(7);
  }
 
 
  /*
   *  Allocate result images
   */
  if(verbose>1) fprintf(stdout, "allocating memory for parametric images\n");
  if(verbose>2) {
    printf("dimz := %d\n", img.dimz);
    printf("dimy := %d\n", img.dimy);
    printf("dimx := %d\n", img.dimx);
  }
  imgInit(&bpimg); imgInit(&r1img); imgInit(&k2img);
  imgInit(&wssimg); imgInit(&errimg);
  ret=imgAllocateWithHeader(&bpimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret==0 && k2file[0])
    ret=imgAllocateWithHeader(&k2img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret==0 && r1file[0])
    ret=imgAllocateWithHeader(&r1img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret==0 && wssfile[0])
    ret=imgAllocateWithHeader(&wssimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret==0 && errfile[0])
    ret=imgAllocateWithHeader(&errimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret) {
    fprintf(stderr, "Error: out of memory.\n");
    imgEmpty(&img); dftEmpty(&tac); dftEmpty(&bf);
    imgEmpty(&bpimg); imgEmpty(&r1img); imgEmpty(&k2img);
    imgEmpty(&wssimg); imgEmpty(&errimg);
    return(6);
  }
  /* Set the rest of image header */
  bpimg.start[0]=0.0; bpimg.end[0]=fittime;
  bpimg.unit=IMGUNIT_UNITLESS;
  if(k2file[0]) k2img.unit=IMGUNIT_PER_MIN;
  if(r1file[0]) r1img.unit=IMGUNIT_UNITLESS;
  if(wssfile[0]) wssimg.unit=IMGUNIT_UNITLESS;
  if(errfile[0]) errimg.unit=IMGUNIT_UNITLESS;
 
 
  /*
   *  Allocate memory for QR
   */
  if(verbose>1) fprintf(stdout, "allocating memory for QR\n");
  N=2; M=fitdimt;
 
  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");
    imgEmpty(&img); dftEmpty(&tac); dftEmpty(&bf);
    imgEmpty(&bpimg); imgEmpty(&r1img); imgEmpty(&k2img);
    imgEmpty(&wssimg); imgEmpty(&errimg);
    return(8);
  }
  for(bi=0; bi<bf.voiNr; bi++) mem[bi]=chain+bi*(M+1)*N;
  for(m=0; m<M; m++) wws[m]=wwschain+m*N;
  ws=wwschain+M*N;
 
  /* Compute QR weights beforehand for faster execution */
  for(m=0; m<M; m++) {
    if(tac.w[m]<=1.0e-20) qrweight[m]=0.0;
    else qrweight[m]=sqrt(tac.w[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");
  for(bi=0; bi<bf.voiNr; bi++) {
 
    /* Define memory site for coefficient matrix and vector tau */
    for(m=0; m<M; m++) A[m]=mem[bi]+m*N;
    tau=mem[bi]+M*N;
 
    /* Initiate matrix  (A = mem[bi]) */
    for(m=0; m<M; m++) {
      A[m][0]=tac.voi[0].y[m];
      A[m][1]=bf.voi[bi].y[m];     
    }
 
    /* Apply data weights */
    for(m=0; m<M; m++) for(n=0; n<N; n++) A[m][n]*=qrweight[m];
 
    /* Compute QR decomposition of the coefficient matrix */
    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); dftEmpty(&tac); dftEmpty(&bf); imgEmpty(&img);
      imgEmpty(&bpimg); imgEmpty(&k2img); imgEmpty(&r1img);
      imgEmpty(&wssimg); imgEmpty(&errimg);
      return(9);
    } 
  } /* next BF */
 
  /*
   *  Compute pixel-by-pixel
   */
  if(verbose>0) fprintf(stdout, "computing QR pixel-by-pixel\n");
  thresholded_nr=0; ct=tac.voi[1].y; cti=tac.voi[1].y2; nosolution_nr=0;
  /* 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(bi=0; bi<bf.voiNr; bi++) bf_opt_nr[bi]=0.0;
  fitStart=clock();
  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++) {
      if(verbose>5 && yi==4*img.dimy/10) printf("  computing row %d\n", yi+1);
      for(xi=0; xi<img.dimx; xi++) {
        if(verbose>5 && yi==4*img.dimy/10 && xi==4*img.dimx/10)
          printf("    computing column %d\n", xi+1);
 
        /* if AUC at the end is less than threshold value, then set values to 0 */
        /* Calculate pixel integral */
        for(m=0; m<tac.frameNr; m++) {ct[m]=img.m[pi][yi][xi][m];}
        ret=petintegral(tac.x1, tac.x2, ct, tac.frameNr, cti, NULL);
        if(verbose>6 && yi==4*img.dimy/10 && xi==4*img.dimx/10)
          printf("Pixel (%d,%d,%d), int= %f, threshold= %g\n",
            pi, yi, xi, cti[M-1], threshold); 
        if(cti[M-1]<threshold) {
          bpimg.m[pi][yi][xi][0]=0.0;
          if(k2file[0]) k2img.m[pi][yi][xi][0]=0.0;
          if(r1file[0]) r1img.m[pi][yi][xi][0]=0.0;
          if(wssfile[0]) wssimg.m[pi][yi][xi][0]=0.0; 
          if(errfile[0]) errimg.m[pi][yi][xi][0]=0.0;
          thresholded_nr++;
          continue;        
        }
        /* Go through all basis functions */
        bi_min=-1; rnorm_min=+1.0E80; p1=p2=p3=0.0;
        for(bi=0; bi<bf.voiNr; bi++) {
 
          /* Define memory site for preset coefficient matrix and vector tau */
          for(m=0; m<M; m++) {A[m]=mem[bi]+ m*N;}
          tau=mem[bi]+M*N;
 
          /* Get data vector */
          for(m=0; m<M; m++) {
            B[m]=img.m[pi][yi][xi][m];
            /* Apply data weights */
            B[m]*=qrweight[m];
          }
 
          /* Compute solution */
          ret=qr_solve(A, M, N, tau, B, X, residual, &RNORM, wws, ws);
          if(ret>0) { /* no solution is possible */
            for(n=0; n<N; n++) X[n]=0.0; 
            RNORM=1.0E80;
          }
          
          /* Check if this was best fit for now */
          if(RNORM<rnorm_min && bf.voi[bi].size!=lambda) {
            rnorm_min=RNORM; bi_min=bi;
            p1=X[0]; p2=X[1]; p3=bf.voi[bi_min].size;
          }
        } /* next basis function */
        if(bi_min<0 || rnorm_min>=1.0E50) {
          bpimg.m[pi][yi][xi][0]=0.0;
          if(k2file[0]) k2img.m[pi][yi][xi][0]=0.0;
          if(r1file[0]) r1img.m[pi][yi][xi][0]=0.0;
          if(wssfile[0]) wssimg.m[pi][yi][xi][0]=0.0; 
          if(errfile[0]) errimg.m[pi][yi][xi][0]=0.0;
          nosolution_nr++;
          continue;
        } else bf_opt_nr[bi_min]+=1;
 
        /* Solve k2 and BP */
        k2=p1*(p3-lambda)+p2; BP=k2/(p3-lambda) - 1.0;
        if(verbose>6 && yi==4*img.dimy/10 && xi==4*img.dimx/10) {
          printf("Pixel (%d,%d,%d), p1=%g p2=%g p3=%g\n",
                 pi, yi, xi, p1, p2, p3);
          if(verbose>10) dftPrint(&tac);
        }
 
        /* Calculate WSS */
        for(m=0, wss=0.0; m<M; m++) {
          f=p1*tac.voi[0].y[m] + p2*bf.voi[bi_min].y[m];
          f-=img.m[pi][yi][xi][m]; wss+=tac.w[m]*f*f;
        }
        /* Put results to output images */
        bpimg.m[pi][yi][xi][0]=BP;
        if(bp_plus_one!=0) bpimg.m[pi][yi][xi][0]+=1.0;
        if(k2file[0]) k2img.m[pi][yi][xi][0]=k2;
        if(r1file[0]) r1img.m[pi][yi][xi][0]=p1;
        if(wssfile[0]) wssimg.m[pi][yi][xi][0]=wss;
        if(errfile[0]) {
          if(bi_min==0) ret=1; else if(bi_min==bf.voiNr-1) ret=2; else ret=0;
          errimg.m[pi][yi][xi][0]=(float)ret;
        }
      } /* next column */
    } /* next row */
  } /* next plane */
  fitFinish=clock();
  if(verbose>0) fprintf(stdout, "done.\n");
  if(verbose>1 || thresholded_nr>0)
    fprintf(stdout, "%d pixels were not fitted due to threshold.\n",
      thresholded_nr);
  if(verbose>1 || nosolution_nr>0)
    fprintf(stdout, "no QR solution for %d pixels.\n", nosolution_nr);
 
  /* No need for dynamic image or ref tac anymore */
  imgEmpty(&img); 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(bffile[0]) {
    for(bi=0; bi<bf.voiNr; bi++)
      sprintf(bf.voi[bi].place, "%d", bf_opt_nr[bi]);
    if(dftWrite(&bf, bffile)) {
      fprintf(stderr, "Error in writing %s: %s\n", bffile, dfterrmsg);
      imgEmpty(&bpimg); imgEmpty(&k2img); imgEmpty(&r1img);
      imgEmpty(&wssimg); imgEmpty(&errimg);
      free(bf_opt_nr); return(11);
    }
    if(verbose>0) fprintf(stdout, "basis functions were written in %s\n", bffile);
  }
 
  /* No need for basis functions anymore */
  dftEmpty(&bf); free(bf_opt_nr);
 
  /*
   *  Set negative BP values to zero, if required
   */
  if(nonnegativity_constraint!=0) {
    if(verbose>1) printf("setting negative BP values to zero\n");
    imgCutoff(&bpimg, 0.0, 1);
  }
 
 
  /*
   *  Save parametric image(s)
   */
  ret=imgWrite(bpfile, &bpimg);
  if(ret) fprintf(stderr, "Error: %s\n", bpimg.statmsg);
  else if(verbose>0) {
    if(bp_plus_one==0) fprintf(stdout, "BP image %s saved.\n", bpfile);
    else fprintf(stdout, "DVR image %s saved.\n", bpfile);
  }
  if(ret==0 && r1file[0]) {
    ret=imgWrite(r1file, &r1img);
    if(ret) fprintf(stderr, "Error: %s\n", r1img.statmsg);
    else if(verbose>0) fprintf(stdout, "R1 image %s saved.\n", r1file);
  }
  if(ret==0 && k2file[0]) {
    ret=imgWrite(k2file, &k2img);
    if(ret) fprintf(stderr, "Error: %s\n", k2img.statmsg);
    else if(verbose>0) fprintf(stdout, "k2 image %s saved.\n", k2file);
  }
  if(ret==0 && wssfile[0]) {
    ret=imgWrite(wssfile, &wssimg);
    if(ret) fprintf(stderr, "Error: %s\n", wssimg.statmsg);
    else if(verbose>0) fprintf(stdout, "WSS image %s saved.\n", wssfile);
  }
  if(ret==0 && errfile[0]) {
    ret=imgWrite(errfile, &errimg);
    if(ret) fprintf(stderr, "Error: %s\n", errimg.statmsg);
    else if(verbose>0) fprintf(stdout, "Error image %s saved.\n", errfile);
  }
 
  imgEmpty(&bpimg); imgEmpty(&r1img); imgEmpty(&k2img);
  imgEmpty(&wssimg); imgEmpty(&errimg);
 
  /* How long did the fitting take */
  if(verbose>1) fprintf(stdout, "parameter estimation time := %.1f [s]\n",
    (double)(fitFinish-fitStart) / CLOCKS_PER_SEC );
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/