imgthbcv.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 "libtpcmodel.h"
#include "libtpcmodext.h"
/*****************************************************************************/
/* Local functions */
double bcvfunc(double p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Based on given absolute threshold, calculate between-class variance from",
  "static PET image file in ECAT, NIfTI, or Analyze format.",
  "If threshold is not given, threshold that gives the highest between-class",
  "variance is estimated.",
  "Between-class variance is calculated using formula by Cao et al. 2018",
  "doi: 10.1109/ACCESS.2018.2889013",
  " ",
  "Usage: @P [Options] imgfile [threshold]",
  " ",
  "Options:",
  " -mask=<filename>",
  "     Save mask file containing pixel values 1 (above threshold) and",
  "     0 (under the threshold).",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imgthrs, imgmask, imgqntls, imginteg, imgslim, imgcutof",
  " ",
  "Keywords: image, threshold",
  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;
/*****************************************************************************/
 
/*****************************************************************************/

double imgBetween2ClassesVariance(
  IMG *img,
  const int fi,
  float threshold,
  int verbose
) {
  if(verbose>0) {printf("%s(IMG, %d, %f, ...)\n", __func__, fi, threshold); fflush(stdout);}
  double bcv=nan("");
  if(img==NULL || img->dimz<1 || img->dimy<1 || img->dimx<1 || img->dimt<1) return(bcv);
  if(fi<1 || fi>img->dimt) return(bcv);
 
  if(verbose>1) fprintf(stdout, "thresholding\n");
  unsigned int p0=0, p1=0;
  double s0=0.0, s1=0.0;
  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(img->m[zi][yi][xi][fi-1]<threshold) {
          p0++; s0+=img->m[zi][yi][xi][fi-1];
        } else if(img->m[zi][yi][xi][fi-1]>=threshold) {
          p1++; s1+=img->m[zi][yi][xi][fi-1];
        } // NaNs are left out with if else if
      }
    }
  }
  if(verbose>1) {
    printf("p0=%u\n", p0);
    printf("p1=%u\n", p1);
  }
  /* If no pixels in either of classes, between-class variance is 0 */
  if(p0==0 || p1==0) {
    bcv=0.0;
  } else { // otherwise calculate means, and the between-class variance
 
    /* Calculate means */
    double mu0, mu1, mu;
    mu0=s0/(double)p0;
    mu1=s1/(double)p1;
    mu=(s0+s1)/(double)(p0+p1);
    if(verbose>1) {
      printf("mu0=%g\n", mu0);
      printf("mu1=%g\n", mu1);
      printf("mu=%g\n", mu);
    }
 
    /* Calculate between-class variance */
    bcv=(double)p0*(double)p1*( (mu0-mu1)*(mu0-mu1) + (mu0-mu)*(mu0-mu) + (mu1-mu)*(mu1-mu) );
  }
 
  if(verbose>1) {printf("BCV := %g\n", bcv); fflush(stdout);}
  return(bcv);
}
/*****************************************************************************/
 
/*****************************************************************************/
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char     petfile[FILENAME_MAX], maskfile[FILENAME_MAX];
  float    threshold=nanf("");
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=maskfile[0]=(char)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, "MASK=", 5)==0) {
      strlcpy(maskfile, cptr+5, FILENAME_MAX); 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;
    if(atof_with_check(argv[ai], &v)!=0) {
      fprintf(stderr, "Error: invalid threshold '%s'.\n", argv[ai]);
      return(1);
    }
    threshold=v;
    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);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    if(!isnan(threshold)) printf("threshold := %g\n", threshold);
    if(maskfile[0]) printf("maskfile := %s\n", maskfile);
    fflush(stdout);
  }
 
 
  /*
   *  Read the contents of the PET file to img data structure
   */
  IMG img; imgInit(&img);
  if(verbose>1) {printf("reading %s\n", petfile); fflush(stdout);}
  if(imgRead(petfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    imgEmpty(&img); return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  /* Check that image is static */
  if(img.dimt>1) {
    fprintf(stderr, "Error: do not use for dynamic image.\n");
    imgEmpty(&img); return(2);
  }
  /* Get the range of pixel values in image */
  float imin, imax;
  if(imgMinMax(&img, &imin, &imax)) {
    fprintf(stderr, "Error: invalid image contents.\n");
    imgEmpty(&img); return(3);
  }
  if(verbose>1) {
    printf("image_min := %g\n", imin);
    printf("image_max := %g\n", imax);
    fflush(stdout);
  }
 
  double bcv=0.0;
  if(!isnan(threshold)) {
    if(threshold<=imin || threshold>=imax) {
      fprintf(stdout, "Image pixel value range: %g - %g\n", imin, imax); fflush(stdout);
      fprintf(stderr, "Error: invalid threshold %g.\n", threshold);
      imgEmpty(&img); return(4);
    }
    /* Calculate between-class variance, print it, and quit */
    bcv=imgBetween2ClassesVariance(&img, 1, threshold, verbose-2);
  } else {
    /* Find the threshold that gives the largest between-class variance */
    if(verbose>1) {printf("finding optimal threshold...\n"); fflush(stdout);}
    double p, fv;
    int ret=nlopt1D(bcvfunc, &img, 0.5*(imin+imax), imin, imax, 0.1*(imax-imin),
                    0.000001*(imax-imin), 1000, &p, &fv, verbose-2);
    if(ret!=0) {
      fprintf(stderr, "Error: invalid image contents.\n");
      imgEmpty(&img); return(4);
    }
    threshold=p;
    bcv=-fv;
    printf("threshold := %g\n", threshold);
    fflush(stdout);
  }
  printf("BCV := %g\n", bcv);
  fflush(stdout);
 
 
  /* If requested, make a mask file based on threshold */
  if(maskfile[0]) {
    IMG mask; imgInit(&mask);
    int ret=imgThresholdMask(&img, threshold, imax+1.0, &mask);
    if(ret) {
      fprintf(stderr, "Error: cannot threshold.\n");
      if(verbose>1) printf("ret=%d\n", ret);
      imgEmpty(&img); imgEmpty(&mask);
      return(6);
    }
    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);
    imgEmpty(&mask);
  }
 
 
  imgEmpty(&img);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double bcvfunc(double p, void *fdata)
{
  return(-imgBetween2ClassesVariance(fdata, 1, p, 0));
}
/*****************************************************************************/
 
/*****************************************************************************/