simimyoc.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.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"
/*****************************************************************************/
#define SBFNR 12
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulate dynamic PET image of heart. Approach is very simplistic and",
  "only applicable to software testing and basic simulations:",
  "only one image plane is produced by default, containing a ring representing",
  "myocardial muscle surrounding LV cavity.",
  " ",
  "Usage: @P [options] tacfile imgfile",
  " ",
  "TAC data file must contain the time-activity curves of myocardial muscle,",
  "blood in LV cavity, and optionally background surrounding the heart.",
  "Sample times in TAC file are used as time frames of the simulated image.",
  " ",
  "Options:",
  " -dim=<Image x and y dimension in pixels>",
  "     Set image dimension; by default 128.",
  " -3D",
  "     Simulation in 3D (note that computation may take time).",
  " -pxlsize=<Voxel size (mm)>",
  "     Set image voxel size; by default 1 mm.",
  " -fwhm=<FWHM (mm)>",
  "     Set image resolution; by default 8 mm.",
  " -thickness=<Myocardial wall thickness (mm)>",
  "     Set (relaxed) myocardial wall thickness; by default 10 mm.",
  " -diameter=<Max LV cavity diameter (mm)>",
  "     Set maximal LV cavity diameter; by default 52 mm.",
  " -diamin=<Min LV cavity diameter (fraction)>",
  "     Set minimum LV cavity diameter as fraction of max diameter;",
  "     by default 0.62.",
  " -bif=<Image filename>",
  "     Save heart beat phases separately as image planes;", 
  "     not applicable with option -3D.",
//  " -hbr=<Heart beat rate (beats/min)",
//  "     Set heart beat rate; by default 70 per min",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example: @P tacs.dat myoc.v",
  " ",
  "See also: simiart, fvar4tac, imgfiltg, sim_3tcm, b2t_h2o, imgadd, simcirc",
  " ",
  "Keywords: simulation, image, software testing, myocardium, LV",
  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      tacfile[FILENAME_MAX], imgfile[FILENAME_MAX], biffile[FILENAME_MAX];
  int       dim=128, dimz=0;
  double    pxlsize=1.0, fwhm=8.0;
  double    thickness=10.0, diameter=52.0, diamin=0.62;
  int       scanner_type=SCANNER_HRPLUS;
  float     zoom=1.41421356;
 
  char     *cptr;
  int       ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=imgfile[0]=biffile[0]=(char)0;
 
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1;
    if(strncasecmp(cptr, "DIM=", 4)==0) {
      dim=atoi(cptr+4); if(dim>4) continue;
    } else if(strcasecmp(cptr, "3D")==0) {
      dimz=1000; continue;
    } else if(strncasecmp(cptr, "PXLSIZE=", 8)==0) {
      if(!atof_with_check(cptr+8, &pxlsize) && pxlsize>0.0) continue;
    } else if(strncasecmp(cptr, "DIAMETER=", 9)==0) {
      if(!atof_with_check(cptr+9, &diameter) && diameter>0.0) continue;
    } else if(strncasecmp(cptr, "DIAMIN=", 7)==0) {
      if(!atof_with_check(cptr+7, &diamin) && diamin>0.0 && diamin<=1.0) continue;
    } else if(strncasecmp(cptr, "FWHM=", 5)==0) {
      if(!atof_with_check(cptr+5, &fwhm) && fwhm>=0.0) continue;
    } else if(strncasecmp(cptr, "THICKNESS=", 10)==0) {
      if(!atof_with_check(cptr+10, &thickness) && thickness>0.0) continue;
    } else if(strncasecmp(cptr, "BIF=", 4)==0) {
      strlcpy(biffile, cptr+4, FILENAME_MAX); if(strlen(biffile)>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 */
  if(ai<argc) {strlcpy(tacfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(imgfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Is something missing? */
  if(!imgfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* Set Z dimension to 1 or to dim */
  if(dimz>0) dimz=dim; else dimz=1;
  if(dimz>1 && biffile[0]) {
    fprintf(stderr, "Warning: 3D beating image can not be saved; ignored.\n");
    strcpy(biffile, "");
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("imgfile := %s\n", imgfile);
    if(biffile[0]) printf("biffile := %s\n", biffile);
    printf("dim := %d\n", dim);
    printf("dimz := %d\n", dimz);
    printf("pxlsize := %g\n", pxlsize);
    printf("diameter := %g\n", diameter);
    printf("diamin := %g\n", diamin);
    printf("fwhm := %g\n", fwhm);
//    printf("hbr := %g\n", hbr);
    printf("thickness := %g\n", thickness);
  }
 
 
  /* 
   *  Read the TACs to be placed in image
   */
  DFT tac; dftInit(&tac);
  if(verbose>1) printf("reading TACs\n");
  if(dftRead(tacfile, &tac)) {
    fprintf(stderr, "Error in reading '%s': %s\n", tacfile, dfterrmsg);
    return(2);
  }
  if(tac.voiNr<2) {
    fprintf(stderr, "Error: TAC file does not contain necessary TACs.\n");
    dftEmpty(&tac); return(2);
  }
  if(tac.voiNr>3) {
    fprintf(stderr, "Error: TAC file contains more than three TACs.\n");
    dftEmpty(&tac); return(2);
  }
 
  /*
   *  Allocate memory for the image
   */
  if(verbose>1) printf("allocating memory for image data\n");
  IMG img; imgInit(&img);
  ret=imgAllocate(&img, dimz, 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>1) printf("setting image headers\n");
  img.type=IMG_TYPE_IMAGE;
  imgSetUnit(&img, tac.unit);
  strlcpy(img.studyNr, tac.studynr, MAX_STUDYNR_LEN);
  img.zoom=zoom; 
  img.decayCorrection=IMG_DC_CORRECTED;
  (void)imgSetScanner(&img, scanner_type);
  img._fileFormat=IMG_UNKNOWN;
  for(int 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;
  for(int zi=0; zi<dimz; zi++) img.planeNumber[zi]=1+zi;
 
 
  /* Allocate also memory for an image with heart beat phases as planes,
     if necessary; not done if simulating 3D heart. */
  IMG bimg; imgInit(&bimg);
  if(biffile[0] && diamin<1.0) {
    if(verbose>1) printf("allocating memory for heart beat phases image\n");
    ret=imgAllocateWithHeader(&bimg, SBFNR, img.dimy, img.dimx, tac.frameNr, &img);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot allocate memory for image data.\n");
      dftEmpty(&tac); imgEmpty(&img); return(4);
    }
    for(int i=0; i<SBFNR; i++) bimg.planeNumber[i]=1+i;
  } else biffile[0]=(char)0;
 
  /* Allocate memory for a single image matrix for temp storage */
  IMG simg; imgInit(&simg);
  if(verbose>1) printf("allocating memory for temp image\n");
  ret=imgAllocateWithHeader(&simg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot allocate memory for image data.\n");
    dftEmpty(&tac); imgEmpty(&img); imgEmpty(&bimg); return(4);
  }
 
 
  /*
   *  Simulate image contents
   */
  if(verbose>0) printf("simulation...\n");
 
  /* Myocardial ring is simulated with 12 different diameters in the range
     [maxdiam,mindiam], assuming that during each PET frame each of these
     radius is present the below specified fraction of time */
  double rtfract[SBFNR];
  if(diamin==1.0) {
    rtfract[0]=1.0; for(int i=1; i<SBFNR; i++) rtfract[i]=0.0;
  } else {
    for(int i=0; i<SBFNR; i++) rtfract[i]=1.0;
  }
 
  /* Set the centre of heart in image */
  double cx, cy, cz;
  cx=0.5*(double)img.dimx*img.sizex;
  cy=0.5*(double)img.dimy*img.sizey;
  cz=0.5*(double)img.dimz*img.sizez;
 
  /* 2D: Calculate muscle wall area in one image plane (excluding Pi) */
  /* 3D: Calculate muscle wall volume (excluding Pi and division by 6) */
  double a;
  {
    double r1inner, r1outer, d1, d2; 
    r1inner=0.5*diameter; r1outer=r1inner+thickness;
    if(simg.dimz==1) {
      a=r1outer*r1outer - r1inner*r1inner;
    } else {
      d1=2.0*r1outer; d2=2.0*r1inner; 
      a=d1*d1*d1 - d2*d2*d2;
    }
  }
  if(verbose>2) printf("a := %g\n", a);
 
  double dI, dO, c_cav, c_myo, c_bkg;
  /* Process frame-by-frame */
  ret=0;
  for(int fi=0; fi<img.dimt; fi++) {
    if(verbose>1) printf("Frame %d/%d\n", 1+fi, img.dimt);
 
    c_myo=tac.voi[0].y[fi]; c_cav=tac.voi[1].y[fi];
    if(tac.voiNr>2) c_bkg=tac.voi[2].y[fi]; else c_bkg=0.0;
 
    /* Initial inner diameter */
    dI=diameter;
 
    /* Simulate image with one radius at a time and sum it to the frame image */
    for(int ri=0; ri<SBFNR; ri++) if(rtfract[ri]>0.0) {
 
      /* Calculate the outer muscle wall diameter, assuming that muscle
         wall area (2D) or volume (3D) does not change when diameter changes */
      double f=0.5*dI; // inner radius
      // outer radius
      double r;
      if(simg.dimz==1) {
        r=sqrt(a + f*f); 
        dO=2.0*r;
      } else {
        dO=cbrt(a + dI*dI*dI); // V=Pi*(Diameter^3)/6
        r=0.5*dO;
      }
      if(verbose>2) {
        printf("Radius %d/%d\n", 1+ri, SBFNR);
        printf("inner diameter := %g\n", dI);
        printf("outer diameter := %g\n", dO);
        if(simg.dimz==1) printf("muscle area := %g\n", M_PI*(r*r-f*f));
        else printf("muscle volume := %g\n", M_PI*(dO*dO*dO-dI*dI*dI)/6.0);
      }
 
      /* Simulate myocardial ring */
      for(long long i=0; i<(long long)simg.dimx*simg.dimy*simg.dimz; i++) simg.pixel[i]=0.0;
      if(simg.dimz==1)
        ret=imgSimulateRing(&simg, 0, 0, cx, cy, f, r, c_myo, c_cav, c_bkg, verbose-3);
      else
        ret=imgSimulateSphere(&simg, 0, cx, cy, cz, f, r, c_myo, c_cav, c_bkg, verbose-3);
      if(ret!=0) break;
 
      /* Add weighted simulated image matrix to the image to be saved */
      int zi, yi, xi;
      for(zi=0; zi<img.dimz; zi++)
        for(yi=0; yi<img.dimy; yi++)
          for(xi=0; xi<img.dimx; xi++)
            img.m[zi][yi][xi][fi]+=rtfract[ri]*simg.m[zi][yi][xi][0];
 
      /* And without weighting to bimg if required */
      if(biffile[0]) {
        for(yi=0; yi<img.dimy; yi++)
          for(xi=0; xi<img.dimx; xi++)
            bimg.m[ri][yi][xi][fi]=simg.m[0][yi][xi][0];
      }
 
      /* Decrease the inner diameter */
      dI-=diameter*(1.0-diamin)/(double)(SBFNR-1);
    } // next diameter
    if(ret!=0) break;
  } // next frame
  imgEmpty(&simg); dftEmpty(&tac);
  /* Check that simulation was successful */
  if(ret!=0) {
    fprintf(stderr, "Error: cannot simulate image (%d).\n", ret);
    imgEmpty(&img); imgEmpty(&bimg);
    return(6);
  }
 
  /* Divide the image frames by the sum of time fractions to get the average during the frame */
  if(verbose>1) printf("averaging heart beat phases\n");
  /* Sum of fractional times */
  double rtsum=0.0;
  for(int i=0; i<SBFNR; i++) if(rtfract[i]>0.0) rtsum+=rtfract[i];
  if(verbose>2) printf("rtsum := %g\n", rtsum);
  ret=imgArithmConst(&img, (float)rtsum, ':', -1.0, verbose-10);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot process simulated image (%d).\n", ret);
    imgEmpty(&img); imgEmpty(&bimg);
    return(7);
  }
 
  /*
   *  Simulate the image resolution effect
   */
  if(fwhm>0.0) {
    if(verbose>0) printf("simulating resolution effect\n");
    char tmp[256];
    double s, d;
    /* Calculate filter SD */
    s=fwhm/2.354820; // 2*sqrt(2*ln(2))
    /* Convert filter SD to pixels */
    d=0.5*(img.sizex+img.sizey); // Z size is assumed to be the same
    s/=d; // printf("filter_sd_in_pixels := %g\n", s);
    /* Apply Gaussian filter */
    if(img.dimz==1)
      ret=imgGaussianFIRFilter(&img, s, s, 0.0, 1.0E-06, verbose-3);
    else
      ret=imgGaussianFIRFilter(&img, s, s, s, 1.0E-06, verbose-3);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot filter simulated image (%d).\n", ret);
      if(verbose>0) fprintf(stderr, "       %s.\n", tmp);
      imgEmpty(&img); imgEmpty(&bimg);
      return(8);
    }
    /* Same for beat image, if necessary (only 2D) */
    if(biffile[0]) {
      ret=imgGaussianFIRFilter(&bimg, s, s, s, 1.0E-06, verbose-3);
      if(ret!=0) {
        fprintf(stderr, "Error: cannot filter simulated image (%d).\n", ret);
        if(verbose>0) fprintf(stderr, "       %s.\n", tmp);
        imgEmpty(&img); imgEmpty(&bimg);
        return(9);
      }
    }
  }
 
 
  /*
   *  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); imgEmpty(&bimg); return(11);
  }
  if(verbose>0) printf("%s written.\n", imgfile);
  imgEmpty(&img);
 
  /*
   *  Save beating image, if necessary
   */
  if(biffile[0]) {
    if(verbose>1) printf("Writing image file %s ...\n", biffile);
    if(imgWrite(biffile, &bimg)) {
      fprintf(stderr, "Error: %s\n", bimg.statmsg); 
      imgEmpty(&bimg); return(13);
    }
    if(verbose>0) printf("%s written.\n", biffile);
    imgEmpty(&bimg);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/