sim_mbf.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 radiowater PET time-radioactivity concentration curves (TACs)",
  "for myocardial muscle and LV cavity based on MBF model by Iida et al (1, 2)",
  "as represented in ref 3.",
  " ",
  "Usage: @P [options] bloodfile MBF PTF Va Beta simfile",
  " ",
  "Options:",
  " -pH2O=<Partition coefficient for water>",
  "     Enter the partition coefficient of water; 0.9464 by default.",
  " -sub | -nosub",
  "     TAC of perfusbale tissue (Ct) is written (-sub)",
  "     or not written (-nosub, default) into the simfile.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "MBF must be given in units mL/(mL*min), and PTF and Va as fractions.",
  "TACs of regional muscle, LV cavity, and optionally perfusable tissue are ",
  "written in the given simfile, overwriting any previous contents.",
  " ",
  "For accurate results, blood input TAC should be noiseless and have very",
  "short sampling intervals. Simulated curves can thereafter",
  "be interpolated to represent PET frames using program simframe.",
  " ",
  "Example:",
  "   @P s3456ab.bld 1.24 0.6 0.3 0.8 s3456sim.tac",
  " ",
  "References:",
  "1. Iida H, Rhodes CG, de Silva R, Yamamoto Y, Araujo LI, Maseri A, Jones T.",
  "   Myocardial tissue fraction - correction for partial volume effects and",
  "   measure of tissue viability. J Nucl Med 1991; 32:2169-2175.",
  "2. Iida H, Rhodes CG, de Silva R, Araujo LI, Bloomfield P, Lammertsma AA,",
  "   Jones T. Use of the left ventricular time-activity curve as a noninvasive",
  "   input function in dynamic oxygen-15-water positron emission tomography.",
  "   J Nucl Med 1992; 33:1669-1677.",
  "3. Oikonen V. Model equations for myocardial perfusion studies with [15O]H2O",
  "   PET. https://www.turkupetcentre.net/reports/tpcmod0005.pdf",
  " ",
  "See also: fitmbf, b2t_h2o, sim_3tcm, tacadd, simframe, simimyoc",
  " ",
  "Keywords: simulation, myocardium, 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];
  double  PTF=-1.0, MBF=-1.0, Beta=-1.0, Va=-1.0, pc=0.9464;
  char   *cptr;
  int     n, ret;
  int     originallySec=0;
  int     save_ct=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  blofile[0]=simfile[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, "P=", 2)==0) {
      pc=atofVerified(cptr+2); if(pc>0.0 && pc<=1.0) continue;
    } else if(strncasecmp(cptr, "PH2O=", 5)==0) {
      pc=atofVerified(cptr+5); if(pc>0.0 && pc<=1.0) continue;
    } else if(strncasecmp(cptr, "NOSUB", 5)==0) {
      save_ct=0; continue;
    } else if(strncasecmp(cptr, "SUB", 3)==0) {
      save_ct=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 */
  if(ai<argc) {strlcpy(blofile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    MBF=atofVerified(argv[ai++]); 
    if(MBF<0.0 || MBF>20.0) {
      fprintf(stderr, "Error: invalid MBF (%g).\n", MBF); return(1);}
  }
  if(ai<argc) {
    PTF=atofVerified(argv[ai++]);
    if(PTF<=0.0 || PTF>1.0) {
      fprintf(stderr, "Error: invalid PTF (%g).\n", PTF); return(1);}
  }
  if(ai<argc) {
    Va=atofVerified(argv[ai++]); 
    if(Va<0.0 || Va>1.0) {
      fprintf(stderr, "Error: invalid Va (%g).\n", Va); return(1);}
  }
  if(ai<argc) {
    Beta=atofVerified(argv[ai++]); 
    if(Beta<=0.0 || Beta>1.0) {
      fprintf(stderr, "Error: invalid Beta (%g).\n", Beta); return(1);}
  }
  if(ai<argc) {strlcpy(simfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    /* we should never get this far */
    fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!simfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("blofile := %s\n", blofile);
    printf("simfile := %s\n", simfile);
    printf("MBF := %g\n", MBF);
    printf("PTF := %g\n", PTF);
    printf("Va := %g\n", Va);
    printf("Beta := %g\n", Beta);
    printf("pH2O := %g\n", pc);
  }
  
 
  /*
   *  Read arterial BTAC
   */
  if(verbose>1) fprintf(stdout, "reading %s\n", blofile);
  TAC input; tacInit(&input);
  ret=tacRead(&input, blofile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), blofile);
    tacFree(&input); return(2);
  }
  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(input.tacNr>1) {
    fprintf(stderr, 
     "Warning: input file contains %d TACs; using the first.\n", input.tacNr);
    input.tacNr=1;
  }
  if(input.sampleNr<2) {
    fprintf(stderr, "Error: too few samples in input TAC.\n");
    tacFree(&input); return(2);
  }
  if(input.sampleNr<10) {
    fprintf(stderr, "Warning: too few samples for reliable simulation.\n");
  }
  /* Guess time units, if necessary */
  if(input.tunit==UNIT_UNKNOWN) {
    if(input.x[input.sampleNr-1]>20.0) {
      if(verbose>1) printf("Note: assuming that times are in seconds.\n");
      input.tunit=UNIT_SEC;
    } else {
      if(verbose>1) printf("Note: assuming that times are in minutes.\n");
      input.tunit=UNIT_MIN;
    }
  }
  /* Convert time to min if necessary */
  if(input.tunit==UNIT_SEC) {
    tacXUnitConvert(&input, UNIT_MIN, &status);
    originallySec=1;
  }
 
 
  /*
   *  Allocate space for simulated data
   */
  TAC sim; tacInit(&sim);
  n=3;
  if(verbose>1) 
    printf("allocating memory for %d samples and %d TACs\n", input.sampleNr, n);
  ret=tacAllocate(&sim, input.sampleNr, n);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(ret));
    tacFree(&sim); tacFree(&input); return(4);
  }
  sim.sampleNr=input.sampleNr; 
  sim.tacNr=n;
  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];
  strcpy(sim.c[0].name, "muscle");
  strcpy(sim.c[1].name, "lvcav");
  strcpy(sim.c[2].name, "Ct");
 
 
  /*
   *  Simulation
   */
  if(verbose>1) printf("simulating\n");
  ret=simC1(input.x, input.c[0].y, input.sampleNr, MBF, MBF/pc, sim.c[2].y);
  if(ret) {
    fprintf(stderr, "Error: cannot simulate TTAC.\n");
    if(verbose>1) printf("ret := %d\n", ret);
    tacFree(&sim); tacFree(&input); return(6);
  }
 
  /*
   *  Calculate regional muscle and LV curves
   */
  for(int i=0; i<sim.sampleNr; i++) {
    /* Calculate muscle region TAC */
    sim.c[0].y[i] = PTF*sim.c[2].y[i] + Va*input.c[0].y[i];
    /* Calculate LV cavity region TAC */
    sim.c[1].y[i] = Beta*input.c[0].y[i] + (1.0-Beta)*sim.c[2].y[i];
  }
 
  /* simulation done */
  tacFree(&input);
 
 
  /* Convert time to sec if necessary */
  if(originallySec) {
    tacXUnitConvert(&sim, UNIT_SEC, &status);
  }
 
 
  /*
   *  Save simulated data 
   */
  if(verbose>1) printf("writing %s\n", simfile);
  if(save_ct==0) sim.tacNr=2;
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/