dftcbv.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
#include "libtpcmodel.h"
#include "libtpcmodext.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 file, calculated",
  "from a [O-15]CO study.",
  " ",
  "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.",
  " -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.0 are assumed to be percentages.",
  "Blood TAC can be given in a separate BTAC file, or as a region id inside",
  "TTAC file.",
  " ",
  "See also: imgcbv, p2blood, interpol, taccalc, tacformat",
  " ",
  "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;
  int     ri, fi, ret;
  int     leave_negat=1, pet_volume=1, add_vb=0;
  char    infile[FILENAME_MAX], blfile[FILENAME_MAX], outfile[FILENAME_MAX],
         *cptr, vbfile[FILENAME_MAX], tmp[256];
  DFT     dft, blood, tdft;
  double  vb=-1.0, t1, t2;
  
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  infile[0]=outfile[0]=blfile[0]=vbfile[0]=(char)0;
  dftInit(&dft); dftInit(&blood); dftInit(&tdft);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "NONEGATIVES", 1)==0) {
      leave_negat=0; continue;
    } else if(strncasecmp(cptr, "PV", 1)==0) {
      pet_volume=1; continue;
    } else if(strncasecmp(cptr, "TV", 1)==0) {
      pet_volume=0; continue;
    } else if(strncasecmp(cptr, "SIMULATE", 3)==0) {
      add_vb=1; continue;
    }
    fprintf(stderr, "Error: unknown option '%s'.\n", argv[ai]);
    return(1);
  } else break;
  
  /* 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 */
  for(; ai<argc; ai++) {
    if(!infile[0]) {strlcpy(infile, argv[ai], FILENAME_MAX); continue;}
    else if(!blfile[0]) {strlcpy(blfile, argv[ai], FILENAME_MAX); continue;}
    else if(!vbfile[0] && vb<0.0) {
      if(access(argv[ai], 0)==0) {
        strlcpy(vbfile, argv[ai], FILENAME_MAX); continue;
      }
      ret=atof_with_check(argv[ai], &vb);
      if(ret==0) {if(vb>=1.0) vb/=100.; if(vb>=0.0) continue;}
      fprintf(stderr, "Error: invalid Vb: '%s'.\n", argv[ai]);
      return(1);
    }
    else if(!outfile[0]) {strlcpy(outfile, argv[ai], FILENAME_MAX); continue;}
    fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!outfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("infile := %s\n", infile);
    printf("blfile := %s\n", blfile);
    printf("vbfile := %s\n", vbfile);
    printf("pet_volume := %d\n", pet_volume);
    printf("vb := %g\n", vb);
    printf("leave_negat := %d\n", leave_negat);
    printf("add_vb := %d\n", add_vb);
  }
 
  if(verbose>1) {
    printf("\napplying formula:\n");
    if(add_vb==0) {
      if(pet_volume==0) printf("Ct=(Cpet-Vb*Cb)/(1-Vb)\n");
      else              printf("Ct=Cpet-Vb*Cb\n");
    } else {
      if(pet_volume==0) printf("Cpet=(1-Vb)*Ct+Vb*Cb\n");
      else              printf("Cpet=Ct+Vb*Cb\n");
    }
    printf("\n");
  }
 
  /*
   *  Read regional TACs
   */
  if(verbose>1) printf("reading %s\n", infile);
  if(dftRead(infile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", infile, dfterrmsg);
    return(2);
  }
 
  /*
   *  Read blood TAC
   */
  if(verbose>1) printf("reading %s\n", blfile);
  ret=dftReadinput(&blood, &dft, blfile, &ri, &t1, &t2, 0, tmp, verbose-3);
  if(ret!=0) {
    fprintf(stderr, "Error in reading '%s': %s\n", blfile, tmp);
    dftEmpty(&dft); return(3);
  }
  if(verbose>2) printf("blood_filetype := %d\n", ri);
  /* Check blood time range */
  if((dft.frameNr>1 && t1>dft.x[1]) || (dft.x[dft.frameNr-1]>1.2*t2)) {
    fprintf(stderr, "Warning: blood TAC is clearly shorter than tissue TAC.\n");
  }
 
  /*
   *  If Vb values are given in a file, then read that
   */
  if(vbfile[0]) {
    if(verbose>1) printf("reading %s\n", vbfile);
    if(dftRead(vbfile, &tdft)) {
      fprintf(stderr, "Error in reading '%s': %s\n", vbfile, dfterrmsg);
      dftEmpty(&dft); dftEmpty(&blood); return(4);
    }
    if(tdft.frameNr>1) {
      fprintf(stderr, "Error: Vb file contains more than one sample time.\n");
      dftEmpty(&dft); dftEmpty(&blood); dftEmpty(&tdft); return(4);
    }
    /* check regions */
    ret=0;
    if(dft.voiNr!=tdft.voiNr) 
      ret=1;
    else {
      for(ri=0; ri<dft.voiNr; ri++) {
        if(strcmp(dft.voi[ri].name, tdft.voi[ri].name)!=0) ret=2;
        if(strcmp(dft.voi[ri].voiname, tdft.voi[ri].voiname)!=0) ret=3;
        if(strcmp(dft.voi[ri].hemisphere, tdft.voi[ri].hemisphere)!=0) ret=4;
        if(strcmp(dft.voi[ri].place, tdft.voi[ri].place)!=0) ret=5;
        if(ret) break;
      }
    }
    if(ret) {
      fprintf(stderr, "Error: Vb file does not contain the same regions.\n");
      dftEmpty(&dft); dftEmpty(&blood); dftEmpty(&tdft); return(4);
    }
    /* check that Vb is fraction */
    for(ri=ret=0; ri<tdft.voiNr; ri++) if(tdft.voi[ri].y[0]>=1.0) ret=1;
    if(ret) {
      for(ri=0; ri<tdft.voiNr; ri++) tdft.voi[ri].y[0]/=100.;
      for(ri=ret=0; ri<tdft.voiNr; ri++) if(tdft.voi[ri].y[0]>=1.0) ret=1;
      if(ret) {
        fprintf(stderr, "Error: invalid contents in Vb file.\n");
        dftEmpty(&dft); dftEmpty(&blood); dftEmpty(&tdft); return(4);
      }
      fprintf(stderr, "Warning: Vb values were converted to fractions.\n");
    }
  }
 
 
  /*
   *  Make subtraction, or add Vb contribution
   */
  if(add_vb==0) {
    if(verbose>1) printf("subtracting blood\n");
    for(ri=0; ri<dft.voiNr; ri++) {
      if(vbfile[0]) {
        vb=tdft.voi[ri].y[0];
        if(verbose>2) printf("Vb=%g for %s\n", vb, dft.voi[ri].name);
      }
      for(fi=0; fi<dft.frameNr; fi++) {
        dft.voi[ri].y2[fi]=dft.voi[ri].y[fi];
        if(!isnan(dft.voi[ri].y2[fi])) {
          dft.voi[ri].y2[fi]-=vb*blood.voi[0].y[fi];
          if(pet_volume==0) dft.voi[ri].y2[fi]/=(1.0-vb);
          if(leave_negat==0 && dft.voi[ri].y2[fi]<0.0) dft.voi[ri].y2[fi]=0.0;
        }
      }
    } /* next region */
  } else {
    if(verbose>1) printf("adding blood\n");
    for(ri=0; ri<dft.voiNr; ri++) {
      if(vbfile[0]) {
        vb=tdft.voi[ri].y[0];
        if(verbose>2) printf("Vb=%g for %s\n", vb, dft.voi[ri].name);
      }
      for(fi=0; fi<dft.frameNr; fi++) {
        dft.voi[ri].y2[fi]=dft.voi[ri].y[fi];
        if(!isnan(dft.voi[ri].y2[fi])) {
          if(pet_volume==0) dft.voi[ri].y2[fi]*=(1.0-vb);
          dft.voi[ri].y2[fi]+=vb*blood.voi[0].y[fi];
          if(leave_negat==0 && dft.voi[ri].y2[fi]<0.0) dft.voi[ri].y2[fi]=0.0;
        }
      }
    } /* next region */
  }
  if(verbose>20) dftPrint(&dft);
 
 
  /*
   *  Save data
   */
  if(verbose>1) printf("writing %s\n", outfile);
  for(ri=0; ri<dft.voiNr; ri++) for(fi=0; fi<dft.frameNr; fi++)
    dft.voi[ri].y[fi]=dft.voi[ri].y2[fi];
  ret=dftWrite(&dft, outfile);
  if(ret) {
    fprintf(stderr, "Error (%d) in writing '%s': %s\n", ret, outfile, dfterrmsg);
    dftEmpty(&dft); dftEmpty(&blood); dftEmpty(&tdft);
    return(11);
  }
  if(verbose>0) printf("  %s written.\n", outfile);
 
  dftEmpty(&dft); dftEmpty(&blood); dftEmpty(&tdft);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/