fitdelay.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 "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
/*****************************************************************************/
#define MAX_DELAY 60
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "For estimation and correction of the delay-time (difference in appearance",
  "times of radioactivity) between PET tissue and input (blood or plasma) TACs.",
  " ",
  "Program is based on the methods used by Meyer et al. (1) and",
  "van den Hoff et al. (2):",
  "The plasma/blood curve is shifted -60 - +60 sec, and a two-tissue",
  "compartment model (with parameters K1, k2, k3, k4 and Vb) in multilinear",
  "form (3) is fitted to the shifted TAC and each regional tissue TAC,",
  "with the nonnegative least squares method (4).",
  "For each region, the delay leading to the lowest sum-of-squares is selected;",
  "the over-all delay value is calculated as a median of the regional delays.",
  "Dispersion is not considered in this application.",
  " ",
  "Usage: @P [options] inputfile tissuefile fittime [inputfile2 [inputfile3 [inputfile4]]]",
  " ",
  "Options:",
  " -o=<Filename>",
  "     Filename for the time delay corrected TAC made from inputfile.",
  " -o2=<Filename>",
  "     Filename for the time delay corrected TAC made from inputfile2.",
  " -o3=<Filename>",
  "     Filename for the time delay corrected TAC made from inputfile3.",
  " -o4=<Filename>",
  "     Filename for the time delay corrected TAC made from inputfile4.",
  " -timeunit=<min|sec>",
  "     If datafile(s) do not contain the unit of sample times, it is",
  "     recommended to specify it with this option. By default, units in data",
  "     files are trusted.",
  " -format=<none|dft|pmod|if>",
  "     Specify the output data format; none means that no title lines are saved.",
  " -fit=<Filename>",
  "     Filename for fitted best TACs; by default these are not saved.",
  " -model=<1|2>",
  "     Select whether 1- or 2-tissue (default) compartment model is applied.",
  " -L[og]",
  "     Time delay and other log information is written as comments in",
  "     the corrected TAC file, if format supports comments.",
  " -matrix=<Filename>",
  "     Filename for saving NNLS matrix in CSV format for testing purposes.",
  "     With this option the tissue file must contain one TAC only.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "As tissue data, the scanner countrate curve is recommended, unless scanned",
  "volume contains heart or large artery or vein where tracer was injected;",
  "It may be possible to use also regional TACs, if datafile contains frame",
  "start and end times. If tissue data contains background, remove it first",
  "with dftrmbkg. The units of sample times should be specified in datafiles;",
  "file format is specified in (5).",
  " ",
  "Fittime must be given in seconds.",
  " ",
  "Estimated tracer appearance times in blood/plasma and tissue curves, and",
  "their differences (time delays and median time delay) are written in stdout.",
  "By default, delay corrected blood/plasma file is written with name *.delay.*",
  "but this can be changed with option -o=<Filename>.",
  "The same correction can be applied to 1-3 additional files, for",
  "example plasma metabolite TACs.",
  " ",
  "Example 1.",
  "Delay correction for C-11 or F-18 labeled tracer data, using",
  "metabolite corrected plasma curve as input and count-rate data as tissue",
  "and correcting also plasma metabolites and total blood for the delay-time:",
  "   fitdelay ut345ap_pure.kbq ut345dy1.img.cr 1800 ut345ap_met.kbq ut345ab.kbq",
  " ",
  "Example 2.",
  "Delay correction for [O-15]water data, using regional tissue",
  "curves as replacement for count-rate data:",
  "   fitdelay ut111ab.kbq ut111dy1.dft 120",
  " ",
  "References:",
  " 1. Meyer. Simultaneous correction for tracer arrival delay and dispersion",
  "    in CBF measurements by the H215O autoradiographic method and dynamic PET.",
  "    J Nucl Med 1989; 30:1069-1078.",
  " 2. van den Hoff et al. Accurate local blood flow measurements with",
  "    dynamic PET: fast determination of input function delay and dispersion",
  "    by multilinear minimization. J Nucl Med 1993; 34:1770-1777.",
  " 3. Blomqvist G. On the construction of functional maps in positron",
  "    emission tomography. J Cereb Blood Flow Metab 1984; 4:629-632.",
  " 4. Lawson CL & Hanson RJ. Solving least squares problems.",
  "    Prentice-Hall, 1974.",
  " 5. https://www.turkupetcentre.net/petanalysis/format_tpc_dft.html",
  " ",
  "See also: tactime, imghead, tacmean, tocr, dftrmbkg, tacunit, fit_h2o",
  " ",
  "Keywords: TAC, modelling, input, blood, time delay",
  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 ri, fi, fj, di, ret;
  int min;
  int make_log=0;
  int time_unit=TUNIT_UNKNOWN;
  int orig_tissue_time_unit, orig_input_time_unit;
  int output_format=DFT_FORMAT_UNKNOWN;
  int model=2;
  double v, f, ss, coeff[5], minv, maxv, length=-1.0, onep;
  double *regional_delay;
  char *cptr, tmp[FILENAME_MAX];
  char pfile[FILENAME_MAX], tfile[FILENAME_MAX], rfile[FILENAME_MAX];
  char p2file[FILENAME_MAX], p3file[FILENAME_MAX], p4file[FILENAME_MAX];
  char r2file[FILENAME_MAX], r3file[FILENAME_MAX], r4file[FILENAME_MAX];
  char ffile[FILENAME_MAX];
  char matfile[FILENAME_MAX];
  DFT pdata, data, ipdata, fdata;
  /* nnls */
  int NNLS_N=5;
  int      n, m, nnls_n, nnls_m, nnls_index[NNLS_N];
  double  *nnls_a[NNLS_N], *nnls_b, *nnls_zz, nnls_x[NNLS_N], *nnls_mat,
           nnls_wp[NNLS_N], *dptr, nnls_rnorm;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  pfile[0]=rfile[0]=tfile[0]=p2file[0]=p3file[0]=r2file[0]=r3file[0]=(char)0;
  p4file[0]=r4file[0]=ffile[0]=(char)0;
  matfile[0]=(char)0;
  dftInit(&pdata); dftInit(&data); dftInit(&ipdata); dftInit(&fdata);
  orig_tissue_time_unit=orig_input_time_unit=TUNIT_UNKNOWN;
 
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1;
    if(strcasecmp(cptr, "L")==0 || strcasecmp(cptr, "LOG")==0) {
      make_log=1; continue;
    } else if(strncasecmp(cptr, "O=", 2)==0) {
      cptr+=2; if(strlen(cptr)>0) {strcpy(rfile, cptr); continue;}
    } else if(strncasecmp(cptr, "O2=", 3)==0) {
      cptr+=3; if(strlen(cptr)>0) {strcpy(r2file, cptr); continue;}
    } else if(strncasecmp(cptr, "O3=", 3)==0) {
      cptr+=3; if(strlen(cptr)>0) {strcpy(r3file, cptr); continue;}
    } else if(strncasecmp(cptr, "O4=", 3)==0) {
      cptr+=3; if(strlen(cptr)>0) {strcpy(r4file, cptr); continue;}
    } else if(strncasecmp(cptr, "MATRIX=", 7)==0) {
      cptr+=7; if(strlen(cptr)>0) {strcpy(matfile, cptr); continue;}
    } else if(strncasecmp(cptr, "TIMEUNIT=", 9)==0) {
      cptr+=9;
      if(strncasecmp(cptr, "S", 1)==0) {time_unit=TUNIT_SEC; continue;}
      if(strncasecmp(cptr, "M", 1)==0) {time_unit=TUNIT_MIN; continue;}
    } else if(strncasecmp(cptr, "FORMAT=", 7)==0) {
      cptr+=7;
      if(strncasecmp(cptr, "DFT", 1)==0) {
        output_format=DFT_FORMAT_STANDARD; continue;}
      if(strncasecmp(cptr, "NONE", 2)==0) {
        output_format=DFT_FORMAT_PLAIN; continue;}
      if(strncasecmp(cptr, "PMOD", 2)==0) {
        output_format=DFT_FORMAT_PMOD; continue;}
      if(strcasecmp(cptr, "IF")==0) {
        output_format=DFT_FORMAT_IF; continue;}
    } else if(strncasecmp(cptr, "MODEL=", 6)==0) {
      model=atoi(cptr+6); if(model==1 || model==2) continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0) {
      cptr+=4; strcpy(ffile, cptr); if(strlen(ffile)>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 */
  for(; ai<argc; ai++) {
    if(!pfile[0]) {strcpy(pfile, argv[ai]); continue;}
    if(!tfile[0]) {strcpy(tfile, argv[ai]); continue;}
    if(length<=0.9) {
      length=atof_dpi(argv[ai]); if(length>0.0) continue;
      fprintf(stderr, "Error: invalid fit time as an argument.\n"); return(1);
    }
    if(!p2file[0]) {strcpy(p2file, argv[ai]); continue;}
    if(!p3file[0]) {strcpy(p3file, argv[ai]); continue;}
    if(!p4file[0]) {strcpy(p4file, argv[ai]); continue;}
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!tfile[0]) {
    fprintf(stderr, "Error: missing tissue file name.\n");
    return(1);
  }
  if(length<60.0) {
    fprintf(stderr, "Error: invalid fit time.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("pfile := %s\n",  pfile);
    printf("p2file := %s\n", p2file);
    printf("p3file := %s\n", p3file);
    printf("p4file := %s\n", p4file);
    printf("tfile := %s\n",  tfile);
    printf("rfile := %s\n",  rfile);
    printf("r2file := %s\n", r2file);
    printf("r3file := %s\n", r3file);
    printf("r4file := %s\n", r4file);
    printf("ffile := %s\n", ffile);
    if(matfile[0]) printf("matfile := %s\n", matfile);
    printf("output_format := %d\n", output_format);
    printf("make_log := %d\n", make_log);
    printf("time_unit := %s\n", petTunit(time_unit));
    printf("length := %g\n", length);
    printf("model := %d\n", model);
  }
 
  /* Check that extra files do exist, if file names were given */
  if(p2file[0] && access(p2file, 0)==-1) {
    fprintf(stderr, "Error: %s not found.\n", p2file); return(1);}
  if(p3file[0] && access(p3file, 0)==-1) {
    fprintf(stderr, "Error: %s not found.\n", p3file); return(1);}
  if(p4file[0] && access(p4file, 0)==-1) {
    fprintf(stderr, "Error: %s not found.\n", p4file); return(1);}
 
  /* Set NNLS_N */
  if(model==1) NNLS_N=3; else NNLS_N=5;
 
 
  /*
   *  Read data
   */
   
  /* Input */
  if(verbose>1) printf("reading %s\n", pfile);
  if(dftRead(pfile, &pdata)) {
    fprintf(stderr, "Error in reading '%s': %s\n", pfile, dfterrmsg);
    dftEmpty(&pdata); return(2);
  }
  if(pdata.voiNr>1) {
    fprintf(stderr, "Warning: %s contains %d TACs; only first is used.\n",
      pfile, pdata.voiNr);
    pdata.voiNr=1;
  }
  orig_input_time_unit=pdata.timeunit;
  if(verbose>13) dftPrint(&pdata);
  
  /* Tissue (or countrate) data */
  if(verbose>1) printf("reading %s\n", tfile);
  if(dftRead(tfile, &data)) {
    fprintf(stderr, "Error in reading '%s': %s\n", tfile, dfterrmsg);
    dftEmpty(&pdata); dftEmpty(&data); return(2);
  }
  orig_tissue_time_unit=data.timeunit;
  if(verbose>2) {
    printf("pdata_time_unit := %s\n", petTunit(pdata.timeunit));
    if(verbose>2) printf("last_t := %g\n", pdata.x[pdata.frameNr-1]);
    printf("data_time_unit := %s\n", petTunit(data.timeunit));
    if(verbose>2) printf("last_t := %g\n", data.x[data.frameNr-1]);
  }
 
  /* If more than one tissue TAC then do not allow saving NNLS matrix */
  if(data.voiNr>1 && matfile[0]) {
    matfile[0]=(char)0;
    fprintf(stderr, "Warning: option -matrix disabled because more than one tissue TAC found.\n");
  }
 
  /*
   *  Check time units and convert times to sec when necessary
   */
  if(verbose>1) printf("checking time units\n");
  /* Make sure that input time unit is set */
  if(pdata.timeunit==TUNIT_UNKNOWN) {
    if(time_unit!=TUNIT_UNKNOWN) {
      /* If user known what time unit is, then use it */
      pdata.timeunit=time_unit;
    } else { /* If not, then we have to start guessing */
      if(data.timeunit!=TUNIT_UNKNOWN) {
        /* If tissue time unit is known, then use it */
        pdata.timeunit=data.timeunit;
      } else {
        /* Guess time unit based on last sample time */
        if(pdata.x[pdata.frameNr-1]>=180) pdata.timeunit=TUNIT_SEC;
        else pdata.timeunit=TUNIT_MIN;
      }
      fprintf(stderr, "Warning: assuming input time_unit := %s\n", petTunit(pdata.timeunit));
    }
  }
  /* Make sure that tissue time unit is set */
  if(data.timeunit==TUNIT_UNKNOWN) {
    if(time_unit!=TUNIT_UNKNOWN) {
      /* If user known what time unit is, then use it */
      data.timeunit=time_unit;
    } else { /* If not, then we have to start guessing */
      if(pdata.timeunit!=TUNIT_UNKNOWN) {
        /* If input time unit was known, then use it */
        data.timeunit=pdata.timeunit;
      } else {
        /* Guess time unit based on last sample time */
        if(data.x[data.frameNr-1]>=180) data.timeunit=TUNIT_SEC;
        else data.timeunit=TUNIT_MIN;
      }
      fprintf(stderr, "Warning: assuming tissue time_unit := %s\n", petTunit(data.timeunit));
    }
  }
  /* Save the original unit, so that units can be converted back */
  orig_input_time_unit=pdata.timeunit;
  orig_tissue_time_unit=data.timeunit;
  if(verbose>2) {
    printf("orig_input_time_unit := %s\n", petTunit(orig_input_time_unit));
    printf("orig_tissue_time_unit := %s\n", petTunit(orig_tissue_time_unit));
  }
  /* Convert data times to sec */
  if(data.timeunit!=TUNIT_SEC) {
    if(verbose>0)
      fprintf(stdout, "Note: tissue sample times are converted to seconds for fitting.\n");
    dftTimeunitConversion(&data, TUNIT_SEC);
  }  
  if(pdata.timeunit!=TUNIT_SEC) {
    if(verbose>0)
      fprintf(stdout, "Note: input sample times are converted to seconds for fitting.\n");
    dftTimeunitConversion(&pdata, TUNIT_SEC);    
  }
  /* Check fit time */
  if(length<data.x[data.frameNr-1]/60.0 || length<pdata.x[pdata.frameNr-1]/60.0)
    fprintf(stderr, "Warning: fit time (%g s) is short, compared with data.\n", length);
 
 
  /* Check if the ascending part of PTAC has been missed */
  if(verbose>1) printf("checking if the ascending part of plasma TAC has been missed\n");
  for(fi=1, n=0, maxv=pdata.voi[0].y[0]; fi<pdata.frameNr; fi++) {
    v=pdata.voi[0].y[fi];
    if(fi>1 && fi<pdata.frameNr-2) { 
      // mean of samples to avoid noise effects
      v+=pdata.voi[0].y[fi-1]+pdata.voi[0].y[fi+1];
      v+=pdata.voi[0].y[fi-2]+pdata.voi[0].y[fi+2];
      v/=5.0;
    } else if(fi>0 && fi<pdata.frameNr-1) { 
      // mean of samples to avoid noise effects
      v+=pdata.voi[0].y[fi-1]+pdata.voi[0].y[fi+1];
      v/=3.0;
    }
    if(v>maxv) {maxv=v; n=fi;}
  }
  if(verbose>2) printf("maxv := %g\nmax_index := %d\nmax_time := %g\n", maxv, n, pdata.x[n]);
  if(n==0) {  /* the first or second value is the highest */
    fprintf(stderr, "Error: missed the ascending phase of plasma/blood data.\n");
    dftEmpty(&pdata); dftEmpty(&data); return(2);
  }
  if(pdata.x[n]>length) {  /* peak after fit end time */
    fprintf(stderr, 
      "Error: missed the plasma/blood peak; check the time unit.\n");
    dftEmpty(&pdata); dftEmpty(&data); return(2);
  }
  /* check tissue data (first TAC) */
  if(verbose>1) printf("checking if the ascending part of tissue TAC has been missed\n");
  for(fi=1, n=0, maxv=data.voi[0].y[0]; fi<data.frameNr; fi++)
    if(data.voi[0].y[fi]>maxv) {maxv=data.voi[0].y[fi]; n=fi;}
  if(verbose>2) printf("maxv := %g\nmax_index := %d\n", maxv, n);
  if(n==0) {
    fprintf(stderr, "Error: missed the ascending phase of tissue data.\n");
    if(verbose>12) dftPrint(&data);
    dftEmpty(&pdata); dftEmpty(&data); return(2);
  }
  if(n<2) fprintf(stderr, "Warning: check the first samples in %s\n", tfile);
 
 
  /* "Remove" tissue data after fit length */
  if(data.timetype==DFT_TIME_STARTEND) {
    if(length>data.x2[data.frameNr-1]) length=data.x2[data.frameNr-1];
  } else {
    if(length>data.x[data.frameNr-1]) length=data.x[data.frameNr-1];
  }
  for(fi=0; fi<data.frameNr; fi++) {
    if(data.timetype==DFT_TIME_STARTEND) {if(data.x1[fi]>=length) break;}
    else {if(data.x[fi]>=length) break;}
  } fi--;
  if(fi<5) {
    fprintf(stderr, "Error: too few time frames included in fit.\n");
    dftEmpty(&data); dftEmpty(&pdata); return(2);
  }
  if(fi<data.frameNr) {
    data.frameNr=fi;
  } else {
    if(data.timetype==3) length=data.x2[data.frameNr-1];
    else length=data.x[data.frameNr-1];
  }
  if(verbose>3) printf("fit_frameNr := %d\nfit_length := %g\n", data.frameNr, length);
  if(verbose>2) {
    if(verbose>10) dftPrint(&data);
    printf("pdata_time_unit := %s\n", petTunit(pdata.timeunit));
    if(verbose>2) printf("last_t := %g\n", pdata.x[pdata.frameNr-1]);
    printf("data_time_unit := %s\n", petTunit(data.timeunit));
    if(verbose>2) printf("last_t := %g\n", data.x[data.frameNr-1]);
  }
 
  /* Check how many plasma/blood samples there are during fit time */
  for(fi=0; fi<pdata.frameNr; fi++) if(pdata.x[fi]>length) break;
  if(verbose>1) printf("nr of input samples in fit range := %d\n", fi);
  if(fi<5 || (fi<10 && fi<data.frameNr)) {
    fprintf(stderr, "Error: too few plasma/blood samples included in fit.\n");
    dftEmpty(&pdata); dftEmpty(&data); return(2);
  }
 
 
  /* Calculate tissue integrals at frame mid times */
  if(verbose>1) printf("integrating tissue TACs\n");
  for(ri=0; ri<data.voiNr; ri++) {
    if(data.timetype==3)
      ret=petintegral(data.x1, data.x2, data.voi[ri].y, data.frameNr,
                      data.voi[ri].y2, data.voi[ri].y3);
      else
        ret=interpolate(data.x, data.voi[ri].y, data.frameNr, data.x,
                      NULL, data.voi[ri].y2, data.voi[ri].y3, data.frameNr);
    if(ret) {
      fprintf(stderr, "Error in integration of tissue data (%d).\n", ret);
      dftEmpty(&data); dftEmpty(&pdata); return(2);
    }
  }
 
  /* Find the time where plasma curve starts to rise by integrating original plasma data 
     and searching the time where integral is >= 1/1000 of final integral */
  ret=integrate(pdata.x, pdata.voi[0].y, pdata.frameNr, pdata.voi[0].y2);
  if(ret) {
    fprintf(stderr, "Error in integration of plasma data (%d).\n", ret);
    dftEmpty(&data); dftEmpty(&pdata); return(2);
  }
  for(fi=0, onep=0.0; fi<pdata.frameNr-1; fi++)
    if(pdata.voi[0].y2[fi]>=0.001*pdata.voi[0].y2[pdata.frameNr-1]) {onep=pdata.x[fi]; break;}
  if(verbose>1) fprintf(stdout, "Plasma curve starts to rise at about %g s.\n", onep);
 
  /*
   *  Interpolate and integrate plasma TAC at frame mid times
   *  Do this also for plasma TACs moved -MAX_DELAY - +MAX_DELAY sec
   */
  if(verbose>1) printf("integrating input TACs\n");
  dftInit(&ipdata);
  if(dftSetmem(&ipdata, data.frameNr, 2*MAX_DELAY+1)) {
    fprintf(stderr, "Error: out of memory.\n");
    dftEmpty(&data); dftEmpty(&pdata); return(3);
  }
  dftCopymainhdr(&pdata, &ipdata); ipdata.timetype=data.timetype;
  ipdata.voiNr=2*MAX_DELAY+1; ipdata.frameNr=data.frameNr;
  for(fi=0; fi<data.frameNr; fi++) {
    ipdata.x[fi]=data.x[fi]; ipdata.x1[fi]=data.x1[fi]; ipdata.x2[fi]=data.x2[fi];}
  /* Original plasma data */
  ri=0;
  if(ipdata.timetype==3) {
    ret=interpolate4pet(pdata.x, pdata.voi[0].y, pdata.frameNr,
        ipdata.x1, ipdata.x2,
        ipdata.voi[ri].y, ipdata.voi[ri].y2, ipdata.voi[ri].y3, ipdata.frameNr);
  } else {
    ret=interpolate(pdata.x, pdata.voi[0].y, pdata.frameNr,
        ipdata.x, ipdata.voi[ri].y, ipdata.voi[ri].y2, ipdata.voi[ri].y3,
        ipdata.frameNr);
  }
  if(ret) {
    fprintf(stderr, "Error (%d) in interpolation of plasma data.\n", ret);
    dftEmpty(&data); dftEmpty(&pdata); dftEmpty(&ipdata); return(4);
  }
  sprintf(ipdata.voi[ri].voiname, "d%d", 0);
  if(verbose>3) printf("%s ip-integrals: %g %g %g\n", ipdata.voi[ri].voiname,
    ipdata.voi[ri].y[ipdata.frameNr-1], ipdata.voi[ri].y2[ipdata.frameNr-1],
    ipdata.voi[ri].y3[ipdata.frameNr-1]);
  /* Subtracted times , skip negative times */
  for(ri=1; ri<=MAX_DELAY; ri++) {
    for(fi=0; fi<pdata.frameNr; fi++) {
      pdata.x[fi]-=1.0; pdata.x1[fi]-=1.0; pdata.x2[fi]-=1.0;}
    if(ipdata.timetype==3)
      ret=interpolate4pet(pdata.x, pdata.voi[0].y, pdata.frameNr,
          ipdata.x1, ipdata.x2, ipdata.voi[ri].y, ipdata.voi[ri].y2,
      ipdata.voi[ri].y3, ipdata.frameNr);
    else
      ret=interpolate(pdata.x, pdata.voi[0].y, pdata.frameNr,
        ipdata.x, ipdata.voi[ri].y, ipdata.voi[ri].y2, ipdata.voi[ri].y3, 
      ipdata.frameNr);
    if(ret) {
      fprintf(stderr, "Error (%d) in interpolation.\n", ret);
      dftEmpty(&data); dftEmpty(&pdata); dftEmpty(&ipdata); return(4);
    }
    sprintf(ipdata.voi[ri].voiname, "d%d", -ri);
    if(verbose>7) printf("%s ip-integrals: %g %g %g\n", ipdata.voi[ri].voiname,
      ipdata.voi[ri].y[ipdata.frameNr-1], ipdata.voi[ri].y2[ipdata.frameNr-1],
      ipdata.voi[ri].y3[ipdata.frameNr-1]);
  }
  for(fi=0; fi<pdata.frameNr; fi++) {
    pdata.x[fi]+=(double)MAX_DELAY;
    pdata.x1[fi]+=(double)MAX_DELAY; pdata.x2[fi]+=(double)MAX_DELAY;
  }
  /* Added times */
  for(ri=MAX_DELAY+1; ri<=2*MAX_DELAY; ri++) {
    for(fi=0; fi<pdata.frameNr; fi++) {
      pdata.x[fi]+=1.0; pdata.x1[fi]+=1.0; pdata.x2[fi]+=1.0;}
    if(ipdata.timetype==DFT_TIME_STARTEND) {
      ret=interpolate4pet(pdata.x, pdata.voi[0].y, pdata.frameNr,
               ipdata.x1, ipdata.x2, ipdata.voi[ri].y, ipdata.voi[ri].y2,
               ipdata.voi[ri].y3, ipdata.frameNr);
      if(ret==3) { /* time ranges do not match any more */
        for(fi=0; fi<ipdata.frameNr; fi++)
          ipdata.voi[ri].y[fi]=ipdata.voi[ri].y2[fi]=ipdata.voi[ri].y3[fi]=0.0;
        ret=0;
      }
    } else {
      ret=interpolate(pdata.x, pdata.voi[0].y, pdata.frameNr,
              ipdata.x, ipdata.voi[ri].y, ipdata.voi[ri].y2, ipdata.voi[ri].y3, ipdata.frameNr);
    }
    if(ret!=0) {
      fprintf(stderr, "Error in interpolation (%d).\n", ret);
      dftEmpty(&data); dftEmpty(&pdata); dftEmpty(&ipdata); return(4);
    }
    sprintf(ipdata.voi[ri].voiname, "d%d", ri-MAX_DELAY);
    if(verbose>7) printf("%s ip-integrals: %g %g %g\n", ipdata.voi[ri].voiname,
      ipdata.voi[ri].y[ipdata.frameNr-1], ipdata.voi[ri].y2[ipdata.frameNr-1],
      ipdata.voi[ri].y3[ipdata.frameNr-1]);
  }
  for(fi=0; fi<pdata.frameNr; fi++) {
    pdata.x[fi]-=(double)MAX_DELAY;
    pdata.x1[fi]-=(double)MAX_DELAY; pdata.x2[fi]-=(double)MAX_DELAY;
  }
  if(verbose>6) dftPrint(&ipdata);
 
 
  /*
   *  Fit K1-k4 & Vb
   */
  if(verbose>1) printf("fitting\n");
  /* Allocate memory required by NNLS */
  nnls_n=NNLS_N; nnls_m=data.frameNr;
  nnls_mat=(double*)malloc(((nnls_n+2)*nnls_m)*sizeof(double));
  if(nnls_mat==NULL) {
    fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
    dftEmpty(&data); dftEmpty(&pdata); dftEmpty(&ipdata);
    return(5);
  }
  for(n=0, dptr=nnls_mat; n<nnls_n; n++) {nnls_a[n]=dptr; dptr+=nnls_m;}
  nnls_b=dptr; dptr+=nnls_m; nnls_zz=dptr;
  /* Allocate memory for delays producing smallest SS */
  regional_delay=(double*)calloc(data.voiNr, sizeof(double));
  if(regional_delay==NULL) {
    fprintf(stderr, "Error: cannot allocate memory for regional delays.\n");
    dftEmpty(&data); dftEmpty(&pdata); dftEmpty(&ipdata); free(nnls_mat);
    return(5);
  }
  /* Allocate memory for fitted (best) curves */
  if(ffile[0]) {
    ret=dftdup(&data, &fdata);
    if(ret) {
      fprintf(stderr, "Error: cannot allocate memory for fitted TAC(s).\n");
      dftEmpty(&data); dftEmpty(&pdata); dftEmpty(&ipdata); 
      free(nnls_mat); free(regional_delay);
      return(5);
    }
  }
 
  /* Fit with different regions */
  if(data.isweight) { /* use quare root of weights */
    for(m=0; m<nnls_m; m++) {
      if(data.w[m]<=1.0e-20) data.w[m]=0.0; else data.w[m]=sqrt(data.w[m]);}
  }
  if(data.voiNr>1) printf("Regional time delays:\n");
  for(ri=0; ri<data.voiNr; ri++) {
 
 
    /* If requested, open file to save NNLS matrix in */
    FILE *mfp=NULL;
    if(matfile[0]) {
      mfp=fopen(matfile, "w");
      if(mfp==NULL) {
        fprintf(stderr, "Error: cannot write file '%s'.\n", matfile);
        dftEmpty(&data); dftEmpty(&pdata); dftEmpty(&ipdata); 
        free(nnls_mat); free(regional_delay);
        return(30);
      }
    }
 
    /* Fit with different delays */
    min=-1; minv=9.9e+99;
    for(di=0; di<=2*MAX_DELAY; di++) {
 
      /* Fill the NNLS data matrix */
      if(model==1) {
        for(m=0; m<nnls_m; m++) {
          nnls_a[0][m]=ipdata.voi[di].y[m];
          nnls_a[1][m]=ipdata.voi[di].y2[m];
          nnls_a[2][m]=-data.voi[ri].y2[m];
          nnls_b[m]=data.voi[ri].y[m];
        }
      } else {
        for(m=0; m<nnls_m; m++) {
          nnls_a[0][m]=ipdata.voi[di].y[m];
          nnls_a[1][m]=ipdata.voi[di].y2[m];
          nnls_a[2][m]=ipdata.voi[di].y3[m];
          nnls_a[3][m]=-data.voi[ri].y2[m];
          nnls_a[4][m]=-data.voi[ri].y3[m];
          nnls_b[m]=data.voi[ri].y[m];
        }
      }
 
      /* Add weights */
      if(data.isweight) for(m=0; m<nnls_m; m++) {
        nnls_b[m]*=data.w[m];
        for(n=0; n<nnls_n; n++) nnls_a[n][m]*=data.w[m];
      }
 
      /* If requested, save NNLS matrix in file */
      if(matfile[0]) {
        int dt;
        if(di<=MAX_DELAY) dt=-di; else dt=di-MAX_DELAY;
        fprintf(mfp, "%d\nB", dt);
        for(int n=0; n<nnls_n; n++) fprintf(mfp, ",A%d", n+1);
        fprintf(mfp, "\n");
        for(int m=0; m<nnls_m; m++) {
          fprintf(mfp, "%g", nnls_b[m]);
          for(int n=0; n<nnls_n; n++) fprintf(mfp, ",%g", nnls_a[n][m]);
          fprintf(mfp, "\n");
        }
      }
 
      /* NNLS */
      ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm, nnls_wp, nnls_zz, nnls_index);
      if(ret>1) { /* no solution is possible */
        continue;
      } else if(ret==1) { /* max iteration count exceeded */ }
      /* Get coefficients */
      for(n=0; n<nnls_n; n++) coeff[n]=nnls_x[n];
      /* Calculate Sum-of-Squares */
      if(model==1) {
        for(m=0; m<nnls_m; m++) {
          nnls_a[0][m]=ipdata.voi[di].y[m];
          nnls_a[1][m]=ipdata.voi[di].y2[m];
          nnls_a[2][m]=-data.voi[ri].y2[m];
          nnls_b[m]=data.voi[ri].y[m];
        }
      } else {
        for(m=0; m<nnls_m; m++) {
          nnls_a[0][m]=ipdata.voi[di].y[m];
          nnls_a[1][m]=ipdata.voi[di].y2[m];
          nnls_a[2][m]=ipdata.voi[di].y3[m];
          nnls_a[3][m]=-data.voi[ri].y2[m];
          nnls_a[4][m]=-data.voi[ri].y3[m];
          nnls_b[m]=data.voi[ri].y[m];
        }
      }
      if(ffile[0]) {
        for(m=0, ss=0.0; m<nnls_m; m++) {
          for(n=0, f=0.0; n<nnls_n; n++) f+=nnls_x[n]*nnls_a[n][m];
          fdata.voi[ri].y2[m]=f;
        }
      }
      if(data.isweight) for(m=0; m<nnls_m; m++) {
        nnls_b[m]*=data.w[m];
        for(n=0; n<nnls_n; n++) nnls_a[n][m]*=data.w[m];
      }
      for(m=0, ss=0.0; m<nnls_m; m++) {
        for(n=0, f=0.0; n<nnls_n; n++) f+=nnls_x[n]*nnls_a[n][m];
        v=f; f-=nnls_b[m]; ss+=f*f;
      }
 
      if(verbose>4) {
        if(model==1) {
          printf("Fit %2d -> SS=%12.3e Vb=%g K1+Vb*k2=%g k2=%g\n", di, ss, coeff[0], coeff[1], coeff[2]);
        } else {
          printf("Fit %2d -> SS=%12.3e Vb=%g K1+Vb*(k2+k3+k4)=%g\n", di, ss, coeff[0], coeff[1]);
          printf("           K1*(k3+k4)=%g k2+k3+k4=%g k2k4=%g\n", coeff[2], coeff[3], coeff[4]);
        }
      }
 
      /* Check if minimum SS */
      if(ss<minv) {
        minv=ss; min=di;
        if(ffile[0]) for(m=0; m<nnls_m; m++) fdata.voi[ri].y[m]=fdata.voi[ri].y2[m];
      }
 
    } /* next delay time */
 
    if(matfile[0]) fclose(mfp);
 
    /* what delay times the minima are? */
    if(verbose>5) printf("  min TAC (%s) := %d\n", data.voi[ri].name, min);
    if(min<=MAX_DELAY) min=-min; else min-=MAX_DELAY;
    regional_delay[ri]=(double)min;
    if(verbose>=0 && data.voiNr>1)
      printf("  time_delay for %s := %g s\n", data.voi[ri].name, regional_delay[ri]);
 
  } /* next region */
 
  /* Free memory of fit data */
  free(nnls_mat);
 
 
  /*
   *  Calculate the median of all delay values
   */
  if(verbose>1) printf("calculating median\n");
  v=dmedian(regional_delay, data.voiNr);
  /* regional delays are not needed anymore */
  free(regional_delay);
  if(verbose>=0) {
    printf("Estimated tracer appearance time in plasma := %.1f s\n", onep);
    printf("Estimated tracer appearance time in tissue := %.1f s\n", onep+(double)v);
    printf("Median time delay := %g s\n", v);
  }
 
  /* Check that plasma is not moved too much to the left */
  /* i.e. it should not start to rise before time 0      */
  /*
  if(-v>onep) v=-onep;
  printf("Applied delay %g s\n", v);
  */
  if(-v>onep) {
    fprintf(stderr, "Warning: possible error in determined tracer appearance time.\n");
    fprintf(stderr, "         Please check the delay corrected input against tissue data.\n");
  }
 
  /* Ok this is the final delay time (in seconds) */
  double delayT=v;
 
 
  /*
   *  Delay correction for input data
   */
  if(verbose>1) printf("correcting input TACs\n");
  for(fi=0; fi<pdata.frameNr; fi++) {
    pdata.x[fi]+=delayT; pdata.x1[fi]+=delayT; pdata.x2[fi]+=delayT;}
  /* Remove samples which have now negative times */
  fi=0; while(fi<pdata.frameNr) {
    if(pdata.x[fi]>=0.0) break;
    for(fj=fi+1; fj<pdata.frameNr; fj++) {
      for(ri=0; ri<pdata.voiNr; ri++) pdata.voi[ri].y[fj-1]=pdata.voi[ri].y[fj];
      pdata.x[fj-1]=pdata.x[fj]; pdata.x1[fj-1]=pdata.x1[fj];
      pdata.x2[fj-1]=pdata.x2[fj];
    }
    pdata.frameNr--;
  }
  if(pdata.frameNr<=0) {
    fprintf(stderr, "Error: no positive sample times left.\n");
    dftEmpty(&data); dftEmpty(&pdata); dftEmpty(&ipdata); dftEmpty(&fdata);
    return(6);
  }
  if(verbose>10) dftPrint(&pdata);
 
 
  /*
   *  Save delay corrected plasma data
   */
  /* If not specified, then create the filename for output */
  if(!rfile[0]) {
    cptr=strrchr(pfile, '.'); if(cptr!=NULL) n=strlen(cptr); else n=0;
    strcpy(rfile, pfile); rfile[strlen(pfile)-n]=(char)0;
    strcat(rfile, ".delay"); if(cptr!=NULL) strcat(rfile, cptr);
    if(verbose>2) printf("rfile := %s\n", rfile);
    if(verbose>=0) printf("  %s -> %s\n", pfile, rfile);
  }
  /* Set output file format */
  if(output_format!=DFT_FORMAT_UNKNOWN) pdata._type=output_format;
  /* Change times back to original unit if necessary */
  if(pdata.timeunit!=orig_input_time_unit) {
    if(verbose>2) printf("converting input %s to %s\n",
      petTunit(pdata.timeunit), petTunit(orig_input_time_unit) );
    if(verbose>3) printf("last_t := %g\n", pdata.x[pdata.frameNr-1]);
    ret=dftTimeunitConversion(&pdata, orig_input_time_unit);
    if(verbose>3) printf("last_t := %g\n", pdata.x[pdata.frameNr-1]);
  }
  /* Write the log information, if required */
  dftSetComments(&pdata);
  if(make_log!=0) {
    //if(strlen(pdata.comments)>0) strcat(pdata.comments, "\n");
    sprintf(tmp, "# delay_fit_time := %g s\n", length);
    strcat(pdata.comments, tmp);
    sprintf(tmp, "# time_delay := %g s\n", delayT);
    strcat(pdata.comments, tmp);
    sprintf(tmp, "# tissue_compartment_nr := %d\n", model);
    strcat(pdata.comments, tmp);
  }
  /* Write the TAC */
  if(verbose>1) printf("writing %s\n", rfile);
  if(dftWrite(&pdata, rfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", rfile, dfterrmsg);
    dftEmpty(&data); dftEmpty(&pdata); dftEmpty(&ipdata); dftEmpty(&fdata);
    return(11);
  }
 
  /*
   *  Free memory
   */
  dftEmpty(&pdata); dftEmpty(&ipdata);
 
  /*
   *  Save fitted TACs, if necessary
   */
  if(ffile[0]) {
    sprintf(fdata.comments, "# Curves fitted from %s\n", tfile);
 
    /* Change times back to original unit if necessary */
    if(data.timeunit!=orig_tissue_time_unit) {
      ret=dftTimeunitConversion(&data, orig_tissue_time_unit);
    }
    dftSetComments(&fdata);
    if(make_log!=0) {
      //if(strlen(pdata.comments)>0) strcat(pdata.comments, "\n");
      sprintf(tmp, "# delay_fit_time := %g s\n", length);
      strcat(pdata.comments, tmp);
      sprintf(tmp, "# time_delay := %g s\n", delayT);
      strcat(pdata.comments, tmp);
      sprintf(tmp, "# tissue_compartment_nr := %d\n", model);
      strcat(pdata.comments, tmp);
    }
    if(verbose>1) printf("writing %s\n", ffile);
    if(dftWrite(&fdata, ffile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", ffile, dfterrmsg);
      dftEmpty(&data); dftEmpty(&fdata);
      return(13);
    }
  }
  dftEmpty(&fdata);
 
 
  /*
   *  Delay correction for other files, if required
   */
  if(p2file[0]) {
    /* Read in */
    if(verbose>1) printf("reading %s\n", p2file);
    if(dftRead(p2file, &pdata)) {
      fprintf(stderr, "Error in reading '%s': %s\n", p2file, dfterrmsg);
      dftEmpty(&data); return(8);
    }
    /* If time unit is unknown, then set it */
    if(pdata.timeunit==TUNIT_UNKNOWN) {
      /* Assume that it is the same as in input data */
      pdata.timeunit=orig_input_time_unit;
    }
    /* Make delay correction, changing delay time to correct units */
    f=delayT; if(pdata.timeunit==TUNIT_MIN) f/=60;
    for(fi=0; fi<pdata.frameNr; fi++) {
      pdata.x[fi]+=f; pdata.x1[fi]+=f; pdata.x2[fi]+=f;}
    /* Remove negative times */
    fi=0; while(fi<pdata.frameNr) {
      if(pdata.x[fi]>=0.0) break;
      for(fj=fi+1; fj<pdata.frameNr; fj++) {
        for(ri=0; ri<pdata.voiNr; ri++) 
          pdata.voi[ri].y[fj-1]=pdata.voi[ri].y[fj];
        pdata.x[fj-1]=pdata.x[fj]; pdata.x1[fj-1]=pdata.x1[fj];
        pdata.x2[fj-1]=pdata.x2[fj];
      }
      pdata.frameNr--;
    }
    if(pdata.frameNr<=0) {
      fprintf(stderr, "Error: no positive sample times left.\n");
      dftEmpty(&data); dftEmpty(&pdata); return(9);
    }
    /* Construct output filename */
    if(!r2file[0]) {
      strcpy(r2file, p2file); cptr=strrchr(p2file, '.');
      if(cptr!=NULL) r2file[strlen(p2file)-strlen(cptr)]=(char)0;
      strcat(r2file, ".delay"); if(cptr!=NULL) strcat(r2file, cptr);
    }
    /* Set output file format */
    if(output_format!=DFT_FORMAT_UNKNOWN) pdata._type=output_format;
    /* Save output file */
    dftSetComments(&pdata);
    if(make_log!=0) {
      //if(strlen(pdata.comments)>0) strcat(pdata.comments, "\n");
      sprintf(tmp, "# delay_fit_time := %g s\n", length);
      strcat(pdata.comments, tmp);
      sprintf(tmp, "# time_delay := %g s\n", delayT);
      strcat(pdata.comments, tmp);
      sprintf(tmp, "# tissue_compartment_nr := %d\n", model);
      strcat(pdata.comments, tmp);
    }
    if(verbose>=0) printf("  %s -> %s\n", p2file, r2file);
    ret=dftWrite(&pdata, r2file); dftEmpty(&pdata);
    if(ret) {
      fprintf(stderr, 
        "Error (%d) in writing '%s': %s\n", ret, r2file, dfterrmsg);
      dftEmpty(&data); dftEmpty(&pdata); return(16);
    }
  }
 
  if(p3file[0]) {
    /* Read in */
    if(verbose>1) printf("reading %s\n", p3file);
    dftEmpty(&pdata);
    if(dftRead(p3file, &pdata)) {
      fprintf(stderr, "Error in reading '%s': %s\n", p3file, dfterrmsg);
      dftEmpty(&data); return(8);
    }
    /* If time unit is unknown, then set it */
    if(pdata.timeunit==TUNIT_UNKNOWN) {
      /* Assume that it is the same as in input data */
      pdata.timeunit=orig_input_time_unit;
    }
    /* Make delay correction, changing delay time to correct units */
    f=delayT; if(pdata.timeunit==TUNIT_MIN) f/=60;
    for(fi=0; fi<pdata.frameNr; fi++) {
      pdata.x[fi]+=f; pdata.x1[fi]+=f; pdata.x2[fi]+=f;}
    /* Remove negative times */
    fi=0; while(fi<pdata.frameNr) {
      if(pdata.x[fi]>=0.0) break;
      for(fj=fi+1; fj<pdata.frameNr; fj++) {
        for(ri=0; ri<pdata.voiNr; ri++) 
          pdata.voi[ri].y[fj-1]=pdata.voi[ri].y[fj];
        pdata.x[fj-1]=pdata.x[fj]; pdata.x1[fj-1]=pdata.x1[fj];
        pdata.x2[fj-1]=pdata.x2[fj];
      }
      pdata.frameNr--;
    }
    if(pdata.frameNr<=0) {
      fprintf(stderr, "Error: no positive sample times left.\n");
      dftEmpty(&data); dftEmpty(&pdata); return(9);
    }
    /* Construct output filename */
    if(!r3file[0]) {
      strcpy(r3file, p3file); cptr=strrchr(p3file, '.');
      if(cptr!=NULL) r3file[strlen(p3file)-strlen(cptr)]=(char)0;
      strcat(r3file, ".delay"); if(cptr!=NULL) strcat(r3file, cptr);
    }
    /* Set output file format */
    if(output_format!=DFT_FORMAT_UNKNOWN) pdata._type=output_format;
    /* Save output file */
    dftSetComments(&pdata);
    if(make_log!=0) {
      //if(strlen(pdata.comments)>0) strcat(pdata.comments, "\n");
      sprintf(tmp, "# delay_fit_time := %g s\n", length);
      strcat(pdata.comments, tmp);
      sprintf(tmp, "# time_delay := %g s\n", delayT);
      strcat(pdata.comments, tmp);
      sprintf(tmp, "# tissue_compartment_nr := %d\n", model);
      strcat(pdata.comments, tmp);
    }
    if(verbose>=0) printf("  %s -> %s\n", p3file, r3file);
    ret=dftWrite(&pdata, r3file); dftEmpty(&pdata);
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s': %s\n", ret, r3file, dfterrmsg);
      dftEmpty(&data); dftEmpty(&pdata); return(17);
    }
  }
 
  if(p4file[0]) {
    /* Read in */
    if(verbose>1) printf("reading %s\n", p4file);
    dftEmpty(&pdata);
    if(dftRead(p4file, &pdata)) {
      fprintf(stderr, "Error in reading '%s': %s\n", p4file, dfterrmsg);
      dftEmpty(&data); return(8);
    }
    /* If time unit is unknown, then set it */
    if(pdata.timeunit==TUNIT_UNKNOWN) {
      /* Assume that it is the same as in input data */
      pdata.timeunit=orig_input_time_unit;
    }
    /* Make delay correction, changing delay time to correct units */
    f=delayT; if(pdata.timeunit==TUNIT_MIN) f/=60;
    for(fi=0; fi<pdata.frameNr; fi++) {
      pdata.x[fi]+=f; pdata.x1[fi]+=f; pdata.x2[fi]+=f;}
    /* Remove negative times */
    fi=0; while(fi<pdata.frameNr) {
      if(pdata.x[fi]>=0.0) break;
      for(fj=fi+1; fj<pdata.frameNr; fj++) {
        for(ri=0; ri<pdata.voiNr; ri++) pdata.voi[ri].y[fj-1]=pdata.voi[ri].y[fj];
        pdata.x[fj-1]=pdata.x[fj]; pdata.x1[fj-1]=pdata.x1[fj];
        pdata.x2[fj-1]=pdata.x2[fj];
      }
      pdata.frameNr--;
    }
    if(pdata.frameNr<=0) {
      fprintf(stderr, "Error: no positive sample times left.\n");
      dftEmpty(&data); dftEmpty(&pdata); return(9);
    }
    /* Construct output filename */
    if(!r4file[0]) {
      strcpy(r4file, p4file); cptr=strrchr(p4file, '.');
      if(cptr!=NULL) r4file[strlen(p3file)-strlen(cptr)]=(char)0;
      strcat(r4file, ".delay"); if(cptr!=NULL) strcat(r4file, cptr);
    }
    /* Set output file format */
    if(output_format!=DFT_FORMAT_UNKNOWN) pdata._type=output_format;
    /* Save output file */
    dftSetComments(&pdata);
    if(make_log!=0) {
      //if(strlen(pdata.comments)>0) strcat(pdata.comments, "\n");
      sprintf(tmp, "# delay_fit_time := %g s\n", length);
      strcat(pdata.comments, tmp);
      sprintf(tmp, "# time_delay := %g s\n", delayT);
      strcat(pdata.comments, tmp);
      sprintf(tmp, "# tissue_compartment_nr := %d\n", model);
      strcat(pdata.comments, tmp);
    }
    if(verbose>=0) printf("  %s -> %s\n", p4file, r4file);
    ret=dftWrite(&pdata, r4file); dftEmpty(&pdata);
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s': %s\n", ret, r4file, dfterrmsg);
      dftEmpty(&data); dftEmpty(&pdata); return(18);
    }
  }
 
  /*
   *  Free memory
   */
  dftEmpty(&data); dftEmpty(&pdata);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/