heartcor.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcidi.h"
/*****************************************************************************/
#ifndef M_2SQRT2LN2
#define M_2SQRT2LN2 2.354820045
#endif
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Corrects regional myocardium and blood time-activity curves (TACs)",
  "measured with PET for recovery and spillover fractions, and contribution",
  "of blood activity in vascular space (10%) (Henze et al., 1983).",
  "This method works properly, if regions-of-interest (ROIs) are drawn",
  "in the middle of cavity and muscle, excluding the border areas.",
  "Method is valid when cavity radius >> wall thickness and FWHM/2.36.",
  " ",
  "Usage: @P tacfile cav_id myo_id diameter thickness FWHM outputfile",
  " ",
  "Tacfile must contain at least the TACs of myocardial cavity and entire",
  "myocardial muscle (or representative region of it). The names or numbers",
  "of these TACs must be given as next command-line parameters.",
  "Corrections based on these TACs are applied to any additional TACs in tacfile.",
  " ",
  "The diameter of cavity, myocardial wall thickness, and PET image resolution",
  "must be given in millimeters.",
  " ",
  "Options:",
  " -Vb=<Vascular volume fraction>",
  "     Set vascular volume fraction (0-0.9) in myocardial tissue;",
  "     by default Vb=0.10",
  " -sim[ulate]",
  "     Recovery and spillover errors are simulated and added to initially",
  "     non-affected data, originating from compartmental model simulations;",
  "     enter blood and muscle tissue TACs instead of cavity and myocardial",
  "     region TACs, respectively.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1: correction of regional PET data",
  "     @P s6368dy.tac lv myo 52 10 10 s6368dy_cor.tac",
  " ",
  "Example 2: simulation of recovery and spillover effects, assuming Vb=8%%",
  "     @P -simulate original.tac 1 2 40 10 8 simulated.tac",
  "Example 3: calculate FMM, FMB, FBM, and FBB, without doing anything else",
  "     @P none 1 2 40 10 8 none",
  " ",
  "Reference:",
  "1. Henze E, Huang S-C, Ratib O, Hoffman E, Phelps ME, Schelbert HR.",
  "   Measurements of regional tissue and blood-pool radiotracer",
  "   concentrations from serial tomographic images of the heart.",
  "   J Nucl Med. 1983;24:987-996.",
  " ",
  "See also: fitmbf, patlak, logan, simimyoc, taccalc, taccbv, spillinp",
  " ",
  "Keywords: TAC, myocardium, cavity, spill-over, recovery, 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     ret, ri, fi;
  double  R=-1.0;        // radius of cavity and circular ROI (mm)
  double  d=-1.0;        // thickness of myocardium (mm)
  double  s=-1.0;        // spatial resolution (mm), s=FWHM/2.354820045
  int     simulation=0;  // 1=add errors, 0=correct errors
  char   *cptr, dfile[FILENAME_MAX], rfile[FILENAME_MAX];
  char   *cav_roi=NULL, *myo_roi=NULL;
  int     lv=-1;         // number of LV TAC in DFT struct
  int     myo=-1;        // number of whole myocardium TAC in DFT struct
  double  Vb=0.1;        // vascular volume fraction in tissue
  DFT     dft;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dfile[0]=rfile[0]=(char)0;
  dftInit(&dft);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(!*cptr) continue;
    if(strncasecmp(cptr, "SIMULATE", 3)==0) {
      simulation=1; continue;
    } else if(strncasecmp(cptr, "VB=", 3)==0) {
      ret=atof_with_check(cptr+3, &Vb);
      if(ret==0 && Vb>=0.0 && Vb<1.0) continue;
    }
    fprintf(stderr, "Error: invalid 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(!dfile[0]) {
      strlcpy(dfile, argv[ai], FILENAME_MAX); continue;
    } else if(cav_roi==NULL) {
      cav_roi=argv[ai]; continue;
    } else if(myo_roi==NULL) {
      myo_roi=argv[ai]; continue;
    } else if(R<0.0) {
      R=0.5*atof_dpi(argv[ai]); if(R>0.0) continue;
      fprintf(stderr, "Error: invalid diameter '%s'\n", argv[ai]); return(1);
    } else if(d<0.0) {
      d=atof_dpi(argv[ai]); if(d>0.0) continue;
      fprintf(stderr, "Error: invalid thickness '%s'\n", argv[ai]); return(1);
    } else if(s<0.0) {
      s=atof_dpi(argv[ai])/M_2SQRT2LN2; if(s>0.0) continue;
      fprintf(stderr, "Error: invalid resolution '%s'\n", argv[ai]); return(1);
    } else if(!rfile[0]) {
      strlcpy(rfile, argv[ai], FILENAME_MAX); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!rfile[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
  if(strcasecmp(dfile, "NONE")==0 || strcasecmp(rfile, "NONE")==0) {
    strcpy(dfile, ""); strcpy(rfile, "");
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dfile := %s\n", dfile);
    printf("rfile := %s\n", rfile);
    printf("cav_roi := %s\n", cav_roi);
    printf("myo_roi := %s\n", myo_roi);
    printf("R := %g\n", R);
    printf("d := %g\n", d);
    printf("s := %.6f\n", s);
    printf("simulation := %d\n", simulation);
    printf("Vb := %g\n", Vb);
  }
 
 
  /* Calculate correction coefficients */
  /* Do this before trying to read datafile, to allow user just to calculate
     these values for some other use */
  double FMM, FMB, FBM, FBB;
  ret=heartRecoverySpilloverCorrectionFactors(R, d, s, Vb,
                                              &FMM, &FMB, &FBM, &FBB);
  if(strlen(dfile)==0 || verbose>0) {
    printf("FMB := %.6f\n", FMB);
    printf("FMM := %.6f\n", FMM);
    printf("FBB := %.6f\n", FBB);
    printf("FBM := %.6f\n", FBM);
  }
  if(strlen(dfile)==0) return(0); // we are done then
 
 
  /* Read datafile */
  if(verbose>1) printf("reading %s\n", dfile);
  if(dftRead(dfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dfile, dfterrmsg);
    return(2);
  }
 
  /* Select the LV and whole myocardium ROIs */
  ret=dftSelectRegions(&dft, cav_roi, 1);
  if(ret<1) {
    fprintf(stderr, "Error: ROI '%s' not found in %s\n", cav_roi, dfile);
    dftEmpty(&dft); return(3);
  } else {
    int len=256, last_len=256;
    for(ri=0; ri<dft.voiNr; ri++)
      if(dft.voi[ri].sw==1) {lv=ri; len=last_len=strlen(dft.voi[ri].name); break;}
    for(;ri<dft.voiNr; ri++)
      if(dft.voi[ri].sw==1 && (len=strlen(dft.voi[ri].name))<last_len) 
        {lv=ri; last_len=len;}
  }
  ret=dftSelectRegions(&dft, myo_roi, 1);
  if(ret<1) {
    fprintf(stderr, "Error: ROI '%s' not found in %s\n", myo_roi, dfile);
    dftEmpty(&dft); return(3);
  } else {
    int len=256, last_len=256;
    for(ri=0; ri<dft.voiNr; ri++)
      if(dft.voi[ri].sw==1) {myo=ri; len=last_len=strlen(dft.voi[ri].name); break;}
    for(;ri<dft.voiNr; ri++)
      if(dft.voi[ri].sw==1 && (len=strlen(dft.voi[ri].name))<last_len)
        {myo=ri; last_len=len;}
  }
  if(verbose>0) {
    printf("cavity := %s\n", dft.voi[lv].name);
    printf("myocardium := %s\n", dft.voi[myo].name);
  }
  if(lv==myo) {
    fprintf(stderr, "Error: same ROI selected for cavity and myocardium.\n");
    dftEmpty(&dft); return(3);
  }
 
  /* Correct the TACs */
  if(verbose>1) printf("correcting TACs\n");
  double CMO, CBO; // imaged concentrations of myocardium and blood pool
  double CM, CB;   // true concentrations in myocardium and in blood
  double CMs, CMOs;// concentrations in smaller myocardial regions
  double K=1.0/(FBB*FMM-FMB*FBM);
  if(verbose>2 && simulation==0) printf("K := %.5f\n", K);
  for(fi=0; fi<dft.frameNr; fi++) {
    if(simulation==1) {
      /* add errors to blood and whole myocardium TACs */
      CM=dft.voi[myo].y[fi];
      CB=dft.voi[lv].y[fi];
      CMO=FMM*CM + FBM*CB; dft.voi[myo].y[fi]=CMO;
      CBO=FMB*CM + FBB*CB; dft.voi[lv].y[fi]=CBO;
      /* add errors to other myocardial curves */
      for(ri=0; ri<dft.voiNr; ri++) if(ri!=lv && ri!=myo) {
        CMs=dft.voi[ri].y[fi];
        CMOs=FMM*CMs + FBM*CB; dft.voi[ri].y[fi]=CMOs;
      }
    } else {
      /* correct the LV and whole myocardium TACs */
      CMO=dft.voi[myo].y[fi];
      CBO=dft.voi[lv].y[fi];
      CM=K*(FBB*CMO-FBM*CBO); dft.voi[myo].y[fi]=CM;
      CB=K*(FMM*CBO-FMB*CMO); dft.voi[lv].y[fi]=CB;
      /* correct other myocardial regions */
      for(ri=0; ri<dft.voiNr; ri++) if(ri!=lv && ri!=myo) {
        CMOs=dft.voi[ri].y[fi];
        //CMs=(CMOs - FBM*((CBO-FMB*CM)/FBB)) / FMM;
        CMs=(CMOs - FBM*CB) / FMM;
        dft.voi[ri].y[fi]=CMs;
      }
    }
  } // next frame time
 
  /* Save data */
  if(dftWrite(&dft, rfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", rfile, dfterrmsg);
    dftEmpty(&dft); return(11);
  }
  if(verbose>0) printf("%s written.\n", rfile);
 
  /* Free memory */
  dftEmpty(&dft);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/