diff --git a/smoothest.cc b/smoothest.cc
index af1abd0d9d9fd63f8e71e6dbe1d5cc2ea6b9b504..1cc8fbaf9355af4c29707b7e62c92909837df691 100644
--- a/smoothest.cc
+++ b/smoothest.cc
@@ -1,37 +1,99 @@
// $Id$
-// COPYRIGHT
-#include <getopt.h>
-
-#include <cstdlib>
+//#include <cstdlib>
#include <cmath>
#include <iostream>
+//#include <fstream>
#include <string>
+#include <map>
+
#include <avwio.h>
-#include <cassert>
+#include "options.h"
+#define _GNU_SOURCE 1
#define POSIX_SOURCE 1
+Option<bool> verbose(string("-V,--verbose"), false,
+ string("switch on diagnostic messages"),
+ false, BaseOption::no_argument);
+Option<bool> help(string("-h,--help"), false,
+ string("display this message"),
+ false, BaseOption::no_argument);
+Option<float> dof(string("-d,--dof"), 100.0,
+ string("number of degrees of freedom"),
+ true, BaseOption::requires_argument);
+Option<string> maskname(string("-m,--mask"), string("mask"),
+ string("brain mask volume"),
+ true, BaseOption::requires_argument);
+Option<string> zstatname(string("-z,--zstat"), "zstat",
+ string("filename of zstat image (not with -d)"),
+ true, BaseOption::requires_argument);
+Option<string> residname(string("-r,--res"), "res4d",
+ string("filename of `residual-fit' image (must use -d)"),
+ true, BaseOption::requires_argument);
+
+class Interpolate {
+public:
+ Interpolate() {
+ lut[5] = 1.5423138; lut[6] = 1.3757105; lut[7] = 1.2842680;
+ lut[8] = 1.2272151; lut[9] = 1.1885232; lut[10] = 1.1606988;
+ lut[11] = 1.1398000; lut[12] = 1.1235677; lut[13] = 1.1106196;
+ lut[14] = 1.1000651; lut[15] = 1.0913060; lut[16] = 1.0839261;
+ lut[17] = 1.0776276; lut[18] = 1.0721920; lut[19] = 1.0674553;
+ lut[20] = 1.0632924; lut[25] = 1.0483053; lut[30] = 1.0390117;
+ lut[40] = 1.0281339; lut[50] = 1.0219834; lut[60] = 1.0180339;
+ lut[70] = 1.0152850; lut[80] = 1.0132621; lut[90] = 1.0117115;
+ lut[100] = 1.0104851; lut[150] = 1.0068808; lut[200] = 1.0051200;
+ lut[300] = 1.0033865; lut[500] = 1.0020191;
+ }
+
+ inline float operator()(float v) {
+
+ float retval = 0;
+
+ map<int, float>::iterator i = lut.lower_bound(v);
+
+ if(i != lut.end()) {
+ if(i != lut.begin()) {
+ map<int, float>::iterator j = i--;
+
+ retval = (j->second - i->second)/(j->first - i->first)*(v - i->first)
+ + j->second;
+ }
+ } else {
+ retval = 1.0321/v + 1;
+ }
+
+ retval = pow(retval, 0.5);
+
+ return retval;
+ }
+
+private:
+ map<int, float> lut;
+};
+
+Interpolate interpolate;
template<class T> class Volume
{
private:
- unsigned int _x, _y, _z, _xy, _xyz;
- T *_data;
+ unsigned int x_, y_, z_, xy_, xyz_;
+ T *data_;
public:
- Volume(): _x(0), _y(0), _z(0),
- _xy(0), _data(NULL) {};
+ Volume(): x_(0), y_(0), z_(0),
+ xy_(0), data_(NULL) {};
Volume(const unsigned int x,
const unsigned int y,
const unsigned int z,
- T *data): _x(x), _y(y), _z(z) {
- _xy = x * y;
- _xyz = _xy * z;
- _data = new T[_xy * z];
+ T *data): x_(x), y_(y), z_(z) {
+ xy_ = x_ * y;
+ xyz_ = xy_ * z;
+ data_ = new T[xyz_];
}
Volume(AVW *ap, const unsigned int v) {
@@ -40,35 +102,35 @@ public:
short x, y, z, M;
AvwGetDim(ap, &x, &y, &z, &M);
- _x = x;
- _xy = x * y;
- _xyz = _xy * z;
- _data = new T[_xyz];
+ x_ = x;
+ xy_ = x_ * y;
+ xyz_ = xy_ * z;
+ data_ = new T[xyz_];
- AvwReadVolumes(ap, (void *)_data, 1);
+ AvwReadVolumes(ap, (void *)data_, 1);
}
unsigned int volsize() {
- return _xyz;
+ return xyz_;
}
inline T& operator () (const unsigned int x,
const unsigned int y,
const unsigned int z) {
- return _data[(z * _xy) + (y * _x) + x];
+ return data_[(z * xy_) + (y * x_) + x];
}
inline T& operator () (const unsigned int n) {
- return _data[n];
+ return data_[n];
}
inline T* buffer () {
- return _data;
+ return data_;
}
~Volume() {
try {
- delete [] _data;
+ delete [] data_;
} catch(...) {
// Er... panic!
}
@@ -81,39 +143,21 @@ inline float sqr(const T& a)
return a * a;
}
-static int verbose = 0;
-
-void usage(void)
-{
- cerr << endl;
- cerr << "usage: smoothest -d <fp-number> -m <filename> -r <filename>" << endl;
- cerr << endl;
- cerr << "Compulsory:" << endl;
- cerr << "-d,--dof" << "\t" << "degrees of freedom" << endl;
- cerr << "-m,--mask" << "\t" << "brain mask volume" << endl;
- cerr << "-r,--res" << "\t" << "4d `residual-of-fit' image" << endl;
- cerr << endl;
- cerr << "Optional:" << endl;
- cerr << "-v,--verbose" << "\t" << "display diagnostic messages" << endl;
- cerr << "-h,--help" << "\t" << "display this message" << endl;
- cerr << endl;
- cerr << endl;
-}
-
+// Standardise the residual field (assuming gaussianity)
unsigned long standardise(Volume<unsigned char> *mask,
Volume<float> **R,
short M)
{
unsigned long count = 0;
- if( M > 2 ) {
- // Standardise the residual field assuming gaussianity
- for ( unsigned int i = 0; i < mask->volsize(); i++) {
+ for ( unsigned int i = 0; i < mask->volsize(); i++) {
- if( (*mask)(i) ) {
+ if( (*mask)(i) ) {
- count ++;
+ count ++;
+ if( M > 2 ) {
+
// For each voxel
// calculate mean and standard deviation...
double Sx = 0.0, SSx = 0.0;
@@ -140,134 +184,98 @@ unsigned long standardise(Volume<unsigned char> *mask,
return count;
}
-int main(int argc, char **argv) {
-
- int hflag= 0, dflag = 1, mflag = 1, rflag = 1;
-
- static struct option long_options[] =
- {
- {"dof", 1, 0, 'd'},
- {"mask", 1, 0, 'm'},
- {"verbose", 0, &verbose, 1},
- {"res", 1, 0, 'r'},
- {"help", 0, &hflag, 1},
- {0, 0, 0, 0}
- };
-
- float v = 0;
- int c; int option_index = 0;
-
- string *maskname = NULL, *residname = NULL;
-
- while( (c = getopt_long(argc, argv, "d:m:r:hv",
- long_options, &option_index)) != -1 )
- {
- switch(c)
- {
- case 0:
- break;
- case 'd':
- v = atof(optarg);
- dflag --;
- break;
- case 'm':
- maskname = new string(optarg);
- mflag --;
- break;
- case 'r':
- residname = new string(optarg);
- rflag --;
- break;
- case 'v':
- verbose = 1;
- break;
- case 'h':
- hflag = 1;
- break;
- case '?':
- break;
- default:
- abort();
- }
- }
+string title = "\
+smoothest (Version 1.0b)\nCopyright(c) 2000, University of Oxford (Dave Flitney)";
- if ( hflag | dflag | mflag | rflag )
- {
- usage();
- return EXIT_FAILURE;
- }
+string examples = "\
+\tsmoothest -d <number> -r <filename> -m <filename>\n\
+\tsmoothest -z <filename> -m <filename>";
+
+int main(unsigned int argc, char **argv) {
+
+ OptionParser options(title, examples);
+
+ options.add(verbose);
+ options.add(help);
+ options.add(dof);
+ options.add(maskname);
+ options.add(residname);
+ options.add(zstatname);
+
+ options.parse_command_line(argc, argv);
+
+// if(verbose.value()) options.check_compulsory_arguments(true);
+
+ if(help.value()) {
+ options.usage();
+ exit(EXIT_SUCCESS);
+ }
+
+ if( !((zstatname.set() && residname.unset()) || (residname.set() && zstatname.unset())) ||
+ (zstatname.set() && (residname.set() || dof.set())) ||
+ (residname.set() && dof.unset()) ||
+ (maskname.unset()) ) {
+ options.usage();
+ cerr << endl;
+ cerr << "***************************************************************************" << endl;
+ cerr << "You must specify either a zstat image OR a 4d residual image." << endl;
+ cerr << "If processing a zstat image then you should not set degrees of freedom." << endl;
+ cerr << "You must specify a mask volume image filename" << endl;
+ cerr << "You must specify the degrees of freedom for processing a 4d residual image." << endl;
+ cerr << "***************************************************************************" << endl;
+ cerr << endl;
+ exit(EXIT_FAILURE);
+ }
- if(verbose) {
- cerr << "maskname = " << *maskname << endl;
- cerr << "residname = " << *residname << endl;
- cerr << "dof = " << v << endl;
+ if(verbose.value()) {
+ cout << "verbose = " << verbose.value() << endl;
+ cout << "help = " << help.value() << endl;
+ cout << "dof = " << dof.value() << endl;
+ cout << "maskname = " << maskname.value() << endl;
+ cout << "residname = " << residname.value() << endl;
+ cout << "zstatname = " << zstatname.value() << endl;
}
// Read the AVW mask image (single volume)
AVW *maskfile;
- if((maskfile=AvwOpen(maskname->c_str(),"r"))==NULL){
- printf("Failed to open mask: %s\n",maskname->c_str());
+ if((maskfile=AvwOpen(maskname.value().c_str(),"r"))==NULL){
+ printf("Failed to open mask: %s\n", maskname.value().c_str());
return(1);
}
- delete maskname;
Volume<unsigned char> *mask = new Volume<unsigned char>(maskfile, 0);
- // Read the AVW residual images (array of volumes)
- AVW *residfile;
- if((residfile=AvwOpen(residname->c_str(),"r"))==NULL){
- printf("Failed to open residuals: %s\n",residname->c_str());
- return(1);
+ AVW *datafile;
+ if(residname.set()) {
+ // Read the AVW residual images (array of volumes)
+ if((datafile = AvwOpen(residname.value().c_str(), "r")) == NULL){
+ printf("Failed to open residuals: %s\n", residname.value().c_str());
+ return(EXIT_FAILURE);
+ }
+ } else {
+ // Read the AVW zstat image (array of one volume)
+ if((datafile = AvwOpen(zstatname.value().c_str(), "r")) == NULL){
+ printf("Failed to open zstat volume: %s\n", zstatname.value().c_str());
+ return(EXIT_FAILURE);
+ }
}
+
short xdim, ydim, zdim, M;
- AvwGetDim(residfile, &xdim, &ydim, &zdim, &M);
- Volume<float> **R = new Volume<float> * [M];
- for ( unsigned short i = 0; i < M; i++ ) {
- R[i] = new Volume<float>(residfile, i);
- }
- delete residname;
-
- unsigned long count = standardise(mask, R, M);
-// if( M > 2 ) {
-// // Standardise the residual field assuming gaussianity
-// for ( unsigned int i = 0; i < mask->volsize(); i++) {
-
-// if( (*mask)(i) ) {
-
-// count ++;
-
-// // For each voxel
-// // calculate mean and standard deviation...
-// double Sx = 0.0, SSx = 0.0;
-
-// for ( unsigned short t = 0; t < M; t++ ) {
-// float& R_it = (*R[t])(i);
-
-// Sx += R_it;
-// SSx += sqr(R_it);
-// }
-
-// float mean = Sx / M;
-// float sd = sqrt( (SSx - (sqr(Sx) / M)) / (M - 1) );
-
-// // ... and use them to standardise to N(0, 1).
-// for ( unsigned short t = 0; t < M; t++ ) {
-// float& R_it = (*R[t])(i);
-
-// R_it = (R_it - mean) / sd;
-// }
-// }
-// }
-// }
-
- if(verbose) cerr << "Masked-in voxels = " << count << endl;
+ AvwGetDim(datafile, &xdim, &ydim, &zdim, &M);
+ Volume<float> **R = new Volume<float> * [M]; // R is an array of pointers to Volume<float> objects
+ for ( unsigned short i = 0; i < M; i++ )
+ R[i] = new Volume<float>(datafile, i);
+
+ unsigned long mask_volume = standardise(mask, R, M);
+
+ if(verbose.value()) cerr << "Masked-in voxels = " << mask_volume << endl;
- count = 0; // Reset count for the sums
+ unsigned long N = 0;
// Estimate the smoothness of the normalised residual field
+ // see TR00DF1 for mathematical description of the algorithm.
enum {X = 0, Y, Z};
- float F[3] = {0, 0, 0};
- float V0 = 0;
+ float SSminus[3] = {0, 0, 0}, S2 = 0;
for ( unsigned short z = 1; z < zdim ; z++ )
for ( unsigned short y = 1; y < ydim ; y++ )
@@ -278,63 +286,76 @@ int main(int argc, char **argv) {
(*mask)(x, y-1, z) &&
(*mask)(x, y, z-1) ) {
- count ++;
+ N ++;
for ( unsigned short t = 0; t < M; t++ ) {
Volume<float>& S_it = *R[t];
// Sum over M
-
- F[X] += S_it(x, y, z) * S_it(x-1, y, z);
- F[Y] += S_it(x, y, z) * S_it(x, y-1, z);
- F[Z] += S_it(x, y, z) * S_it(x, y, z-1);
- V0 += sqr(S_it(x, y, z));
+ SSminus[X] += S_it(x, y, z) * S_it(x-1, y, z);
+ SSminus[Y] += S_it(x, y, z) * S_it(x, y-1, z);
+ SSminus[Z] += S_it(x, y, z) * S_it(x, y, z-1);
+
+ S2 += sqr(S_it(x, y, z));
}
}
-
- if(verbose) {
- cerr << "Non-edge voxels = " << count << endl;
- cerr << "(v - 2)/(v - 1) = " << (v - 2)/(v - 1) << endl;
- }
- float norm = (v - 2) / ((v - 1) * count * M);
- F[X] *= norm;
- F[Y] *= norm;
- F[Z] *= norm;
+
+ float norm = 1.0/(float) N;
+ float v = dof.value(); // v - degrees of freedom (nu)
+ if(M > 1) {
+ if(verbose.value()) {
+ cerr << "Non-edge voxels = " << N << endl;
+ cerr << "(v - 2)/(v - 1) = " << (v - 2)/(v - 1) << endl;
+ }
+ norm = (v - 2) / ((v - 1) * N * M);
+ }
+
+ SSminus[X] *= norm;
+ SSminus[Y] *= norm;
+ SSminus[Z] *= norm;
- V0 *= norm;
+ S2 *= norm;
// Convert to sigma squared
- float W[3];
- W[X] = -1.0 / (4 * log(fabs(F[X]/V0)));
- W[Y] = -1.0 / (4 * log(fabs(F[Y]/V0)));
- W[Z] = -1.0 / (4 * log(fabs(F[Z]/V0)));
-
- // Convert to full width half maximum
- float FWHMx, FWHMy, FWHMz;
- FWHMx = sqrt(8 * log(2) * W[X]);
- FWHMy = sqrt(8 * log(2) * W[Y]);
- FWHMz = sqrt(8 * log(2) * W[Z]);
+ float sigmasq[3];
+ sigmasq[X] = -1.0 / (4 * log(fabs(SSminus[X]/S2)));
+ sigmasq[Y] = -1.0 / (4 * log(fabs(SSminus[Y]/S2)));
+ sigmasq[Z] = -1.0 / (4 * log(fabs(SSminus[Z]/S2)));
- cout << "FWHMx = " << FWHMx << " voxels, "
- << "FWHMy = " << FWHMy << " voxels, "
- << "FWHMz = " << FWHMz << " voxels" << endl;
-
- float Wc = pow(FWHMx * FWHMy * FWHMz, 1.0/3.0) / sqrt(4 * log(2));
+ float dLh = pow(sigmasq[X] * sigmasq[Y] * sigmasq[Z], -0.5) * pow(8, -0.5);
+ if(M > 1) dLh *= interpolate(v);
+ // Convert to full width half maximum
+ float FWHM[3];
+ FWHM[X] = sqrt(8 * log(2) * sigmasq[X]);
+ FWHM[Y] = sqrt(8 * log(2) * sigmasq[Y]);
+ FWHM[Z] = sqrt(8 * log(2) * sigmasq[Z]);
+ float resels = FWHM[X] * FWHM[Y] * FWHM[Z];
+
+ if(verbose.value()) {
+ cerr << "FWHMx = " << FWHM[X] << " voxels, "
+ << "FWHMy = " << FWHM[Y] << " voxels, "
+ << "FWHMz = " << FWHM[Z] << " voxels" << endl;
+ }
+
// Scale results from voxels to real world dimensions, e.g., mm.
float dx, dy, dz, dtr;
- AvwGetVoxDim(residfile, &dx, &dy, &dz, &dtr);
- FWHMx *= dx; FWHMy *= dy; FWHMz *= dz;
-
- cout << "FWHMx = " << FWHMx << " mm, "
- << "FWHMy = " << FWHMy << " mm, "
- << "FWHMz = " << FWHMz << " mm" << endl;
-
- float Fc = pow(FWHMx * FWHMy * FWHMz, 1.0/3.0);
-
- cout << "Fc = " << Fc << endl;
- cout << "Wc = " << Wc << endl;
+ AvwGetVoxDim(datafile, &dx, &dy, &dz, &dtr);
+ FWHM[X] *= dx; FWHM[Y] *= dy; FWHM[Z] *= dz;
+
+ if(verbose.value()) {
+ cerr << "FWHMx = " << FWHM[X] << " mm, "
+ << "FWHMy = " << FWHM[Y] << " mm, "
+ << "FWHMz = " << FWHM[Z] << " mm" << endl;
+ cerr << "DLH " << dLh << " voxels^-3" << endl;
+ cerr << "VOLUME " << mask_volume << " voxels" << endl;
+ cerr << "RESELS " << resels << " voxels per resel" << endl;
+ }
+ cout << "DLH " << dLh << endl;
+ cout << "VOLUME " << mask_volume << endl;
+ cout << "RESELS " << resels << endl;
+
return EXIT_SUCCESS;
}