o2_p2w.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"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Converts arterial plasma TAC from a [O-15]O2 PET study",
  "to blood [O-15]H2O TAC (1).",
  " ",
  "Usage: @P [Options] plasmafile bloodfile",
  " ",
  "Options:",
  " -HCT=<Hematocrit>",
  "     Specify the hematocrit; when specified, the plasma-to-blood ratio",
  "     is calculated using it, assuming that water contents of RBCs and",
  "     plasma are 0.73 and 0.94, respectively.",
  " -PBR=<Plasma-to-blood ratio>",
  "     Specify the plasma-to-blood ratio; by default 1.12; effective only",
  "     when HCT is not given.",
  " -SCT=<Sample collect time>",
  "     Plasma sample collection time (usually 20 s) should be given for",
  "     accurate fit with fit_o2bl.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "References:",
  "1. Lubberink M, Wong YY, Raijmakers PGHM, Schuit RC, Luurtsema G,",
  "   Boellaard R, Knaapen P, Vonk-Noordegraaf A, Lammertsma AA. Myocardial",
  "   oxygen extraction fraction measured using bolus inhalation of 15O-oxygen",
  "   gas and dynamic PET. J Nucl Med. 2011;52(1):60-66.",
  " ",
  "See also: fit_o2bl, o2metab, sim_o2bl, b2t_mo2, fit_mo2, tac2svg, taccalc",
  " ",
  "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    pfile[FILENAME_MAX], bfile[FILENAME_MAX], *cptr;
  double  pb_ratio=1.12, bp_ratio, HCT=-1.0, coll_time=-1.0;
  TAC     tac;
  double  w_p=0.94, w_rbc=0.73;
  int     times_converted=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  pfile[0]=bfile[0]=(char)0;
  tacInit(&tac);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1;
    if(strncasecmp(cptr, "PBR=", 4)==0) {
      if(!atofCheck(cptr+4, &pb_ratio) && pb_ratio>=1.0 && pb_ratio<2.0)
        continue;
    } else if(strncasecmp(cptr, "HCT=", 4)==0) {
      if(!atofCheck(cptr+4, &HCT) && HCT>0.0 && HCT<1.0) {
        pb_ratio=-1.0; continue;
      }
    } else if(strncasecmp(cptr, "SCT=", 4)==0) {
      if(!atofCheck(cptr+4, &coll_time) && coll_time>=0.0 && coll_time<60.0)
        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);}
  
  /* Arguments */
  for(; ai<argc; ai++) {
    if(!pfile[0]) {
      strlcpy(pfile, argv[ai], FILENAME_MAX); continue;
    } else if(!bfile[0]) {
      strlcpy(bfile, argv[ai], FILENAME_MAX); continue;
    }
    fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!bfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    for(ai=0; ai<argc; ai++)
      printf("%s ", argv[ai]);
    printf("\n");
    printf("bfile := %s\n", bfile);
    printf("pfile := %s\n", pfile);
    if(HCT<=0.0) printf("pb_ratio := %g\n", pb_ratio);
    else printf("HCT := %g\n", HCT);
    if(coll_time>=0.0) printf("coll_time := %g [s]\n", coll_time);
  }
 
  /* Calculate plasma-to-blood ratio from HCT, if necessary */
  /* And calculate blood-to-plasma ratio */
  if(HCT>0.0) {
    if(verbose>1) {
      printf("w_rbc := %g\n", w_rbc);
      printf("w_plasma := %g\n", w_p);
    }
    bp_ratio= 1.0 - HCT*(1.0 - w_rbc/w_p);
    pb_ratio=1.0/bp_ratio;
    if(verbose>1) printf("pb_ratio := %g\n", pb_ratio);
  } else {
    bp_ratio=1.0/pb_ratio;
  }
  if(verbose>1) printf("bp_ratio := %g\n", bp_ratio);
 
 
  /*
   *  Read plasma TAC
   */
  if(verbose>1) printf("reading %s\n", pfile);
  ret=tacRead(&tac, pfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  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: only first TAC in plasma file is used.\n");
    tac.tacNr=1;
  }
  /* Sort by sample times */
  tacSortByTime(&tac, NULL);
  /* Check for missing values */
  ret=tacNaNs(&tac);
  if(ret>0) {
    if(verbose>1) printf("missing concentrations.\n");
    /* Try to fix missing concentrations */
    ret=tacFixNaNs(&tac);
    if(ret!=0) {
      fprintf(stderr, "Error: missing concentrations in %s.\n", pfile);
      tacFree(&tac); return(2);
    }
  }
  /* Convert sample times from min to sec, if necessary */
  if(tac.tunit==UNIT_UNKNOWN) {
    /* If sample collect time is longer than TAC end time, then
       TAC is probably in minutes */
    if(coll_time<tac.x[tac.sampleNr-1]) {
      tac.tunit=UNIT_MIN;
      fprintf(stderr, "Warning: sample times are assumed to be in minutes.\n");
    } else {
      fprintf(stderr, "Warning: sample times are assumed to be in seconds.\n");
      tac.tunit=UNIT_SEC;
    }
  }
  if(tac.tunit==UNIT_MIN) {
    ret=tacXUnitConvert(&tac, UNIT_SEC, &status);
    if(!ret) {
      if(verbose>1) printf("Sample times are converted to sec\n");
      times_converted=1;
    }
  }
 
  /*
   *  Multiply plasma data with bp_ratio
   */
  for(int i=0; i<tac.sampleNr; i++) tac.c[0].y[i]*=bp_ratio;
 
  /*
   *  Set time "frames" based on the collection time
   */
  if(coll_time>0.0) {
    if(verbose>1) printf("setting collection time\n");
    for(int i=0; i<tac.sampleNr; i++) {
      /* Check that sample times are not overlapping */
      if(i>0 && ((tac.x[i]-tac.x[i-1])<coll_time)) {
        fprintf(stderr, "Error: sample and collection times do not match.\n");
        tacFree(&tac); return(4);
      }
      /* Set frame start and end times based on collection time */
      tac.x1[i]=tac.x[i]-0.5*coll_time;
      tac.x2[i]=tac.x[i]+0.5*coll_time;
    }
    /* Set TAC type accordingly */
    tac.isframe=1; /* start and end times */
    if(tac.format==TAC_FORMAT_SIMPLE) tac.format=TAC_FORMAT_PMOD;
  }
 
  strcpy(tac.c[0].name, "Blood_H2O");
 
  /*
   *  Convert sample times back to minutes, if necessary
   */
  if(times_converted) {
    if(verbose>1) printf("Sample times are converted back to min\n");
    tacXUnitConvert(&tac, UNIT_MIN, &status);
    times_converted=0;
  }
 
 
  /*
   *  Write blood file
   */
  if(verbose>1) printf("writing %s\n", bfile);
  FILE *fp; fp=fopen(bfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", bfile);
    tacFree(&tac); return(11);
  }
  ret=tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
  fclose(fp); tacFree(&tac);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    return(12);
  }
  if(verbose>=0) printf("Blood curve written in %s\n", bfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/