liverinp.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[] = {
  "Input for liver models consists of two parts: portal vein and hepatic artery.",
  "This program calculates a weighted average of TACs of portal vein and",
  "hepatic artery.",
  " ",
  "Usage: @P [Options] portalfile arteryfile portal_fraction outputfile",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "  @P ue65pv.blo ue65ha.blo 0.90 ue65hi.blo",
  " ",
  "See also: liverpv, taccalc, interpol, taccat, tacformat, b2plasma",
  " ",
  "Keywords: input, blood, TAC, modelling, simulation, liver",
  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    fi, ret, dataNr=0;
  char   pfile[FILENAME_MAX], afile[FILENAME_MAX], wfile[FILENAME_MAX];
  double portalf=-1.0;
  TAC    pdata, adata;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  pfile[0]=afile[0]=wfile[0]=(char)0;
  tacInit(&pdata); tacInit(&adata);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==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(!afile[0]) {
      strlcpy(afile, argv[ai], FILENAME_MAX); continue;
    } else if(portalf<0.0) {
      if(atofCheck(argv[ai], &portalf)==0) {
        if(portalf>1.0) portalf/=100.0;
        if(portalf>=0.0 && portalf<=1.0) continue;
      }
      fprintf(stderr, "Error: invalid proportion of portal vein input.\n");
      return(1);
    } else if(!wfile[0]) {
      strlcpy(wfile, argv[ai], FILENAME_MAX); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!wfile[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("afile := %s\n", afile);
    printf("wfile := %s\n", wfile);
    printf("portalf := %g\n", portalf);
  }
 
 
  /*
   *  Read TAC files
   */
 
  /* Read portal vein data */
  if(verbose>0) printf("reading %s\n", pfile);
  ret=tacRead(&pdata, pfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&pdata); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(pdata.format));
    printf("tacNr := %d\n", pdata.tacNr);
    printf("sampleNr := %d\n", pdata.sampleNr);
    printf("xunit := %s\n", unitName(pdata.tunit));
    printf("yunit := %s\n", unitName(pdata.cunit));
  }
  if(pdata.tacNr>1) {
    fprintf(stderr, "Warning: only first TAC in portal data is used.\n");
    pdata.tacNr=1;
  }
 
  /* Read hepatic artery data */
  if(verbose>0) printf("reading %s\n", afile);
  ret=tacRead(&adata, afile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&pdata); tacFree(&adata); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(adata.format));
    printf("tacNr := %d\n", adata.tacNr);
    printf("sampleNr := %d\n", adata.sampleNr);
    printf("xunit := %s\n", unitName(adata.tunit));
    printf("yunit := %s\n", unitName(adata.cunit));
  }
  if(adata.tacNr>1) {
    fprintf(stderr, "Warning: only first TAC in hepatic arterial data is used.\n");
    adata.tacNr=1;
  }
 
  /* Sort by sample times */
  tacSortByTime(&pdata, NULL);
  tacSortByTime(&adata, NULL);
 
  /* Convert units */
  ret=tacXUnitConvert(&adata, pdata.tunit, &status);
  if(ret!=TPCERROR_OK && verbose>0) {
    fprintf(stderr, "Warning: different or unknown time units.\n");
  }
  ret=tacYUnitConvert(&adata, pdata.cunit, &status);
  if(ret!=TPCERROR_OK && verbose>0) {
    fprintf(stderr, "Warning: different or unknown concentration units.\n");
  }
 
  /* Check that hepatic artery TAC has higher peak */
  {
    double pmax, amax;
    ret=tacYRange(&pdata, 0, NULL, &pmax, NULL, NULL, NULL, NULL);
    if(ret==0) ret=tacYRange(&adata, 0, NULL, &amax, NULL, NULL, NULL, NULL);
    if(ret) {
      fprintf(stderr, "Error: invalid data.\n");
      tacFree(&pdata); tacFree(&adata); return(2);
    }
    if(amax<pmax) {
      fprintf(stderr, "Warning: unexpectedly portal vein has higher peak than hepatic artery.\n");
    }
  }
 
  /* Interpolate hepatic artery data to same sample times as portal vein data */
  ret=tacAllocateMore(&pdata, 1);
  if(!ret) 
    ret=liInterpolate(adata.x, adata.c[0].y, adata.sampleNr,
                      pdata.x, pdata.c[1].y, NULL, NULL, pdata.sampleNr, 4, 1, verbose-2);
  if(ret) {
    fprintf(stderr, "Error: cannot interpolate plasma TAC.\n");
    tacFree(&pdata); tacFree(&adata); return(4);
  }
 
 
  /*
   *  Calculate weighted average TAC
   */
  if(verbose>1) printf("calculating weighted average TAC\n");
  for(fi=0, dataNr=0; fi<pdata.sampleNr; fi++) {
    if(isnan(pdata.c[0].y[fi]) || isnan(pdata.c[1].y[fi])) {
      pdata.c[0].y[fi]=nan(""); continue;
    }
    if(pdata.x[fi]<=(1.1*adata.x[adata.sampleNr-1])) {
      pdata.c[0].y[fi]= portalf*pdata.c[0].y[fi] + (1.0-portalf)*pdata.c[1].y[fi];
      dataNr++;
    } // else keep all data as portal TAC
  }
  if(dataNr<2) {
    fprintf(stderr, "Error: check the TAC sample times.\n");
    tacFree(&pdata); tacFree(&adata);
    return(6);
  }
  if(dataNr+1<pdata.sampleNr) {
    fprintf(stderr,
      "Warning: last %d samples are taken from hepatic vein alone.\n",
      pdata.sampleNr-dataNr
    );
  }
 
 
  /*
   *  Save average TAC
   */
  if(verbose>1) printf("writing %s\n", wfile);
  {
    char tmp[256];
    tpcProgramName(argv[0], 1, 0, tmp, 128);
    iftPut(&pdata.h, "program", tmp, 1, NULL);
    filenameRmPath(pfile);
    iftPut(&pdata.h, "portal_vein_file", pfile, 1, NULL);
    filenameRmPath(afile);
    iftPut(&pdata.h, "hepatic_artery_file", afile, 1, NULL);
    sprintf(tmp, "%g", portalf);
    iftPut(&pdata.h, "portal_vein_fraction", tmp, 1, NULL);
  }
  pdata.tacNr=1; // just to be sure
  FILE *fp; fp=fopen(wfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", wfile);
    tacFree(&pdata); tacFree(&adata); return(11);
  }
  ret=tacWrite(&pdata, fp, TAC_FORMAT_UNKNOWN, 1, &status);
  fclose(fp); tacFree(&pdata); tacFree(&adata);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    return(12);
  }
  if(verbose>0) printf("Liver input TAC written in %s\n", wfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/