ecatmrp.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 <ctype.h>
#include <dirent.h>
#include <sys/stat.h>
#include <sys/file.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
#include "libtpcrec.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Reconstruct PET image from attenuation and normalization corrected",
  "2D PET sinogram using median root prior (MRP) [1,2,3].",
  "This program does not correct image data for decay.",
  "Image is corrected for the dead-time based on factors in the sinogram.",
  "Image data is written in units ECAT counts per pixel per second.",
  " ",
  "Not for clinical use! Program is not thoroughly tested!",
  " ",
  "Usage: @P [Options] scnfile imgfile",
  " ",
  "Options:",
  " -zoom=<value>",
  "     Set zoom factor; by default 1.0 (no zoom).",
  " -dim=<value>",
  "     Set image x and y dimensions; by default 128.",
  " -rot[ation]=<value>",
  "     Set image rotation in degrees, -180 - 180; by default 0 (no rotation).",
  " -x=<value>",
  "     Set image shift in x direction (in cm); by default 0.",
  " -y=<value>",
  "     Set image shift in y direction (in cm); by default 0.",
  " -perv",
  "     Divide image voxel values by voxel volume (mL); by default the voxel",
  "     volume is corrected later during image calibration.",
  " -beta=<value>",
  "     Set beta value, usually 0.1 - 0.9; by default 0.3.",
  " -mask=<3|5>",
  "     Set mask dimension for median filter, either 3 (3x3) or 5 (5x5 without",
  "     corener pixels); by default 3.",
  " -iter=<value>",
  "     Set maximum number of iterations; by default 150.",
  " -skip=<value>",
  "     Set the number of iterations without prior; by default 1.",
  " -os=<1|2|4|8|16>",
  "     Set the number of OS sets (acceleration); by default 1.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "  @P -dim=256 a2345dy1.scn a2345dy1.img",
  " ",
  "References:",
  "1. Alenius S, Ruotsalainen U. Bayesian image reconstruction for emission",
  "   tomography based on median root prior.",
  "   Eur J Nucl Med. 1997;24(3): 258-265.",
  "2. Alenius S, Ruotsalainen U. Generalization of median root prior",
  "   reconstruction. IEEE Trans Med Imaging 2002; 21(11): 1413-1420.",
  "3. Katoh C, Ruotsalainen U, Laine H, Alenius S, Iida H, Nuutila P, Knuuti J.",
  "   Iterative reconstruction based on median root prior in quantification of",
  "   myocardial blood flow and oxygen metabolism.",
  "   J Nucl Med. 1999;40(5):862-867.",
  " ",
  "See also: ecatnorm, edecay, ecalibr, atnmake, eframe, ecatfbp, img2scn",
  " ",
  "Keywords: ECAT, sinogram, image, reconstruction, MRP",
  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         scnfile[FILENAME_MAX], imgfile[FILENAME_MAX];
  int          maxIterNr=150;
  int          skipPriorNr=1;
  float        beta=0.3;
  float        rotate=0.0;
  float        zoom=1.0;
  float        shiftX=0.0, shiftY=0.0;
  int          pervol=0;
  int          imgDim=128;
  int          maskDim=3;
  int          osSetNr=1;
  char        *cptr;
  int          ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  scnfile[0]=imgfile[0]=(char)0;
  /* 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, "ITER=", 5)==0) {
      maxIterNr=atoi(cptr+5); if(maxIterNr>1) continue;
    } else if(strncasecmp(cptr, "SKIP=", 5)==0) {
      skipPriorNr=atoi(cptr+5); if(skipPriorNr>0) continue;
    } else if(strncasecmp(cptr, "OS=", 3)==0) {
      osSetNr=atoi(cptr+3); 
      if(osSetNr==0 || osSetNr==1 || osSetNr==2) continue;
      if(osSetNr==4 || osSetNr==8 || osSetNr==16) continue;
      if(osSetNr==32 || osSetNr==64 || osSetNr==128) continue;
    } else if(strncasecmp(cptr, "DIM=", 4)==0) {
      imgDim=atoi(cptr+4); if(imgDim>1 && imgDim<=2048) continue;
    } else if(strncasecmp(cptr, "MASK=", 5)==0) {
      maskDim=atoi(cptr+5); if(maskDim==3 || maskDim==5) continue;
    } else if(strncasecmp(cptr, "ZOOM=", 5)==0) {
      zoom=atof_dpi(cptr+5); if(zoom>0.1 && zoom<100.0) continue;
    } else if(strncasecmp(cptr, "BETA=", 5)==0) {
      beta=atof_dpi(cptr+5); if(beta>0.05 && beta<100.0) continue;
    } else if(strncasecmp(cptr, "ROTATION=", 9)==0) {
      double v;
      ret=atof_with_check(cptr+9, &v);
      if(v>=-180.0 && v<=180.0) {rotate=v; continue;}
    } else if(strncasecmp(cptr, "ROT=", 4)==0) {
      double v;
      ret=atof_with_check(cptr+4, &v);
      if(v>-360.0 && v<360.0) {rotate=v; continue;}
    } else if(strncasecmp(cptr, "X=", 2)==0) {
      double v;
      ret=atof_with_check(cptr+2, &v);
      if(v>-50.0 && v<50.0) {shiftX=10.0*v; continue;}
    } else if(strncasecmp(cptr, "Y=", 2)==0) {
      double v;
      ret=atof_with_check(cptr+2, &v);
      if(v>-50.0 && v<50.0) {shiftY=10.0*v; 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(scnfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(imgfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: too many arguments.\n");
    return(1);
  }
 
  /* Is something missing? */
  if(!imgfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("scnfile := %s\n", scnfile);
    printf("imgfile := %s\n", imgfile);
    printf("Beta := %g\n", beta);
    printf("Zoom := %g\n", zoom);
    printf("Dimension := %d\n", imgDim);
    printf("Mask := %d\n", maskDim);
    printf("Rotation := %g\n", rotate);
    printf("X_shift_mm := %g\n", shiftX);
    printf("Y_shift_mm := %g\n", shiftY);
    printf("volume_correction := %d\n", pervol);
    printf("maxIterNr := %d\n", maxIterNr);
    printf("skipPriorNr := %d\n", skipPriorNr);
    printf("osSetNr := %d\n", osSetNr);
  }
 
  /* Check that input and output files do not have the same filename */
  if(strcasecmp(imgfile, scnfile)==0) {
    fprintf(stderr, "Error: original sinogram must not be overwritten.\n");
    return(1);
  }
  /* Check that image dimension is even */
  if(imgDim%2) {
    fprintf(stderr, "Error: invalid image dimension.\n");
    return(1);
  }
 
 
  /*
   *  Read sinogram
   */
  if(verbose>0) printf("reading sinogram %s\n", scnfile);
  IMG scn; imgInit(&scn);
  ret=imgRead(scnfile, &scn);
  if(ret) {
    fprintf(stderr, "Error: %s\n", scn.statmsg);
    if(verbose>1) printf("ret := %d\n", ret);
    return(2);
  }
  if(verbose>1) { 
    printf("dimx := %d\n", scn.dimx);
    printf("dimy := %d\n", scn.dimy);
    printf("dimz := %d\n", scn.dimz);
    printf("dimt := %d\n", scn.dimt);
  }
  /* Check that we have a sinogram */
  if(scn.type!=IMG_TYPE_RAW) {
    fprintf(stderr, "Error: file is not a sinogram.\n");
    imgEmpty(&scn);
    return(2);
  }
 
  /*
   *  Reconstruction
   */
  IMG img; imgInit(&img);
  if(verbose>0) printf("MRP reconstruction... (may take a long time)\n");
  ret=imgMRP(&scn, &img, imgDim, zoom, shiftX, shiftY, rotate,
             maxIterNr, skipPriorNr, beta, maskDim, osSetNr,
             verbose);
  if(ret) {
    fprintf(stderr, "Error in reconstruction.\n");
    if(verbose>1) printf("ret := %d\n", ret);
    imgEmpty(&scn);
    return(3);
  }
  /* Set header info in the output data */
  img.zoom=zoom;
 
  /* Pixel volume correction, if required */
  if(pervol) {
    if(verbose>0) printf("correcting for voxel volume\n");
    float vol=0.001*img.sizex*img.sizey*img.sizez;
    if(vol>0.0) {
      int zi, yi, xi, ti;
      for(zi=0; zi<img.dimz; zi++)
        for(yi=0; yi<img.dimy; yi++) 
          for(xi=0; xi<img.dimx; xi++)
            for(ti=0; ti<img.dimt; ti++)
              img.m[zi][yi][xi][ti]/=vol;
    } else {
      fprintf(stderr, "Warning: cannot correct for pixel volume.\n");
    }
  }
 
 
  /*
   *  Save the image */
  if(verbose>0) printf("writing image %s\n", imgfile);
  if(imgWrite(imgfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    imgEmpty(&img); imgEmpty(&scn);
    return(11);
  }
 
  imgEmpty(&img); imgEmpty(&scn);
  if(verbose>0) fprintf(stdout, "done.\n");
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/