imgcbv.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 "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Subtracts the contribution of vascular radioactivity from dynamic PET",
  "image. Vascular volume fraction Vb can be given as a value that is common",
  "to all image voxels, or as separate value for each image voxel in a",
  "(preferably smoothed) Vb image, usually calculated from a [O-15]CO study.",
  "Vb must be given as a volume fraction, not percentage, whether it is given",
  "directly or inside the Vb image.",
  " ",
  "Usage: @P [Options] imgfile btacfile Vb outputfile",
  " ",
  "Options:",
  " -noneg",
  "     Negative voxel 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: Vascular volume fraction in all tissue is assumed to be 0.04:",
  "     imgcbv ua2918dy1.img ua2918ab.kbq 0.04 ua2918cbv.img",
  "Example 2: Vascular volume fraction is given in specific image file:",
  "     imgcbv ua2918dy1.v ua2918ab.kbq ua2918vb.v ua2918cbv.v",
  " ",
  "See also: taccbv, imgcalc, p2blood, imgunit, tacunit",
  " ",
  "Keywords: image, 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 ecat7CopyHeadersNoQuant(
  const char *file1,
  const char *file2,
  int verbose
) {
  int ii, mi, ret;
  char tmp[1024];
  ECAT_HEADERS ehdr;
 
  if(verbose>1) printf("ecatCopyHeaders(%s, %s)\n", file1, file2);  
  if(file1==NULL || file2==NULL) return(1);
  
  
  /* Read headers */
  ehdrInitiate(&ehdr);
  ret=ecat7ReadHeaders(file1, &ehdr, verbose-1); if(ret!=0) return(10);
  /* Remove main header information that should NOT be copied */
  strcpy(tmp, "user_process_code");
  ii=iftGet(&ehdr.mh, tmp, 0); iftDeleteItem(&ehdr.mh, ii, 0);
  strcpy(tmp, "ecat_calibration_factor");
  ii=iftGet(&ehdr.mh, tmp, 0); iftDeleteItem(&ehdr.mh, ii, 0);
 
  /* Remove subheader information that should NOT be copied */
  for(mi=0; mi<ehdr.nr; mi++) {
    strcpy(tmp, "scale_factor");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "image_min");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "image_max");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "scan_min");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "scan_max");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
  }
  
  /* Write headers */
  ret=ecat7WriteHeaders(file2, &ehdr, verbose-1);
  ehdrEmpty(&ehdr);
  if(ret!=0) return(10+ret);
 
  return 0;
}

/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  int     ret;
  int     leave_negat=1, 
          pet_volume=1, 
          add_vb=0;
  char    infile[FILENAME_MAX], blfile[FILENAME_MAX], outfile[FILENAME_MAX],
         *cptr, tmp[FILENAME_MAX+128],
          vbfile[FILENAME_MAX];
  double  vb=0.0;
  double  endtime;
  float   f;
  DFT     blood;
  IMG     img, vbimg;
  int     isvalue=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  infile[0]=outfile[0]=blfile[0]=vbfile[0]=(char)0;
  imgInit(&img); imgInit(&vbimg); dftInit(&blood);
  /* 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(isvalue==0 && !vbfile[0]) {
      double v;
      ret=atof_with_check(argv[ai], &v);
      if(ret==0) {
        isvalue=1; vb=v; 
        if(vb>=1.0) {
          fprintf(stderr, "Warning: Vb fraction is set to %.3f\n", vb);
        }
      } else { 
        strlcpy(vbfile, argv[ai], FILENAME_MAX);
      }
      continue;
    } 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");
  }
 
  /* Check that output file is not the same as the input file */
  if(strcasecmp(infile, outfile)==0) {
    fprintf(stderr, "Error: original file must not be overwritten.\n");
    return(1);
  }
 
 
  /*
   *  Read PET image and blood TAC
   */
  if(verbose>1) printf("reading data files\n");
  endtime=-1.0;
  ret=imgReadModelingData(
    infile, NULL, blfile, NULL, NULL, &endtime, &ai, &img,
    NULL, &blood, 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(verbose>2) printf("endtime := %g min\n", endtime);
  /* Convert blood sample times to min */
  dftTimeunitConversion(&blood, TUNIT_MIN);
  if(verbose>3) dftPrint(&blood);
 
 
  /*
   *  If Vb values are given in an image file, then read that
   */
  if(vbfile[0]) {
    if(verbose>1) printf("reading %s\n", vbfile);
    ret=imgRead(vbfile, &vbimg);
    if(ret) {
      fprintf(stderr, "Error: %s\n", vbimg.statmsg);
      imgEmpty(&img); dftEmpty(&blood); return(8);
    }
    /* Check if data is raw or image */
    if(vbimg.type!=IMG_TYPE_IMAGE) {
      fprintf(stderr, "Error: %s is not an image.\n", vbfile);
      imgEmpty(&img); dftEmpty(&blood); imgEmpty(&vbimg); return(8);
    }
    if(vbimg.dimt>1) {
      fprintf(stderr, "Error: Vb image contains > 1 frame.\n");
      imgEmpty(&img); dftEmpty(&blood); imgEmpty(&vbimg); return(8);
    }
    /* check dimensions */
    ret=0;
    if(img.dimz!=vbimg.dimz) ret=1;
    if(img.dimx!=vbimg.dimx) ret+=2;
    if(img.dimy!=vbimg.dimy) ret+=4;
    if(ret) {
      fprintf(stderr, "Error: Vb image has different dimensions.\n");
      imgEmpty(&img); dftEmpty(&blood); imgEmpty(&vbimg); return(8);
    }
    /* check that Vb is a fraction */
    ret=imgMax(&vbimg, &f);
    if(verbose>2) printf("Maximum value in Vb image: %g\n", f);
    if(f>1.0) {
      ret=imgArithmConst(&vbimg, 0.01, '*', 1.0, verbose-4);
      fprintf(stderr, "Warning: Vb values were converted to fractions.\n");
    }
  }
 
  /*
   *  Make subtraction
   */
  if(verbose>1) printf("subtracting\n");
  {
  int fi, zi, xi, yi;
  for(zi=0; zi<img.dimz; zi++)
    for(yi=0; yi<img.dimy; yi++)
      for(xi=0; xi<img.dimx; xi++) {
        if(vbfile[0]) vb=vbimg.m[zi][yi][xi][0];
        for(fi=0; fi<img.dimt; fi++) {
          img.m[zi][yi][xi][fi]-=vb*blood.voi[0].y[fi];
          if(pet_volume==0) img.m[zi][yi][xi][fi]/=(1.0-vb);
          if(leave_negat==0 && img.m[zi][yi][xi][fi]<0.0)
            img.m[zi][yi][xi][fi]=0.0;
        }
      }
  }
  dftEmpty(&blood); imgEmpty(&vbimg);
 
 
  /*
   *  Save new image
   */
  if(verbose) fprintf(stdout, "writing image %s\n", outfile);
  ret=imgWrite(outfile, &img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    imgEmpty(&img); dftEmpty(&blood); imgEmpty(&vbimg); return(13);
  }
  if(verbose>1) fprintf(stdout, "Vb corrected image %s saved.\n", outfile);
 
  if(img._fileFormat==IMG_E7 || img._fileFormat==IMG_E7_2D) {
    imgEmpty(&img); 
    if(verbose>1) printf("copying header information\n");
    ret=ecat7CopyHeadersNoQuant(infile, outfile, verbose-1);
    if(verbose>0 && ret!=0) 
      printf("copying headers not successful (%d).\n", ret);
  }
 
  imgEmpty(&img); 
  if(verbose>0) printf("done.\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/