fslstats.cc

/*  fslstats.cc

    Mark Jenkinson and Matthew Webster, FMRIB Image Analysis Group

    Copyright (C) 2003-2007 University of Oxford  */

/*  CCOPYRIGHT  */

#include "miscmaths/miscmaths.h"
#include "newimage/newimageall.h"
#include "newimage/costfns.h"

using namespace NEWIMAGE;
using namespace MISCMATHS;

bool masks_used=false;
bool lthr_used=false;
bool uthr_used=false;

#if defined ( __CYGWIN__ ) ||  defined (__sun)
extern "C" { 
#include <ieeefp.h> //for finite
}
#endif

void print_usage(const string& progname) {
  cout << "Usage: fslstats <input> [options]" << endl << endl;
  cout << endl;
  cout << "-l <lthresh> : set lower threshold" << endl;
  cout << "-u <uthresh> : set upper threshold" << endl;
  cout << "-r           : output <robust min intensity> <robust max intensity>" << endl;
  cout << "-R           : output <min intensity> <max intensity>" << endl;
  cout << "-e           : output mean entropy ; mean(-i*ln(i))" << endl;
  cout << "-E           : output mean entropy (of nonzero voxels)" << endl;
  cout << "-v           : output <voxels> <volume>" << endl;
  cout << "-V           : output <voxels> <volume> (for nonzero voxels)" << endl;
  cout << "-m           : output mean" << endl;
  cout << "-M           : output mean (for nonzero voxels)" << endl;
  cout << "-s           : output standard deviation" << endl;
  cout << "-S           : output standard deviation (for nonzero voxels)" << endl;
  cout << "-w           : output smallest ROI <xmin> <xsize> <ymin> <ysize> <zmin> <zsize> <tmin> <tsize> containing nonzero voxels" << endl;
  cout << "-x           : output co-ordinates of maximum voxel" << endl;
  cout << "-X           : output co-ordinates of minimum voxel" << endl;
  cout << "-c           : output centre-of-gravity (cog) in mm coordinates" << endl;
  cout << "-C           : output centre-of-gravity (cog) in voxel coordinates" << endl;
  cout << "-p <n>       : output nth percentile (n between 0 and 100)" << endl;
  cout << "-P <n>       : output nth percentile (for nonzero voxels)" << endl;
  cout << "-a           : use absolute values of all image intensities"<< endl;
  cout << "-n           : treat NaN or Inf as zero for subsequent stats" << endl;
  cout << "-k <mask>    : use the specified image (filename) for masking - overrides lower and upper thresholds" << endl;
  cout << "-h <nbins>   : output a histogram (for the thresholded/masked voxels only) with nbins" << endl; 
  cout << "-H <nbins> <min> <max>   : output a histogram (for the thresholded/masked voxels only) with nbins and histogram limits of min and max" << endl; 
  cout << endl;
  cout << "Note - thresholds are not inclusive ie lthresh<allowed<uthresh" << endl;
}


// Some specialised nonzero functions just for speedup
//  (it avoids generating masks when not absolutely necessary)

long int nonzerocount(const volume4D<float>& vol)
{
  long int totn=0;
  for (int t=vol.mint(); t<=vol.maxt(); t++) {
    for (int z=vol.minz(); z<=vol.maxz(); z++) {
      for (int y=vol.miny(); y<=vol.maxy(); y++) {
	for (int x=vol.minx(); x<=vol.maxx(); x++) {
	  if (vol(x,y,z,t)!=0.0) {
	    totn++;
	  }
	}
      }
    }
  }
  return totn;
}

double nonzeromean(const volume4D<float>& vol)
{
  double totv=0.0;
  long int totn=0;
  for (int t=vol.mint(); t<=vol.maxt(); t++) {
    for (int z=vol.minz(); z<=vol.maxz(); z++) {
      for (int y=vol.miny(); y<=vol.maxy(); y++) {
	for (int x=vol.minx(); x<=vol.maxx(); x++) {
	  if (vol(x,y,z,t)!=0.0) {
	    totv+=(double) vol(x,y,z,t);
	    totn++;
	  }
	}
      }
    }
  }
  double meanval=0.0;
  if (totn>0) {
    meanval=totv/totn;
  }
  return meanval;
}

double nonzerostddev(const volume4D<float>& vol)
{
  double totv=0.0, totvv=0.0;
  long int totn=0;
  for (int t=vol.mint(); t<=vol.maxt(); t++) {
    for (int z=vol.minz(); z<=vol.maxz(); z++) {
      for (int y=vol.miny(); y<=vol.maxy(); y++) {
	for (int x=vol.minx(); x<=vol.maxx(); x++) {
	  if (vol(x,y,z,t)!=0.0) {
	    float v=vol(x,y,z,t);
	    totvv+=(double) v*v;
	    totv+=(double) v;
	    totn++;
	  }
	}
      }
    }
  }
  double var=0.0;
  if (totn>1) {
    double meanval = totv/totn;
    var = (totvv - totn*meanval*meanval)/(totn-1);
  }
  return std::sqrt(var);
}

int generate_masks(volume4D<float> &mask, volume4D<float> &masknz, const volume4D<float> &vin,
		   float& lthr, float& uthr) 
{
  if (!lthr_used) { lthr=vin.min()-1; }
  if (!uthr_used) { uthr=vin.max()+1; }
  mask = binarise(vin,lthr,uthr,exclusive);
  if (mask.tsize()!=1) masknz = (binarise(vin,0.0f, 0.0f)-1.0f)*-1.0f*mask;
  else masknz = (binarise(vin,0.0f, 0.0f)-1.0f)*-1.0f*mask[0];
  return 0;
}

int generate_masks(const volume4D<float> &mask, volume4D<float> &masknz, const volume4D<float> &vin)
{
  if (mask.tsize()!=1) masknz = (binarise(vin,0.0f, 0.0f)-1.0f)*-1.0f*mask;
  else masknz = (binarise(vin,0.0f, 0.0f)-1.0f)*-1.0f*mask[0];
  return 0;
}


int fmrib_main_float(int argc, char* argv[]) 
{

  cout.setf(ios::dec); 
  cout.setf(ios::fixed, ios::floatfield); 
  cout.setf(ios::left, ios::adjustfield); 
  cout.precision(6);  

  volume4D<float> vin, vol, mask, masknz;
  read_volume4D(vol,argv[1]);

  float lthr=0, uthr=0;  // these initial values are not used
  if (masks_used) {
    vin = vol;
    generate_masks(mask,masknz,vin,lthr,uthr);
    vol = vin * mask;
  }

  int narg=2;
  string sarg;
 
  while (narg<argc) {
    sarg=argv[narg];
    if (sarg=="-n") {
      for (int t=vol.mint(); t<=vol.maxt(); t++)
        for (int z=vol.minz(); z<=vol.maxz(); z++)
          for (int y=vol.miny(); y<=vol.maxy(); y++)
            for (int x=vol.minx(); x<=vol.maxx(); x++)
              if (!finite((double)vol(x,y,z,t)))
	        vol(x,y,z,t)=0;
    } else if (sarg=="-m") {
      if (masks_used) cout <<  vol.mean(mask) << " ";
      else cout << vol.mean() << " ";
    } else if (sarg=="-M") {
      if (masks_used) cout << vol.mean(masknz) << " ";
      else {
	double nzmean=0;
	nzmean = nonzeromean(vol);
	cout << nzmean << " ";
      }
    } else if (sarg=="-X") {
      ColumnVector coord(4);
      coord(4)=1.0;
      if (masks_used) {
	coord(1) = vol.mincoordx(mask);
	coord(2) = vol.mincoordy(mask);
	coord(3) = vol.mincoordz(mask);
      } else {
	coord(1) = vol.mincoordx();
	coord(2) = vol.mincoordy();
	coord(3) = vol.mincoordz();
      }
      coord = (vol[0].niftivox2newimagevox_mat()).i() * coord;
      cout << MISCMATHS::round(coord(1)) << " " << 
	MISCMATHS::round(coord(2)) << " " << MISCMATHS::round(coord(3)) << " ";
    } else if (sarg=="-x") { 
      ColumnVector coord(4);
      coord(4)=1.0;
      if (masks_used) {
	coord(1) = vol.maxcoordx(mask);
	coord(2) = vol.maxcoordy(mask);
	coord(3) = vol.maxcoordz(mask);
      } else {
	coord(1) = vol.maxcoordx();
	coord(2) = vol.maxcoordy();
	coord(3) = vol.maxcoordz();
      }
      coord = (vol[0].niftivox2newimagevox_mat()).i() * coord;
      cout << MISCMATHS::round(coord(1)) << " " << 
	MISCMATHS::round(coord(2)) << " " << MISCMATHS::round(coord(3)) << " ";
    } else if (sarg=="-w") {
      if (!masks_used) { 
	if (vin.nvoxels()<1) { vin = vol; }
	masks_used=true;
	generate_masks(mask,masknz,vin,lthr,uthr); 
	vol = vin * mask; 
      }
      int xmin=masknz.maxx(),xmax=masknz.minx(),ymin=masknz.maxy(),ymax=masknz.miny(),zmin=masknz.maxz(),zmax=masknz.minz(),tmin=masknz.maxt(),tmax=masknz.mint();
      
      for(int t=masknz.mint();t<=masknz.maxt();t++) {
	for(int z=masknz.minz();z<=masknz.maxz();z++) {
	  for(int y=masknz.miny();y<=masknz.maxy();y++) {
	    for(int x=masknz.minx();x<=masknz.maxx();x++) {
	      if (masknz(x,y,z,t)>0.5) {
		// if (masknz(x,y,z)>0.5) {
		if (x<xmin) xmin=x;
		if (x>xmax) xmax=x;
		if (y<ymin) ymin=y;
		if (y>ymax) ymax=y;
		if (z<zmin) zmin=z;
		if (z>zmax) zmax=z;
		if (t<tmin) tmin=t;
		if (t>tmax) tmax=t;
	      }
	    }
	  }
	}
      }
      cout << xmin << " " << 1+xmax-xmin << " " << ymin << " " << 1+ymax-ymin << " " << zmin << " " << 1+zmax-zmin << " " << tmin << " " << 1+tmax-tmin << " ";
    } else if (sarg=="-e") {
      if (!masks_used) { 
	if (vin.nvoxels()<1) { vin = vol; }
	masks_used=true;
	generate_masks(mask,masknz,vin,lthr,uthr); 
	vol = vin * mask; 
      }
      ColumnVector hist;
      int nbins=1000;
      double entropy=0;
      hist = vol.histogram(nbins,mask);
      double ntot = hist.Sum();
      for (int j=1; j<=nbins; j++) {
	if (hist(j)>0) {
	  entropy -= (hist(j)/ntot) * log(hist(j)/ntot);	
	}
      }
      entropy /= log((double) nbins);
      cout << entropy << " ";
    } else if (sarg=="-E") { 
      ColumnVector hist;
      int nbins=1000;
      double entropy=0;
      hist = vol.histogram(nbins,masknz);
      double ntot = hist.Sum();
      for (int j=1; j<=nbins; j++) {
	if (hist(j)>0) {
	  entropy -= (hist(j)/ntot) * log(hist(j)/ntot);	
	}
      }
      entropy /= log((double) nbins);
      cout << entropy << " ";
    } else if (sarg=="-k") {
      narg++;
      if (narg>=argc) {
	cerr << "Must specify an argument to -k" << endl;
	exit(2);
      }
      read_volume4D(mask,argv[narg]);
      if (!samesize(mask[0],vol[0])) {
	cerr << "Mask and image must be the same size" << endl;
	exit(3);
      }
      if ( mask.tsize() > vol.tsize() ) {
	cerr << "Mask and image must be the same size" << endl;
	exit(3);
      }
      if ( mask.tsize() != vol.tsize() && mask.tsize() != 1) {
	// copy the last max volume until the correct size is reached
	while (mask.tsize() < vol.tsize() ) {
   	  mask.addvolume(mask[mask.maxt()]);
        }
      }
      if (!masks_used) {
	masks_used=true;
	vin = vol;
      }
      float th= 0.5;
      if (th!=0) {
        mask.binarise(th);
      } else {
        mask.binarise(1);
      }
      generate_masks(mask,masknz,vin);
        if (mask.tsize()!=1) vol=vin*mask; 
	else vol=vin*mask[0];
    } else if (sarg=="-l") {
      narg++;
      if (narg<argc) {
        lthr = atof(argv[narg]);
      } else {
	cerr << "Must specify an argument to -l" << endl;
	exit(2);
      }
      lthr_used=true;
      if (!masks_used) {
	masks_used=true;
	vin = vol;
      }
      generate_masks(mask,masknz,vin,lthr,uthr);
      vol = vin * mask;
    } else if (sarg=="-u") {
      narg++;
      if (narg<argc) {
        uthr = atof(argv[narg]);
      } else {
	cerr << "Must specify an argument to -u" << endl;
	exit(2);
      }
      uthr_used=true;
      if (!masks_used) {
	masks_used=true;
	vin = vol;
      }
      generate_masks(mask,masknz,vin,lthr,uthr);
      vol = vin * mask;
    } else if (sarg=="-a") {
      vol = abs(vol);
    } else if (sarg=="-v") {
      if (masks_used) {
	cout << (long int) mask.sum() << " " 
	     << mask.sum() * vol.xdim() * vol.ydim() * vol.zdim() << " ";
      } else {
	cout << (long int) vol.nvoxels() << " "
	     << vol.nvoxels() * vol.xdim() * vol.ydim() * vol.zdim() << " ";
      }
    } else if (sarg=="-V") {
      if (masks_used) {
	cout << (long int) masknz.sum() << " " 
	     << masknz.sum() * vol.xdim() * vol.ydim() * vol.zdim() << " ";
      } else {
	long int nzvox;
	nzvox = nonzerocount(vol);
	cout << nzvox << " " 
	     << nzvox * vol.xdim() * vol.ydim() * vol.zdim() << " ";
      }
    } else if (sarg=="-d") {
	// hidden debug option!
      cout << vol.sum() << " ";
    } else if (sarg=="-s") {
	if (masks_used) cout << vol.stddev(mask) << " ";
	else cout << vol.stddev() << " ";
    } else if (sarg=="-S") {
      if (masks_used) {
	cout << vol.stddev(masknz) << " ";
      } else {
	cout << nonzerostddev(vol) << " ";
      }
    } else if (sarg=="-r") {
        if (masks_used) cout << vol.robustmin(mask) << " " << vol.robustmax(mask) << " ";
        else cout << vol.robustmin() << " " << vol.robustmax() << " ";
    } else if (sarg=="-R") {
	if (masks_used) cout << vol.min(mask) << " " << vol.max(mask) << " ";
	else cout << vol.min() << " " << vol.max() << " ";
    } else if (sarg=="-c") {
	ColumnVector cog(4);
	// convert from fsl mm to voxel to nifti sform coord
	cog.SubMatrix(1,3,1,1) = vol[0].cog();
	cog(4) = 1.0;
	cog = vol[0].newimagevox2mm_mat() * cog; 
	cout << cog(1) << " " << cog(2) << " " << cog(3) << " " ;
    } else if (sarg=="-C") {
    ColumnVector cog(4);
	// convert from fsl mm to fsl voxel coord to nifti voxel coord
	cog.SubMatrix(1,3,1,1) = vol[0].cog();
	cog(4) = 1.0;
	cog = (vol[0].niftivox2newimagevox_mat()).i() * cog;
	cout << cog(1) << " " << cog(2) << " " << cog(3) << " " ;
    } else if (sarg=="-p") {
      float n;
      narg++;
      if (narg<argc) {
        n = atof(argv[narg]);
      } else {
	cerr << "Must specify an argument to -p" << endl;
	exit(2);
      }
      if ( (n<0) || (n>100) ) {
    	cerr << "Percentile must be between 0 and 100" << endl;
    	exit(1);
      }
      if (masks_used) cout << vol.percentile((float) n/100.0, mask) << " ";
      else cout << vol.percentile((float) n/100.0) << " ";
    } else if (sarg=="-P") { 
      float n;
      narg++;
      if (narg<argc) {
        n = atof(argv[narg]);
      } else {
	cerr << "Must specify an argument to -P" << endl;
	exit(2);
      }
      if ( (n<0) || (n>100) ) {
    	cerr << "Percentile must be between 0 and 100" << endl;
    	exit(1);
      }
      if (!masks_used) {
	if (vin.nvoxels()<1) { vin = vol; }
	masks_used=true;
	generate_masks(mask,masknz,vin,lthr,uthr); 
	vol = vin * mask; 
      }
      cout << vol.percentile((float) n/100.0,masknz) << " ";
    } else if (sarg=="-h") {
      float n;
      narg++;
      if (narg<argc) {
        n = atof(argv[narg]);
      } else {
	cerr << "Must specify the number of bins" << endl;
	exit(2);
      }
      int nbins = (int) n;
      if (nbins<1) {
    	cerr << "Must specify at least 1 bin" << endl;
    	exit(1);
      }
      if (masks_used) {
	cout << vol.histogram(nbins,vol.min(),vol.max(),mask) << " ";
      } else {
	cout << vol.histogram(nbins,vol.min(),vol.max()) << " ";
      }
   } else if (sarg=="-H") {
      float n;
      narg++;
      if (narg<argc) {
        n = atof(argv[narg]);
      } else {
	cerr << "Must specify the number of bins" << endl;
	exit(2);
      }
      int nbins = (int) n;
      if (nbins<1) {
    	cerr << "Must specify at least 1 bin" << endl;
    	exit(1);
      }
      float min=0;
      narg++;
      if (narg<argc) {
        min = atof(argv[narg]);
      } else {
	cerr << "Must specify the histogram minimum intensity" << endl;
	exit(2);
      }
      float max=0;
      narg++;
      if (narg<argc) {
        max = atof(argv[narg]);
      } else {
	cerr << "Must specify the histogram maximum intensity" << endl;
	exit(2);
      }
      if (masks_used) {
	cout << vol.histogram(nbins,min,max,mask) << " ";
      } else {
	cout << vol.histogram(nbins,min,max) << " ";
      }
    } else {
	cerr << "Unrecognised option: " << sarg << endl;
	exit(3);
    }
  
    narg++;
  }

  cout << endl;
  return 0;
}



int main(int argc,char *argv[])
{

  Tracer tr("main");

  string progname=argv[0];
  if (argc<3) { 
    print_usage(progname);
    return 1; 
  }

  return fmrib_main_float(argc,argv);

}