imgthrs.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[] = {
  "Threshold PET image file in ECAT, NIfTI, or Analyze format.",
  "Program calculates an integral (AUC) over time frames of a dynamic image,",
  "finds the maximal integral value of all planes, and cuts off (sets to zero)",
  "the pixels in dynamic image which have lower or higher integral than",
  "the specified threshold-%.",
  "Original image file is not modified.",
  " ",
  "Usage: @P [Options] imgfile lowerthreshold upperthreshold [outputfile]",
  " ",
  "Options:",
  " -mask=<filename>",
  "     Save cutoff mask file containing pixels values 1 (between threshold",
  "     limits) and 0 (under lower or over upper threshold).",
  " -abs",
  "     Thresholds are given as absolute values instead of percentages",
  "     of maximum in the static image; not applicable to dynamic images.",
  " -start=<time (min)>",
  "     AUC calculation starts from given time. By default from scan start.",
  " -end=<time (min)>",
  "     AUC calculation ends to given time. By default until scan end.",
  " --force",
  "     Program does not mind if none of pixels are thresholded;",
  "     otherwise that is considered as an error.",
  " -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.",
  " ",
  "Example: threshold the background and cerebrospinal fluid:",
  "     @P b123dy1.v 30 100 b123thres.v",
  " ",
  "See also: img2dft, imgmask, imgqntls, imginteg, imgposv, imgslim, imgcutof",
  " ",
  "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], maskfile[FILENAME_MAX], outfile[FILENAME_MAX];
  float    lowerThreshold=-1.0, upperThreshold=-1.0;
  int      absThreshold=0; // 0=percentages; 1=values
  int      forceMode=0; // 0=off; 1=on
  char    *cptr;
  int      ret=0;
  IMG      img, sum, mask;
  float    tstart, tstop;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=outfile[0]=maskfile[0]=(char)0;
  tstart=tstop=-1.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, "MASK=", 5)==0) {
      strlcpy(maskfile, cptr+5, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "ABSOLUTE", 3)==0) {
      absThreshold=1; continue;
    } else if(strncasecmp(cptr, "START=", 6)==0) {
      double v;
      if(atof_with_check(cptr+6, &v)==0) {
        if(v>0.0) tstart=60.0*v; else tstart=0.0;
        continue;
      }
    } else if(strncasecmp(cptr, "END=", 4)==0) {
      double v;
      if(atof_with_check(cptr+4, &v)==0) {
        if(v>0.0) {tstop=60.0*v; continue;}
      }
    } else if(strncasecmp(cptr, "STOP=", 5)==0) {
      double v;
      if(atof_with_check(cptr+5, &v)==0) {
        if(v>0.0) {tstop=60.0*v; 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;
    if(absThreshold==0) lowerThreshold*=0.01;
    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;
    if(absThreshold==0) {
      if(upperThreshold<=0.0) {
        fprintf(stderr, "Error: invalid upper threshold '%s'.\n", argv[ai]);
        return(1);
      }
      upperThreshold*=0.01;
    }
    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 threshold range.\n");
    if(verbose>0) {
      if(absThreshold==0) 
        printf("Threshold range: %g - %g\n", 100.*lowerThreshold, 100.*upperThreshold);
      else
        printf("Threshold range: %g - %g\n", lowerThreshold, upperThreshold);
    }
    return(1);
  }
  if(lowerThreshold==upperThreshold) {
    /* Image data is stored with limited precision, therefore limits should not be equal */
    float f;
    if(absThreshold==0) {
      f=0.001*lowerThreshold; lowerThreshold-=f; upperThreshold+=f;
      if(verbose>1) printf("Threshold range: %g - %g\n", 100.*lowerThreshold, 100.*upperThreshold);
    } else {
      f=0.001*lowerThreshold; if(f>0.1) f=0.1;
      lowerThreshold-=f; upperThreshold+=f;
      if(verbose>1) printf("Threshold range: %g - %g\n", lowerThreshold, upperThreshold);
      if(verbose>10) printf("Threshold range: %.2f - %.2f\n", lowerThreshold, upperThreshold);
    }
  }
 
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    if(outfile[0]) printf("outfile := %s\n", outfile);
    if(maskfile[0]) printf("maskfile := %s\n", maskfile);
    if(absThreshold==0) {
      printf("upper_threshold_percentage := %g\n", 100.*upperThreshold);
      printf("lower_threshold_percentage := %g\n", 100.*lowerThreshold);
    } else {
      printf("upper_threshold := %g\n", upperThreshold);
      printf("lower_threshold := %g\n", lowerThreshold);
    }
    if(tstart>=0) printf("start_time := %g\n", tstart/60.0);
    if(tstop>=0) printf("end_time := %g\n", tstop/60.0);
  }
 
 
  /*
   *  Read the contents of the PET file to img data structure
   */
  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);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  /* Absolute thresholds can be used for static image only */
  if(absThreshold!=0 && img.dimt>1) {
    fprintf(stderr, "Error: do not use absolute thresholds for dynamic image.\n");
    imgEmpty(&img); return(2);
  }
  /* Check that frame times are available */
  if(!imgExistentTimes(&img)) {
    if(tstart>0.0 || tstop>0.0) {
      fprintf(stderr, "Error: image does not contain frame times.\n");
      imgEmpty(&img); return(2);
    }
    if(img.dimt>1)
      fprintf(stderr, "Warning: image does not contain frame times.\n");
  }
  /* For raw data, divide counts by frame duration */
  ret=imgRawCountsPerTime(&img, 1);
  if(ret) {
    fprintf(stderr, "Error: cannot correct sinogram counts.\n");
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); return(2);
  }
 
 
  /* 
   *  Calculate the integral image from dynamic image
   */
  if(img.dimt>1) {
    if(verbose>0) fprintf(stdout, "calculating integral image\n");
    imgInit(&sum);
    if(tstart<=0.0 && tstop<=0.0) {
      /* Time range not specified, thus just sum over all frames */
      ret=imgFrameIntegral(&img, 0, img.dimt-1, &sum, verbose-4);
    } else {
      /* Check that reasonable time range was given by user */
      if(tstart<img.start[0]) tstart=img.start[0];
      if(tstop<0.01 || tstop>img.end[img.dimt-1]) tstop=img.end[img.dimt-1];
      if(verbose>1) {
        printf("tstart := %g\n", tstart);
        printf("tstop := %g\n", tstop);
      }
      ret=imgTimeIntegral(&img, tstart, tstop, &sum, 0, NULL, verbose-4);
    }
    if(ret) {
      fprintf(stderr, "Error: cannot integrate.\n");
      if(verbose>1) printf("ret=%d\n", ret);
      imgEmpty(&img); imgEmpty(&sum); return(4);
    }
  } else {
    /* If image is static, then just copy it. Do _not_ multiply with
       frame duration, because it would hamper thresholding with absolute
       values, would not work if frame times are zero, and 
       would not help in thresholding with percentages of max. */
    if(verbose>2) fprintf(stdout, "static image\n");
    imgInit(&sum);
    if(imgDup(&img, &sum)!=0) {
      fprintf(stderr, "Error: cannot duplicate static image.\n");
      imgEmpty(&img); imgEmpty(&sum); return(4);
    }
  }
 
  /*
   *  Thresholding
   */
  if(verbose>0) fprintf(stdout, "thresholding\n");
 
  if(absThreshold==0) {
    /* Search for maximum pixel value */
    if(verbose>1) printf("searching maximum pixel value...\n");
    float max;
    ret=imgMax(&sum, &max);
    if(ret) {
      fprintf(stderr, "Error: cannot find maximum.\n");
      if(verbose>1) printf("ret=%d\n", ret);
      imgEmpty(&img); imgEmpty(&sum);
      return(5);
    }
    if(verbose>1) printf("max := %g\n", max);
    lowerThreshold*=max; upperThreshold*=max;
    if(verbose>1)
      printf("lowerThresholdValue := %g\nupperThresholdValue:=%g\n",
        lowerThreshold, upperThreshold);
  }
  
  /* Make a mask image */
  imgInit(&mask);
  ret=imgThresholdMask(&sum, lowerThreshold, upperThreshold, &mask);
  if(ret) {
    fprintf(stderr, "Error: cannot threshold.\n");
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); imgEmpty(&sum); imgEmpty(&mask);
    return(6);
  }
  /* Sum image is not needed any more */
  imgEmpty(&sum);
 
  /* Report the number of retained pixels, in verbose mode in each plane */
  {
    int zi, xi, yi;
    long long nr, total_nr=0, thrs_nr;
    if(verbose>1 && mask.dimz>1) printf("Nr of retained pixels:\n");
    for(zi=0; zi<mask.dimz; zi++) {
      nr=0;
      for(yi=0; yi<mask.dimy; yi++)
        for(xi=0; xi<mask.dimx; xi++)
          if(fabs(mask.m[zi][yi][xi][0])>0.01) nr++;
      if(verbose>1 && mask.dimz>1) printf("  plane %d: %lld\n", zi+1, nr);
      total_nr+=nr;
    }
    if(verbose>0 || !outfile[0]) printf("retained_pixels := %lld\n", total_nr);
    thrs_nr=mask.dimz*mask.dimy*mask.dimx-total_nr;
    if(verbose>0 || !outfile[0]) printf("thresholded_pixels := %lld\n", thrs_nr);
    /* If no pixels were thresholded, and not in force mode, then do no more */
    if(thrs_nr==0) {
      if(forceMode) {
        fprintf(stderr, "Warning: no pixels thresholded.\n");
      } else {
        fprintf(stderr, "Error: no pixels thresholded.\n");
        imgEmpty(&img); imgEmpty(&mask);
        return(7);
      }
    }
  }
 
  /* Write the mask image, if requested */
  if(maskfile[0]) {
    if(verbose>1) fprintf(stdout, "writing %s\n", maskfile);
    ret=imgWrite(maskfile, &mask);
    if(ret) {
      fprintf(stderr, "Error: %s\n", img.statmsg); 
      if(verbose>1) printf("ret=%d\n", ret);
      imgEmpty(&img); imgEmpty(&mask);
      return(11);
    }
    if(verbose>0) printf("Written %s\n", maskfile);
  }
  
  /* If thresholded image was not requested, then quit */
  if(!outfile[0]) {
    if(verbose>0) printf("done.\n");
    imgEmpty(&img); imgEmpty(&mask);
    return(0);
  }
 
  /* According to the mask, threshold the original image */
  ret=imgThresholdByMask(&img, &mask, 0.0);
  if(ret) {
    fprintf(stderr, "Error: cannot threshold by the mask.\n");
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); imgEmpty(&mask);
    return(8);
  }
  imgEmpty(&mask);
 
 
  /*
   *  Write the thresholded image
   */
  /* For raw data, multiply counts by frame duration */
  ret=imgRawCountsPerTime(&img, 0);
  if(ret) {
    fprintf(stderr, "Error: cannot correct sinogram counts.\n");
    imgEmpty(&img); return(9);
  }
  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(13);
  }
  if(verbose>0) printf("Written %s\n", outfile);
  imgEmpty(&img);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/