b2rbc.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 "tpcli.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculates the TAC of blood cells (erythrocytes, RBC) using",
  "haematocrit (HCT), blood TAC, and plasma TAC, based on equation",
  " ",
  "  Blood = HCT*RBC + (1-HCT)*Plasma ",
  " ",
  "HCT is normally between 0.40-0.51 in men and 0.36-0.47 in women.",
  " ",
  "Usage: @P [Options] bloodfile plasmafile HCT rbcfile",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: bpr2cpr, p2blood, b2plasma, taccalc, fit_bpr",
  " ",
  "Keywords: input, modelling, simulation, RBC, plasma, blood",
  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], rfile[FILENAME_MAX];
  TAC    plasma, blood;
  double HCT=-1.0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  pfile[0]=bfile[0]=rfile[0]=(char)0;
  tacInit(&blood); tacInit(&plasma);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    //cptr=argv[ai]+1;
    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(!bfile[0]) {
      strlcpy(bfile, argv[ai], FILENAME_MAX); continue;
    } else if(!pfile[0]) {
      strlcpy(pfile, argv[ai], FILENAME_MAX); continue;
    } else if(HCT<0.0) {
      if(!atofCheck(argv[ai], &HCT)) {
        if(HCT>1.0) HCT/=100.0;
        if(HCT>0.0 && HCT<=1.0) continue;
      }
      fprintf(stderr, "Error: invalid HCT.\n"); 
      return(1);
    } else if(!rfile[0]) {
      strlcpy(rfile, argv[ai], FILENAME_MAX); continue;
    }
    fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!rfile[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("pfile := %s\n", pfile);
    printf("bfile := %s\n", bfile);
    printf("rfile := %s\n", rfile);
    printf("HCT := %g\n", HCT);
  }
 
 
  /*
   *  Read blood TAC
   */
  if(verbose>1) printf("reading %s\n", bfile);
  ret=tacRead(&blood, bfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&blood); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(blood.format));
    printf("tacNr := %d\n", blood.tacNr);
    printf("sampleNr := %d\n", blood.sampleNr);
    printf("xunit := %s\n", unitName(blood.tunit));
    printf("yunit := %s\n", unitName(blood.cunit));
  }
  if(blood.tacNr>1) {
    fprintf(stderr, "Warning: only first TAC in blood file is used.\n");
    blood.tacNr=1;
  }
 
  /*
   *  Read plasma TAC
   */
  if(verbose>1) printf("reading %s\n", pfile);
  ret=tacRead(&plasma, pfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&blood); tacFree(&plasma); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(plasma.format));
    printf("tacNr := %d\n", plasma.tacNr);
    printf("sampleNr := %d\n", plasma.sampleNr);
    printf("xunit := %s\n", unitName(plasma.tunit));
    printf("yunit := %s\n", unitName(plasma.cunit));
  }
  if(plasma.tacNr>1) {
    fprintf(stderr, "Warning: only first TAC in plasma file is used.\n");
    plasma.tacNr=1;
  }
 
  /* Sort by sample times */
  tacSortByTime(&blood, NULL); tacSortByTime(&plasma, NULL);
 
  /* Check for missing values */
  ret=tacNaNs(&blood);
  if(ret>0) {
    if(verbose>1) printf("missing blood concentrations.\n");
    /* Try to fix missing concentrations */
    ret=tacFixNaNs(&blood);
    if(ret!=0) {
      fprintf(stderr, "Error: missing concentrations in %s.\n", bfile);
      tacFree(&blood); tacFree(&plasma); return(2);
    }
  }
  ret=tacNaNs(&plasma);
  if(ret>0) {
    if(verbose>1) printf("missing plasma concentrations.\n");
    /* Try to fix missing concentrations */
    ret=tacFixNaNs(&plasma);
    if(ret!=0) {
      fprintf(stderr, "Error: missing concentrations in %s.\n", pfile);
      tacFree(&blood); tacFree(&plasma); return(2);
    }
  }
 
  /* Convert units */
  ret=tacXUnitConvert(&plasma, blood.tunit, &status);
  if(ret!=TPCERROR_OK && verbose>0) {
    fprintf(stderr, "Warning: different or unknown time units.\n");
  }
  ret=tacYUnitConvert(&plasma, blood.cunit, &status);
  if(ret!=TPCERROR_OK && verbose>0) {
    fprintf(stderr, "Warning: different or unknown concentration units.\n");
  }
 
  /* Check the sample nr */
  if(plasma.sampleNr!=blood.sampleNr ||
     plasma.x[plasma.sampleNr-1] != blood.x[blood.sampleNr-1])
  {
    fprintf(stderr, "Warning: different sample times.\n");
  }
 
  /* Interpolate plasma data to the same sample times as blood data */
  ret=tacAllocateMore(&blood, 1);
  if(!ret) 
    ret=liInterpolate(plasma.x, plasma.c[0].y, plasma.sampleNr,
                      blood.x, blood.c[1].y, NULL, NULL, blood.sampleNr, 
                      4, 1, verbose-2);
  if(ret) {
    fprintf(stderr, "Error: cannot interpolate plasma TAC.\n");
    tacFree(&blood); tacFree(&plasma); return(4);
  }
 
  /* Calculate RBC curve */
  for(int i=0; i<blood.sampleNr; i++)
    blood.c[0].y[i] = (blood.c[0].y[i] - (1.0-HCT)*blood.c[1].y[i]) / HCT;
  tacFree(&plasma);
 
  /* 
   *  Save RBC TAC 
   */
  if(verbose>1) printf("writing %s\n", rfile);
  blood.tacNr=1; // just to be sure
  FILE *fp; fp=fopen(rfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", rfile);
    tacFree(&blood); return(11);
  }
  ret=tacWrite(&blood, fp, TAC_FORMAT_UNKNOWN, 1, &status);
  fclose(fp); tacFree(&blood);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    return(12);
  }
  if(verbose>=0) printf("RBC curve written in %s\n", rfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/