sim_h2o.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 "tpcfileutil.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpctacmod.h"
#include "tpccm.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulation of PET tissue time-radioactivity concentration curves (TACs)",
  "from arterial blood (Ca) TACs, based on compartmental model for radiowater.",
  "Model has two parameters, perfusion (blood flow, f) and partition",
  "coefficient of water (p).",
  "Delay (deltat) and dispersion (tau) of arterial input function (AIF),",
  "perfusable tissue fraction (PTF), and arterial volume fraction (Va) can be",
  "included in the simulation.",
  " ",
  "  dAIF(t)/dt = (1/tau)*Ca(T-deltat) - (1/tau)*AIF(T-deltat) ",
  "  dCt(t)/dt = f*AIF(T) - (f/p)*Ct(T) ",
  "  Cpet(T)= PTF*Ct(T) + Va*AIF(T) ",
  " ",
  "Usage: @P [options] parfile [bloodfile simfile]",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "To create a template parameter file, do not enter names for input and",
  "simulated TACs.",
  "If parameter file does not contain units, then per min and per mL units",
  "are assumed for f, p, PTF and Va, and seconds for deltat and tau.",
  "For accurate results, blood TAC should have very short sampling intervals.",
  " ",
  "See also: b2t_h2o, sim_3tcm, fit_h2o, fit_wrlv, tacadd, simframe",
  " ",
  "Keywords: TAC, simulation, compartmental model, perfusion, radiowater",
  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   blofile[FILENAME_MAX], simfile[FILENAME_MAX], parfile[FILENAME_MAX];
  PAR    par;
  int    ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  blofile[0]=simfile[0]=parfile[0]=(char)0;
  parInit(&par);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
  
  /* 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(parfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(blofile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(simfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!parfile[0]) {
    fprintf(stderr, "Error: missing parameter file; use option --help\n");
    return(1);
  }
  if(blofile[0] && !simfile[0]) {
    fprintf(stderr, "Error: missing file name.\n");
    return(1);
  }
  if(fileExist(parfile) && !simfile[0]) {
      fprintf(stderr, "Error: missing file names.\n");
      parFree(&par); return(1);
    }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("parfile := %s\n", parfile);
    printf("blofile := %s\n", blofile);
    printf("simfile := %s\n", simfile);
    fflush(stdout);
  }
 
 
  /*
   *  Make template parameter file, if no other file names were given, and then quit
   */
  if(!blofile[0]) {
    /* We have already checked that parameter file does not exist */
 
     /* Allocate space for 6 parameters, 4 TTACs */
    ret=parAllocate(&par, 6, 4);
    if(ret) {
      fprintf(stderr, "Error: cannot allocate memory for parameters.\n");
      parFree(&par); return(1);
    }
    par.parNr=6; par.tacNr=4;
    strcpy(par.n[0].name, "f"); par.n[0].unit=UNIT_ML_PER_ML_MIN;
    strcpy(par.n[1].name, "p"); par.n[1].unit=UNIT_ML_PER_ML;
    strcpy(par.n[2].name, "PTF"); par.n[2].unit=UNIT_ML_PER_ML;
    strcpy(par.n[3].name, "Va"); par.n[3].unit=UNIT_ML_PER_ML;
    strcpy(par.n[4].name, "deltaT"); par.n[4].unit=UNIT_SEC;
    strcpy(par.n[5].name, "tau"); par.n[5].unit=UNIT_SEC;
    for(int i=0; i<par.tacNr; i++) {
      sprintf(par.r[i].name, "tac%d", 1+i);
      par.r[i].model=modelCodeIndex("radiowater2");
    }
    int ri=0; // first TTAC
    par.r[ri].p[0]=1.0;
    par.r[ri].p[1]=1.0;
    par.r[ri].p[2]=0.6;
    par.r[ri].p[3]=0.2;
    par.r[ri].p[4]=0.0;
    par.r[ri].p[5]=0.0;
    ri++;
    par.r[ri].p[0]=4.0;
    par.r[ri].p[1]=1.0;
    par.r[ri].p[2]=0.6;
    par.r[ri].p[3]=0.2;
    par.r[ri].p[4]=0.0;
    par.r[ri].p[5]=0.0;
    ri++;
    par.r[ri].p[0]=0.2;
    par.r[ri].p[1]=0.8;
    par.r[ri].p[2]=0.97;
    par.r[ri].p[3]=0.03;
    par.r[ri].p[4]=5.0;
    par.r[ri].p[5]=0.0;
    ri++;
    par.r[ri].p[0]=0.2;
    par.r[ri].p[1]=0.8;
    par.r[ri].p[2]=0.97;
    par.r[ri].p[3]=0.03;
    par.r[ri].p[4]=5.0;
    par.r[ri].p[5]=5.0;
 
    if(verbose>1) printf("writing %s\n", parfile);
    FILE *fp; fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", parfile);
      parFree(&par); return(11);
    }
    ret=parWrite(&par, fp, PAR_FORMAT_TSV_UK, 1, &status);
    fclose(fp); parFree(&par); 
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot write %s\n", parfile);
      return(12);
    }
    if(verbose>=0) {printf("%s saved.\n", parfile); fflush(stdout);}
    return(0);
  }
 
  
  /*
   *  Read parameter file to get parameters for the simulation
   */
  if(verbose>1) fprintf(stdout, "reading %s\n", parfile);
  if(parRead(&par, parfile, &status)) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), parfile);
    parFree(&par);
    return(3);
  }
  if(verbose>5) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  /* Check the validity of model and its parameters */
  for(int i=0; i<par.tacNr; i++) {
    if(verbose>3) printf("model %s for tac %s\n", modelCode(par.r[i].model), par.r[i].name);
    if(par.r[i].model!=modelCodeIndex("radiowater2")) {
      fprintf(stderr, "Error: parameter file has invalid model '%s'.\n", modelDesc(par.r[i].model));
      parFree(&par); return(3);
    }
    /* check that we got at least as many parameters as the model has */
    if(par.parNr<(int)modelParNr(par.r[i].model)) {
      fprintf(stderr, "Error: invalid parameters for selected model.\n");
      parFree(&par); return(3);
    }
  }
  /* check that we can find the obligatory model parameters */
  int i_deltaT=-1, i_tau=-1, i_f=-1, i_p=-1, i_PTF=-1, i_Va=-1;
  {
    ret=0;
    if((i_f=parFindParameter(&par, "f"))<0) ret++;
    if((i_p=parFindParameter(&par, "p"))<0) ret++;
    if((i_PTF=parFindParameter(&par, "PTF"))<0) ret++;
    if((i_Va=parFindParameter(&par, "Va"))<0) ret++;
    if((i_deltaT=parFindParameter(&par, "deltaT"))<0 && 
       (i_deltaT=parFindParameter(&par, "delayT"))<0) ret++;
    if((i_tau=parFindParameter(&par, "tau"))<0) ret++;
    if(ret) {
      fprintf(stderr, "Error: required parameters not available.\n");
      parFree(&par); return(3);
    }
    if(verbose>5) {
      printf("parameter indices:\n");
      printf("  f=p[%d]\n", i_f);
      printf("  p=p[%d]\n", i_p);
      printf("  PTF=p[%d]\n", i_PTF);
      printf("  Va=p[%d]\n", i_Va);
      printf("  deltaT=p[%d]\n", i_deltaT);
      printf("  tau=p[%d]\n", i_tau);
    }
  }
 
 
  /*
   *  Read input TAC
   */
  if(verbose>1) fprintf(stdout, "reading input TAC\n");
  TAC input; tacInit(&input);
  ret=tacReadModelingInput(blofile, NULL, NULL, &input, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s (input files)\n", errorMsg(status.error));
    tacFree(&input); parFree(&par); return(4);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(input.format));
    printf("tacNr := %d\n", input.tacNr);
    printf("sampleNr := %d\n", input.sampleNr);
    printf("xunit := %s\n", unitName(input.tunit));
    printf("yunit := %s\n", unitName(input.cunit));
  }
  if(verbose>10) tacWrite(&input, stdout, TAC_FORMAT_PMOD, 1, NULL);
  if(input.sampleNr<3) {
    fprintf(stderr, "Error: too few samples in input TAC.\n");
    tacFree(&input); parFree(&par); return(4);
  }
  if(input.sampleNr<10) {
    fprintf(stderr, "Warning: too few samples for reliable simulation.\n"); fflush(stderr);
  }
  if(tacXUnitConvert(&input, UNIT_SEC, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&input); parFree(&par); return(4);
  }
 
 
  /*
   *  Allocate space for simulated data
   */
  if(verbose>1) fprintf(stdout, "allocating space for simulated TACs\n");
  TAC sim; tacInit(&sim);
  ret=tacAllocateWithPAR(&sim, &par, input.sampleNr, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&input); tacFree(&sim); parFree(&par); return(4);
  }
  sim.format=TAC_FORMAT_PMOD;
  sim.weighting=WEIGHTING_OFF;
  sim.cunit=input.cunit;
  sim.tunit=input.tunit;
  sim.format=input.format;
  sim.isframe=0;
  for(int j=0; j<sim.sampleNr; j++) sim.x[j]=input.x[j];
 
 
  /*
   *  Simulation
   */
  if(verbose>1) {printf("simulating...\n"); fflush(stdout);}
  for(int ri=0; ri<sim.tacNr; ri++) {
    if(verbose>2) {printf("  %s\n", sim.c[ri].name); fflush(stdout);}
    double cf;
 
    /* Make delayed and dispersed input TAC */
    double deltaT=par.r[ri].p[i_deltaT];
    cf=unitConversionFactor(par.n[i_deltaT].unit, UNIT_SEC); if(isfinite(cf)) deltaT*=cf;
    if(verbose>4) printf("    deltaT := %g s\n", deltaT);
 
    double tau=par.r[ri].p[i_tau];
    cf=unitConversionFactor(par.n[i_tau].unit, UNIT_SEC); if(isfinite(cf)) tau*=cf;
    if(verbose>4) printf("    tau := %g s\n", tau);
 
    double buf1[3*sim.sampleNr]; double *buf2=buf1+sim.sampleNr; double *ca=buf2+sim.sampleNr;
    doubleCopy(buf1, input.c[0].y, sim.sampleNr);
    if(simDispersion(sim.x, buf1, sim.sampleNr, tau, 0.0, buf2)) {
      fprintf(stderr, "Error: cannot add dispersion for '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
    for(int i=0; i<sim.sampleNr; i++) buf2[i]=input.x[i]+deltaT;
    if(liInterpolate(buf2, buf1, sim.sampleNr, sim.x, ca, NULL, NULL, sim.sampleNr, 3, 1, 0)!=0) {
      fprintf(stderr, "Error: cannot add delay time for '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
 
    /* Simulate pure tissue TAC */
    double f=par.r[ri].p[i_f];
    cf=unitConversionFactor(par.n[i_f].unit, UNIT_ML_PER_ML_SEC); 
    if(isfinite(cf)) f*=cf; else f/=60.0;
    if(verbose>4) printf("    f := %g mL/(s*mL)\n", f);
    if(!(f>=0.0)) {
      fprintf(stderr, "Error: invalid f for '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
 
    double p=par.r[ri].p[i_p];
    cf=unitConversionFactor(par.n[i_p].unit, UNIT_ML_PER_ML); if(isfinite(cf)) p*=cf;
    if(verbose>4) printf("    p := %g mL/mL\n", p);
    double k2=f/p; if(!(k2>=0.0)) {
      fprintf(stderr, "Error: invalid k2=f/p for '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
 
    if(simC1(sim.x, ca, sim.sampleNr, f, k2, sim.c[ri].y)) {
      fprintf(stderr, "Error: cannot simulate '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
 
    /* Multiply tissue with PTF */
    double PTF=par.r[ri].p[i_PTF];
    cf=unitConversionFactor(par.n[i_PTF].unit, UNIT_ML_PER_ML); if(isfinite(cf)) PTF*=cf;
    if(verbose>4) printf("    PTF := %g mL/mL\n", PTF);
    if(!(PTF>=0.0)) {
      fprintf(stderr, "Error: invalid PTF for '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
    for(int i=0; i<sim.sampleNr; i++) sim.c[ri].y[i]*=PTF;
 
    /* Add the contribution of arterial blood in tissue vasculature */
    double Va=par.r[ri].p[i_Va];
    cf=unitConversionFactor(par.n[i_Va].unit, UNIT_ML_PER_ML); if(isfinite(cf)) Va*=cf;
    if(verbose>4) printf("    Va := %g mL/mL\n", Va);
    if(!(Va>=0.0)) {
      fprintf(stderr, "Error: invalid Va for '%s'.\n", sim.c[ri].name);
      tacFree(&input); tacFree(&sim); parFree(&par); return(5);
    }
    for(int i=0; i<sim.sampleNr; i++) sim.c[ri].y[i]+=Va*ca[i];
 
  }
  /* simulation done */
  tacFree(&input); parFree(&par);
 
 
  /*
   *  Save simulated data 
   */
  if(verbose>1) printf("writing %s\n", simfile);
  FILE *fp; fp=fopen(simfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", simfile);
    tacFree(&sim); return(11);
  }
  ret=tacWrite(&sim, fp, TAC_FORMAT_PMOD, 1, &status);
  fclose(fp); tacFree(&sim);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    return(12);
  }
  if(verbose>=0) {printf("%s saved.\n", simfile); fflush(stdout);}
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/