fvar4img.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Program for adding Gaussian noise to dynamic PET image",
  "using equations (1, 2):",
  "  SD(t) = TAC(t) * sqrt(Pc/(TAC(t)*exp(-lambda*t)*deltat)",
  "  TAC_noisy(t) = TAC(t) + SD(t)*G(0,1) ,", 
  "where Pc is the proportionality constant that determines the level of noise,",
  "TAC(t) is the mean activity concentration in the image frame,",
  "Deltat is the scan frame length, and G(0,1) is a pseudo random number from",
  "a Gaussian distribution with zero mean and SD of one.",
  " ",
  "Image must be in ECAT (6.3 or 7) format, or, Analyze 7.5 or NIfTI-1 format",
  "with SIF.",
  "Image is assumed to be decay corrected to zero time.",
  " ",
  "Usage: @P [Options] image Pc outputimage ",
  " ",
  "Options:",
  " -seed=<seed for random number generator>",
  "     Computer clock is used by default.",
  " -i=<O-15|N-13|C-11|F-18|Ge-68|Ga-68|Br-76|Rb-82|Cu-62>",
  "     Specifies the isotope.",
  " -minsd=<SD>",
  "     Set a minimum SD to add a minimum level of noise also to the image",
  "     frames with no activity.",
  " -efile=<TAC filename>",
  "     Save image mean TAC and SD that were used in noise simulation.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P -i=C-11 simulated.v 5.0 noisy.v",
  " ",
  "References:",
  "1. Chen K, Huang SC, Yu DC. Phys Med Biol 1991;36:1183-1200.",
  "2. Varga J, Szabo Z. J Cereb Blood Flow Metab 2002;22:240-244.",
  " ",
  "See also: dft2img, flat2img, eframe, imgdecay, fvar4tac, img2tif",
  " ",
  "Keywords: image, simulation, noise",
  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=0, fi, pi, ri, ci;
  unsigned long int seed=0L;
  char      imgfile[FILENAME_MAX], outfile[FILENAME_MAX], efile[FILENAME_MAX];
  char     *cptr, tmp[FILENAME_MAX];
  double    pc=-1.0, minsd=-1.0, halflife=0.0;
  IMG       img;
  DFT       tac;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=outfile[0]=efile[0]=(char)0;
  imgInit(&img); dftInit(&tac);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "S=", 2)==0) {
      seed=atol(cptr+2); if(seed>0L) continue;
    } else if(strncasecmp(cptr, "SEED=", 5)==0) {
      seed=atol(cptr+5); if(seed>0L) continue;
    } else if(strncasecmp(cptr, "I=", 2)==0) {
      cptr=hlCorrectIsotopeCode(cptr+2); if(cptr!=NULL) {
        halflife=60.0*hlFromIsotope(cptr); if(halflife>0.0) continue;}
      fprintf(stderr, "Error: invalid isotope with option -i\n"); return(1);
    } else if(strncasecmp(cptr, "MINSD=", 6)==0) {
      minsd=atof_dpi(cptr+6); if(minsd>0.0) continue;
    } else if(strncasecmp(cptr, "EFILE=", 6)==0) {
      strcpy(efile, cptr+6); if(strlen(efile)>1) 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(!imgfile[0]) {
      strcpy(imgfile, argv[ai]); continue;
    } else if(pc<0.0) {
      pc=atof_dpi(argv[ai]); if(pc>0.0) continue;
    } else if(!outfile[0]) {
      strcpy(outfile, argv[ai]); continue;
    }
    /* we should never get this far */
    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; try %s --help\n",
      argv[0]);
    return(1);
  }
  if(strcasecmp(imgfile, outfile)==0) {
    fprintf(stderr, "Error: must not overwrite input image.\n"); return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("imgfile := %s\n", imgfile);
    printf("outfile := %s\n", outfile);
    printf("proportionality_constant := %g\n", pc);
    if(halflife>0.0) printf("halflife := %g\n", halflife);
    if(minsd>0.0) printf("minsd := %g\n", minsd);
    printf("efile := %s\n", efile);
  }
 
 
  /*
   *  Read the contents of the ECAT file to img data structure
   */
  if(verbose>1) printf("reading image %s\n", imgfile);
  if(imgRead(imgfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    if(verbose>1) imgInfo(&img);
    imgEmpty(&img);
    return(2);
  }
  if(verbose>9) imgInfo(&img);
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  /* Check that image is image */
  if(img.type!=IMG_TYPE_IMAGE) {
    fprintf(stderr, "Error: %s is not an image.\n", imgfile);
    imgEmpty(&img); return(1);
  }
  /* Check that frame times are available */
  if(!imgExistentTimes(&img)) {
    fprintf(stderr, "Error: image does not contain frame times.\n");
    imgEmpty(&img); return(2);
  }
  /* Check that halflife is known */
  if(img.isotopeHalflife<=0.0 && halflife<=0.0) {
    fprintf(stderr, "Error: set isotope with option -i\n");
    imgEmpty(&img); return(3);
  }
  /* If user specified the isotope, then use it */
  if(halflife>0.0) img.isotopeHalflife=halflife;
  /* Calculate decay correction factors */
  ret=imgSetDecayCorrFactors(&img, 1);
  if(ret) {
    fprintf(stderr, "Error %d in setting decay correction factors.\n", ret);
    imgEmpty(&img); return(3);
  }
 
 
  /* Check if outputfile exists; rename it as backup file, if necessary */
  ret=backupExistingFile(outfile, NULL, tmp);
  if(ret!=0) {
    fprintf(stderr, "Error: %s\n", tmp);
    imgEmpty(&img); return(11);
  }
 
 
  /*
   *  Calculate an average TAC of all image pixels
   */
  if(verbose>1) printf("calculating average TAC of all image pixels\n");
  /* Allocate memory for TAC */
  ret=dftSetmem(&tac, img.dimt, 2);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    imgEmpty(&img); return(5);
  }
  /* Set TAC information */
  strcpy(tac.isotope, imgIsotope(&img));
  tac._type=1; tac.frameNr=img.dimt; tac.voiNr=2;
  strcpy(tac.unit, imgUnit(img.unit));
  tac.timetype=DFT_TIME_STARTEND; tac.timeunit=TUNIT_SEC;
  strcpy(tac.voi[0].voiname,"Mean"); strcpy(tac.voi[0].name,tac.voi[0].voiname);
  strcpy(tac.voi[1].voiname,"SD");   strcpy(tac.voi[1].name,tac.voi[1].voiname);
  for(fi=0; fi<img.dimt; fi++) {
    tac.x1[fi]=img.start[fi]; tac.x2[fi]=img.end[fi];
    tac.x[fi]=img.mid[fi];
  }
  tac.decayCorrected=DFT_DECAY_CORRECTED;
  /* Calculate average of pixel values */
  ret=imgAverageTAC(&img, img.sd);
  if(ret) {
    fprintf(stderr, "Error: cannot calculate average TAC (%d).\n", ret);
    imgEmpty(&img); dftEmpty(&tac);
    return(6);
  }
  for(fi=0; fi<img.dimt; fi++) tac.voi[0].y[fi]=img.sd[fi];
  if(verbose>3) dftPrint(&tac);
 
 
  /*
   *  Compute SD values for each frame
   */
  if(verbose>1) printf("calculating SD for each frame\n");
  ret=noiseSD4SimulationFromDFT(&tac, 0, pc, tac.voi[1].y, tmp, verbose-3);
  if(ret!=0) {
    fprintf(stderr,
            "Error: cannot calculate SD for noise simulation (%d).\n", ret);
    if(verbose>1) printf("     %s\n", tmp);
    imgEmpty(&img); dftEmpty(&tac);
    return(7);
  }
  /* Set min SD, if specified by user */
  if(minsd>0.0)
   for(fi=0; fi<img.dimt; fi++)
     if(minsd>tac.voi[1].y[fi]) tac.voi[1].y[fi]=minsd;
  if(verbose>1) {
    printf("Frame\t\t\tMean\tSD\n");
    for(fi=0; fi<img.dimt; fi++)
      printf("%d\t%g\t%g\t%g\t%g\n",
             1+fi, tac.x1[fi], tac.x2[fi], tac.voi[0].y[fi], tac.voi[1].y[fi]);
  }
  if(verbose>8) dftPrint(&tac);
  if(efile[0]) {  /* Save image mean TAC and SD, if required */
    if(verbose>1) printf("writing %s\n", efile);
    DFT_NR_OF_DECIMALS=6;
    if(dftWrite(&tac, efile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", efile, dfterrmsg);
      imgEmpty(&img); dftEmpty(&tac); return(11);
    }
    if(verbose>0) printf("Mean TAC and noise SD written in %s\n", efile);
  }
 
 
  /*
   *  Add noise to the image pixels
   */
  if(verbose>1) printf("adding noise to image\n");
  if(seed>0) GAUSSDEV_SEED=seed;
  else GAUSSDEV_SEED=drandSeed(0);
  for(fi=0; fi<img.dimt; fi++) {
    if(verbose>0 && img.dimt>1) {fprintf(stdout, "."); fflush(stdout);}
    for(pi=0; pi<img.dimz; pi++)
      for(ri=0; ri<img.dimy; ri++)
        for(ci=0; ci<img.dimx; ci++)
          img.m[pi][ri][ci][fi] += tac.voi[1].y[fi]*gaussdev();
  }
  if(verbose>0 && img.dimt>1) {fprintf(stdout, "\n"); fflush(stdout);}
 
 
  /*
   *  Write the noisy image
   */
  if(verbose>1) printf("writing noisy image %s\n", outfile);
  if(imgWrite(outfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    imgEmpty(&img); dftEmpty(&tac);
    return(21);
  }
 
 
  /*
   *  Free up memory
   */
  imgEmpty(&img);
  dftEmpty(&tac);
  if(verbose>0) printf("done.\n");
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/