fslstats.cc

/*  avwstats++.cc

    Mark Jenkinson, FMRIB Image Analysis Group

    Copyright (C) 2003 University of Oxford  */

/*  CCOPYRIGHT  */


// Like avwstats but better!  :)

#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;

void print_usage(const string& progname) {
  cout << "Usage: avwstats++ <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 << "-k <mask>    : use the specified image (filename) for masking - overrides lower and upper thresholds" << 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);
  masknz = mask * (1.0f - binarise(vin,0.0f, 0.0f));
  return 0;
}

int generate_masks(const volume4D<float> &mask, volume4D<float> &masknz, const volume4D<float> &vin)
{
  masknz = mask * (1.0f - binarise(vin,0.0f, 0.0f));
  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=="-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") {
      if (masks_used) {
	cout << vol.mincoordx(mask) << " " << vol.mincoordy(mask) << " " 
	   << vol.mincoordz(mask) << " ";
      } else {
	cout << vol.mincoordx() << " " << vol.mincoordy() << " " 
	   << vol.mincoordz() << " ";
      }
    } else if (sarg=="-x") { 
      if (masks_used) {
	cout << vol.maxcoordx(mask) << " " << vol.maxcoordy(mask) << " " 
	     << vol.maxcoordz(mask) << " ";
      } else {
	cout << vol.maxcoordx() << " " << vol.maxcoordy() << " " 
	     << vol.maxcoordz() << " ";
      }
    } 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() ) {
	// 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);
      vol = vin * mask;
    } 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(vin);
    } 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) {
	float rmin=vol.robustmin(mask);
	float rmax=vol.robustmax(mask);
	if (rmin>rmax) { 
	  float tmp = rmax;
	  rmax = rmin;
	  rmin = tmp;
	}
	cout << rmin << " " << rmax << " ";
      } 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 sform coord
	cog.SubMatrix(1,3,1,1) = vol[0].cog();
	cog(4) = 1.0;
	if (vol[0].sform_code()!=NIFTI_XFORM_UNKNOWN) {
	  cog = vol[0].sform_mat() * (vol[0].sampling_mat()).i() * cog; 
	}
	cout << cog(1) << " " << cog(2) << " " << cog(3) << " " ;
    } else if (sarg=="-C") {
	ColumnVector cog(4);
	// convert from fsl mm to voxel coord
	cog.SubMatrix(1,3,1,1) = vol[0].cog();
	cog(4) = 1.0;
	cog = (vol[0].sampling_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 {
	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);

}