imgqntls.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 "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcsvg.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Divide the voxel values in static PET image into quantiles.",
  "PET image file must be in ECAT, NIfTI, or Analyze format.",
  "User must specify the number of groups, and program writes the group",
  "cut points (quantiles, nr of groups - 1) into stdout.",
  " ",
  "Usage: @P [Options] petfile groupnr",
  " ",
  "Options:",
  " -pos[itive]",
  "     Include only voxels with value > 0.",
  " -abs[olute]",
  "     Use absolute voxel values.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imghist, imginteg, imgpeak, imgthrs, imgcutof, imgmask",
  " ",
  "Keywords: image, histogram, 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;
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
static int statFloatCompAsc(const void *i, const void *j)
{
  const float *di = (const float *)i;
  const float *dj = (const float *)j;
  return(*di > *dj) - (*di < *dj);
}
static int statFloatCompDesc(const void *i, const void *j)
{
  const float *di = (const float *)i;
  const float *dj = (const float *)j;
  return(*di < *dj) - (*di > *dj);
}
 
void statSortFloat(
  float *data,
  unsigned int n,
  int order
) {
  if(n<2 || data==NULL) return;
  if(order==0) qsort(data, n, sizeof(float), statFloatCompAsc);
  else qsort(data, n, sizeof(float), statFloatCompDesc);
}
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char     petfile[FILENAME_MAX];
  int      groupNr=0;
  char    *cptr;
  int      ret=0;
  int      positives=0;
  int      absolutes=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[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, "POSITIVES", 3)==0) {
      positives=1; continue;
    } else if(strncasecmp(cptr, "ABSOLUTES", 3)==0) {
      absolutes=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);
  if(ai<argc) groupNr=atoi(argv[ai++]);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!petfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(groupNr<2) {
    fprintf(stderr, "Error: invalid nr of groups.\n");
    return(1);
  }
  if(absolutes && positives) {
    fprintf(stderr, "Error: invalid combination of options -abs and -pos.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("groupNr := %d\n", groupNr);
    printf("positives_only := %d\n", positives);
    printf("absolute_values := %d\n", absolutes);
    fflush(stdout);
  }
 
 
 
  /*
   *  Read PET data
   */
  if(verbose>0) fprintf(stdout, "reading %s\n", petfile);
  IMG img; imgInit(&img);
  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(verbose>1) {
    printf("image_size := %d x %d x %d\n", img.dimx, img.dimy, img.dimz);
    if(img.dimt>1) printf("frame_nr := %d\n", img.dimt);
  }
  if(img.dimt>1) {
    fprintf(stderr, "Error: cannot process dynamic images.\n");
    imgEmpty(&img); return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
 
  /* Copy the pixel values into an array */
  if(verbose>1) fprintf(stdout, "copying voxel values\n");
  long long pxlNr=img.dimx*img.dimy*img.dimz;
  float *pxls=(float*)malloc(pxlNr*sizeof(float));
  if(pxls==NULL) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    imgEmpty(&img); return(3);
  }
  pxlNr=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(positives && img.m[zi][yi][xi][0]<1.0E-010) continue;
        if(absolutes==0) pxls[pxlNr]=img.m[zi][yi][xi][0];
        else pxls[pxlNr]=fabsf(img.m[zi][yi][xi][0]);
        pxlNr++;
      }
  if(verbose>1) printf("included_pixel_nr := %lld\n", pxlNr);
  if(pxlNr<=groupNr) {
    if(positives) fprintf(stderr, "Error: no positive voxel values.\n");
    else fprintf(stderr, "Error: more groups than voxels.\n");
    imgEmpty(&img); free(pxls);
    return(4);
  }
  /* No need for original image */
  imgEmpty(&img);
 
  /* Sort pixel values */
  if(verbose>0) fprintf(stdout, "sorting voxel values\n");
  statSortFloat(pxls, pxlNr, 1);
 
  /* Determine the quantiles */
  int nrPerGroup=pxlNr/groupNr;
  if(verbose>1) printf("pixels_per_group := %d\n", nrPerGroup);
  float pxlQuantiles[groupNr-1];
  for(int i=1; i<groupNr; i++) {
    pxlQuantiles[i-1]=0.5*(pxls[i*nrPerGroup]+pxls[i*nrPerGroup - 1]);
    if(verbose>10) printf("quantile[%d] := %g\n", i, pxlQuantiles[i-1]);
  }
  /* No need for pixel value list */
  free(pxls);
 
  /* Print quantiles in stdout */
  fflush(stdout);
  for(int i=1; i<groupNr; i++) {
    fprintf(stdout, "quantile[%d] := %g\n", i, pxlQuantiles[i-1]);
  }
  fflush(stdout);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/