b2t_mo2.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulation of PET tissue time-radioactivity concentration curve (TTAC)",
  "in skeletal muscle [O-15]O2 PET studies from decay corrected arterial blood",
  "[O-15]O2 and [O-15]H2O curves (oxygen and water BTACs).",
  "The compartmental model for [O-15]O2 in skeletal muscle is:",
  " ",
  "      Blood         Muscle                 ",
  "  _______________________________________  ",
  " |           K1 |                        | ",
  " |    O    ---------->  O   +  O -Mb     | ",
  " |     2   <----------   2      2        | ",
  " |             k2                        | ",
  " |              |        | k3            | ",
  " |              |        V               | ",
  " |           K1 |                        | ",
  " |    H O  ---------->  H O              | ",
  " |     2   <----------   2               | ",
  " |              | K1/pH2O                | ",
  " |______________|________________________| ",
  " ",
  "Model definitions:",
  "  K1=perfusion (f, mL/min/dL), OER=k3/(k2+k3),",
  "  Ki=f*OER (mL/min/dL), and",
  "  Metabolic rate of oxygen MRO2=Ki*[O2]a (mmol/min/dL)",
  " ",
  "Usage: @P [options] obtacfile wbtacfile f OER ttacfile",
  " ",
  "Options:",
  " -sub | -nosub",
  "     TACs of sub-compartments are written (-sub)",
  "     or not written (-nosub, default) into the output file.",
  " -add",
  "     Simulated TACs are added to an existing tissue data file.",
  "     By default, existing file is overwritten.",
  " -Vb=<Blood volume (%)>",
  "     Set the simulated blood volume; default is 3.5%.",
  " -Af=<Arterial proportion (%)>",
  "     Set the simulated arterial proportion of total blood volume;",
  "     default is 30%.",
  " -pH2O=<value>",
  "     Set the partition coefficient of water; default is 0.99.",
  " -K1k2=<value>",
  "     Set the K1/k2 for oxygen; by default the K1/k2 is estimated from",
  "     OER and saturation curves for hemoglobin and myoglobin.",
  " -Sao2=<value>",
  "     Saturation of arterial blood hemoglobin; default is 0.97.",
  " -p50Hb=<value>",
  "     Half-saturation pressure for hemoglobin; default is 3.6 kPa.",
  " -p50Mb=<value>",
  "     Half-saturation pressure for myoglobin; default is 0.319 kPa.",
  " -nHb=<value>",
  "     Hill coefficient for hemoglobin; default is 2.7.",
  " -Mb=<value>",
  "     Concentration of myoglobin in muscle; default is 4.7 mg/g.",
  " -Hb=<value>",
  "     Concentration of hemoglobin in blood; default is 150 mg/g.",
  " -fpt",
  "     Blood flow (perfusion) is assumed to be given per perfusable tissue",
  "     volume excluding vascular volume. TTAC will still be simulated per",
  "     regional PET volume including vascular volume.",
  " -voiname=<text>",
  "     Enter a name (1-6 chars without spaces) for the simulated TTAC.", 
  " -venao=<filename>",
  "     Save the simulated venous [O-15]O2 BTAC",
  " -venaw=<filename>",
  "     Save the simulated venous [O-15]H2O BTAC",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "For accurate results, input BTACs should be noiseless and have very",
  "short sampling intervals. Simulated curves can thereafter be interpolated",
  "to represent PET frames using program simframe.",
  "Calculated tissue activities are written in the specified file with",
  "these data columns. Columns 2-6 will be saved optionally (-sub):",
  "  0) Time",
  "  1) Total regional radioactivity concentration (2+3+4+5+6)",
  "  2) Labelled [O2] component of TTAC",
  "  3) Labelled [H2O] component of TTAC",
  "  4) Arterial BTAC component of TTAC",
  "  5) Venous labelled [O2] component of TTAC",
  "  6) Venous labelled [H2O] component of TTAC",
  " ",
  "References:",
  "1. Nuutila P, Peltoniemi P, Oikonen V, Larmola K, Kemppainen J, Takala T,",
  "   Sipila H, Oksanen A, Ruotsalainen U, Bolli GB, Yki-Jarvinen H.",
  "   Enhanced stimulation of glucose uptake by insulin increases",
  "   exercise-stimulated glucose uptake in skeletal muscle in humans: studies",
  "   using [15O]O2, [15O]H2O, [18F]fluoro-deoxy-glucose, and positron emission",
  "   tomography. Diabetes 2000; 49:1084-1091.",
  "2. Oikonen V, Nuutila P, Sipila H, Tolvanen T, Peltoniemi P, Ruotsalainen U.",
  "   Quantification of oxygen consumption in skeletal muscle with PET and",
  "   oxygen-15 bolus. Eur J Nucl Med. 1998; 25: 1151.",
  "3. Oikonen V. Modelling of low oxygen consumption. In: J. Knuuti, J. Rinne,",
  "   P.Tenhonen (ed.), Medical Applications of Cyclotrons VIII. Abstracts of",
  "   the VIII Symposium on the Medical Applications of Cyclotrons.",
  "   Annales Universitatis Turkuensis D346:16, 1999.",
  " ",
  "See also: fit_mo2, o2metab, o2_p2w, fit_o2bl, sim_o2bl, tacadd, simframe",
  " ",
  "Keywords: TAC, simulation, modelling, oxygen, skeletal muscle",
  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     ofile[FILENAME_MAX], wfile[FILENAME_MAX], dfile[FILENAME_MAX];
  char     vofile[FILENAME_MAX], vwfile[FILENAME_MAX];
  char    *cptr, voiname[MAX_REGIONNAME_LEN+1];
  int      save_only_total=1;
  int      flow_per_tissue=0;
  int      add_to_previous=0;
  int      K1k2_def_changed=0;
  double   fK1k2; 
  double   Af=0.30;
  double   Vb=0.035;
  double   pH2O=0.99;
  double   Flow;
  double   OER;
  double   SaO2=DEFAULT_SAO2;
  double   p50Hb=DEFAULT_P50HB;
  double   p50Mb=DEFAULT_P50MB;
  double   nHb=DEFAULT_NHB;
  double   cHb=DEFAULT_CHB;
  double   cMb=DEFAULT_CMB;
 
  DFT      data, input;
  int      n, ret, times_changed=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ofile[0]=wfile[0]=dfile[0]=vofile[0]=vwfile[0]=(char)0;
  strcpy(voiname, "Muscle");
  fK1k2=Flow=OER=nan("");
  dftInit(&input); dftInit(&data);
  /* 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;
    if(strncasecmp(cptr, "NOSUB", 5)==0) {
      save_only_total=1; continue;
    } else if(strncasecmp(cptr, "SUB", 3)==0) {
      save_only_total=0; continue;
    } else if(strcasecmp(cptr, "FPT")==0) {
      flow_per_tissue=1; continue;
    } else if(strcasecmp(cptr, "ADD")==0) {
      add_to_previous=1; continue;
    } else if(strncasecmp(cptr, "VOINAME=", 8)==0) {
      n=strlcpy(voiname, cptr+8, MAX_REGIONNAME_LEN);
      if(n<1 || n>=MAX_REGIONNAME_LEN) {
        fprintf(stderr, "Error: invalid VOI name '%s'.\n", cptr+8); return(1);}
      continue;
    } else if(strncasecmp(cptr, "VENAO=", 6)==0) {
      n=strlcpy(vofile, cptr+6, FILENAME_MAX);
      if(n>0 && n<FILENAME_MAX) continue;
    } else if(strncasecmp(cptr, "VENAW=", 6)==0) {
      n=strlcpy(vwfile, cptr+6, FILENAME_MAX);
      if(n>0 && n<FILENAME_MAX) continue;
    } else if(strncasecmp(cptr, "Vb=", 3)==0) {
      if(!atof_with_check(cptr+3, &Vb)) {Vb*=0.01; continue;}
    } else if(strncasecmp(cptr, "Af=", 3)==0) {
      if(!atof_with_check(cptr+3, &Af)) {Af*=0.01; continue;}
    } else if(strncasecmp(cptr, "pH2O=", 5)==0) {
      if(!atof_with_check(cptr+5, &pH2O)) continue;
    } else if(strncasecmp(cptr, "K1k2=", 5)==0) {
      if(!atof_with_check(cptr+5, &fK1k2)) continue;
    } else if(strncasecmp(cptr, "sao2=", 5)==0) {
      if(!atof_with_check(cptr+5, &SaO2)) {K1k2_def_changed=1; continue;}
    } else if(strncasecmp(cptr, "p50hb=", 6)==0) {
      if(!atof_with_check(cptr+5, &p50Hb)) {K1k2_def_changed=1; continue;}
    } else if(strncasecmp(cptr, "p50mb=", 6)==0) {
      if(!atof_with_check(cptr+5, &p50Mb)) {K1k2_def_changed=1; continue;}
    } else if(strncasecmp(cptr, "nhb=", 4)==0) {
      if(!atof_with_check(cptr+5, &nHb)) {K1k2_def_changed=1; continue;}
    } else if(strncasecmp(cptr, "hb=", 3)==0) {
      if(!atof_with_check(cptr+5, &cHb)) {K1k2_def_changed=1; continue;}
    } else if(strncasecmp(cptr, "mb=", 3)==0) {
      if(!atof_with_check(cptr+5, &cMb)) {K1k2_def_changed=1; 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(!ofile[0]) {strlcpy(ofile, argv[ai], FILENAME_MAX); continue;}
    else if(!wfile[0]) {strlcpy(wfile, argv[ai], FILENAME_MAX); continue;}
    else if(isnan(Flow)) {if(!atof_with_check(argv[ai], &Flow)) continue;}
    else if(isnan(OER)) {if(!atof_with_check(argv[ai], &OER)) continue;}
    else if(!dfile[0]) {strlcpy(dfile, argv[ai], FILENAME_MAX); continue;}
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  Flow/=6000.;
 
  /* Did we get all the information that we need? */
  if(!dfile[0]) {
    fprintf(stderr, "Error: missing command-line argument.\n\n");
    tpcPrintUsage(argv[0], info, stdout); return(1);
  }
 
  /* Check parameters */
  if(Af<0.0 || Af>1.0) {
    fprintf(stderr, "Error: invalid arterial fraction.\n"); return(1);
  }
  if(Af>0.0 && Af<0.01) { 
    // check if user accidentally gave fraction instead of percentage
    fprintf(stderr, "Warning: arterial fraction was set to %g%%.\n", 100.0*Af);
  }
  if(Vb<0.0 || Vb>=1.0) {
    fprintf(stderr, "Error: invalid blood volume.\n"); return(1);
  }
  if(Vb>0.0 && Vb<0.01) { 
    // check if user accidentally gave fraction instead of percentage
    fprintf(stderr, "Warning: arterial fraction was set to %g%%.\n", 100.0*Vb);
  }
  if(pH2O>1.0 || pH2O<=0.) {
    fprintf(stderr, "Error: invalid partition constant for water.\n");
    return(1);
  }
  if(OER>10.0 && OER<95.0) OER/=100.;
  if(OER<=0.0 || OER>=1.0) {
    fprintf(stderr, "Error: invalid OER.\n"); return(1);
  }
  if(K1k2_def_changed!=0 && fK1k2>0.0)
    fprintf(stderr, "Warning: option(s) superseded by -K1k2.\n");
  if(SaO2<=0.0 || SaO2>1.0) {
    fprintf(stderr, "Error: invalid SaO2.\n"); return(1);
  } 
  if(p50Hb<=0.0) {fprintf(stderr, "Error: invalid p50Hb.\n"); return(1);} 
  if(p50Mb<=0.0) {fprintf(stderr, "Error: invalid p50Mb.\n"); return(1);} 
  if(nHb<=0.0) {fprintf(stderr, "Error: invalid nHb.\n"); return(1);} 
  if(cMb<=0.0) {fprintf(stderr, "Error: invalid [Mb].\n"); return(1);} 
  if(cHb<=0.0) {fprintf(stderr, "Error: invalid [Hb].\n"); return(1);}
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("save_only_total := %d\n", save_only_total);
    printf("flow_per_tissue := %d\n", flow_per_tissue);
    printf("add_to_previous := %d\n", add_to_previous);
    printf("ofile := %s\n", ofile);
    printf("wfile := %s\n", wfile);
    printf("dfile := %s\n", dfile);
    if(vofile[0]) printf("vofile := %s\n", vofile);
    if(vwfile[0]) printf("vwfile := %s\n", vwfile);
    printf("Vb := %g%%\n", 100.*Vb);
    printf("Af := %g%%\n", 100.*Af);
    printf("pH2O := %g\n", pH2O);
    printf("Flow := %g mL/dL/min\n", 6000.*Flow);
    printf("OER := %g\n", OER);
    printf("SaO2 := %g\n", SaO2);
    printf("p50Hb := %g\n", p50Hb);
    printf("p50Mb := %g\n", p50Mb);
    printf("nHb := %g\n", nHb);
    printf("cMb := %g\n", cMb);
    printf("cHb := %g\n", cHb);
    if(!isnan(fK1k2)) printf("fK1k2 := %g\n", fK1k2);
    if(voiname[0]) printf("voiname := %s\n", voiname);
  }
 
 
  /*
   *  Read input data
   */
  if(verbose>1) printf("reading '%s'\n", ofile);
  if(dftRead(ofile, &input)) {
    fprintf(stderr, "Error in reading '%s': %s\n", ofile, dfterrmsg);
    return(2);
  }
  if(input.frameNr<3) {
    fprintf(stderr, "Error: not enough input samples for decent simulation.\n");
    dftEmpty(&input); return(2);
  }
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&input);
  /* Check for NA's */
  if(dft_nr_of_NA(&input) > 0) {
    fprintf(stderr, "Error: missing sample(s) in %s\n", ofile);
    dftEmpty(&input); return(2);
  }
  if(input.voiNr>1) {
    fprintf(stderr, "Warning: [O-15]O2 BTAC file may not be valid.\n");
    input.voiNr=1;
  }
  /* If sample times are minutes, change them to sec */
  if(input.timeunit!=TUNIT_SEC && input.timeunit!=TUNIT_MIN) {
    if(input.x[input.frameNr-1]<20.) {
      if(verbose>0) fprintf(stdout, "Note: sample times assumed to be in min.\n");
      input.timeunit=TUNIT_MIN;
    } else {
      if(verbose>0) fprintf(stdout, "Note: sample times assumed to be in sec.\n");
      input.timeunit=TUNIT_SEC;
    }
  }
  if(input.timeunit==TUNIT_MIN) {
    if(verbose>0) 
      fprintf(stderr, "Note: Oxygen BTAC sample times are converted to sec.\n");
    dftMin2sec(&input);
    times_changed=1;
  }
  /*
   *  Read H2O BTAC
   */
  if(verbose>1) printf("reading '%s'\n", wfile);
  if(dftRead(wfile, &data)) {
    fprintf(stderr, "Error in reading '%s': %s\n", wfile, dfterrmsg);
    dftEmpty(&input); dftEmpty(&data); return(3);
  }
  if(data.frameNr<input.frameNr) {
    fprintf(stderr, "Error: not enough samples for decent simulation.\n");
    dftEmpty(&input); dftEmpty(&data); return(3);
  }
  /* Check for NA's */
  if(dft_nr_of_NA(&data) > 0) {
    fprintf(stderr, "Error: missing sample(s) in %s\n", wfile);
    dftEmpty(&input); dftEmpty(&data); return(3);
  }
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&data);
  if(data.voiNr>1) {
    fprintf(stderr, "Warning: [O-15]H2O BTAC file may not be valid.\n");
    data.voiNr=1;
  }
  /* If sample times are minutes, change them to sec */
  if(data.timeunit!=TUNIT_SEC && data.timeunit!=TUNIT_MIN) {
    if(data.x[data.frameNr-1]<20.) {
      if(verbose>0) fprintf(stdout, "Note: sample times assumed to be in min.\n");
      data.timeunit=TUNIT_MIN;
    } else {
      if(verbose>0) fprintf(stdout, "Note: sample times assumed to be in sec.\n");
      data.timeunit=TUNIT_SEC;
    }
  }
  if(data.timeunit==TUNIT_MIN) {
    if(verbose>0) 
      fprintf(stderr, "Note: Radiowater BTAC sample times are converted to sec.\n");
    dftMin2sec(&data);
  }
  /* Check that water TAC is not much shorter than oxygen TAC */
  if(data.x[data.frameNr-1]<0.75*input.x[input.frameNr-1]) {
    fprintf(stderr, "Error: [O-15]H2O BTAC is shorter than [O-15]O2 BTAC.\n");
    dftEmpty(&input); dftEmpty(&data); return(4);
  }
  if(data.x[data.frameNr-1]<0.95*input.x[input.frameNr-1])
    fprintf(stderr, "Warning: [O-15]H2O BTAC is shorter than [O-15]O2 BTAC.\n");
  /* Copy and interpolate water TAC to input TAC structure */
  if(dftAddmem(&input, 1)) {
    fprintf(stderr, "Error: cannot allocate more memory.\n");
    dftEmpty(&input); dftEmpty(&data); return(4);
  }
  dftCopyvoihdr(&data, 0, &input, 1);
  if(interpolate(data.x, data.voi[0].y, data.frameNr,
                 input.x, input.voi[1].y, NULL, NULL, input.frameNr)) 
  {
    fprintf(stderr, "Error: cannot interpolate [O-15]H2O BTAC.\n");
    dftEmpty(&input); dftEmpty(&data); return(4);
  }
  input.voiNr=2;
  /* Remove the originally read water data; pointer is re-used later */
  dftEmpty(&data);
  if(verbose>6) dftPrint(&input);
 
 
  /* Integrate BTACs; integral will be more correct in case of frames */
  ret=0;
  for(int ri=0; ri<input.voiNr && !ret; ri++) {
    if(input.timetype==DFT_TIME_STARTEND)
      ret=petintegral(input.x1, input.x2, input.voi[ri].y,
                  input.frameNr, input.voi[ri].y2, NULL);
    else
      ret=integrate(input.x, input.voi[ri].y, input.frameNr, input.voi[ri].y2);
  }
  if(ret) {
    fprintf(stderr, "Error: Cannot integrate input BTACs.\n");
    dftEmpty(&input); return(5);
  }
 
 
  /*
   *  Allocate memory for tissue data
   */
  if(verbose>1) printf("allocating memory for simulated data\n");
  if(dftSetmem(&data, input.frameNr, 6)) {
    fprintf(stderr, "Error: cannot allocate memory for simulated TTACs.\n");
    dftEmpty(&input); return(6);
  }
  data.frameNr=input.frameNr; data.voiNr=6;
  dftCopymainhdr(&input, &data);
  data._type=DFT_FORMAT_PMOD; // fileformat
  strcpy(data.isotope, "O-15");
  strcpy(data.radiopharmaceutical, "[O-15]O2");
  strcpy(data.injectionTime, input.injectionTime);
  strcpy(data.scanStartTime, input.scanStartTime);
  data.decayCorrected=DFT_DECAY_CORRECTED;
  /* Copy times */
  for(int fi=0; fi<input.frameNr; fi++) {
    data.x[fi]=input.x[fi];
    data.x1[fi]=input.x1[fi]; data.x2[fi]=input.x2[fi];
  }
  /* Set TTAC names */
  strcpy(data.voi[0].hemisphere, "Total");
  strcpy(data.voi[1].hemisphere, "O2");
  strcpy(data.voi[2].hemisphere, "H2O");
  strcpy(data.voi[3].hemisphere, "artery");
  strcpy(data.voi[4].hemisphere, "vena"); strcpy(data.voi[4].place, "O2");
  strcpy(data.voi[5].hemisphere, "vena"); strcpy(data.voi[5].place, "H2O");
  n=rnameSplit(voiname, data.voi[0].voiname, data.voi[0].hemisphere, 
               data.voi[0].place, MAX_REGIONSUBNAME_LEN);
  for(int ri=1; ri<data.voiNr; ri++) {
    strcpy(data.voi[ri].voiname, data.voi[0].voiname);
    if(n>1) strcpy(data.voi[ri].hemisphere, data.voi[0].hemisphere);
    if(n>2) strcpy(data.voi[ri].place, data.voi[0].place);
  }
  for(int ri=0; ri<data.voiNr; ri++)
    rnameCatenate(data.voi[ri].name, MAX_REGIONNAME_LEN, data.voi[ri].voiname,
                    data.voi[ri].hemisphere, data.voi[ri].place, '_');
 
 
  /*
   *  Simulate the TTACs
   */
  if(verbose>1) printf("simulating TTACs\n");
  double k1=Flow; if(flow_per_tissue==0) k1/=(1.0-Vb);
  double k2w=k1/pH2O;
  double k1k2; 
  if(!isnan(fK1k2)) k1k2=fK1k2; 
  else k1k2=mo2k1k2(OER, SaO2, p50Hb, p50Mb, nHb, cHb, cMb, verbose-1); 
  double k2o=k1/k1k2;
  double k3=k2o*OER/(1.0-OER);
  if(verbose>0) {
    double v;
    if((k2o+k3)>0.0) {
      printf("Vt_o := %g\n", k1/(k2o+k3));
      v=k1*k3/(k2o+k3);
      printf("K1*k3/(k2+k3) := %g\n", v);
    } else {
      v=0.0;
    }
    if(k2w>0.0) printf("Vt_w := %g\n", (v+k1)/k2w);
  }
  ret=simOxygen(
    input.x, input.voi[0].y, input.voi[1].y, input.voi[0].y2, input.voi[1].y2, 
    input.frameNr, k1, k2o, k3, k1, k2w, Vb, Af, 
    data.voi[0].y, data.voi[1].y, data.voi[2].y, data.voi[3].y,
    data.voi[4].y, data.voi[5].y, data.voi[4].y2, data.voi[5].y2,
    verbose-1 
  );
  dftEmpty(&input);
  if(ret!=0) {
    fprintf(stderr, "Error in simulation (%d).\n", ret);
    dftEmpty(&data); return(9);
  }
 
 
 
 
  /*
   *  Save simulated tissue TAC
   */
  if(verbose>1) printf("saving simulated TTAC.\n");
  if(times_changed) dftSec2min(&data);
  DFT_NR_OF_DECIMALS=6;
  if(save_only_total!=0) data.voiNr=1;
  dftSetComments(&data);
  if(add_to_previous!=0 && access(dfile, 0)!=-1) {
    if(verbose>0) printf("adding to existing file\n");
    DFT prevdft;
    /* read previous file */
    dftInit(&prevdft);
    if(dftRead(dfile, &prevdft)) {
      fprintf(stderr, "Error in reading '%s': %s\n", dfile, dfterrmsg);
      dftEmpty(&prevdft); dftEmpty(&data); return(2);
    }
    if(verbose>2) 
      printf("existing %d TACs with %d frames\n",prevdft.voiNr,prevdft.frameNr);
    /* add present simulations */
    for(int ri=0; ri<data.voiNr; ri++) {
      ret=dftAdd(&prevdft, &data, ri);
      if(ret) {
        fprintf(stderr, "Error: cannot add simulated TAC in %s\n", dfile);
        dftEmpty(&prevdft); dftEmpty(&data); return(13);
      }
    }
    /* save that */
    ret=dftWrite(&prevdft, dfile);
    /* Clear it from memory */
    dftEmpty(&prevdft);
  } else {
    /* Set file format to PMOD, if extension is .tac */
    if(!strcasecmp(filenameGetExtension(dfile), ".tac")) data._type=DFT_FORMAT_PMOD;
    /* save only present simulations */
    ret=dftWrite(&data, dfile);
  }
  if(ret) {
    fprintf(stderr, "Error in writing '%s': %s\n", dfile, dfterrmsg);
    dftEmpty(&data);
    return(12);
  }
  if(verbose>=0) fprintf(stdout, "Simulated TTAC(s) written in %s\n", dfile);
 
  
  /*
   *  Save simulated BTACs if requested
   */
  data.voi[0].voiname[0]=data.voi[0].hemisphere[0]=data.voi[0].place[0]=(char)0;
  data.voiNr=1;
  if(vofile[0]) {
    for(int i=0; i<data.frameNr; i++) data.voi[0].y[i]=data.voi[4].y2[i];
    strcpy(data.voi[0].name, "Venous_O2");
    if(dftWrite(&data, vofile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", vofile, dfterrmsg);
      dftEmpty(&data); return(21);
    }
    if(verbose>0) 
      fprintf(stdout, "Simulated venous [O-15]O2 BTAC written in %s\n", vofile);
  }
  if(vwfile[0]) {
    for(int i=0; i<data.frameNr; i++) data.voi[0].y[i]=data.voi[5].y2[i];
    strcpy(data.voi[0].name, "Venous_H2O");
    if(dftWrite(&data, vwfile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", vwfile, dfterrmsg);
      dftEmpty(&data); return(22);
    }
    if(verbose>0) 
      fprintf(stdout, "Simulated venous [O-15]H2O BTAC written in %s\n", vwfile);
  }
 
  dftEmpty(&data);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/