b2t_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 "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulation of PET tissue time-radioactivity concentration curve (TTAC)",
  "in [O-15]H2O PET studies from decay corrected arterial blood TAC based on",
  "the one-tissue compartment model.",
  " ",
  "Usage: @P [options] bloodfile f p E Vb fA simfile",
  " ",
  "Perfusion (f) must be given in units mL/(min*100 mL), partition coefficient",
  "(p) in units mL/mL, and vascular volume (Vb) and the arterial fraction of it",
  "(fA) as percentages.",
  "By default, blood flow (perfusion) and tissue concentratios are",
  "assumed to represented per PET volume (including vascular volume).",
  "Extraction coefficient, E=1-exp(-PS/f), can usually be set to 1.0.",
  "Blood sample times must be in seconds, unless specified in the file.",
  "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 fr4sim.",
  " ",
  "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.",
  " -fpt",
  "     Blood flow (perfusion) is assumed to be given per perfusable tissue",
  "     volume excluding vascular volume. TAC will still be simulated per",
  "     regional PET volume including vascular volume.",
  " -voiname=<text>",
  "     Enter a name (1-6 chars without spaces) for the simulated TAC", 
  " -vena=<filename>",
  "     Save the simulated venous blood TAC",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Simulated TACs are written in ASCII format with columns:",
  "  1) Sample time (s)",
  "  2) Total tissue activity concentration (Cpet)",
  "  3) Activity concentration in tissue, (1-Vb)*Ct (optional)",
  "  4) Arterial contribution to tissue activity, Vb*fA*Cab (optional)",
  "  5) Venous contribution to tissue activity, Vb*(1-fA)*Cvb (optional)",
  " ",
  "See also: fit_h2o, sim_mbf, sim_3tcm, avgbolus, tacadd, simframe, tacunit",
  " ",
  "Keywords: TAC, simulation, modelling, 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;
  int        ri, fi, ret;
  int        save_only_total=1;
  int        flow_per_tissue=0;
  int        add_to_previous=0;
  int        tunit;
  char       bfile[FILENAME_MAX], tfile[FILENAME_MAX], vfile[FILENAME_MAX];
  char      *cptr, voiname[MAX_REGIONNAME_LEN+1];
  double     Flow, pH2O, E, Vb, fA, Va, Vv;
  double     K1, k2;
  DFT        blood, dft;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  bfile[0]=tfile[0]=vfile[0]=voiname[0]=(char)0;
  Flow=pH2O=E=Vb=fA=Va=Vv=-1.0;
  dftInit(&blood); dftInit(&dft);
  /* 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) {
      ret=strlcpy(voiname, cptr+8, MAX_REGIONNAME_LEN);
      if(ret<1 || ret>=MAX_REGIONNAME_LEN) {
        fprintf(stderr, "Error: invalid VOI name '%s'.\n", cptr+8); return(1);}
      continue;
    } else if(strncasecmp(cptr, "VENA=", 5)==0) {
      ret=strlcpy(vfile, cptr+5, FILENAME_MAX);
      if(ret<1 || ret>=FILENAME_MAX) {
        fprintf(stderr, "Error: invalid VOI name '%s'.\n", cptr+5); return(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 */
  if(ai<argc) {strlcpy(bfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    Flow=atof_dpi(argv[ai])/6000.;
    if(!(Flow>=0.0)) {fprintf(stderr, "invalid f (%s).\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    pH2O=atof_dpi(argv[ai]);
    if(!(pH2O>=0.0)) {fprintf(stderr, "invalid pH2O (%s).\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    E=atof_dpi(argv[ai]);
    if(!(E>=0.0) || E>1.0) {fprintf(stderr, "invalid E (%s).\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    Vb=atof_dpi(argv[ai])/100.;
    if(!(Vb>=0.0) || Vb>1.0) {fprintf(stderr, "invalid Vb (%s).\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    fA=atof_dpi(argv[ai])/100.;
    if(!(fA>=0.0) || fA>1.0) {fprintf(stderr, "invalid fA (%s).\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {strlcpy(tfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]); return(1);}
 
  /* Did we get all the information that we need? */
  if(!tfile[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
  
  /* Calculate Va and Vb */
  Va=fA*Vb; Vv=(1.0-fA)*Vb;
  if(E<=0. || E>1. || Va<0. || Va>=1. || Vv<0. || Vv>=1. || Va+Vv>=1.) {
    fprintf(stderr, "Erroneous parameter values.\n");
    return(1);
  }
  
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("bfile := %s\n", bfile);
    printf("tfile := %s\n", tfile);
    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);
    if(vfile[0]) printf("vfile := %s\n", vfile);
    if(voiname[0]) printf("voiname := %s\n", voiname);
    printf("Flow := %g\n", Flow);
    printf("pH2O := %g\n", pH2O);
    printf("E := %g\n", E);
    printf("Vb := %g\n", Vb);
    printf("fA := %g\n", fA);
    printf("Va := %g\n", Va);
    printf("Vv := %g\n", Vv);
  }
 
 
  /*
   *  Read blood data
   */
  if(verbose>1) printf("reading blood TAC in %s\n", bfile);
  if(dftRead(bfile, &blood)) {
    fprintf(stderr, "Error in reading '%s': %s\n", bfile, dfterrmsg);
    dftEmpty(&blood); return(2);
  }
  if(blood.frameNr<2) {
    fprintf(stderr, "Error: too few samples in '%s'.\n", bfile);
    dftEmpty(&blood); return(2);
  }
  if(blood.frameNr<4) {
    fprintf(stderr, "Warning: few samples in '%s'.\n", bfile);
  }
  if(blood.voiNr>1) {
    fprintf(stderr, "Error: blood data contains more than one curve.\n");
    dftEmpty(&blood); return(2);
  }
  if(blood.timeunit==TUNIT_UNKNOWN) {
    blood.timeunit=TUNIT_SEC;
    if(verbose>=0) printf("assuming that blood time unit is sec\n");
  }
  tunit=blood.timeunit; // save original time units
  if(blood.timeunit!=TUNIT_SEC) {
    ret=dftTimeunitConversion(&blood, TUNIT_SEC);
    if(ret) {
      fprintf(stderr, "Error: cannot convert blood data sample times to sec.\n");
      dftEmpty(&blood); return(2);
    }
  }
  strcpy(blood.voi[0].voiname, "Blood");
  strcpy(blood.voi[0].name, blood.voi[0].voiname);
 
 
  /*
   *  Allocate memory for simulated data
   */
  if(verbose>1) printf("allocating memory for simulated TACs\n");
  ret=dftSetmem(&dft, blood.frameNr, 4);
  if(ret) {
    fprintf(stderr, "Error in allocating memory for simulated data.\n");
    dftEmpty(&blood); dftEmpty(&dft); return(3);
  }
  dft.frameNr=blood.frameNr; dft.voiNr=4;
  /* Copy times & header info */
  dftCopymainhdr(&blood, &dft);
  for(fi=0; fi<dft.frameNr; fi++) {
    dft.x[fi]=blood.x[fi]; dft.x1[fi]=blood.x1[fi]; dft.x2[fi]=blood.x2[fi];
  }
  strcpy(dft.isotope, "O-15");
  strcpy(dft.radiopharmaceutical, "[O-15]H2O");
  strcpy(dft.injectionTime, blood.injectionTime);
  strcpy(dft.scanStartTime, blood.scanStartTime);
  dft.timeunit=blood.timeunit;
  dft.decayCorrected=DFT_DECAY_CORRECTED;
  for(ri=0; ri<dft.voiNr; ri++) {
    strcpy(dft.voi[ri].place, "");
    dft.voi[ri].size=100.;
  }
  if(strlen(voiname)==0) {
    strcpy(dft.voi[0].voiname, "Total"); strcpy(dft.voi[0].hemisphere, "Cpet");
  } else {
    strcpy(dft.voi[0].name, voiname);
    rnameSplit(voiname, dft.voi[0].voiname, dft.voi[0].hemisphere, 
               dft.voi[0].place, MAX_REGIONSUBNAME_LEN);
  }
  strcpy(dft.voi[1].voiname, "Tissue"); strcpy(dft.voi[1].hemisphere, "Ct");
  strcpy(dft.voi[2].voiname, "Artery"); strcpy(dft.voi[2].hemisphere, "Ca");
  dft.voi[2].size=100.*Va;
  strcpy(dft.voi[3].voiname, "Vena"); strcpy(dft.voi[3].hemisphere, "Cv");
  dft.voi[3].size=100.*Vv;
  for(ri=1; ri<dft.voiNr; ri++)
    strcpy(dft.voi[ri].name, dft.voi[ri].voiname);
 
 
  /*
   *  Simulate TACs
   */
  if(verbose>1) printf("simulating tissue data\n");
  K1=Flow*E; if(flow_per_tissue==0) K1/=(1.0-Va-Vv);
  k2=K1/pH2O;
  if(verbose>1) printf("K1 := %g\nk2 := %g\n", K1, k2);
  ret=simC3s(blood.x, blood.voi[0].y, blood.frameNr, K1, k2, 0., 0., 0., 0.,
             dft.voi[1].y, NULL, NULL, NULL);
  if(ret) {
    fprintf(stderr, "Error %d in simulating tissue data.\n", ret);
    dftEmpty(&blood); dftEmpty(&dft); return(5);
  }
  if(verbose>2) printf("simulating PET curves\n");
  for(fi=0; fi<dft.frameNr; fi++) {
    dft.voi[2].y[fi]=blood.voi[0].y[fi];
    dft.voi[3].y[fi]=(1.0-E)*blood.voi[0].y[fi] + (E/pH2O)*dft.voi[1].y[fi];
  }
 
 
  /* 
   *  Save vena TAC if required    
   */
  if(vfile[0]) {
    if(verbose>0) printf("saving vena data in %s\n", vfile);
    DFT_NR_OF_DECIMALS=6;
    DFT temp;
    dftInit(&temp); dftdup(&dft, &temp); dftCopyvoi(&temp, 3, 0); temp.voiNr=1;
    dftTimeunitConversion(&temp, tunit);
    dftSetComments(&temp);
    if(ret!=0 || dftWrite(&temp, vfile)!=0) {
      fprintf(stderr, "Error: cannot save vena data in %s\n", vfile);
      if(verbose>1) fprintf(stderr, "Error: %s\n", dfterrmsg);
      dftEmpty(&temp); dftEmpty(&blood); dftEmpty(&dft); return(11);
    }
    dftEmpty(&temp);
  }
 
 
  /* Calculate the activities per PET volume */
  for(fi=0; fi<dft.frameNr; fi++) {
    dft.voi[1].y[fi]*=(1.0-Va-Vv);
    dft.voi[2].y[fi]*=Va;
    dft.voi[3].y[fi]*=Vv;
    dft.voi[0].y[fi]=dft.voi[1].y[fi]+dft.voi[2].y[fi]+dft.voi[3].y[fi];
  }
  /* Set data info */
  sprintf(dft.comments,
    "# Flow := %g\n# pH2O :=%g\n# E := %g\n# Va := %g\n# Vv := %g\n",
    Flow, pH2O, E, Va, Vv);
  /* Blood is not needed anymore */
  dftEmpty(&blood);
  /* Convert sample times to original units */
  dftTimeunitConversion(&dft, tunit);
 
 
  /*
   *  Write simulated TACs
   */
  if(verbose>1) printf("saving PET curves\n");
  DFT_NR_OF_DECIMALS=6;
  if(save_only_total) dft.voiNr=1;
  /* Set file format to PMOD, if extension is .tac */
  if(add_to_previous==0 && !strcasecmp(filenameGetExtension(tfile), ".tac"))
    dft._type=DFT_FORMAT_PMOD;
  /* Some format has to be set anyway, and simple format would lose information */
  if(dft._type==DFT_FORMAT_PLAIN || dft._type==DFT_FORMAT_UNKNOWN)
    dft._type=DFT_FORMAT_STANDARD;
  dftSetComments(&dft);
  if(add_to_previous!=0 && access(tfile, 0)!=-1) {
    if(verbose>0) printf("adding to existing file\n");
    DFT prevdft;
    /* read previous file */
    dftInit(&prevdft);
    if(dftRead(tfile, &prevdft)) {
      fprintf(stderr, "Error in reading '%s': %s\n", tfile, dfterrmsg);
      dftEmpty(&prevdft); dftEmpty(&dft); return(2);
    }
    if(verbose>2) printf("existing %d TACs with %d frames\n",
      prevdft.voiNr, prevdft.frameNr);
    /* add present simulations */
    for(ri=0; ri<dft.voiNr; ri++) {
      ret=dftAdd(&prevdft, &dft, ri);
      if(ret) {
        fprintf(stderr, "Error: cannot add simulated TAC in %s\n", tfile);
        dftEmpty(&prevdft); dftEmpty(&dft); return(13);
      }
    }
    /* save that */
    ret=dftWrite(&prevdft, tfile);
    /* Clear it from memory */
    dftEmpty(&prevdft);
  } else {
    /* save only present simulations */
    ret=dftWrite(&dft, tfile);
  }
  if(ret) {
    fprintf(stderr, "Error in writing '%s': %s\n", tfile, dfterrmsg);
    dftEmpty(&dft);
    return(12);
  }
  if(verbose>=0) fprintf(stdout, "simulated TAC(s) written in %s\n", tfile);
 
  /* Free memory */
  dftEmpty(&dft);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/