imgcutof.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[] = {
  "Cut off PET image (in ECAT, NIfTI, or Analyze format) pixel values that are",
  "outside given range of values.",
  "Original image file is modified, if file name for output is not given.",
  " ",
  "Usage: @P [Options] imgfile lowerlimit upperlimit [outputfile]",
  " ",
  "Options:",
  " -zero",
  "     Set pixel value to zero if original value is outside of the range;",
  "     by default value is set to the limit.",
  " -mask",
  "     Set pixel value to one or zero, depending whether original value is",
  "     inside or outside of the range.",
  " --force",
  "     Optional output file is written even in the case that no pixel values",
  "     fell outside of limits; by default image is saved only if contents were",
  "     modified.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Lower and upper limit can be given directly as values on the command line,",
  "or as names of files containing only the value, or value with key 'lower' or",
  "'upper', respectively.",
  " ",
  "See also: imgthrs, imgmask, imgqntls, imginteg, imgunit",
  " ",
  "Keywords: image, threshold, mask",
  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];
  float    lowerThreshold=+1.0, upperThreshold=-1.0;
  int      zeroMode=0; // 0=off; 1=on
  int      maskMode=0; // 0=off; 1=on
  int      forceMode=0; // 0=off; 1=on
  char    *cptr;
  int      ret=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=outfile[0]=(char)0;
  /* Get options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strncasecmp(cptr, "ZERO", 1)==0) {
      zeroMode=1; continue;
    } else if(strncasecmp(cptr, "MASK", 1)==0) {
      maskMode=1; continue;
    } else if(strcasecmp(cptr, "F")==0 || strcasecmp(cptr, "FORCE")==0) {
      forceMode=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 */
  if(ai<argc) {strlcpy(petfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    double v;
    char *s=iftReadValue(argv[ai], "lower", 0);
    if(s!=NULL) {ret=atof_with_check(s, &v); free(s);}
    else ret=atof_with_check(argv[ai], &v);
    if(ret!=0) {
      fprintf(stderr, "Error: invalid lower limit '%s'.\n", argv[ai]);
      return(1);
    }
    lowerThreshold=v;
    ai++;
  } else {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n", argv[0]);
    return(1);
  }
  if(ai<argc) {
    double v;
    char *s=iftReadValue(argv[ai], "upper", 0);
    if(s!=NULL) {ret=atof_with_check(s, &v); free(s);}
    else ret=atof_with_check(argv[ai], &v);
    if(ret!=0) {
      fprintf(stderr, "Error: invalid lower limit '%s'.\n", argv[ai]);
      return(1);
    }
    upperThreshold=v;
    ai++;
  } else {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n", argv[0]);
    return(1);
  }
  if(ai<argc) {strlcpy(outfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
 
  /* Did we get all the information that we need? */
  if(!petfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n", argv[0]);
    return(1);
  }
  if(lowerThreshold>upperThreshold) {
    fprintf(stderr, "Error: invalid value range.\n");
    return(1);
  }
  if(lowerThreshold==upperThreshold) {
    /* Image data is stored with limited precision, therefore limits should not be equal */
    float f;
    f=0.001*lowerThreshold; lowerThreshold-=f; upperThreshold+=f;
    if(verbose>1) printf("Threshold range: %g - %g\n", lowerThreshold, upperThreshold);
  }
  /* Overwrite original file, if output file name was not given */
  if(!outfile[0]) strcpy(outfile, petfile);
  /* Mask-mode overrides zero-mode */
  if(maskMode) zeroMode=0;
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    if(outfile[0]) printf("outfile := %s\n", outfile);
    printf("upper_limit := %g\n", upperThreshold);
    printf("lower_limit := %g\n", lowerThreshold);
    printf("zeroMode := %d\n", zeroMode);
    printf("maskMode := %d\n", maskMode);
  }
 
 
  /*
   *  Read the contents of the PET file to img data structure
   */
  IMG img; imgInit(&img);
  if(verbose>0) printf("reading %s\n", petfile);
  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);
  }
 
 
  /*
   *  Cut-off
   */
  if(verbose>0) fprintf(stdout, "verifying pixel values...\n");
  if(maskMode==0) {
    /* Report the number of cut-off pixels, in verbose mode in each plane and frame */
    int fi, zi, xi, yi;
    long long nr, total_nr=0;
    if(verbose>1) printf("Nr of cut-off pixels:\n");
    for(fi=0; fi<img.dimt; fi++) {
      for(zi=0; zi<img.dimz; zi++) {
        nr=0;
        for(yi=0; yi<img.dimy; yi++) {
          for(xi=0; xi<img.dimx; xi++) {
            if(!isfinite(img.m[zi][yi][xi][fi]) || img.m[zi][yi][xi][fi]<lowerThreshold) {
              if(zeroMode==0) img.m[zi][yi][xi][fi]=lowerThreshold; else img.m[zi][yi][xi][fi]=0.0;
              nr++; continue;
            }
            if(!isfinite(img.m[zi][yi][xi][fi]) || img.m[zi][yi][xi][fi]>upperThreshold) {
              if(zeroMode==0) img.m[zi][yi][xi][fi]=upperThreshold; else img.m[zi][yi][xi][fi]=0.0;
              nr++; continue;
            }
          }
        }
        if(verbose>1 && img.dimz>1 && img.dimt>1) 
          printf("  frame %d plane %d: %lld\n", fi+1, zi+1, nr);
        total_nr+=nr;
      }
    }
    if(verbose>0) printf("cutoff_pixels := %lld\n", total_nr);
 
    /* If no pixels were cut-off, and not in force mode, then do no more */
    if(total_nr==0) {
      if(forceMode) {
        fprintf(stderr, "Warning: no pixels thresholded.\n");
      } else {
        fprintf(stderr, "Note: no pixels thresholded, image not saved.\n");
        imgEmpty(&img);
        return(0);
      }
    }
 
  } else { // Mask mode
    long long nr=0;
    for(int fi=0; fi<img.dimt; fi++) {
      for(int zi=0; zi<img.dimz; zi++) {
        for(int yi=0; yi<img.dimy; yi++) {
          for(int xi=0; xi<img.dimx; xi++) {
            if(!isfinite(img.m[zi][yi][xi][fi]) || img.m[zi][yi][xi][fi]<lowerThreshold) {
              img.m[zi][yi][xi][fi]=0.0; continue;}
            if(!isfinite(img.m[zi][yi][xi][fi]) || img.m[zi][yi][xi][fi]>upperThreshold) {
              img.m[zi][yi][xi][fi]=0.0; continue;}
            img.m[zi][yi][xi][fi]=1.0; nr++;
          }
        }
      }
    }
    if(verbose>0) printf("retained_pixels := %lld\n", nr);
 
    /* If no pixels were retained, and not in force mode, then do no more */
    if(nr==0) {
      if(forceMode) {
        fprintf(stderr, "Warning: no pixels retained.\n");
      } else {
        fprintf(stderr, "Note: no pixels retained, image not saved.\n");
        imgEmpty(&img);
        return(0);
      }
    }
 
  }
 
 
  /*
   *  Write the image
   */
  if(verbose>1) fprintf(stdout, "writing %s\n", outfile);
  if(imgWrite(outfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg); 
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img);
    return(11);
  }
  if(verbose>0) printf("Written %s\n", outfile);
  imgEmpty(&img);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/