imgfiltg.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 <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Applying Gaussian filter for a static or dynamic PET image in ECAT, NIfTI,",
  "or Analyze format. Filtering can be applied to x,y,z-volume (3D) by",
  "specifying the SDs of Gaussian filter for x,y- and z-dimensions, or only",
  "to the x,y-planes by giving the SD only for x,y-dimensions.",
  "SDs are assumed to be in pixels, if units are not given.",
  "SD=FWHM/2.355.",
  " ",
  "Usage: @P [Options] imgfile outputfile SDxy [SDz]",
  " ",
  "Options:",
  " -FWHM",
  "     FWHM value(s) are given instead of SDs.",
  " -method=<<FIR>|<AM>>",
  "     Specify the Gaussian filtering method; default is the finite impulse",
  "     response method (FIR); Alvarez-Mazorra approximate Gaussian image",
  "     filtering method (AM) is less accurate but may be faster.",
  " -tol=<tolerance>",
  "     Default tolerance can be changed with this option.",
  " -steps=<nr>",
  "     More steps provides better accuracy but longer execution time;",
  "     default is 4. Applies only to AM method.",
  " -xysize=<pixel size>",
  "     Set image pixel x,y size (mm), if missing or wrong in image data.",
  " -zsize=<pixel size>",
  "     Set image pixel z size (mm), if missing or wrong in image data.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example #1: 3D filtering with Gaussian SDs given in mm",
  "     @P s5998dy1.v s5998dy1_filt.v 3.4mm 4.1mm",
  " ",
  "Example #2: 2D filtering with Gaussian SDs given in pixels",
  "     @P s5998dy1.v s5998dy1_filt.v 4pxl",
  " ",
  "Example #3: 3D filtering with FWHMs given in mm",
  "     @P -fwhm s5998dy1.v s5998dy1_filt.v 8mm 9.66mm",
  " ",
  "See also: imgdysmo, imgfsegm, imgthrs, imgbkgrm, fvar4img, imgprofi",
  " ",
  "Keywords: image, smoothing, gaussian, FWHM",
  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     petfile[FILENAME_MAX], outfile[FILENAME_MAX];
  int      gaussian_method=1; // 1= FIR; 2=Alvarez-Mazorra; 3=EBox
  double   xsize=nan(""), ysize=nan(""), zsize=nan("");
  int      useSD=1; // 0=FWHM; 1=SD
  double   SDx=nan(""), SDy=nan(""), SDz=nan("");
  int      SDxUnit, SDyUnit, SDzUnit; // 0=pixels, 1=mm, 2=cm
  double   tolerance=1.0E-04;
  int      stepNr=4;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=outfile[0]=(char)0;
  SDxUnit=SDyUnit=SDzUnit=0; 
  /* Get options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strncasecmp(cptr, "METHOD=", 7)==0 && strlen(cptr)>7) {
      cptr+=7;
      if(strcasecmp(cptr, "FIR")==0) {gaussian_method=1; continue;}
      if(strcasecmp(cptr, "AM")==0) {gaussian_method=2; continue;}
      if(strcasecmp(cptr, "EBOX")==0) {gaussian_method=3; continue;}
    } else if(strcasecmp(cptr, "SD")==0) {
      useSD=1; continue;
    } else if(strcasecmp(cptr, "FWHM")==0) {
      useSD=0; continue;
    } else if(strncasecmp(cptr, "STEPS=", 6)==0) {
      stepNr=atoi(cptr+6); if(stepNr>0) continue;
    } else if(strncasecmp(cptr, "XYSIZE=", 7)==0) {
      xsize=ysize=atof_dpi(cptr+7); if(xsize>0.0) continue;
    } else if(strncasecmp(cptr, "ZSIZE=", 6)==0) {
      zsize=atof_dpi(cptr+6); if(zsize>0.0) continue;
    } else if(strncasecmp(cptr, "TOL=", 4)==0) {
      tolerance=atof_dpi(cptr+4); if(tolerance>0.0) continue;
    } else if(strncasecmp(cptr, "TOLERANCE=", 10)==0) {
      tolerance=atof_dpi(cptr+10); if(tolerance>0.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(petfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    strlcpy(outfile, argv[ai++], FILENAME_MAX);
    double v; // check that user did not forget file name
    if(atof_with_check(outfile, &v)==0) {
      fprintf(stderr, "Error: invalid output file name '%s'\n", outfile);
      return(1);
    }
  }
  if(ai<argc) {
    double v;
    int ret=atof_with_check(argv[ai], &v);
    if(ret==0 && v>=0.0) {
      SDx=SDy=(float)v;
      /* Check if unit was specified too */
      if(strcasestr(argv[ai], "pxl")!=NULL) SDxUnit=SDyUnit=0;
      if(strcasestr(argv[ai], "mm")!=NULL) SDxUnit=SDyUnit=1;
      if(strcasestr(argv[ai], "cm")!=NULL) SDxUnit=SDyUnit=2;
      ai++;
    } else {
      fprintf(stderr, "Error: invalid filter S.D.\n");
      return(1);
    }
  }
  if(ai<argc) {
    double v;
    int ret=atof_with_check(argv[ai], &v);
    if(ret==0 && v>=0.0) {
      SDz=(float)v;
      /* Check if unit was specified too */
      if(strcasestr(argv[ai], "pxl")!=NULL) SDzUnit=0;
      if(strcasestr(argv[ai], "mm")!=NULL) SDzUnit=1;
      if(strcasestr(argv[ai], "cm")!=NULL) SDzUnit=2;
      ai++;
    } else {
      fprintf(stderr, "Error: invalid filter S.D.\n");
      return(1);
    }
  }
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Did we get all the information that we need? */
  if(!outfile[0] || !(SDx>=0.0)) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n", argv[0]);
    return(1);
  }
  if(strcasecmp(outfile, petfile)==0) {
    fprintf(stderr, "Error: check the output file name.\n");
    return(1);
  }
  if(!(SDx>0.0) && !(SDy>0.0) && !(SDz>0.0)) {
    fprintf(stderr, "Warning: with given SD no filtering is applied.\n");
    fflush(stderr);
  }
 
 
  /* Convert FWHM to SD, if necessary */
  if(!useSD) {SDx/=2.355; SDy/=2.355; if(SDz>0.0) SDz/=2.355;}
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("outfile := %s\n", outfile);
    if(!isnan(xsize)) printf("pixel_xysize := %g mm\n", xsize);
    if(!isnan(zsize)) printf("pixel_zsize := %g mm\n", zsize);
    printf("gaussian_method := %d\n", gaussian_method);
    printf("stepNr := %d\n", stepNr);
    fflush(stdout);
  }
 
 
  /*
   *  Read the contents of the PET file to img data structure
   */
  IMG img; imgInit(&img);
  if(verbose>0) printf("reading %s\n", petfile);
  {
    int ret=imgRead(petfile, &img);
    if(ret) {
      fprintf(stderr, "Error: %s\n", img.statmsg);
      if(verbose>1) printf("ret=%d\n", ret);
      imgEmpty(&img); return(2);
    }
    /* Check if PET data is raw or image */
    if(img.type!=IMG_TYPE_IMAGE) {
      fprintf(stderr, "Error: %s is not an image.\n", petfile);
      imgEmpty(&img); return(2);
    }
    if(imgNaNs(&img, 1)>0)
      if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  }
 
  /* Set or check pixel sizes */
  if(xsize>0.0) img.sizex=xsize;
  if(ysize>0.0) img.sizey=ysize;
  if(zsize>0.0) img.sizez=zsize;
  if(verbose>1) printf("image_xyz_sizes := %g %g %g\n", img.sizex, img.sizey, img.sizez);
  if(!isnan(xsize) || !isnan(ysize)) {
    if(img.sizex!=img.sizey) {
      fprintf(stderr, "Warning: different pixel x and y size.\n");
      fflush(stderr);
    }
  }
  if(verbose>1)
    printf("image_xyzf_dimensions := %d %d %d %d\n", img.dimx, img.dimy, img.dimz, img.dimt);
 
  /* Check that image size in smoothing dimension is at least 4 */
  if(SDz>0 && img.dimz<4) {
    fprintf(stderr, "Error: invalid Z dimension for 3D filtering.\n");
    imgEmpty(&img); return(3);
  }
  if((SDx>0 && img.dimx<4) || (SDy>0 && img.dimy<4)) {
    fprintf(stderr, "Error: invalid x,y dimensions for filtering.\n");
    imgEmpty(&img); return(3);
  }
 
 
  /*
   *  Set Gaussian SD in units of pixels
   */
  if((SDxUnit!=0 && !(img.sizex>0.0)) || (SDyUnit!=0 && !(img.sizey>0.0)) ||
     (SDzUnit!=0 && !(img.sizez>0.0)))
  {
    fprintf(stderr, "Error: unknown image pixel size.\n");
    imgEmpty(&img); return(3);
  }
  /* Convert cm -> mm */
  if(SDxUnit==2) {SDx*=10.0; SDxUnit=1;}
  if(SDyUnit==2) {SDy*=10.0; SDyUnit=1;}
  if(SDzUnit==2) {SDz*=10.0; SDzUnit=1;}
  /* Convert mm -> pixels */
  if(SDxUnit==1) {SDx/=img.sizex; SDxUnit=0;}
  if(SDyUnit==1) {SDy/=img.sizey; SDyUnit=0;}
  if(SDzUnit==1) {SDz/=img.sizez; SDzUnit=0;}
  if(verbose>1) {
    printf("gaussian_sdx := %g pxl\n", SDx);
    printf("gaussian_sdy := %g pxl\n", SDy);
    printf("gaussian_sdz := %g pxl\n", SDz);
  }
 
 
  /*
   *  Gaussian smoothing
   */
  int ret=0;
  if(gaussian_method==1) {
    if(verbose>0) fprintf(stdout, "Gaussian FIR filtering\n");
    ret=imgGaussianFIRFilter(&img, SDx, SDy, SDz, tolerance, verbose-2);
  } else if(gaussian_method==2) {
    if(verbose>0) fprintf(stdout, "Gaussian AM filtering\n");
    ret=imgGaussianAMFilter(&img, SDx, SDy, SDz, stepNr, tolerance, verbose);
  } else { // EBox
    if(verbose>0) fprintf(stdout, "Gaussian Ebox filtering\n");
    ret=imgGaussianEBoxFilter(&img, SDx, SDy, SDz, stepNr, verbose);
  }
  if(ret!=0) {
    fprintf(stderr, "Error: cannot filter the image.\n");
    if(verbose>1) fprintf(stderr, "ret := %d\n", ret);
    imgEmpty(&img); return(7);
  }
 
 
  /*
   *  Save filtered image
   */
  if(verbose>2) fprintf(stdout, "writing %s\n", outfile);
  if(imgWrite(outfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg); 
    imgEmpty(&img);
    return(11);
  }
  if(verbose>0) printf("Written %s\n", outfile);
  imgEmpty(&img);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/