simiart.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
//#include <unistd.h>
#include <string.h>
#include <math.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodext.h"
#include "libtpcidi.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulate dynamic PET image plane with background activity surrounding",
  "a circular object that is assumed to be extending across image planes",
  "such as abdominal aorta. Spill-over and spill-in effects can be simulated",
  "by Gaussian 2D smoothing, or by using method based on Germano et al ",
  "JNM 1992; 33: 613-620, which is only valid for vessels with very large",
  "diameter. To simulate the circular vessel correctly in 2D image matrix use",
  "the equations in Brix et al Nuklearmedizin 2002;41:184-190 instead;",
  "however, that would require numerical solution to double integrals",
  " ",
  "Usage: @P [Options] bloodfile bkgfile image",
  " ",
  "Options:",
  " -dim=<Image x and y dimension in pixels>",
  "     Set image dimension; by default 128.",
  " -pxlsize=<Voxel size (mm)>",
  "     Set image voxel size; by default 1 mm.",
  " -diameter=<Vessel diameter (mm)>",
  "     Set vessel diameter; by default 25 mm; vessel border is simulated",
  "     by simple cutoff, thus diameter should be several pixels.",
  " -fwhm=<FWHM (mm)>",
  "     Set image resolution; by default 8 mm.",
  " -xpos=<vessel position (mm)>",
  "     Set vessel distance from centre of image; by default 0 mm;",
  "     negative distance moves vessel to the left, positive to the right.",
  " -ypos=<vessel position (mm)>",
  "     Set vessel distance from centre of image; by default 0 mm",
  "     negative distance moves vessel upwards, positive downwards.",
  " -gaussian | -bar | -nopve",
  "     Select the spill-over and spill-in simulation method; by default",
  "     no PVE simulation. ",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1:",
  "  @P blood.tac background.tac vessel.v",
  "Example 2:",
  "  @P -fwhm=4.3 -diameter=7 -pxlsize=2.3 blood.tac background.tac vessel.v",
  " ",
  "See also: fvar4tac, imgfiltg, imgadd, imgcalc, simcirc, simimyoc, flat2img",
  " ",
  "Keywords: simulation, image, software testing, vessel, input",
  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      fi, ret;
  char     blofile[FILENAME_MAX], bkgfile[FILENAME_MAX], imgfile[FILENAME_MAX];
  char    *cptr, tmp[256];
  int      pve_method=-1; // -1=No PVE, 0=Germano, 1=Gaussian
  int      dim=128;
  double   pxlsize=1.0, diameter=25.0, fwhm=8.0, xpos=0.0, ypos=0.0;
  int      scanner_type=SCANNER_HRPLUS;
  float    zoom=1.41421356;
  DFT      tac;
  IMG      img;
  double   cx, cy;
 
  
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  blofile[0]=bkgfile[0]=imgfile[0]=(char)0;
  /* Get 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=='-') cptr++; if(cptr==NULL) continue;
    if(strncasecmp(cptr, "DIM=", 4)==0) {
      if(atoi_with_check(cptr+4, &dim)==0 && dim>4) continue;
    } else if(strncasecmp(cptr, "PXLSIZE=", 8)==0) {
      if(atof_with_check(cptr+8, &pxlsize)==0 && pxlsize>0.0) continue;
    } else if(strncasecmp(cptr, "DIAMETER=", 9)==0) {
      if(atof_with_check(cptr+9, &diameter)==0 && diameter>0.0) continue;
    } else if(strncasecmp(cptr, "FWHM=", 5)==0) {
      if(atof_with_check(cptr+5, &fwhm)==0 && fwhm>0.0) continue;
    } else if(strncasecmp(cptr, "XPOS=", 5)==0) {
      if(atof_with_check(cptr+5, &xpos)==0) continue;
    } else if(strncasecmp(cptr, "YPOS=", 5)==0) {
      if(atof_with_check(cptr+5, &ypos)==0) continue;
    } else if(strncasecmp(cptr, "NOPVE", 2)==0) {
      if(pve_method>=0) {
        fprintf(stderr, "Error: invalid combination of simulation options.\n");
        return(1);
      }
      pve_method=-1; continue;
    } else if(strcasecmp(cptr, "BAR")==0) {
      if(pve_method>=0) {
        fprintf(stderr, "Error: invalid combination of simulation options.\n");
        return(1);
      }
      pve_method=0; continue;
    } else if(strcasecmp(cptr, "GAUSSIAN")==0) {
      if(pve_method>=0) {
        fprintf(stderr, "Error: invalid combination of simulation options.\n");
        return(1);
      }
      pve_method=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 */
  if(ai<argc) {strlcpy(blofile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(bkgfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(imgfile, 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(!imgfile[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("blofile := %s\n", blofile);
    printf("bkgfile := %s\n", bkgfile);
    printf("imgfile := %s\n", imgfile);
    printf("dim := %d\n", dim);
    printf("pxlsize := %g\n", pxlsize);
    printf("diameter := %g\n", diameter);
    printf("fwhm := %g\n", fwhm);
    printf("xpos := %g\n", xpos);
    printf("ypos := %g\n", ypos);
    printf("pve_method := %d\n", pve_method);
  }
 
 
  /* 
   *  Read blood and background TACs
   */
  dftInit(&tac);
  if(verbose>1) printf("reading blood TAC\n");
  if(dftRead(blofile, &tac)) {
    fprintf(stderr, "Error in reading '%s': %s\n", blofile, dfterrmsg);
    return(2);
  }
  if(tac.voiNr>1) {
    fprintf(stderr, "Error: blood data contains more than one curve.\n");
    dftEmpty(&tac); return(2);
  }
  if(verbose>1) printf("reading background TAC\n");
  ret=dftReadReference(&tac, bkgfile, NULL, NULL, tmp, verbose-4);
  if(ret<1) {
    fprintf(stderr, "Error in reading background TAC: %s\n", tmp);
    dftEmpty(&tac); return(2);
  }
  if(ret>1) {
    fprintf(stderr, "Error: background data contains more than one curve.\n");
    dftEmpty(&tac); return(2);
  }
 
 
  /*
   *  Allocate memory for the image
   */
  if(verbose>1) printf("allocating memory for image data\n");
  imgInit(&img);
  ret=imgAllocate(&img, 1, dim, dim, tac.frameNr);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot allocate memory for image data.\n");
    dftEmpty(&tac); return(4);
  }
  /* Set header information */
  if(verbose>2) printf("setting image headers\n");
  img.type=IMG_TYPE_IMAGE;
  imgSetUnit(&img, tac.unit);
  strcpy(img.studyNr, tac.studynr);
  img.zoom=zoom; 
  img.decayCorrection=IMG_DC_CORRECTED;
  (void)imgSetScanner(&img, scanner_type);
  img._fileFormat=IMG_UNKNOWN;
  img.planeNumber[0]=1;
  for(fi=0; fi<tac.frameNr; fi++) {
    img.start[fi]=tac.x1[fi]; img.end[fi]=tac.x2[fi]; img.mid[fi]=tac.x[fi];
    if(tac.timeunit==TUNIT_MIN) {
      img.start[fi]*=60.; img.end[fi]*=60.; img.mid[fi]*=60.;}
  }
  img.sizex=img.sizey=img.sizez=pxlsize;
 
 
  /*
   *  Simulate image contents one frame at a time
   */
  cx=xpos+0.5*(double)img.dimx*img.sizex;
  cy=ypos+0.5*(double)img.dimy*img.sizey;
  if(verbose>2) printf("vessel_pos := %g,%g\n", cx, cy);
 
  ret=idiSimulateTubeImgPlane(pve_method, &img, 0, cx, cy, 0.5*diameter, fwhm,
                              tac.voi[1].y, tac.voi[0].y);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot simulate image (%d).\n", ret);
    imgEmpty(&img); dftEmpty(&tac);
    return(8);
  }
  if(pve_method==1) {
    double s=fwhm/2.355;
    s/=img.sizex;
    ret=imgGaussianFIRFilter(&img, s, s, 0, 1.0E-06, verbose-3);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot simulate PVE (%d).\n", ret);
      imgEmpty(&img); dftEmpty(&tac);
      return(9);
    }
  }
  dftEmpty(&tac);
 
 
  /*
   *  Save image
   */
  if(verbose>1) printf("Writing image file %s ...\n", imgfile);
  if(imgWrite(imgfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg); 
    imgEmpty(&img); return(11);
  }
  if(verbose>0) printf("%s written\n", imgfile);
  imgEmpty(&img);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/