o2metab.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculates the blood time-activity curves (BTACs) of labelled water and",
  "oxygen from arterial blood TAC measured after [O-15]O2 bolus using",
  "a dedicated [O-15]O2 metabolite model (1, 2, 3) to be used as input in",
  "[O-15]O2 model.",
  " ",
  "Usage: @P [Options] btacfile parameterfile obtacfile wtacfile",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Model parameters for k3 model (k1, k1+k3, delay) or k4 model (k1, k3,",
  "k3/k4, delay) can be given in a RES or IFT file made by program fit_o2bl.",
  "Parameters and data units are assumed to be in seconds, in case units",
  "are not given in the files.",
  " ",
  "References:",
  "1. Huang S-C et al. Modelling approach for separating blood time-activity",
  "   curves in positron emission tomographic studies.",
  "   Phys Med Biol. 1991;36(6):749-761.",
  "2. Iida H et al. Modelling approach to eliminate the need to separate",
  "   arterial plasma in oxygen-15 inhalation positron emission tomography.",
  "   J Nucl Med. 1993;34:1333-1340.",
  "3. Kudomi N et al. A physiologic model for recirculation water correction",
  "   in CMRO2 assessment with 15O2 inhalation PET.", 
  "   J Cereb Blood Flow Metab. 2009;29(2):355-364.",
  " ",
  "See also: fit_o2bl, sim_o2bl, o2_p2w, fit_mo2, b2t_mo2, taccalc, tac2svg",
  " ",
  "Keywords: input, oxygen, blood, metabolite correction",
  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      ret;
  char     blofile[FILENAME_MAX], parfile[FILENAME_MAX],
           oxyfile[FILENAME_MAX], watfile[FILENAME_MAX];
  double   k1, k2, k3, k4, delay, k3k4, k1k3;
  DFT      blood, sim;
  RES      res;
  IFT      ift;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dftInit(&blood); dftInit(&sim); resInit(&res); iftInit(&ift);
  blofile[0]=parfile[0]=oxyfile[0]=watfile[0]=(char)0;
  k1=k2=k3=k4=k3k4=k1k3=nan(""); delay=0.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==NULL) 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(!blofile[0]) {strlcpy(blofile, argv[ai], FILENAME_MAX); continue;}
    if(!parfile[0]) {strlcpy(parfile, argv[ai], FILENAME_MAX); continue;}
    if(!oxyfile[0]) {strlcpy(oxyfile, argv[ai], FILENAME_MAX); continue;}
    if(!watfile[0]) {strlcpy(watfile, argv[ai], FILENAME_MAX); continue;}
    fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!watfile[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("parfile := %s\n", parfile);
    printf("oxyfile := %s\n", oxyfile);
    printf("watfile := %s\n", watfile);
  }
 
  /* Read parameter file */
  if(verbose>1) printf("reading parameters in %s\n", parfile);
  /* First, try RES format */
  if(resRead(parfile, &res, verbose-2)==0) {
    if(verbose>2) printf("result format\n");
    ret=res2ift(&res, &ift, verbose-3);
    resEmpty(&res);
    if(ret) {
      fprintf(stderr, "Error in converting results.\n");
      dftEmpty(&blood); dftEmpty(&sim); iftEmpty(&ift);
      return(3);
    }
  /* If not RES then try IFT */
  } else if(iftRead(&ift, parfile, 1, 0)==0) {
    if(verbose>2) printf("ift format\n");
  } else {
    fprintf(stderr, "Error: cannot read '%s'\n", parfile);
    dftEmpty(&blood); dftEmpty(&sim); iftEmpty(&ift);
    return(2);
  }
  /* Now try to find the required parameters */
  int i, n; char key[12]; double v; int u; char buf1[128], buf2[128];
  /* k1 */
  strcpy(key, "k1"); i=iftGet(&ift, key, 0);
  if(i<0) {
    fprintf(stderr, "Error: k1 not found in '%s'.\n", parfile);
    dftEmpty(&blood); dftEmpty(&sim); iftEmpty(&ift);
    return(4);
  }
  n=sscanf(ift.item[i].value, "%128s %128s", buf1, buf2);
  if(n==0 || atof_with_check(buf1, &v)) {
    fprintf(stderr, "Error: valid k1 not found in '%s'.\n", parfile);
    dftEmpty(&blood); dftEmpty(&sim); iftEmpty(&ift);
    return(4);
  }
  if(n>1) {u=petCunitId(buf2); if(u==CUNIT_PER_MIN) v/=60.0;}
  k1=v; if(verbose>2) printf("k1 := %g\n", k1);
  /* delay; if not found, keep it at zero */
  strcpy(key, "delay"); i=iftGet(&ift, key, 0);
  if(i<0) strcpy(key, "delayt");
  i=iftGet(&ift, key, 0);
  if(i<0) strcpy(key, "timedelay");
  i=iftGet(&ift, key, 0);
  if(i>=0) {
    n=sscanf(ift.item[i].value, "%128s %128s", buf1, buf2);
    if(n==0 || atof_with_check(buf1, &v)) {
      fprintf(stderr, "Error: valid delay not found in '%s'.\n", parfile);
      dftEmpty(&blood); dftEmpty(&sim); iftEmpty(&ift);
      return(4);
    }
    if(n>1) {u=petTunitId(buf2); if(u==TUNIT_MIN) v*=60.0;}
    delay=v; if(verbose>2) printf("delay := %g\n", k1);
  }
  /* k3 */
  strcpy(key, "k3"); i=iftGet(&ift, key, 0);
  if(i>=0) {
    n=sscanf(ift.item[i].value, "%128s %128s", buf1, buf2);
    if(n>0 && !atof_with_check(buf1, &v)) {
      if(n>1) {u=petCunitId(buf2); if(u==CUNIT_PER_MIN) v/=60.0;}
      k3=v; if(verbose>2) printf("k3 := %g\n", k3);
    }
  }
  /* k3/k4 */
  strcpy(key, "k3/k4"); i=iftGet(&ift, key, 0);
  if(i>=0) {
    n=sscanf(ift.item[i].value, "%128s %128s", buf1, buf2);
    if(n>0 && !atof_with_check(buf1, &v)) {
      k3k4=v; if(verbose>2) printf("k3/k4 := %g\n", k3k4);
    }
  }
  /* k1+k3 */
  strcpy(key, "k1+k3"); i=iftGet(&ift, key, 0);
  if(i>=0) {
    n=sscanf(ift.item[i].value, "%128s %128s", buf1, buf2);
    if(n>0 && !atof_with_check(buf1, &v)) {
      if(n>1) {u=petCunitId(buf2); if(u==CUNIT_PER_MIN) v/=60.0;}
      k1k3=v; if(verbose>2) printf("k1+k3 := %g\n", k1k3);
    }
  }
  /* check what we got */
  k2=k1;
  if(isnan(k3) && isnan(k1k3)) {
    fprintf(stderr, "Error: missing k3 or k1+k3 in '%s'.\n", parfile);
    dftEmpty(&blood); dftEmpty(&sim); iftEmpty(&ift);
    return(4);
  }
  if(!isnan(k3)) {
    if(isnan(k3k4)) {
      fprintf(stderr, "Error: missing k3/k4 in '%s'.\n", parfile);
      dftEmpty(&blood); dftEmpty(&sim); iftEmpty(&ift);
      return(4);
    }
    if(k3k4<=1.0E-100) k4=0.0; else k4=k3/k3k4;
  }
  if(isnan(k3)) {k3=k1k3-k1; k4=0.0;}
  if(verbose>1) {
    printf("k2 := %g\n", k2);
    printf("k3 := %g\n", k3);
    printf("k4 := %g\n", k4);
  }
  if(k1<0.0 || k3<0.0 || k4<0.0) {
    fprintf(stderr, "Error: invalid parameter value(s) in '%s'.\n", parfile);
    dftEmpty(&blood); dftEmpty(&sim); iftEmpty(&ift);
    return(4);
  }
  iftEmpty(&ift);
 
 
  /*
   *  Read blood file
   */
  if(verbose>1) printf("reading %s\n", blofile);
  if(dftRead(blofile, &blood)) {
    fprintf(stderr, "Error in reading '%s': %s\n", blofile, dfterrmsg);
    return(5);
  }
  if(blood.voiNr>1) {
    blood.voiNr=1;
    fprintf(stderr, "Warning: input file contains extra column(s).\n");
  }
  if(blood.frameNr<3) {
    fprintf(stderr, "Error: blood data not valid.\n");
    dftEmpty(&blood); dftEmpty(&sim);
    return(5);
  }
  /* Make sure that sample times are in seconds */
  int times_changed=0;
  if(blood.timeunit==TUNIT_UNKNOWN && blood.x[blood.frameNr-1]<15.0) {
    fprintf(stderr, "Warning: assuming that BTAC sample times are in minutes.\n");
    blood.timeunit=TUNIT_MIN;
  }
  if(blood.timeunit==TUNIT_MIN) {
    dftMin2sec(&blood); times_changed=1;
  }
 
 
  /*
   *  Allocate memory for simulated water and oxygen TACs
   */
  if(dftdup(&blood, &sim)) {
    fprintf(stderr, "Error: cannot allocate memory for simulated TACs.\n");
    dftEmpty(&blood); dftEmpty(&sim);
    return(6);
  }
  strcpy(sim.voi[0].hemisphere, ""); strcpy(sim.voi[0].place, "");
 
  /* 
   *  Simulate the water TAC
   */
  if(verbose>1) printf("simulating\n");
  ret=simC3s(blood.x, blood.voi[0].y, blood.frameNr, k1, k2, k3, k4, 0, 0,
    blood.voi[0].y2, blood.voi[0].y3, NULL, NULL);
  if(ret) {
    fprintf(stderr, "Error in simulation of metabolism (%d).\n", ret);
    dftEmpty(&blood); dftEmpty(&sim); return(7);
  }
  /* Correct simulated water TAC for delay */
  for(int i=0; i<blood.frameNr; i++) blood.x1[i]=blood.x[i]+delay;
  blood.timetype=DFT_TIME_MIDDLE;
  ret=interpolate(blood.x1, blood.voi[0].y3, blood.frameNr, sim.x,
    sim.voi[0].y, sim.voi[0].y2, NULL, sim.frameNr);
  if(ret) {
    fprintf(stderr, "Error in simulation of delay (%d).\n", ret);
    dftEmpty(&blood); dftEmpty(&sim); return(8);
  }
  /* Write the file */
  if(verbose>1) printf("writing O15-water TAC in %s\n", watfile);
  if(times_changed) dftSec2min(&sim);
  strcpy(sim.voi[0].voiname, "Water");
  strcpy(sim.voi[0].name, sim.voi[0].voiname);
  if(dftWrite(&sim, watfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", watfile, dfterrmsg);
    dftEmpty(&blood); dftEmpty(&sim); return(11);
  }
  if(verbose==1) printf("O15-water BTAC written in %s\n", watfile);
 
  /*
   *  Calculate the oxygen-15 blood TAC
   */
  /* Subtract delayed water TAC from blood curve; this equals O15-oxygen  */
  for(int i=0; i<sim.frameNr; i++)
    sim.voi[0].y[i]=blood.voi[0].y[i]-sim.voi[0].y[i];
  /* Write the file */
  if(verbose>1) printf("writing O15-oxygen TAC in %s\n", oxyfile);
  strcpy(sim.voi[0].voiname, "Oxygen");
  strcpy(sim.voi[0].name, sim.voi[0].voiname);
  if(dftWrite(&sim, oxyfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", oxyfile, dfterrmsg);
    dftEmpty(&blood); dftEmpty(&sim); return(12);
  }
  if(verbose==1) printf("O15-oxygen BTAC written in %s\n", oxyfile);
 
 
  dftEmpty(&blood); dftEmpty(&sim); 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/