imgcbvp.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodel.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Estimate Vb and (1-Vb)*Ct based on arterial blood curve (BTAC), assuming",
  "that at BTAC peak tissue concentration is zero, thus Vb=Cpet/Cb",
  "Vb is then reduced to assure that other TACs remain non-negative.",
  " ",
  "Do not use this oversimplified method for quantitative analyses!",
  " ",
  "Usage: @P [Options] imgfile btacfile ctimgfile",
  " ",
  "Data in image and BTAC file must be stored in the same units.",
  " ",
  "Options:",
  " -vbfile=<filename>",
  "     Save image file containing Vb as pixel values.",
  " -LV | -RV",
  "     Correct either arterial Vb or BTAC in injection veins, or",
  "     by default, both.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imgcbv, taccbvp, imgpeak, imgcalc, imgthrs",
  " ",
  "Keywords: image, 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     petfile[FILENAME_MAX], btacfile[FILENAME_MAX], vbfile[FILENAME_MAX], ctfile[FILENAME_MAX];
  int      corrmode=3; // 0=no correction, 1=LV, 2=RV, or 3=both
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=btacfile[0]=ctfile[0]=vbfile[0]=(char)0;
  /* Get 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==NULL) continue;
    if(strncasecmp(cptr, "VB=", 3)==0) {
      strlcpy(vbfile, cptr+3, FILENAME_MAX); continue;
    } else if(strcasecmp(cptr, "LV")==0) {
      corrmode=1; continue;
    } else if(strcasecmp(cptr, "RV")==0) {
      corrmode=2; 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 */
  if(ai<argc) {strlcpy(petfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(btacfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(ctfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
 
  /* Did we get all the information that we need? */
  if(!ctfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n", argv[0]);
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("btacfile := %s\n", btacfile);
    printf("ctfile := %s\n", ctfile);
    if(vbfile[0]) printf("vbfile := %s\n", vbfile);
    printf("corrmode := %d\n", corrmode);
    fflush(stdout);
  }
 
 
  /*
   *  Read PET image and blood TAC
   */
  if(verbose>1) printf("reading data files\n");
  DFT btac;
  IMG img, vbimg;
  imgInit(&img); imgInit(&vbimg); dftInit(&btac);
  double endtime=-1.0;
  char tmp[256];
  int ret=imgReadModelingData(
    petfile, NULL, btacfile, NULL, NULL, &endtime, &ai, &img,
    &btac, NULL, 0, stdout, verbose-1, tmp);
  if(ret!=0) {
    fprintf(stderr, "Error in reading data: %s.\n", tmp);
    if(verbose>1) printf("  ret := %d\n", ret);
    return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  if(verbose>2) printf("endtime := %g min\n", endtime);
  if(verbose>10) dftPrint(&btac);
 
  /* Get the range of pixel values in image */
  float imgmin, imgmax;
  if(imgMinMax(&img, &imgmin, &imgmax)) {
    fprintf(stderr, "Error: invalid image contents.\n");
    imgEmpty(&img); return(3);
  }
  if(verbose>1) {
    printf("image_min := %g\n", imgmin);
    printf("image_max := %g\n", imgmax);
    fflush(stdout);
  }
 
 
  /*
   *  Search BTAC peak
   */
  if(verbose>1) printf("finding BTAC peak\n");
  double bpeak=btac.voi[0].y[0];
  int ipeak=0;
  for(int i=1; i<btac.frameNr; i++)
    if(btac.voi[0].y[i]>bpeak) {bpeak=btac.voi[0].y[i]; ipeak=i;}
  if(verbose>1) {
    printf("  btac_maxv := %g\n", bpeak);
    printf("  btac_maxx := %g\n", btac.x[ipeak]);
  }
  /* If peak is not the first sample, then add previous sample */
  if(ipeak>0) bpeak+=btac.voi[0].y[ipeak-1];
  if(verbose>2) {
    printf("  bpeak := %g\n", bpeak);
  }
 
 
  /*
   *  Allocate memory for Vb image
   */
  ret=imgAllocateWithHeader(&vbimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    dftEmpty(&btac); imgEmpty(&img); return(4);
  }
  vbimg.unit=CUNIT_UNITLESS;
 
 
  /* Check if correction for the venous injection activity is needed */
  double rvtac[img.dimt];
  if(corrmode==3 || corrmode==2) {
 
    /* Find pixels where TTAC peak is somewhat higher than BTAC peak, 
       and calculate mean TAC of those */
    int nrv=0;
    for(int zi=0; zi<img.dimz; zi++) {
      for(int yi=0; yi<img.dimy; yi++) {
        for(int xi=0; xi<img.dimx; xi++) {
          float tpeak=img.m[zi][yi][xi][ipeak];
          if(ipeak>0) tpeak+=img.m[zi][yi][xi][ipeak-1];
          if(tpeak/bpeak>2.0) {
            for(int ti=0; ti<img.dimt; ti++) rvtac[ti]+=img.m[zi][yi][xi][ti];
            nrv++;
          }
        }
      }
    }
    if(verbose>2) printf("nrv := %d\n", nrv);
    if(nrv>0) for(int ti=0; ti<img.dimt; ti++) rvtac[ti]/=(double)nrv;
    else for(int ti=0; ti<img.dimt; ti++) rvtac[ti]=0.0;
    
    if(nrv==0) { // No injection vena evident
      if(verbose>0) printf("RV correction not applied.\n");
      if(corrmode==3) corrmode=1; // from both to only arterial Vb correction
      else if(corrmode==2) corrmode=0; // from venous to none
    }
  }
 
  if(corrmode==0) {
    fprintf(stderr, "Warning: no corrections applied, no files written.\n");
    dftEmpty(&btac); imgEmpty(&img); imgEmpty(&vbimg);
    return(0);
  }
 
 
  /* Correct for the venous injection activity */
  if(corrmode==0 || corrmode==2) {
 
    /* Find peak of injection venous curve */
    double ymaxrv, xmaxrv; int imaxrv;
    ymaxrv=rvtac[0];
    xmaxrv=img.mid[0];
    imaxrv=0;
    for(int ti=1; ti<img.dimt; ti++) if(rvtac[ti]>ymaxrv) {
      ymaxrv=rvtac[ti];
      xmaxrv=img.mid[ti];
      imaxrv=ti;
    }
    if(verbose>1) {
      printf("  rv_maxv := %g\n", ymaxrv);
      printf("  rv_maxx := %g\n", xmaxrv);
    }
 
    /* Correct the image data for RV TAC */
    for(int zi=0; zi<img.dimz; zi++) {
      for(int yi=0; yi<img.dimy; yi++) {
        for(int xi=0; xi<img.dimx; xi++) {
 
          /* Calculate the maximum of Vb that is possible in this pixel from T/B ratio */
          float vb=img.m[zi][yi][xi][imaxrv]/ymaxrv;
          if(vb>=1.0) { // this pixel must be mostly blood
            for(int ti=0; ti<img.dimt; ti++) img.m[zi][yi][xi][ti]=0.0;
            vbimg.m[zi][yi][xi][0]=vb;
            continue;
          }
 
          /* Make initial BV corrected TTAC */
          float ct[img.dimt], ctmin=0.0; int ictmin=0;
          for(int ti=0; ti<img.dimt; ti++) {
            if(img.m[zi][yi][xi][ti]>0.0) {
              ct[ti]=img.m[zi][yi][xi][ti]-vb*rvtac[ti];
              if(ct[ti]<ctmin) {ctmin=ct[ti]; ictmin=ti;}
            } else ct[ti]=0.0;
          }
          /* Make Vb smaller to prevent negative TAC values and correct again */
          if(ctmin<0.0) {
            vb+=ctmin/rvtac[ictmin];
            for(int ti=0; ti<img.dimt; ti++) {
              if(img.m[zi][yi][xi][ti]>0.0) {
                ct[ti]=img.m[zi][yi][xi][ti]-vb*rvtac[ti];
              }
              if(ct[ti]<0.0) ct[ti]=0.0;
            }
          }
          /* Put correct TAC in place */
          for(int ti=0; ti<img.dimt; ti++) img.m[zi][yi][xi][ti]=ct[ti];
          /* Store the Vb */
          vbimg.m[zi][yi][xi][0]=vb;
        }
      }
    }
  }
 
 
  /* Correct for the arterial BTAC */
  if(corrmode==3 || corrmode==1) {
 
    for(int zi=0; zi<img.dimz; zi++) {
      for(int yi=0; yi<img.dimy; yi++) {
        for(int xi=0; xi<img.dimx; xi++) {
 
          /* Calculate the maximum of Vb that is possible in this pixel from T/B ratio */
          float tpeak=img.m[zi][yi][xi][ipeak];
          if(ipeak>0) tpeak+=img.m[zi][yi][xi][ipeak-1];
          float vb=tpeak/bpeak;
          if(vb>=1.0) { // this pixel must be mostly blood
            for(int ti=0; ti<img.dimt; ti++) img.m[zi][yi][xi][ti]=0.0;
            vbimg.m[zi][yi][xi][0]+=vb; //1.0;
            continue;
          }
 
          /* Make initial BV corrected TTAC */
          float ct[img.dimt], ctmin=0.0; int ictmin=0;
          for(int ti=0; ti<img.dimt; ti++) {
            if(img.m[zi][yi][xi][ti]>0.0) {
              ct[ti]=img.m[zi][yi][xi][ti]-vb*btac.voi[0].y[ti];
              if(ct[ti]<ctmin) {ctmin=ct[ti]; ictmin=ti;}
            } else ct[ti]=0.0;
          }
          /* Make Vb smaller to prevent negative TAC values and correct again */
          if(ctmin<0.0) {
            vb+=ctmin/btac.voi[0].y[ictmin];
            for(int ti=0; ti<img.dimt; ti++) {
              if(img.m[zi][yi][xi][ti]>0.0) {
                ct[ti]=img.m[zi][yi][xi][ti]-vb*btac.voi[0].y[ti];
              }
              if(ct[ti]<0.0) ct[ti]=0.0;
            }
          }
          /* Put correct TAC in place */
          for(int ti=0; ti<img.dimt; ti++) img.m[zi][yi][xi][ti]=ct[ti];
          /* Store the Vb */
          vbimg.m[zi][yi][xi][0]+=vb;
        }
      }
    }
  }
 
 
  /* If requested, save (1-Vb)Ct image */
  if(ctfile[0]) {
    if(verbose>1) fprintf(stdout, "writing %s\n", ctfile);
    ret=imgWrite(ctfile, &img);
    if(ret) {
      fprintf(stderr, "Error: %s\n", img.statmsg); 
      if(verbose>1) printf("ret=%d\n", ret);
      dftEmpty(&btac); imgEmpty(&img); imgEmpty(&vbimg);
      return(11);
    }
    if(verbose>0) printf("Written %s\n", ctfile);
  }
 
  /* If requested, save Vb image */
  if(vbfile[0]) {
    if(verbose>1) fprintf(stdout, "writing %s\n", vbfile);
    ret=imgWrite(vbfile, &vbimg);
    if(ret) {
      fprintf(stderr, "Error: %s\n", vbimg.statmsg); 
      if(verbose>1) printf("ret=%d\n", ret);
      dftEmpty(&btac); imgEmpty(&img); imgEmpty(&vbimg);
      return(11);
    }
    if(verbose>0) printf("Written %s\n", vbfile);
  }
 
 
  dftEmpty(&btac); imgEmpty(&img); imgEmpty(&vbimg);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/