ecatfbp.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 filtered back-projection (FBP).",
  "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:",
  " -filter=<None|Ramp|Butter|Hann|Hamm|Parzen|Shepp>",
  "     Set the reconstruction filter; by default Hann.",
  " -cutoff=<value>",
  "     Set cut-off value; by default 0.3.",
  " -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.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "  @P -dim=256 a2345dy1.scn a2345dy1.img",
  " ",
  "See also: ecatnorm, edecay, ecalibr, atnmake, eframe, ecatmrp, img2scn",
  " ",
  "Keywords: ECAT, sinogram, image, reconstruction, FBP",
  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          filter=FBP_FILTER_HANN;
  float        cutoff=0.3;
  float        rotate=0.0;
  float        zoom=1.0;
  float        shiftX=0.0, shiftY=0.0;
  int          pervol=0;
  int          imgDim=128;
  char        *cptr;
  int          ret;
  static char *fbpFilterName[]={
                 "None","Ramp","Butter","Hann","Hamm","Parzen","Shepp", 0
               };
 
 
  /*
   *  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, "FILTER=", 7)==0) {
      cptr+=7; filter=-1; 
      int i=0;
      while(fbpFilterName[i]!=0) {
        if(!strcasecmp(cptr, fbpFilterName[i])) {
          filter=i; break;
        } else i++;
      }
      if(filter>=0) continue;
    } else if(strncasecmp(cptr, "DIM=", 4)==0) {
      imgDim=atoi(cptr+4); if(imgDim>1 && imgDim<=2048) 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, "CUTOFF=", 7)==0) {
      cutoff=atof_dpi(cptr+7); if(cutoff>0.0 && cutoff<1.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("Filter := %s\n", fbpFilterName[filter]);
    printf("Cut-off := %g\n", cutoff);
    printf("Zoom := %g\n", zoom);
    printf("Dimension := %d\n", imgDim);
    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);
  }
 
  /* 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("FBP reconstruction...\n");
  ret=imgFBP(&scn, &img, imgDim, zoom, filter, cutoff, shiftX, shiftY, rotate,
             verbose-2);
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/