imgdysmo.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Smoothing of dynamic PET image in ECAT 6.3 and 7, Analyze 7.5, and NIfTI",
  "format.",
  "Replaces pixel TACs with an average TAC of up to 5x5x5 neighbouring pixels",
  "which are selected based on the MRL (maximum run length) and AUC difference.",
  "Note that before application in data analysis, the quantitativity of",
  "the regional results must be validated.",
  " ",
  "Usage: @P [Options] dynamic_image smoothed_image",
  " ",
  "Options:",
  " -m=<3|5>",
  "     Set filter mask to 3x3x3 or 5x5x5 pixels (default).",
  " -s=<1|2|3|4|5>",
  "     Set smoothing strength: 1=min, 2=default, 5=max.",
  " -a=<nr>",
  "     Set the nr of pixels to average inside smoothing mask;",
  "     by default this is set automatically based on the mask size.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1.",
  "  @P b123dy1.v b123dy1smoothed.v",
  " ",
  "See also: imgthrs, imgbkgrm, imgfsegm, imgfiltg, fvar4img",
  " ",
  "Keywords: image, smoothing, modelling",
  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      ret=0, maskSize=5, smoothing=2, avgNr=0;
  char     imgfile[FILENAME_MAX], smofile[FILENAME_MAX], tmp[FILENAME_MAX];
  char    *cptr;
  IMG      img, smoimg;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=smofile[0]=(char)0;
  imgInit(&img); imgInit(&smoimg);
  /* 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;
    cptr=argv[ai]+1;
    if(strncasecmp(cptr, "M=", 2)==0) {
      cptr+=2; maskSize=atoi(cptr); if(maskSize==3 || maskSize==5) continue;
    } else if(strncasecmp(cptr, "S=", 2)==0) {
      cptr+=2; smoothing=atoi(cptr); if(smoothing>0 && smoothing<6) continue;
    } else if(strncasecmp(cptr, "A=", 2)==0) {
      cptr+=2; avgNr=atoi(cptr); if(avgNr>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 */
  for(; ai<argc; ai++) {
    if(!imgfile[0]) {
      strcpy(imgfile, argv[ai]); continue;
    } else if(!smofile[0]) {
      strcpy(smofile, argv[ai]); continue;
    }
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!smofile[0]) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n",
      argv[0]);
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("imgfile := %s\n", imgfile);
    printf("smofile := %s\n", smofile);
    printf("maskSize := %d\n", maskSize);
    printf("smoothing=%d\n", smoothing);
    printf("avgNr := %d\n", avgNr);
  }
 
 
  /*
   *  Read dynamic image
   */
  if(verbose>=0) fprintf(stdout, "reading dynamic image %s\n", imgfile);
  if(imgRead(imgfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    if(verbose>1) imgInfo(&img);
    imgEmpty(&img);
    return(2);
  }
  if(verbose>10) imgInfo(&img);
  /* Check that image is image */
  if(img.type!=IMG_TYPE_IMAGE) {
    fprintf(stderr, "Error: raw PET data cannot be smoothed with this program.\n");
    imgEmpty(&img); return(1);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  /* Check that MRL and frame nr are ok */
  if(img.dimt<2) {
    fprintf(stderr, "Error: image has only one time frame.\n");
    imgEmpty(&img); return(1);
  }
  /* Check that plane nr is ok */
  if(img.dimz<3) {
    fprintf(stderr, "Warning: smoothing would work better for images with more planes.\n");
  }
 
 
  /*
   *  Calculate the smoothed image
   */
  /* Set avgNr if necessary */
  if(avgNr<2) switch(maskSize) {
    case 3:
      switch(smoothing) {
        case 1: avgNr=3; break;
        case 2: avgNr=5; break;
        case 3: avgNr=7; break;
        case 4: avgNr=9; break;
        case 5: avgNr=11; break;
        default: avgNr=5; break;
      }
      break;
    case 5:
    default:
      switch(smoothing) {
        case 1: avgNr=5; break;
        case 2: avgNr=9; break;
        case 3: avgNr=15; break;
        case 4: avgNr=25; break;
        case 5: avgNr=35; break;
        default: avgNr=9; break;
      }
      break;
  } else { /* Check that avgNr is not set to stupid values by the user */
    ret=maskSize*maskSize*maskSize/2; if(avgNr>ret) avgNr=ret;
  }
  if(verbose>1) {
    printf("Average of %d best pixels inside %dx%dx%d matrix.\n",
      avgNr, maskSize, maskSize, maskSize);
  }
  if(verbose>=0) printf("Calculating smoothed image\n");
  ret=imgsegmSimilar(&img, maskSize, avgNr, &smoimg, verbose-1);
  if(ret) {
    fprintf(stderr, "Error (%d): cannot smooth.\n", ret);
    imgEmpty(&img); imgEmpty(&smoimg); return(4);
  }
 
  /* Save the smoothed image */
  if(verbose>=0) printf("writing smoothed image %s\n", smofile);
  backupExistingFile(smofile, NULL, tmp);
  if(imgWrite(smofile, &smoimg)) {
    fprintf(stderr, "Error: %s\n", smoimg.statmsg);
    imgEmpty(&img); imgEmpty(&smoimg);
    return(11);
  }
 
  /*
   *  Free up memory
   */
  imgEmpty(&img); imgEmpty(&smoimg);
 
  if(verbose>=0) fprintf(stdout, "done.\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/