imghist.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[] = {
  "Make histogram of the voxel values in PET image.",
  "PET image file must be in ECAT, NIfTI, or Analyze format.",
  "Resulting histogram file contains the bin as the two first columns,",
  "and the proportion (0-1) of values in the bin in the next column.",
  "In case of dynamic PET data a separate histogram is made from each frame",
  "and saved as its own column in the histogram file.",
  " ",
  "Usage: @P [Options] petfile [histogram]",
  " ",
  "Options:",
  " -bin=<Binsize as absolute value> or -per=<Binsize as percentage>",
  "     Set the bin size using these options; by default, bin size is 10%",
  "     of maximum-minimum voxel value in the whole image.",
  " -nr",
  "     Nr of voxels in each bin is written in histogram file; proportion",
  "     of them is written by default.",
  " -pos[itive]",
  "     Include only voxels with value > 0.",
  " -svg=<filename>",
  "     Histogram is plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imginteg, imgmax, imgposv, imglkup, hist4dat, imgqntls, imgthrs",
  " ",
  "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;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char     petfile[FILENAME_MAX], svgfile[FILENAME_MAX], hisfile[FILENAME_MAX];
  char    *cptr;
  int      ret=0;
  int      report_proportions=1;
  int      positives=0;
  double   bin_size=-1.0, bin_perc=-1.0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=hisfile[0]=svgfile[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, "SVG=", 4)==0) {
      if(strlcpy(svgfile, cptr+4, FILENAME_MAX)>1) continue;
    } else if(strncasecmp(cptr, "BIN=", 4)==0) {
      ret=atof_with_check(cptr+4, &bin_size);
      if(ret==0 && bin_size>0.0) continue;
    } else if(strncasecmp(cptr, "PER=", 4)==0) {
      ret=atof_with_check(cptr+4, &bin_perc);
      if(ret==0 && bin_perc>0.0) continue;
    } else if(strcasecmp(cptr, "NR")==0) {
      report_proportions=0; continue;
    } else if(strncasecmp(cptr, "POSITIVES", 3)==0) {
      positives=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) strlcpy(hisfile, argv[ai++], FILENAME_MAX);
  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(bin_size>0.0 && bin_perc>0.0) {
    fprintf(stderr, "Error: options -bin and -per cannot be used together.\n");
    return(1);
  }
  if(bin_size<=0.0 && bin_perc<=0.0) bin_perc=10.0;
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    if(hisfile[0]) printf("hisfile := %s\n", hisfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    if(bin_size>0.0) printf("bin_size := %g\n", bin_size);
    if(bin_perc>0.0) printf("bin_percentage := %g\n", bin_perc);
    printf("positives_only := %d\n", positives);
  }
 
 
 
  /*
   *  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(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
 
  /* Find the minimum and maximum voxel value in the whole PET data */
  float maxv, minv;
  imgMinMax(&img, &minv, &maxv);
  if(verbose>1) printf("min := %g\nmax := %g\n", minv, maxv);
  /* Refine the minimum */
  if(positives) {
    if(minv<0.0) minv=0.0;
  } else {
    /* If minimum is just a little over zero, then set it to zero */
    if(minv>0.0 && minv<0.05*maxv)
      minv=0.0;
  }
  if(verbose>1) printf("refined_min := %g\n", minv);
  /* Calculate the value range */
  double rangev=(double)maxv-(double)minv;
  if(verbose>2) printf("max-min := %g\n", rangev);
 
  /* Set bin size based on bin percentage, if necessary */
  if(bin_size<=0.0) {
    bin_size=0.01*bin_perc*rangev;
    if(verbose>1) printf("bin_size := %g\n", bin_size);
  }
  if(bin_size<1.0E-020) {
    fprintf(stderr, "Error: invalid bin size.\n");
    imgEmpty(&img); return(2);
  }
 
  /* How many bins we may need? */
  int bin_nr=(int)ceil(rangev/bin_size);
  if(verbose>1) printf("estimated_bin_nr := %d\n", bin_nr);
  if(bin_nr<2) {
    fprintf(stderr, "Error: invalid number of bins.\n");
    imgEmpty(&img); return(2);
  } else if(bin_nr>10000) {
    fprintf(stderr, "Error: too many bins.\n");
    imgEmpty(&img); return(2);
  }
 
  /* Allocate memory for histogram data */
  DFT hist; dftInit(&hist);
  ret=dftSetmem(&hist, bin_nr+1, img.dimt);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot allocate memory for the histogram.\n");
    imgEmpty(&img); return(3);
  }
  hist.frameNr=bin_nr+1; 
  hist.voiNr=img.dimt;
  hist.timetype=DFT_TIME_STARTEND;
  hist._type=DFT_FORMAT_PMOD;
  /* Set bins as TAC frame times */
  /* Bins should be in calculated from zero, even if bins closest to zero
     are not actually used */
  {
    int n=0;
    double x=0.0;
    if(minv<0.0) {while(x>minv) x-=bin_size;}
    else if(minv>0.0) {while(x+bin_size<minv) x+=bin_size;}
    //printf("n=%d x=%g\n", n, x);
    hist.x1[n]=x; hist.x2[n]=hist.x1[n]+bin_size;
    n++; x+=bin_size; 
    while(x<=maxv && n<hist.frameNr) {
      //printf("n=%d x=%g\n", n, x);
      hist.x1[n]=x; hist.x2[n]=hist.x1[n]+bin_size;
      x+=bin_size; 
      n++;
    }
    hist.frameNr=n; 
  }
  for(int i=0; i<hist.frameNr; i++) hist.x[i]=0.5*(hist.x1[i]+hist.x2[i]);
 
  /* Set TAC column name(s) */
  if(hist.voiNr==1) {
    strcpy(hist.voi[0].name, "static");
  } else {
    for(int i=0; i<hist.voiNr; i++) sprintf(hist.voi[i].name, "fr%d", 1+i);
  }
 
 
  /*
   *  Go through the PET data and bin voxel values
   */
  if(verbose>1) printf("computing histogram...\n");
  long long voxNr=0;
  int zi, yi, xi, ti, bi;
  for(zi=0; zi<img.dimz; zi++)
    for(yi=0; yi<img.dimy; yi++)
      for(xi=0; xi<img.dimx; xi++)
        for(ti=0; ti<img.dimt; ti++) {
          if(positives && img.m[zi][yi][xi][ti]<1.0E-010) continue;
          voxNr++;
          for(bi=bin_nr-1; bi>0; bi--) {
            if(img.m[zi][yi][xi][ti]>hist.x1[bi]) {
              hist.voi[ti].y[bi]+=1.0; break;
            }
          }
          if(bi==0) hist.voi[ti].y[0]+=1.0;
        }
  if(verbose>1) printf("included_pixel_nr := %lld\n", voxNr);
  if(voxNr<1) {
    fprintf(stderr, "Error: no positive voxel values.\n");
    imgEmpty(&img); dftEmpty(&hist); 
    return(4);
  }
  /* Calculate proportion of pixels in each bin, if requested */
  if(report_proportions) {
    double c=1.0/(double)voxNr;
    for(ti=0; ti<hist.voiNr; ti++)
      for(bi=0; bi<hist.frameNr; bi++)
        hist.voi[ti].y[bi]*=c;
  }
 
  /* PET data is needed no more */
  imgEmpty(&img);
 
 
  /*
   *  Save or print the histogram, if required
   */
  if(!hisfile[0] && !svgfile[0]) {
    printf("\nHistogram:\n"); fflush(stdout);
    dftPrint(&hist);
    dftEmpty(&hist); return(0);
  }
  if(hisfile[0]) {
    if(verbose>1) printf("writing %s\n", hisfile);
    if(dftWrite(&hist, hisfile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", hisfile, dfterrmsg);
      dftEmpty(&hist); return(11);
    }
    if(verbose>0) printf("Histogram written in %s\n", hisfile);
  }
 
  /*
   *  Saving and/or plotting of fitted TACs
   */
  if(svgfile[0]) {
    if(verbose>1) printf("saving SVG plot\n");
 
    /* Create a second TAC struct for plotting lines */
    DFT hist2; dftInit(&hist2); ret=dftdup(&hist, &hist2);
    if(ret) {
      fprintf(stderr, "Error: cannot plot histogram.\n");
      dftEmpty(&hist); dftEmpty(&hist2);
      return(21);
    }
 
    /* Save SVG plot */
    ret=plot_fit_svg(&hist, &hist2, "", svgfile, verbose-8);
    if(ret) {
      fprintf(stderr, "Error: cannot plot the histogram.\n");
      dftEmpty(&hist); dftEmpty(&hist2);
      return(22);
    }
    if(verbose>0) printf("Histogram plotted in %s\n", svgfile);
 
    dftEmpty(&hist2);
  }
 
  dftEmpty(&hist);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/