taccbv.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpcfileutil.h"
#include "tpcli.h"
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Subtracts the contribution of vascular radioactivity from regional",
  "PET TTACs. Vascular volume fraction Vb can be given as a value that is",
  "common to all regions, or as regional Vb values in a TAC or PAR file,",
  "calculated from a [O-15]CO study or estimated by model fitting.",
  " ",
  "Usage: @P ttacfile btacfile Vb [outputfile]",
  " ",
  "Options:",
  " -noneg",
  "     Negative TAC values are set to 0.",
  " -pv | -tv",
  "     Equation Ct=Cpet-Vb*Cb is applied by default or with option -pv;",
  "     with option -tv equation Ct=(Cpet-Vb*Cb)/(1-Vb) is applied.",
  " -sim",
  "     Simulate the contribution of vascular radioactivity instead of",
  "     correcting for it, calculating Cpet from Ct using equations above.",
  " --force",
  "     Program does not mind if the time or calibration units",
  "     cannot be converted to match, or if TAC names do not match.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1:",
  "     @P uo372.tac uo372ab.bld 0.045 uo372cbv.tac",
  "Example 2:",
  "     @P uo372.dft uo372ab.kbq uo372vb.dft uo372cbv.dft",
  " ",
  "Vb values that are >=1 are assumed to be percentages.",
  "Blood TAC can be given in a separate BTAC file, or as a region id inside",
  "TTAC file. Original TTAC file is modified, if output file is not given.",
  " ",
  "See also: imgcbv, p2blood, fitvb, tacadd, interpol, taccalc, tacunit",
  " ",
  "Keywords: TAC, modelling, vascular fraction, simulation",
  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;
  char    ttacfile[FILENAME_MAX], btacfile[FILENAME_MAX], outfile[FILENAME_MAX],
          parfile[FILENAME_MAX];
  int     leaveNegatives=1;
  int     petVolume=1;
  int     addVb=0;
  double  Vb=nan("");
  int     forceMode=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ttacfile[0]=btacfile[0]=outfile[0]=parfile[0]=(char)0;
 
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "NONEGATIVES", 1)==0) {
      leaveNegatives=0; continue;
    } else if(strncasecmp(cptr, "PV", 1)==0) {
      petVolume=1; continue;
    } else if(strncasecmp(cptr, "TV", 1)==0) {
      petVolume=0; continue;
    } else if(strncasecmp(cptr, "SIMULATE", 3)==0) {
      addVb=1; continue;
    } else if(strcasecmp(cptr, "F")==0 || strcasecmp(cptr, "FORCE")==0) {
      forceMode=1; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'\n", argv[ai]);
    return(1);
  } else break;
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-3;
  
  /* 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(ttacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(btacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(fileExist(argv[ai])) strlcpy(parfile, argv[ai++], FILENAME_MAX);
    else if(!atofCheck(argv[ai], &Vb) && Vb>=0.0) ai++;
    else {fprintf(stderr, "Error: invalid Vb: '%s'.\n", argv[ai]); return(1);}
  }
  if(ai<argc) strlcpy(outfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(isnan(Vb) && !parfile[0]) { // note that output file is optional
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(!outfile[0]) strcpy(outfile, ttacfile);
  if(Vb>=1.0) Vb/=100.;
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ttacfile := %s\n", ttacfile);
    printf("btacfile := %s\n", btacfile);
    printf("outfile := %s\n", outfile);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(!isnan(Vb)) printf("Vb := %g %%\n", 100.0*Vb);
    printf("petVolume := %d\n", petVolume);
    printf("leaveNegatives := %d\n", leaveNegatives);
    printf("addVb := %d\n", addVb);
    printf("forceMode := %d\n", forceMode);
    fflush(stdout);
  }
  if(verbose>1) {
    printf("\nApplying formula:\n");
    if(addVb==0) {
      if(petVolume==0) printf("Ct=(Cpet-Vb*Cb)/(1-Vb)\n");
      else             printf("Ct=Cpet-Vb*Cb\n");
    } else {
      if(petVolume==0) printf("Cpet=(1-Vb)*Ct+Vb*Cb\n");
      else             printf("Cpet=Ct+Vb*Cb\n");
    }
    printf("\n"); fflush(stdout);
  }
 
 
  /*
   *  Read the TTAC file
   */
  if(verbose>1) printf("reading %s\n", ttacfile);
  TAC ttac; tacInit(&ttac);
  if(tacRead(&ttac, ttacfile, &status) || tacSortByTime(&ttac, &status)) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); return(2);
  }
  if(verbose>2) {
    printf("ttac.fileformat := %s\n", tacFormattxt(ttac.format));
    printf("ttac.tacNr := %d\n", ttac.tacNr);
    printf("ttac.sampleNr := %d\n", ttac.sampleNr);
    printf("ttac.xunit := %s\n", unitName(ttac.tunit));
    printf("ttac.yunit := %s\n", unitName(ttac.cunit));
  }
 
 
  /*
   *  Read the BTAC file
   */
  if(verbose>1) printf("reading %s\n", btacfile);
  TAC btac; tacInit(&btac);
  if(tacRead(&btac, btacfile, &status) || tacSortByTime(&btac, &status)) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); return(3);
  }
  if(verbose>2) {
    printf("btac.fileformat := %s\n", tacFormattxt(btac.format));
    printf("btac.tacNr := %d\n", btac.tacNr);
    printf("btac.sampleNr := %d\n", btac.sampleNr);
    printf("btac.xunit := %s\n", unitName(btac.tunit));
    printf("btac.yunit := %s\n", unitName(btac.cunit));
  }
  /* Check the number of TACs */
  if(btac.tacNr>1) {
    fprintf(stderr, "Warning: BTAC file contains %d TACs.\n", btac.tacNr);
    fprintf(stdout, "Note: using the first TAC in BTAC file.\n");
    fflush(stderr); fflush(stdout);
    btac.tacNr=1;
  }
  /* Try to convert time units */
  {
    int ret=tacXUnitConvert(&btac, ttac.tunit, &status);
    if(ret!=TPCERROR_OK && verbose>1) {
      if(verbose>0) fprintf(stdout, "Note: TAC files have different or unknown time units.\n");
      if(verbose>2) fprintf(stderr, "Status: %s\n", errorMsg(status.error));
    }
    /* Error does not matter, if user told so */
    if(forceMode) ret=TPCERROR_OK;
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); return(3);
    }
  }
  /* Try to convert concentration units */
  {
    int ret=tacYUnitConvert(&btac, ttac.cunit, &status);
    if(ret!=TPCERROR_OK && verbose>1) {
      if(verbose>0) fprintf(stdout, "Note: TAC files have different or unknown concentration units.\n");
      if(verbose>2) fprintf(stderr, "Status: %s\n", errorMsg(status.error));
    }
    /* Error does not matter, if user told so */
    if(forceMode) ret=TPCERROR_OK;
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); return(3);
    }  
  }
  /* Interpolate BTAC to TTAC sample times */
  TAC bitac; tacInit(&bitac);
  if(tacInterpolate(&btac, &ttac, &bitac, NULL, NULL, &status) != TPCERROR_OK) {
    fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); tacFree(&bitac); return(3);
  }
  tacFree(&btac);
 
 
  /*
   *  If Vb values are given in a file, then read that
   */
  double VbArray[ttac.tacNr];
  if(parfile[0]) {
    if(verbose>1) printf("reading %s\n", parfile);
    PAR par; parInit(&par);
    int ret=parRead(&par, parfile, &status);
    if(ret!=TPCERROR_OK) {
      TAC tmp; tacInit(&tmp);
      ret=tacRead(&tmp, parfile, &status);
      if(ret==TPCERROR_OK) ret=tacToPAR(&tmp, &par, &status);
      tacFree(&tmp);
    }
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
      fprintf(stderr, "Error: cannot read Vb values from %s\n", parfile);
      tacFree(&ttac); tacFree(&bitac); parFree(&par);
      return(4);
    }
    if(verbose>2) {
      printf("\n Vb file contents:\n");
      parWrite(&par, stdout, PAR_FORMAT_CSV_UK, 0, NULL); fflush(stdout);
    }
    if(par.tacNr!=ttac.tacNr) {
      fprintf(stderr, "Error: different number of TACs in Vb file.\n");
      tacFree(&ttac); tacFree(&bitac); parFree(&par);
      return(4);
    }
    int VbIndex; if(par.parNr>1) VbIndex=parFindParameter(&par, "Vb"); else VbIndex=0;
    if(VbIndex<0) {
      fprintf(stderr, "Error: parameter Vb not identified in %s\n", parfile);
      tacFree(&ttac); tacFree(&bitac); parFree(&par);
      return(4);
    }
    /* check TAC names */
    ret=0;
    for(int ri=0; ri<par.tacNr; ri++) {
      if(strcasecmp(par.r[ri].name, ttac.c[ri].name)!=0) {
        if(verbose>1) {
          printf("  not matching: '%s' vs '%s'\n", par.r[ri].name, ttac.c[ri].name);
          fflush(stdout);
        }
        ret++;
      }
    }
    if(ret>0) {
      if(forceMode) fprintf(stderr, "Warning: Vb file does not contain the same TACs.\n");
      else {
        fprintf(stderr, "Error: Vb file does not contain the same TACs.\n");
        tacFree(&ttac); tacFree(&bitac); parFree(&par);
        return(4);
      }
    }
    /* Copy Vb values */
    for(int ri=0; ri<par.tacNr; ri++) VbArray[ri]=par.r[ri].p[VbIndex];
    parFree(&par);
    /* Check that Vb is fraction */
    double VbMax, VbMin; doubleRange(VbArray, ttac.tacNr, &VbMin, &VbMax);
    if(verbose>1) {printf("Vb range: %g - %g\n", VbMin, VbMax); fflush(stdout);}
    if(VbMin<0.0 || VbMax>100.0) {
      fprintf(stderr, "Error: invalid range of Vb values.\n");
      tacFree(&ttac); tacFree(&bitac);
      return(4);
    }
    if(VbMax>1.0) {
      for(int ri=0; ri<ttac.tacNr; ri++) VbArray[ri]*=0.01;
      fprintf(stderr, "Warning: Vb values were converted to fractions.\n"); fflush(stderr);
    }
  }
 
 
  /*
   *  Make subtraction, or add Vb contribution
   */
  {
    int ret=0;
    if(!isnan(Vb)) {
      ret=tacVb(&ttac, -1, &bitac, Vb, addVb, petVolume, &status);
    } else {
      for(int ri=0; ri<ttac.tacNr; ri++)
        if((ret=tacVb(&ttac, ri, &bitac, VbArray[ri], addVb, petVolume, &status))) break;
    }
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&bitac);
      return(8);
    }
    if(!leaveNegatives) {
      for(int ri=0; ri<ttac.tacNr; ri++) for(int fi=0; fi<ttac.sampleNr; fi++)
        if(ttac.c[ri].y[fi]<0.0) ttac.c[ri].y[fi]=0.0;
    }
  }
  tacFree(&bitac);
 
 
  /*
   *  Save data
   */
  {
    if(verbose>1) printf("writing %s\n", outfile);
    FILE *fp; fp=fopen(outfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", outfile);
      tacFree(&ttac); return(11);
    }
    int ret=tacWrite(&ttac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); tacFree(&ttac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      return(12);
    }
  }
 
  if(verbose>=0) printf("%s saved.\n", outfile);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/