sim_rtcm.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 "tpccm.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulation of PET tissue time-radioactivity concentration curves (TAC) from",
  "reference tissue TAC based on a selected reference tissue input compartment",
  "model.",
  "Available reference tissue input models, and model parameters:",
  "  FRTM; (Full) reference tissue model; R1, k2, k3, BP [1]",
  "  SRTM: Simplified reference tissue model; R1, k2, BP [2]",
  "  RRTM: Reduced reference tissue model; R1, k2, k3 (k4=0)",
  "  TRTM: Transport-limited reference tissue model; R1, k2, k3 [3]",
  " ",
  "Usage: @P [options] parfile [reffile simfile]",
  " ",
  "Options:",
  " -model=<FRTM|SRTM|RRTM|TRTM>",
  "     Select the model used to simulate the TACs; by default the model is",
  "     read from the parameter file.",
  " -sub | -nosub",
  "     TACs of sub-compartments (Cf and Cb) are also written (-sub) into",
  "     simfile, or not written (-nosub, default); effective only with FRTM.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "To create a template parameter file, do not enter names for reference TAC",
  "and simulated TACs, but give the model with option -model.",
  " ",
  "References:",
  "1. Cunningham VJ, Hume SP, Price GR, Ahier RG, Cremer JE, Jones AKP.",
  "   Compartmental analysis of diprenorphine binding to opiate receptors",
  "   in the rat in vivo and its comparison with equilibrium data in vitro.",
  "   J Cereb Blood Flow Metab 1991;11:1-9.",
  "2. Lammertsma AA, Hume SP. Simplified reference tissue model for PET",
  "   receptor studies. NeuroImage 1996;4:153-158.",
  "3. Herholz K, Lercher M, Wienhard K, Bauer B, Lenz O, Heiss W-D.",
  "   PET measurement of cerebral acetylcholine esterase activity without",
  "   blood sampling. Eur J Nucl Med 2001;28:472-477.",
  " ",
  "See also: sim_3tcm, fit_frtm, fit_srtm, fit_rrtm, fit_trtm, tacadd, simframe",
  " ",
  "Keywords: TAC, simulation, reference tissue, 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;
/*****************************************************************************/
 
/*****************************************************************************/
//enum {MODEL_UNKNOWN, MODEL_FRTM, MODEL_RRTM, MODEL_SRTM, MODEL_TRTM};
//static char *model_str[] = {"Unknown","FRTM","RRTM","SRTM","TRTM",0};
//static int model_parNr[6]={0,4,3,3,3,0};
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int    ai, help=0, version=0, verbose=1;
  char   reffile[FILENAME_MAX], simfile[FILENAME_MAX], parfile[FILENAME_MAX];
  char  *cptr;
  unsigned int model=0, parNr=0;
  int    save_only_total=1;
  TAC    tac, sim;
  PAR    par;
  int    n, ret;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacInit(&tac); tacInit(&sim); parInit(&par);
  reffile[0]=simfile[0]=parfile[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(strncasecmp(cptr, "MODEL=", 6)==0 && strlen(cptr)>7) {
      cptr+=6;
      model=modelCodeIndex(cptr);
      //for(int i=1; i<5; i++) if(strcasecmp(cptr, model_str[i])==0) {
      //  model=i; break;}
      if(model==0 && strcasecmp(cptr, "RTCM")==0) model=modelCodeIndex("FRTM");
      if(model!=0) continue;
      fprintf(stderr, "Error: invalid model selection.\n"); return(1);
    } else if(strncasecmp(cptr, "SUB", 3)==0) {
      save_only_total=0; continue;
    } else if(strncasecmp(cptr, "NOSUB", 3)==0) {
      save_only_total=1; 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 */
  for(; ai<argc; ai++) {
    if(!parfile[0]) {
      strlcpy(parfile, argv[ai], FILENAME_MAX); continue;
    } else if(!reffile[0]) {
      strlcpy(reffile, argv[ai], FILENAME_MAX); continue;
    } else if(!simfile[0]) {
      strlcpy(simfile, argv[ai], FILENAME_MAX); continue;
    }
    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(reffile[0] && !simfile[0]) {
    fprintf(stderr, "Error: missing filename for simulated TACs.\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    if(model!=0) printf("model := %s\n", modelCode(model));
    printf("parfile := %s\n", parfile);
    printf("reffile := %s\n", reffile);
    printf("simfile := %s\n", simfile);
    printf("save_only_total := %d\n", save_only_total);
  }
  
  /*
   *  Make template parameter file, if no other filenames were given,
   *  and then quit
   */
  if(!reffile[0]) {
    parNr=modelParNr(model);
    if(model==0 || parNr==0) {
      fprintf(stderr, "Error: unknown simulation model; use option -model.\n");
      return(1);
    }
    ret=parAllocate(&par, parNr, 3);
    if(ret) {
      fprintf(stderr, "Error: cannot allocate memory for parameters.\n");
      parFree(&par); return(1);
    }
    par.parNr=parNr; par.tacNr=3;
    for(int i=0; i<par.tacNr; i++) {
      sprintf(par.r[i].name, "tac%d", 1+i);
      par.r[i].model=model;
      //par.r[i].start=0.0;
      //par.r[i].end=90.0;
      par.r[i].p[0]=1.0;
      par.r[i].p[1]=0.2;
    }
    strcpy(par.n[0].name, "R1"); par.n[0].unit=UNIT_UNITLESS;
    strcpy(par.n[1].name, "k2"); par.n[1].unit=UNIT_PER_MIN;
    if(!strcmp(modelCode(model), "FRTM")) {
      strcpy(par.n[2].name, "k3"); par.n[2].unit=UNIT_PER_MIN;
      strcpy(par.n[3].name, "BPnd"); par.n[3].unit=UNIT_UNITLESS;
      par.r[0].p[2]=0.05;
      par.r[0].p[3]=1.0;
      for(int i=1; i<par.tacNr; i++) {
        par.r[i].p[2]=2.0*par.r[i-1].p[2];
        par.r[i].p[3]=2.0*par.r[i-1].p[3];
      }
    } else if(!strcmp(modelCode(model), "SRTM")) {
      strcpy(par.n[2].name, "BPnd"); par.n[2].unit=UNIT_UNITLESS;
      par.r[0].p[2]=1.0;
      for(int i=1; i<par.tacNr; i++) {
        par.r[i].p[2]=2.0*par.r[i-1].p[2];
      }
    } else if(!strcmp(modelCode(model), "RRTM") || 
              !strcmp(modelCode(model), "TRTM")) 
    {
      strcpy(par.n[2].name, "k3"); par.n[2].unit=UNIT_PER_MIN;
      par.r[0].p[2]=0.05;
      for(int i=1; i<par.tacNr; i++) {
        par.r[i].p[2]=2.0*par.r[i-1].p[2];
      }
    }
    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_CSV_UK, 1, &status);
    fclose(fp);
    parFree(&par); return(0);
  }
  
  /*
   *  Read parameter file to get parameters for the simulation, and 
   *  set model in case it was given with command-line option.
   */
  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);
 
  if(model!=0) for(int i=0; i<par.tacNr; i++) par.r[i].model=model;
  /* 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);
    /* check that we got at least as many parameters as the model has */
    if(modelParNr(par.r[i].model)>(unsigned int)par.parNr) {
      fprintf(stderr, "Error: invalid parameters for selected model.\n");
      parFree(&par); return(3);
    }
    /* check that we can find the obligatory model parameters */
    ret=0;
    if(parFindParameter(&par, "R1")<0) ret++;
    if(parFindParameter(&par, "k2")<0) ret++;
    if(!strcmp(modelCode(par.r[i].model), "FRTM")) {
      if(parFindParameter(&par, "k3")<0) ret++;
      if(parFindParameter(&par, "BPnd")<0) ret++;
    } else if(!strcmp(modelCode(par.r[i].model), "SRTM")) {
      if(parFindParameter(&par, "BPnd")<0) ret++;
    } else if(!strcmp(modelCode(par.r[i].model), "RRTM") || 
              !strcmp(modelCode(par.r[i].model), "TRTM")) 
    {
      if(parFindParameter(&par, "k3")<0) ret++;
    } else {
      fprintf(stderr, "Error: invalid model for this simulation.\n");
      if(verbose>1) printf("model := %s\n", modelCode(par.r[i].model));
      parFree(&par); return(3);
    }
    if(ret) {
      fprintf(stderr, "Error: required parameters not available.\n");
      parFree(&par); return(3);
    }
  }
 
  /*
   *  Read reference input TAC
   */
  if(verbose>1) fprintf(stdout, "reading %s\n", reffile);
  ret=tacRead(&tac, reffile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), reffile);
    tacFree(&tac); parFree(&par); return(4);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
  }
  if(tac.tacNr>1) {
    fprintf(stderr, 
     "Warning: reference file contains %d TACs; using the first.\n", tac.tacNr);
    tac.tacNr=1;
  }
  if(tac.sampleNr<2) {
    fprintf(stderr, "Error: too few samples in reference TAC.\n");
    tacFree(&tac); parFree(&par); return(4);
  }
  if(tac.sampleNr<10) {
    fprintf(stderr, "Warning: too few samples for reliable simulation.\n");
  }
 
  /*
   *  Allocate space for simulated data
   */
  n=par.tacNr; if(!save_only_total) n*=3;
  if(verbose>1) 
    printf("allocating memory for %d samples and %d TACs\n", tac.sampleNr, n);
  ret=tacAllocate(&sim, tac.sampleNr, n);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(ret));
    tacFree(&sim); tacFree(&tac); parFree(&par); return(5);
  }
  sim.sampleNr=tac.sampleNr; 
  sim.tacNr=n;
  sim.weighting=WEIGHTING_OFF;
  sim.cunit=tac.cunit;
  sim.tunit=tac.tunit;
  sim.format=tac.format;
  sim.isframe=0;
  for(int j=0; j<sim.sampleNr; j++) sim.x[j]=tac.x[j];
 
  /*
   *  Simulation
   */
  if(verbose>1) printf("simulating\n");
  n=0; // sim tac struct index
  double *px=tac.x;
  double *py=tac.c[0].y;
  double *pf, *pb, *pt;
  double R1, k2, k3, k4, BP;
  int i, ii;
  for(i=0; i<par.tacNr; i++) {
    if(verbose>2) printf("simulating %s\n", sim.c[i].name);
    ii=parFindParameter(&par, "R1"); if(ii>=0) R1=par.r[i].p[ii]; else R1=0.0;
    ii=parFindParameter(&par, "k2"); if(ii>=0) k2=par.r[i].p[ii]; else k2=0.0;
    ii=parFindParameter(&par, "k3"); if(ii>=0) k3=par.r[i].p[ii]; else k3=0.0;
    ii=parFindParameter(&par, "k4"); if(ii>=0) k4=par.r[i].p[ii]; else k4=0.0;
    ii=parFindParameter(&par, "BPnd"); if(ii>=0) BP=par.r[i].p[ii]; else BP=0.0;
    if(fabs(BP)>1.0E-50) k4=k3/BP; else if(fabs(k4)>1.0E-50) BP=k3/k4;
    if(verbose>3) {printf("R1 := %g\n", R1); printf("k2 := %g\n", k2);}
    strcpy(sim.c[n].name, par.r[i].name);
    pt=sim.c[n].y; pf=pb=NULL;
    ret=0;
    if(!strcmp(modelCode(par.r[i].model), "FRTM")) {
      if(verbose>3) {printf("k3 := %g\n", k3); printf("k4 := %g\n", k4);}
      if(!save_only_total) {
        strcpy(sim.c[n+1].name, par.r[i].name); strcat(sim.c[n+1].name, "_free");
        pf=sim.c[n+1].y; 
        strcpy(sim.c[n+2].name, par.r[i].name); strcat(sim.c[n+2].name, "_bound");
        pb=sim.c[n+2].y;
      }
      ret=simRTCM(px, py, tac.sampleNr, R1, k2, k3, k4, pt, pf, pb);
      if(save_only_total) n++; else n+=3;
    } else if(!strcmp(modelCode(par.r[i].model), "SRTM")) {
      if(verbose>3) {printf("BP := %g\n", BP);}
      ret=simSRTM(px, py, tac.sampleNr, R1, k2, BP, pt);
      n++;
    } else if(!strcmp(modelCode(par.r[i].model), "RRTM")) {
      if(verbose>3) {printf("k3 := %g\n", k3);}
      ret=simRTCM(px, py, tac.sampleNr, R1, k2, k3, 0.0, pt, NULL, NULL);
      n++;
    } else if(!strcmp(modelCode(par.r[i].model), "TRTM")) {
      if(verbose>3) {printf("k3 := %g\n", k3);}
      ret=simTRTM(px, py, tac.sampleNr, R1, k2, k3, pt);
      n++;
    }
    if(ret) {
      fprintf(stderr, "Error: invalid data for simulation.\n");
      if(verbose>1) printf("sim_return_code := %d\n", ret);
      tacFree(&sim); tacFree(&tac); parFree(&par); return(8);
    }
  }
  sim.tacNr=n;
  /* simulation done */
  tacFree(&tac); parFree(&par);
 
 
 
  /*
   *  Write useful header information
   */
  ii=iftFindKey(&par.h, "studynr", 0);
  if(ii<0) ii=iftFindKey(&par.h, "study", 0);
  if(ii>=0) {
    iftPut(&sim.h, par.h.item[ii].key, par.h.item[ii].value, 1, NULL);
  }
 
 
  /*
   *  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_UNKNOWN, 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);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/