bfmsrtm.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[] = {
  "Estimate the binding potential (BPnd) from simplified reference tissue",
  "model (SRTM) (Lammertsma & Hume, 1996). The model is solved using",
  "the basis function method (BFM; Gunn et al., 1997).",
  " ",
  "Radioactivity concentration in region(s)-of-interest is given in TTAC file;",
  "data must be corrected for physical decay.",
  "Radioactivity concentration in reference region (RTAC) is given in separate",
  "file, or as name or number of the reference region inside the TTAC file.",
  "Sample times must be in minutes in all data files, unless specified inside",
  "the files. TTAC file should include weights.",
  " ",
  "Usage: @P [options] TTAC RTAC results",
  " ",
  "Options:",
  " -i=<isotope>",
  "     To apply BFM to non-decay-corrected data, as proposed in the original",
  "     method publication, enter the isotope with this option, using codes",
  "     C-11, F-18, Ga-68, Cu-64, ...",
  " -t3min=<value (1/min)>",
  "     Set minimum value for theta3; it must be >= k2min/(1+BPmax)+lambda;",
  "     it has to be > lambda. Default is lambda+0.001 min-1.",
  "     Lambda for F-18 is 0.0063 and for C-11 0.034;",
  "     If isotope is not given, then lambda=0.",
  " -t3max=<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 5000.",
  " -end=<Fit end time (min)>",
  "     Use data from 0 to end time; by default, model is fitted to all frames.",
  " -DVR",
  "     Instead of BPnd, program saves the DVR (=BPnd+1) values.",
  " -w1",
  "     All weights are set to 1.0 (no weighting); by default, weights in",
  "     data file are used, if available.",
  " -wf",
  "     Weight by sampling interval.",
  " -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.",
  " -bf=<filename>",
  "     Basis functions are written in specified DFT file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P a4567.tac cer C-11 a4567bp.res",
  " ",
  "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: tacweigh, tacdecay, logan, imgbfbp, fit_srtm, lhsrtm, fitk2",
  " ",
  "Keywords: TAC, 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;
  char    rtacfile[FILENAME_MAX], ttacfile[FILENAME_MAX], resfile[FILENAME_MAX],
          fitfile[FILENAME_MAX], svgfile[FILENAME_MAX], bffile[FILENAME_MAX];
  int     bfNr=5000;
  int     bp_plus_one=0; // 0=BPnd, 1=DVR
  double  t3min=nan(""), t3max=nan(""); // 1/min
  int     isotope=ISOTOPE_UNKNOWN;
  double  fitdur=nan("");
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  char   *cptr;
  int     ret;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  rtacfile[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, "DVR")==0) {
      bp_plus_one=1; continue;
    } else if(strncasecmp(cptr, "BPND", 2)==0) {
      bp_plus_one=0; 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, "min=", 4)==0) {
      if(atofCheck(cptr+4, &t3min)==0 && t3min>=0.0) continue;
    } else if(strncasecmp(cptr, "max=", 4)==0) {
      if(atofCheck(cptr+4, &t3max)==0 && t3max>=0.0) continue;
    } else if(strncasecmp(cptr, "t3min=", 6)==0) {
      if(atofCheck(cptr+6, &t3min)==0 && t3min>=0.0) continue;
    } else if(strncasecmp(cptr, "t3max=", 6)==0) {
      if(atofCheck(cptr+6, &t3max)==0 && t3max>=0.0) continue;
    } else if(strncasecmp(cptr, "end=", 4)==0) {
      if(atofCheck(cptr+4, &fitdur)==0) {
        if(fitdur<=0.0) fitdur=1.0E+99;
        continue;
      } else {
        fprintf(stderr, "Error: invalid fit time '%s'.\n", argv[ai]);
        return(1);
      }
    } else if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "I=", 2)==0) {
      isotope=isotopeIdentify(cptr+2);
      if(isotope==ISOTOPE_UNKNOWN) {
        fprintf(stderr, "Error: invalid isotope code '%s'\n", cptr+2);
        return(1);
      }
      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(ttacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(rtacfile, argv[ai++], FILENAME_MAX);
  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);
  }
  if(isnan(fitdur)) fitdur=1.0E+99;
  /* Set or check theta3 limits */
  double lambda=0.0;
  if(isotope!=ISOTOPE_UNKNOWN)
    lambda=lambdaFromHalflife(isotopeHalflife(isotope));
  if(verbose>2) printf("lambda := %g\n", lambda);
  if(isnan(t3min)) t3min=lambda+0.001;
  if(isnan(t3max)) t3max=0.60;
  if(t3min>=t3max || t3min<=lambda) {
    fprintf(stderr, "Error: invalid theta3 bounds (%g - %g).\n", t3min, t3max); 
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ttacfile := %s\n", ttacfile);
    printf("rtacfile := %s\n", rtacfile);
    printf("resfile := %s\n", resfile);
    printf("fitfile := %s\n", fitfile);
    printf("svgfile := %s\n", svgfile);
    printf("bffile := %s\n", bffile);
    printf("bfNr := %d\n", bfNr);
    printf("bp_plus_one := %d\n", bp_plus_one);
    /*if(!isnan(t3min))*/ printf("t3min := %g min-1\n", t3min);
    /*if(!isnan(t3max))*/ printf("t3max := %g min-1\n", t3max);
    if(isotope!=ISOTOPE_UNKNOWN) 
      printf("isotope := %s\n", isotopeName(isotope));
    printf("required_fittime := %g min\n", fitdur);
    printf("weights := %d\n", weights);
  }
 
 
  /*
   *  Read TAC files
   */
  TAC ttac, rtac; tacInit(&ttac); tacInit(&rtac);
 
  /* Read tissue TAC file */
  if(verbose>1) printf("reading %s\n", ttacfile);
  ret=tacRead(&ttac, ttacfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(ttac.format));
    printf("tacNr := %d\n", ttac.tacNr);
    printf("sampleNr := %d\n", ttac.sampleNr);
    printf("xunit := %s\n", unitName(ttac.tunit));
    printf("yunit := %s\n", unitName(ttac.cunit));
  }
  /* Set isotope into TTAC header */
  if(isotope!=ISOTOPE_UNKNOWN) tacSetIsotope(&ttac, isotope);
 
  /* Read reference tissue TAC file */
  if(verbose>1) printf("reading %s\n", rtacfile);
  int refindex;
  ret=tacReadReference(&ttac, rtacfile, &rtac, &refindex, &status);
  if(ret==0) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&rtac); return(2);
  }
  if(rtac.tacNr>1) {
    fprintf(stderr, "Warning: several reference regions found: %s selected.\n",
            rtac.c[refindex].name);
  } else if(verbose>1) {
    printf("reference_region := %s\n", rtac.c[refindex].name);
  }
  if(verbose>2) {
    printf("ref_fileformat := %s\n", tacFormattxt(rtac.format));
    printf("ref_tacNr := %d\n", rtac.tacNr);
    printf("ref_sampleNr := %d\n", rtac.sampleNr);
    printf("ref_xunit := %s\n", unitName(rtac.tunit));
    printf("ref_yunit := %s\n", unitName(rtac.cunit));
  }
 
  /* Try to fix missing values, if any */
  if(tacFixNaNs(&ttac) || tacFixNaNs(&rtac)) {
    fprintf(stderr, "Error: missing sample(s) in data.\n");
    tacFree(&ttac); tacFree(&rtac); return(2);
  }
 
  /* Convert sample times to minutes */
  if(tacXUnitConvert(&ttac, UNIT_MIN, &status) ||
     tacXUnitConvert(&rtac, UNIT_MIN, &status))
  {
    fprintf(stderr, "Warning: %s\n", errorMsg(status.error));
    // Do not stop in error, assume that user has set units correctly
  }
 
  /* Remove small frame overlaps and gaps */
  if(tacCorrectFrameOverlap(&ttac, &status) || tacCorrectFrameOverlap(&rtac, &status)) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&rtac); return(2);
  }
 
  /* Add data weights, if requested */
  if(weights==1) {
    ttac.weighting=WEIGHTING_OFF; 
    for(int i=0; i<ttac.sampleNr; i++) ttac.w[i]=1.0;
  } else if(weights==2) {
    if(tacWByFreq(&ttac, ISOTOPE_UNKNOWN, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&rtac); return(2);
    }
  } else if(!tacIsWeighted(&ttac)) {
    fprintf(stderr, "Warning: data is not weighted.\n");
  }
 
  /* Set fit duration */
  double starttime=0.0;
  double endtime=fitdur;
  int fitSampleNr;
  fitSampleNr=tacFittime(&ttac, &starttime, &endtime, NULL, NULL, &status);
  if(fitSampleNr<0) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&rtac); return(2);
  } else if(fitSampleNr<4) {
    fprintf(stderr, "Error: too few data points for a decent fit.\n");
    tacFree(&ttac); tacFree(&rtac); return(2);
  }
  if(verbose>1) {
    printf("starttime := %g\n", starttime);
    printf("endtime := %g\n", endtime);
    //printf("first := %d\n", first);
    //printf("last := %d\n", last);
    printf("fitSampleNr := %d\n", fitSampleNr);
  }
  fitdur=endtime;
 
  
  /*
   *  Remove decay correction from TTACs and RTAC(s), if requested
   */
  if(isotope!=ISOTOPE_UNKNOWN) {
    if(verbose>1) printf("removing decay correction.\n");
    ret=tacDecayCorrection(&ttac, isotope, 0, &status);
    if(ret==TPCERROR_OK) ret=tacDecayCorrection(&rtac, isotope, 0, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&rtac); return(2);
    }
  }
 
  /*
   *  Calculate the basis functions
   */
  if(verbose>1) printf("integrate reference region TAC for BF calculation\n");
  TAC irtac; tacInit(&irtac);
  ret=tacInterpolate(&rtac, &ttac, NULL, &irtac, NULL, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&rtac); tacFree(&irtac); return(3);
  }
 
  if(verbose>1) printf("calculating basis functions\n");
  TAC bf; tacInit(&bf);
  ret=bfmSRTM(ttac.x, irtac.c[refindex].y, fitSampleNr, bfNr, 
              t3min, t3max, &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(&rtac); tacFree(&irtac); return(3);
  }
  tacFree(&irtac); // not needed anymore
  bf.tunit=ttac.tunit;
 
  /* 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(&rtac); tacFree(&bf); 
      return(3);
    }
    ret=tacWrite(&bf, fp, TAC_FORMAT_PMOD, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&rtac); 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);
  if(verbose<2) ret=parAllocateWithTAC(&par, &ttac, 3, &status);
  else ret=parAllocateWithTAC(&par, &ttac, 5, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&rtac); 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, "datafile", ttacfile, 0, NULL);
    iftPut(&par.h, "refname", rtac.c[refindex].name, 0, NULL);
    /* set isotope, if given */
    if(isotope!=ISOTOPE_UNKNOWN) 
      iftPut(&par.h, "isotope", isotopeName(isotope), 0, NULL);
    /* Fit method */
    iftPut(&par.h, "fitmethod", "BFM", 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=starttime;
      par.r[i].end=fitdur;
      /* and nr of fitted parameters */
      par.r[i].fitNr=3;
    }
    /* Set the parameter names and units */
    i=0; strcpy(par.n[i].name, "R1"); par.n[i].unit=UNIT_UNITLESS;
    i++; strcpy(par.n[i].name, "k2"); par.n[i].unit=UNIT_PER_MIN;
    i++; par.n[i].unit=UNIT_UNITLESS;
    if(bp_plus_one==0) strcpy(par.n[i].name, "BPnd");
    else strcpy(par.n[i].name, "DVR");
    if(par.parNr==5) {
      i++; strcpy(par.n[i].name, "theta2"); par.n[i].unit=UNIT_UNITLESS;
      i++; strcpy(par.n[i].name, "theta3"); par.n[i].unit=UNIT_UNITLESS;
    }
  }
 
 
 
  /*
   *  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(&rtac); 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;
  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(&rtac); 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]=rtac.c[refindex].y[j];
        mat[j][1]=bf.c[bi].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(&rtac); 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) && bf.c[bi].size!=lambda) {
        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]=rtac.c[refindex].y[j];
      mat[j][1]=bf.c[bi].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(&rtac); 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<rowNr; j++) ftac.c[i].y[j]=rhs[j];
    }
 
    /* Solve final parameters */
    double R1, k2, BPnd, DVR, theta2, theta3;
    R1=sol[0]; theta2=sol[1]; theta3=bf.c[bi_min].size;
    k2=R1*(theta3-lambda) + theta2;
    DVR= R1 + theta2/(theta3-lambda); BPnd=DVR-1.0;
 
    /* Put the final parameters in place */
    par.r[i].p[0]=R1;
    par.r[i].p[1]=k2;
    if(bp_plus_one==0) par.r[i].p[2]=BPnd; else par.r[i].p[2]=DVR;
    if(par.parNr==5) {
      par.r[i].p[3]=theta2;
      par.r[i].p[4]=theta3;
    }
 
  } /* 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(&rtac); 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(&rtac); 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, "SRTM-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(&rtac); tacFree(&bf);
      tacFree(&ftac); parFree(&par);
      return(21);
    }
    if(verbose>0) printf("Plots written in %s.\n", svgfile);
  }
 
  /*
   *  Save fitted TTACs
   */
  if(fitfile[0]) {
    /*
     *  Return decay correction, if needed
     */
    if(isotope!=ISOTOPE_UNKNOWN) {
      if(verbose>1) printf("decay correction for fitted TTACs.\n");
      ret=tacDecayCorrection(&ftac, isotope, 1, &status);
      if(ret==TPCERROR_OK) ret=tacDecayCorrection(&rtac, isotope, 0, &status);
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&ttac); tacFree(&rtac); return(31);
      }
    }
    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(&rtac); tacFree(&bf);
      tacFree(&ftac); parFree(&par);
      return(32);
    }
    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(&rtac); tacFree(&bf);
      tacFree(&ftac); parFree(&par);
      return(33);
    }
    if(verbose>0) printf("fitted TACs saved in %s.\n", fitfile);
  }
 
 
  tacFree(&ttac); tacFree(&rtac); tacFree(&bf); tacFree(&ftac); 
  parFree(&par);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/