bfmh2o.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpcbfm.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Basis function method (BFM) fitting of one-tissue compartmental water model",
  "to regional PET [O-15]H2O study data (Koeppe et al., 1985; Boellaard et al.,",
  "2005). The model parameters are blood flow (F), partition coefficient",
  "(pWater), and arterial blood volume (Va). Venous radioactivity is assumed to",
  "be the same as the tissue concentration.",
  "Extraction factor of water is assumed to be 1.0.",
  " ",
  "The blood flow obtained using PET [O-15]H2O techniques reflects tissue",
  "perfusion, since the method is dependent on the tissue exchange of",
  "labelled water. Non-nutritive flow (blood flowing through arteriovenous",
  "shunts) is not measured (Lammertsma & Jones, 1983).",
  " ",
  "Radioactivity concentration in arterial blood is given in BTAC file, and",
  "data must be corrected for physical decay and time delay.",
  "Radioactivity concentration in tissue region(s) is given in TTAC file, and",
  "data must be corrected for physical decay.",
  "Sample times must be in seconds in all data files, unless specified inside",
  "the files. Fitting time must be given in sec.",
  " ",
  "Usage: @P [options] BTAC TTAC fittime results",
  " ",
  "Options:",
  " -ml or -dl",
  "     In results the units of F and Va will be given per mL or per dL,",
  "     respectively. By default, units are per dL.",
  " -fpt",
  "     Report blood flow (perfusion) per perfusable tissue volume",
  "     (PET volume minus vascular volume). By default perfusion is reported",
  "     per PET volume.",
  " -Va=<Va(%)>",
  "     Enter a fixed Va; fitted by default.",
  " -pH2O=<Partition coefficient for water>",
  "     Enter the true partition coefficient for water to report perfusion as",
  "     F=k2*pH2O instead of default F=K1; apparent pH2O is fitted in any case.",
  " -svg=<Filename>",
  "     Plots of original and fitted TTACs are written in specified file in",
  "     SVG format.",
  " -fit=<Filename>",
  "     Fitted regional TTACs are written in specified file.",
  " -k2min=<Min k2> and -k2max=<Max k2>",
  "     Enter the minimum and maximum k2 for BFM in units 1/min (for example",
  "     0.01 and 2.5).",
  " -fmin=<Min K1> and -fmax=<Max K1>",
  "     Enter the minimum and maximum perfusion value for BFM; defaults are",
  "     0.5 and 400 ml/(dl*min), respectively.",
  " -pmin=<Min p> and -max=<pmax>",
  "     Enter the minimum and maximum value for apparent partition coefficient",
  "     for water; defaults are 0.3 and 1.0 ml/ml, respectively.",
  " -nr=<value>",
  "     Set number of basis functions; default is 5000.",
  " -bf=<filename>",
  "     Basis function curves are written in specified file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P uo1234bl.bld uo1234.tac 999 uo1234f.res",
  " ",
  "References:",
  "1. Koeppe RA et al. J Cereb Blood Flow metab. 1985;5:224-234.",
  "2. Boellaard R et al. Mol Imaging Biol. 2005;7:273-285.",
  "3. Lammertsma AA, Jones T. J Cereb Blood Flow Metab. 1983;3:416-424.",
  " ",
  "See also: fitdelay, b2t_h2o, fit_h2o, imgbfh2o, arlkup, fitk2",
  " ",
  "Keywords: TAC, modelling, perfusion, blood flow, radiowater, 1TCM",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  char    btacfile[FILENAME_MAX], ttacfile[FILENAME_MAX], resfile[FILENAME_MAX],
          fitfile[FILENAME_MAX], svgfile[FILENAME_MAX], bffile[FILENAME_MAX];
  char   *cptr;
  int     per_dl=1; // 0 or 1
  double  fVaFixed=nan("");
  double  fitdur=nan("");
  double  pWaterTrue=nan("");
  double  flowMin=0.005, flowMax=4.0; // mL/(min*mL)
  double  pWaterMin=0.3, pWaterMax=1.0; // mL/mL
  double  k2min=nan(""), k2max=nan("");; // 1/min
  int     bfNr=5000;
  int     flow_per_perfusable_tissue=0; // 0 or 1
  int     ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  btacfile[0]=ttacfile[0]=resfile[0]=fitfile[0]=svgfile[0]=bffile[0]=(char)0;
  /* 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) continue;
    if(strcasecmp(cptr, "DL")==0) {per_dl=1; continue;}
    else if(strcasecmp(cptr, "ML")==0) {per_dl=0; continue;} 
    else if(strcasecmp(cptr, "FPT")==0) {
      flow_per_perfusable_tissue=1; continue;
    } else if(strncasecmp(cptr, "Va=", 3)==0) {
      fVaFixed=0.01*atofVerified(cptr+3);
      if(!isnan(fVaFixed) && fVaFixed>=0.0 && fVaFixed<1.0) {
        if(fVaFixed>0.0 && fVaFixed<0.01)
          fprintf(stderr, "Warning: Va was set to %g%%\n", 100.*fVaFixed);
        continue;
      }
    } else if(strncasecmp(cptr, "pH2O=", 5)==0 && strlen(cptr)>5) {
      pWaterTrue=atofVerified(cptr+5);
      if(!isnan(pWaterTrue) && pWaterTrue>0.0 && pWaterTrue<=1.25) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); if(strlen(svgfile)>0) continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0) {
      strlcpy(fitfile, cptr+4, FILENAME_MAX); if(strlen(fitfile)>0) continue;
    } else if(strncasecmp(cptr, "NR=", 3)==0) {
      if(atoiCheck(cptr+3, &bfNr)==0 && bfNr>5) continue;
    } else if(strncasecmp(cptr, "BF=", 3)==0) {
      strlcpy(bffile, cptr+3, FILENAME_MAX); if(strlen(bffile)>0) continue;
    } else if(strncasecmp(cptr, "k2min=", 6)==0) {
      if(atofCheck(cptr+6, &k2min)==0 && k2min>=0.0) continue;
    } else if(strncasecmp(cptr, "k2max=", 6)==0) {
      if(atofCheck(cptr+6, &k2max)==0 && k2max>=0.0) continue;
    } else if(strncasecmp(cptr, "fmin=", 5)==0) {
      if(atofCheck(cptr+5, &flowMin)==0 && flowMin>=0.0) {flowMin*=0.01; continue;}
    } else if(strncasecmp(cptr, "fmax=", 5)==0) {
      if(atofCheck(cptr+5, &flowMax)==0 && flowMax>=0.0) {flowMax*=0.01; continue;}
    } else if(strncasecmp(cptr, "pmin=", 5)==0) {
      if(atofCheck(cptr+5, &pWaterMin)==0 && pWaterMin>0.0) continue;
    } else if(strncasecmp(cptr, "pmax=", 5)==0) {
      if(atofCheck(cptr+5, &pWaterMax)==0 && pWaterMax<1.25) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-3;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) strlcpy(btacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(ttacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(atofCheck(argv[ai], &fitdur)) {
      fprintf(stderr, "Error: invalid fit time '%s'.\n", argv[ai]);
      return(1);
    }
    if(fitdur<=0.0) fitdur=1.0E+99;
    ai++;
  }
  if(ai<argc) strlcpy(resfile, 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(!resfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("btacfile := %s\n", btacfile);
    printf("ttacfile := %s\n", ttacfile);
    printf("resfile := %s\n", resfile);
    printf("fitfile := %s\n", fitfile);
    printf("svgfile := %s\n", svgfile);
    printf("bffile := %s\n", bffile);
    printf("per_dl := %d\n", per_dl);
    printf("flow_per_perfusable_tissue := %d\n", flow_per_perfusable_tissue);
    printf("required_fitdur := %g s\n", fitdur);
    if(!isnan(fVaFixed)) printf("fVaFixed := %g%%\n", 100.*fVaFixed);
    if(!isnan(pWaterTrue)) printf("pWaterTrue := %g mL/mL\n", pWaterTrue);
    if(!isnan(k2min)) printf("k2min := %g min-1\n", k2min);
    if(!isnan(k2max)) printf("k2max := %g min-1\n", k2max);
    printf("flowMin := %g mL/(min*mL)\n", flowMin);
    printf("flowMax := %g mL/(min*mL)\n", flowMax);
    printf("pWaterMin := %g mL/mL\n", pWaterMin);
    printf("pWaterMax := %g mL/mL\n", pWaterMax);
    printf("bfNr := %d\n", bfNr);
  }
 
  /* Check user-defined parameter ranges and calculate range of k2 */
  if(flowMin>=flowMax) {
    fprintf(stderr, "Error: invalid range for perfusion (%g - %g).\n",
    100.*flowMin, 100.*flowMax);
    return(1);
  }
  if(pWaterMin>=pWaterMax) {
    fprintf(stderr, "Error: invalid range for p (%g - %g).\n",
    pWaterMin, pWaterMax);
    return(1);
  }
  if(isnan(k2min)) {
    k2min=flowMin/pWaterMax;
    if(verbose>1) printf("k2min := %g min-1\n", k2min);
  }
  if(isnan(k2max)) {
    k2max=flowMax/pWaterMin;
    if(verbose>1) printf("k2max := %g min-1\n", k2max);
  }
  if(k2max<=k2min || k2min<=0.) {
    fprintf(stderr, "Error: invalid range for k2 (%g - %g).\n", k2min, k2max);
    return(1);
  }
  /* Convert k2 from 1/min to 1/sec */
  k2min/=60.; k2max/=60.0;
  if(verbose>2) {
    printf("k2min := %g sec-1\n", k2min);
    printf("k2max := %g sec-1\n", k2max);
  }
 
 
  /*
   *  Read tissue and input data
   */
  if(verbose>1) printf("reading tissue and input data\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC btac, ttac; tacInit(&btac); tacInit(&ttac);
  int fitSampleNr;
  fitdur/=60.0;
  ret=tacReadModelingData(ttacfile, btacfile, NULL, NULL, &fitdur, 0,
                          &fitSampleNr, &ttac, &btac, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  fitdur*=60.0;
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(ttac.format));
    printf("tacNr := %d\n", ttac.tacNr);
    printf("ttac.sampleNr := %d\n", ttac.sampleNr);
    printf("btac.sampleNr := %d\n", btac.sampleNr);
    printf("fitSampleNr := %d\n", fitSampleNr);
    printf("xunit := %s\n", unitName(ttac.tunit));
    printf("yunit := %s\n", unitName(ttac.cunit));
    printf("fitdur := %g s\n", fitdur);
    int i=iftFindKey(&ttac.h, "studynr", 0);
    if(i<0) i=iftFindKey(&ttac.h, "study_number", 0);
    if(i>=0) printf("%s := %s\n", ttac.h.item[i].key, ttac.h.item[i].value);
 
  }
  if(fitSampleNr<4 || btac.sampleNr<4) {
    fprintf(stderr, "Error: too few samples in specified fit duration.\n");
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  /* Sample times should now be in minutes, change them to sec */
  tacXUnitConvert(&ttac, UNIT_SEC, NULL);
  tacXUnitConvert(&btac, UNIT_SEC, NULL);
  /* Check weighting */
  if(!tacIsWeighted(&ttac) && verbose>0) fprintf(stderr, "Note: data is not weighted.\n");
 
  /* Interpolate BTAC to TTAC sample times; needed for Va correction */
  if(verbose>1) printf("interpolating BTAC to TTAC sample times\n");
  TAC tbtac; tacInit(&tbtac);
  ret=tacInterpolate(&btac, &ttac, &tbtac, NULL, NULL, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); return(2);
  }
 
  /* If fixed Va was given by user, then subtract it from the tissue data */
  if(!isnan(fVaFixed) && fVaFixed>0.0) {
    if(verbose>1) printf("subtracting Va*BTAC from TTACs\n");
    for(int i=0; i<fitSampleNr; i++) {
      if(isnan(tbtac.c[0].y[i])) continue;
      for(int j=0; j<ttac.tacNr; j++) {
        if(isnan(ttac.c[j].y[i])) continue;
        ttac.c[j].y[i]-=fVaFixed*tbtac.c[0].y[i];
      }
    }
  }
 
 
  /*
   *  Calculate the basis functions
   */
  if(verbose>1) printf("calculating basis functions\n");
  TAC bf; tacInit(&bf);
  ret=bfm1TCM(&btac, &ttac, bfNr, k2min, k2max, 0, &bf, &status);
  if(ret!=TPCERROR_OK) {
    if(verbose>1) fprintf(stderr, "Error: cannot calculate basis functions.\n");
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); return(3);
  }
  /* Save basis functions if required */
  if(bffile[0]) {
    if(verbose>1) printf("writing %s\n", bffile);
    FILE *fp; fp=fopen(bffile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", bffile);
      tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf); 
      return(3);
    }
    ret=tacWrite(&bf, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf); 
      return(3);
    }
    if(verbose>0) printf("basis functions saved in %s.\n", bffile);
  }
 
 
 
  /*
   *  Prepare the room for results
   */
  if(verbose>1) printf("initializing result data\n");
  PAR par; parInit(&par);
  ret=parAllocateWithTAC(&par, &ttac, 3, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf); 
    return(4);
  }
  /* Copy titles & filenames */
  {
    int i;
    char buf[256];
    time_t t=time(NULL);
    /* set program name */
    tpcProgramName(argv[0], 1, 1, buf, 256);
    iftPut(&par.h, "program", buf, 0, NULL);
    /* set file names */
    iftPut(&par.h, "bloodfile", btacfile, 0, NULL);
    iftPut(&par.h, "datafile", ttacfile, 0, NULL);
    /* Fit method */
    iftPut(&par.h, "fitmethod", "BFM", 0, NULL);
    /* Constants */
    if(fVaFixed>=0.0) {
      sprintf(buf, "%g %%", 100.0*fVaFixed);
      iftPut(&par.h, "Vb", buf, 0, NULL);
    }
    if(pWaterTrue>0) {
      sprintf(buf, "%g", pWaterTrue);
      iftPut(&par.h, "pH2O", buf, 0, NULL);
    }
    /* Set current time to results */
    iftPut(&par.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
    /* Set fit times for each TAC */
    for(i=0; i<par.tacNr; i++) {
      par.r[i].dataNr=fitSampleNr;
      par.r[i].start=0.0;
      par.r[i].end=fitdur/60.0;
      /* and nr of fitted parameters */
      par.r[i].fitNr=3; if(fVaFixed>0.0) par.r[i].fitNr--;
    }
    /* Set the parameter names and units */
    i=0; strcpy(par.n[i].name, "Flow"); 
    if(per_dl==0) par.n[i].unit=UNIT_ML_PER_ML_MIN;
    else par.n[i].unit=UNIT_ML_PER_DL_MIN;
    i++; strcpy(par.n[i].name, "pWater"); par.n[i].unit=UNIT_ML_PER_ML;
    i++; strcpy(par.n[i].name, "Va"); 
    if(per_dl==0) par.n[i].unit=UNIT_ML_PER_ML;
    else par.n[i].unit=UNIT_PERCENTAGE;
  }
 
 
  /*
   *  Prepare the room for fitted TTACs, if requested
   */
  TAC ftac; tacInit(&ftac);
  if(fitfile[0] || svgfile[0]) {
    if(verbose>1) printf("allocating space for fitted TTACs\n");
    ret=tacDuplicate(&ttac, &ftac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot allocate space for fitted TACs.\n");
      tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf);
      parFree(&par); 
      return(4);
    }
    ftac.sampleNr=fitSampleNr;
  }
 
  
  /*
   *  Allocate memory for QR
   */
  if(verbose>1) printf("allocating memory for QR\n");
  int colNr, rowNr;
  colNr=2; if(fVaFixed>0.0) colNr--;
  rowNr=fitSampleNr;
  if(verbose>2) {
    printf("QR_colNr := %d\n", colNr);
    printf("QR_rowNr := %d\n", rowNr);
  }
  double *buf, **mat, *rhs, *sol, r2;
  buf=(double*)calloc(colNr*rowNr+rowNr+colNr, sizeof(double));
  mat=(double**)calloc(rowNr, sizeof(double*));
  if(buf==NULL || mat==NULL) {
    fprintf(stderr, "Error: cannot allocate memory for QR\n");
    tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf); 
    tacFree(&ftac); parFree(&par);
    return(6);
  }
  for(int i=0; i<rowNr; i++) mat[i]=buf+(i*colNr);
  rhs=buf+(rowNr*colNr); sol=buf+(rowNr*colNr+rowNr);
 
  /*
   *  BF fitting to each regional TAC
   */
  if(verbose>1) printf("BFM fitting to TACs.\n");
  double r2_min;
  int bi, bi_min;
  for(int i=0; i<ttac.tacNr; i++) {
 
    if(verbose>1 && ttac.tacNr>1) printf("Region %d %s\n", i+1, ttac.c[i].name);
 
    /* Go through all basis functions */
    bi_min=-1; r2_min=nan("");
    for(bi=0; bi<bf.tacNr; bi++) {
 
      if(verbose>5) printf("bi=%d\n", bi);
 
      /* Initiate matrix */
      for(int j=0; j<rowNr; j++) {
        mat[j][0]=bf.c[bi].y[j];
        if(colNr>1) mat[j][1]=tbtac.c[0].y[j];
        rhs[j]=ttac.c[i].y[j];
      }
      /* Apply data weights */
      if(verbose>5) printf("  weighting\n");
      /* Preallocate temp memory at some point */
      if(tacIsWeighted(&ttac)) qrWeight(colNr, rowNr, mat, rhs, ttac.w, NULL);
 
      /* Compute QR */
      if(verbose>5) printf("  QR\n");
      ret=qrLSQ(mat, rhs, sol, rowNr, colNr, &r2);
      if(ret!=0) {
        fprintf(stderr, "Error: no QR solution for BFM\n");
        tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf);
        tacFree(&ftac); parFree(&par);
        free(buf); free(mat);
        return(6);
      }
      if(verbose>5) {
        printf("solution (%d):", bi); 
        for(int j=0; j<colNr; j++) printf(" %g", sol[j]);
        printf("\nR^2= %g  R=%g\n", r2, sqrt(r2));
      }
      /* Check if this was best fit for now */
      if(isnan(r2_min) || r2_min>r2) {
        r2_min=r2; bi_min=bi;
      }
    } /* next basis function */
    if(verbose>2) printf("Min basis function nr %d with R2=%g\n", bi_min+1, r2_min);
 
    /* Compute the best BF again to retrieve the solution */
    bi=bi_min;
    /* Initiate matrix */
    for(int j=0; j<rowNr; j++) {
      mat[j][0]=bf.c[bi].y[j];
      if(colNr>1) mat[j][1]=tbtac.c[0].y[j];
      rhs[j]=ttac.c[i].y[j];
    }
    /* Apply data weights */
    if(verbose>5) printf("  weighting\n");
    /* Preallocate temp memory at some point */
    if(tacIsWeighted(&ttac)) qrWeight(colNr, rowNr, mat, rhs, ttac.w, NULL);
    /* Compute QR */
    if(verbose>5) printf("  QR\n");
    ret=qrLSQ(mat, rhs, sol, rowNr, colNr, &r2);
    if(ret!=0) {
      fprintf(stderr, "Error: no QR solution for BFM\n");
      tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf);
      tacFree(&ftac); parFree(&par);
      free(buf); free(mat);
      return(6);
    }
    par.r[i].wss=r2;
    if(verbose>5) {
      printf("best_solution (%d):", bi); 
      for(int j=0; j<colNr; j++) printf(" %g", sol[j]); 
      printf("\nR^2= %g  R=%g\n", r2, sqrt(r2));
    }
    /* Remove weights */
    if(tacIsWeighted(&ttac)) qrWeightRm(colNr, rowNr, mat, rhs, ttac.w, NULL);
    /* Fitted TTAC is in rhs[] */
    if(fitfile[0] || svgfile[0]) {
      for(int j=0; j<ftac.sampleNr; j++) ftac.c[i].y[j]=rhs[j];
    }
    /* Solve final parameters Flow, pH2O, and Va */
    double Flow, pH2O, Va, k2;
    k2=bf.c[bi_min].size; //printf("k2=%g\n", k2);
    if(isnan(fVaFixed)) Va=sol[1]; else Va=fVaFixed; //printf("Va=%g\n", Va);
    if(pWaterTrue>0) {
      /* True p for radiowater is given, therefore lets calculate flow from k2 */
      Flow=k2*pWaterTrue;
    } else {
      /* True p for radiowater is not given, therefore lets assume that flow=K1 */
      Flow=sol[0]/(1.0-Va);
    }
    // printf("Flow=%g\n", Flow);
    /* Apparent p for radiowater is K1/k2 */
    pH2O=sol[0]/((1.0-Va)*k2); //printf("pH2O=%g\n", pH2O);
 
    /* Put the final parameters in place */
    if(flow_per_perfusable_tissue==0) Flow*=(1.0-Va);
    if(per_dl!=0) {Flow*=100.0; Va*=100;}
    par.r[i].p[0]=60.0*Flow; // 1/sec -> 1/min
    par.r[i].p[1]=pH2O;
    par.r[i].p[2]=Va;
 
  } /* next region */
 
  /* Free the memory allocated for QR */
  free(buf); free(mat);
 
 
  /*
   *  Print results on screen
   */
  if(verbose>0 && par.tacNr<50) parWrite(&par, stdout, PAR_FORMAT_RES, 0, &status);
 
  /*
   *  Save results
   */
  if(resfile[0]) {
    if(verbose>1) printf("writing %s\n", resfile);
    par.format=parFormatFromExtension(resfile);
    if(verbose>2)
      printf("result file format := %s\n", parFormattxt(par.format));
    if(par.format==PAR_FORMAT_UNKNOWN) par.format=PAR_FORMAT_TSV_UK;
    FILE *fp; fp=fopen(resfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", resfile);
      tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf);
      tacFree(&ftac); parFree(&par);
      return(11);
    }
    ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    //parFree(&par);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf);
      tacFree(&ftac); parFree(&par);
      return(12);
    }
    if(verbose>0) printf("Results saved in %s.\n", resfile);
  }
 
 
  /*
   *  SVG plot of fitted and original data
   */
  if(svgfile[0]) {
 
    if(verbose>1) printf("saving SVG plot\n");
    int i;
    char buf[128];
    sprintf(buf, "Radiowater BFM fit");
    i=iftFindKey(&ttac.h, "studynr", 0);
    if(i<0) i=iftFindKey(&ttac.h, "study_number", 0);
    if(i>=0) {strcat(buf, ": "); strcat(buf, ttac.h.item[i].value);}
    ret=tacPlotFitSVG(&ttac, &ftac, buf, 0.0, nan(""), 0.0, nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf);
      tacFree(&ftac); parFree(&par);
      return(21);
    }
    if(verbose>0) printf("Plots written in %s.\n", svgfile);
  }
 
  /*
   *  Save fitted TTACs
   */
  if(fitfile[0]) {
    if(verbose>1) printf("writing %s\n", fitfile);
    FILE *fp; fp=fopen(fitfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", fitfile);
      tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf);
      tacFree(&ftac); parFree(&par);
      return(31);
    }
    ret=tacWrite(&ftac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf);
      tacFree(&ftac); parFree(&par);
      return(32);
    }
    if(verbose>0) printf("fitted TACs saved in %s.\n", fitfile);
  }
 
 
  tacFree(&ttac); tacFree(&btac); tacFree(&tbtac); tacFree(&bf);
  tacFree(&ftac); parFree(&par);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/