taccbvp.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 "tpcmodels.h"
#include "tpcpar.h"
#include "tpcfunc.h"
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Estimate Vb and Ct based on blood curve (BTAC), assuming that at BTAC peak",
  "tissue concentration is zero, thus Vb=Ct/Cb.",
  "Vb is then reduced to assure that other TACs remain non-negative.",
  "The name or number of BTAC inside TAC file must be given.",
  " ",
  "Do not use this oversimplified method for quantitative analyses!",
  " ",
  "Usage: @P [options] tacfile BTAC [outputfile]",
  " ",
  "Options:",
  " -par=<filename>",
  "     Estimated Vb is written in file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imgcbvp, taccbv, tacpeak, taccalc",
  " ",
  "Keywords: TAC, peak, vascular fraction",
  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   tacfile[FILENAME_MAX], parfile[FILENAME_MAX], corfile[FILENAME_MAX];
  char   bname[MAX_TACNAME_LEN+1];
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=parfile[0]=corfile[0]=(char)0;
  bname[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, "PAR=", 4)==0) {
      strlcpy(parfile, cptr+4, FILENAME_MAX); if(strlen(parfile)>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 */
  if(ai<argc) strlcpy(tacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(bname, argv[ai++], MAX_TACNAME_LEN+1);
  if(ai<argc) strlcpy(corfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing? */
  if(!bname[0]) {tpcPrintUsage(argv[0], info, stdout); 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("tacfile := %s\n", tacfile);
    printf("bname := %s\n", bname);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(corfile[0]) printf("corfile := %s\n", corfile);
  }
 
 
  /*
   *  Read the file
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", 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, "Error: invalid file.\n");
    tacFree(&tac); return(2);
  } else if(tac.tacNr<2) {
    if(verbose>0) fprintf(stderr, "Error: file contains just one TAC.\n");
    tacFree(&tac); return(2);
  }
  /* Check NaNs */
  if(tacNaNs(&tac)>0) {
    fprintf(stderr, "Error: data contains missing values.\n");
    tacFree(&tac); return(2);
  }
  /* Sort the data by sample times (x) */
  if(tacSortByTime(&tac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  /* Get x range */
  double xmin, xmax;
  if(tacXRange(&tac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&tac); return(2);
  }
  if(verbose>1) {
    printf("xmin := %g\n", xmin);
    printf("xmax := %g\n", xmax);
  }
  /* Take average of any duplicate samples */
  if(tacMultipleSamples(&tac, 1, &tac, verbose-2)!=0) {
    fprintf(stderr, "Error: cannot process duplicate samples.\n");
    tacFree(&tac); return(2);
  }
  if(tac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&tac); return(2);
  }
 
  /* Select the LV BTAC */
  int ib=-1;
  {
    int n=tacSelectTACs(&tac, bname, 1, NULL);
    if(n<1) {
      fprintf(stderr, "Error: specified BTAC not found in %s.\n", tacfile);
      tacFree(&tac); return(3);
    } else if(n>1) {
      fprintf(stderr, "Error: %d TACs match '%s' in %s.\n", n, bname, tacfile);
      tacFree(&tac); return(3);
    }
    ib=tacFirstSelected(&tac);
    if(verbose>1) {
      printf("BTAC name := %s\n", tac.c[ib].name);
      printf("BTAC index := %d\n", ib);
    }
  }
 
 
  /*
   *  Prepare space for results
   */
  if(verbose>1) printf("preparing space for parameters\n");
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &tac, 1, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(3);
  }
  par.tacNr=tac.tacNr;
  par.parNr=1;
  par.format=PAR_FORMAT_TSV_UK;
  for(int i=0; i<tac.tacNr; i++) {
    par.r[i].model=0;
    par.r[i].dataNr=tacWSampleNr(&tac);
    par.r[i].start=xmin; 
    par.r[i].end=xmax;
  }
  /* Set parameter names */
  strcpy(par.n[0].name, "Vb");
 
 
  /*
   *  Search BTAC peak.
   *  Note that duplicate samples were removed above, and data is sorted by x.
   */
  if(verbose>1) printf("searching peak\n");
  double bpeak=tac.c[ib].y[0];
  int ipeak=0;
  for(int i=1; i<tac.sampleNr; i++)
    if(tac.c[ib].y[i]>bpeak) {bpeak=tac.c[ib].y[i]; ipeak=i;}
  if(verbose>1) {
    printf("  maxv := %g\n", bpeak);
    printf("  maxx := %g\n", tac.x[ipeak]);
  }
  /* If peak is not the first sample, then add previous sample */
  if(ipeak>0) bpeak+=tac.c[ib].y[ipeak-1];
  if(verbose>2) {
    printf("  bpeak := %g\n", bpeak);
  }
 
  /* 
   *  Estimate Vb and correct TACs
   */
  for(int ci=0; ci<tac.tacNr; ci++) if(ci!=ib) {
 
    /* Set any negative values to zero */
    for(int i=0; i<tac.sampleNr; i++) if(tac.c[ci].y[i]<0.0) tac.c[ci].y[i]=0.0;
 
    /* Calculate the maximum of Vb that is possible in this pixel from T/B ratio */
    double tpeak=tac.c[ci].y[ipeak];
    if(ipeak>0) tpeak+=tac.c[ci].y[ipeak-1];
    double vb=tpeak/bpeak;
    if(verbose>2) printf("  tpeak[%d] := %g\n", ci, tpeak);
 
    if(vb>=1.0) { // this pixel must be mostly blood
      for(int i=0; i<tac.sampleNr; i++) tac.c[ci].y[i]=0.0;
      par.r[ci].p[0]=vb;
      continue;
    }
 
    /* Make initial BV corrected TTAC */
    double ct[tac.sampleNr], ctmin=0.0; 
    int ictmin=0;
    for(int i=0; i<tac.sampleNr; i++) {
      ct[i]=tac.c[ci].y[i]-vb*tac.c[ib].y[i];
      if(ct[i]<ctmin) {ctmin=ct[i]; ictmin=i;}
    }
    /* Make Vb smaller to prevent negative TAC values and correct again */
    if(ctmin<0.0) {
      vb+=ctmin/tac.c[ib].y[ictmin];
      for(int i=0; i<tac.sampleNr; i++) {
        ct[i]=tac.c[ci].y[i]-vb*tac.c[ib].y[i];
        if(ct[i]<0.0) ct[i]=0.0;
      }
    }
    /* Put corrected TAC in place */
    for(int i=0; i<tac.sampleNr; i++) tac.c[ci].y[i]=ct[i];
    /* Store the Vb */
    par.r[ci].p[0]=vb;
  }
 
  /* Set Vb=1 for the BTAC */
  par.r[ib].p[0]=1.0;
 
 
  /* 
   *  Print and save parameters
   */
  if(verbose>0 || !parfile[0]) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  if(parfile[0]) {
    /* Save file */
    if(verbose>1) printf("  saving %s\n", parfile);
    FILE *fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&tac); parFree(&par); return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&par); return(12);
    }
    if(verbose>0) printf("parameters saved in %s\n", parfile);
  }
 
  /*
   *  If requested, save Vb-corrected tissue data.
   */
  if(corfile[0]) {
    int ret=TPCERROR_OK;
    FILE *fp; fp=fopen(corfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing '%s'.\n", corfile);
      ret=TPCERROR_FAIL;
    } else {
      ret=tacWrite(&tac, fp, TAC_FORMAT_PMOD, 1, &status);
      fclose(fp);
      if(ret!=TPCERROR_OK) fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    }
    if(ret!=TPCERROR_OK) {
      tacFree(&tac); parFree(&par);
      return(21);
    }
    if(verbose>0) printf("Corrected TACs saved in %s.\n", corfile);
  }
 
  tacFree(&tac); parFree(&par);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/