diff --git a/cluster.cc b/cluster.cc
index 621ff8d0598c156b552da456cf2713eb2a296dbf..ada28bece47b230ac8278dca46826471de49568f 100644
--- a/cluster.cc
+++ b/cluster.cc
@@ -8,56 +8,155 @@
 
 #include <vector>
 #include <algorithm>
+#include <iomanip>
 #include "fmribmain.h"
 #include "newimageall.h"
 #include "infer.h"
-#include <iomanip>
+#include "options.h"
 
 using namespace NEWIMAGE;
+using std::vector;
+
+
+string title="cluster (Version 1.0)\nCopyright(c) 2000, University of Oxford (Mark Jenkinson)";
+string examples="cluster --in <filename> --thresh <value> [options]";
+
+Option<bool> verbose(string("-v,--verbose"), false, 
+		     string("switch on diagnostic messages"), 
+		     false, no_argument);
+Option<bool> help(string("-h,--help"), false,
+		  string("display this message"),
+		  false, no_argument);
+Option<bool> medxcoords(string("-m,--medx"), false,
+		  string("use MEDx coordinate conventions"),
+		  false, no_argument);
+Option<bool> mm(string("--mm"), false,
+		  string("use mm, not voxel, coordinates"),
+		  false, no_argument);
+Option<bool> minv(string("--min"), false,
+		  string("find minima instead of maxima"),
+		  false, no_argument);
+Option<bool> fractional(string("--fractional"), false,
+			string("interprets the threshold as a fraction of the robust range"),
+			false, no_argument);
+Option<bool> no_table(string("--no_table"), false,
+		      string("suppresses printing of the table info"),
+		      false, no_argument);
+Option<int> voxvol(string("--volume"), 0,
+		   string("number of voxels in the mask"),
+		   false, requires_argument);
+Option<int> numconnected(string("--connectivity"), 26,
+		   string("the connectivity of voxels (default 26)"),
+		   false, requires_argument);
+Option<float> dLh(string("-d,--dlh"), 4.0,
+		  string("smoothness estimate = sqrt(det(Lambda))"),
+		  false, requires_argument);
+Option<float> thresh(string("-t,--thresh,--zthresh"), 2.3,
+		     string("threshold for input volume"),
+		     true, requires_argument);
+Option<float> pthresh(string("-p,--pthresh"), 0.01,
+		      string("p-threshold for clusters"),
+		      false, requires_argument);
+Option<string> inputname(string("-i,--in,-z,--zstat"), string(""),
+			 string("filename of input volume"),
+			 true, requires_argument);
+Option<string> copename(string("-c,--cope"), string(""),
+			string("filename of input cope volume"),
+			false, requires_argument);
+Option<string> outindex(string("-o,--oindex"), string(""),
+			string("filename for output of cluster index (in size order)"),
+			false, requires_argument);
+Option<string> outthresh(string("--othresh"), string(""),
+			 string("filename for output of thresholded image"),
+			 false, requires_argument);
+Option<string> outsize(string("--osize"), string(""),
+		       string("filename for output of size image"),
+		       false, requires_argument);
+Option<string> outmax(string("--omax"), string(""),
+		      string("filename for output of max image"),
+		      false, requires_argument);
+Option<string> outmean(string("--omean"), string(""),
+		       string("filename for output of mean image"),
+		       false, requires_argument);
+
+
+template <class T>
+struct triple { T x; T y; T z; };
+
+template <class T, class S>
+void copyconvert(const vector<triple<T> >& oldcoords, 
+		 vector<triple<S> >& newcoords)
+{
+  newcoords.erase(newcoords.begin(),newcoords.end());
+  newcoords.resize(oldcoords.size());
+  for (unsigned int n=0; n<oldcoords.size(); n++) {
+    newcoords[n].x = (S) oldcoords[n].x;
+    newcoords[n].y = (S) oldcoords[n].y;
+    newcoords[n].z = (S) oldcoords[n].z;
+  }
+}
+
+template <class T, class S>
+void vox2medxcoord(vector<triple<T> >& coords, const volume<S>& vol)
+{
+  int ymax = vol.ysize();
+  for (unsigned int n=0; n<coords.size(); n++) {
+    coords[n].y = (T) (ymax - 1) - coords[n].y;
+  }
+}
 
-void print_usage(char *argv[])
+void vox2mmcoord(vector<triple<float> >& coords, 
+				   float vx, float vy, float vz)
 {
-  cerr << "Usage: " << argv[0] << " <input file> <output file> <cope-image>"  
-       << " <z-threshold> <p-threshold> <output_format> <DetLambdaHalf>"
-       << " <No. Voxels in Mask> [options]" << endl;
-  cerr << "[-f] : threshold is fractional (of robust range) instead of "
-       << "absolute"<<endl;
-  cerr << "[-6,18,26] : number of connections at each point (default is 26)"
-       <<endl;
-  cerr << "If threshold is -ve then binarization is inverted"<<endl;
-  cerr << "Output formats: the output intensity signifies:"<<endl;
-  cerr << "0 - size"<<endl;
-  cerr << "1 - size order (cluster index)"<<endl;
-  cerr << "2 - highest / lowest original value"<<endl;
-  cerr << "3 - mean original value"<<endl;
-  cerr << "4 - save thresholded map and print out stats"<<endl;
-  cerr << "5 - save p-thresholded cluster index map and print out stats"<<endl;
-  //  cerr << "6 - unprocessed unique integer label"<<endl;
+  for (unsigned int n=0; n<coords.size(); n++) {
+    coords[n].x = coords[n].x * vx;
+    coords[n].y = coords[n].y * vy;
+    coords[n].z = coords[n].z * vz;
+  }
+}
+
+
+template <class S>
+void vox2mmcoord(vector<triple<float> >& coords, 
+				   const volume<S>& vol)
+{
+  float vx = vol.xdim(), vy = vol.ydim(), vz = vol.zdim();
+  vox2mmcoord(coords,vx,vy,vz);
+}
+
+
+template <class T, class S>
+void addoriginoffset(vector<triple<T> >& coords, const volume<S>& vol)
+{
+  ColumnVector orig;
+  orig = vol.getorigin();
+  if (orig.SumAbsoluteValue() < 0.1)  return;  // zero origin -> no work
+  for (unsigned int n=0; n<coords.size(); n++) {
+    coords[n].x -= (T) (orig(1) - 1.0);   // -1.0 convert from SPM convention
+    coords[n].y -= (T) (orig(2) - 1.0);   // -1.0 convert from SPM convention
+    coords[n].z -= (T) (orig(3) - 1.0);   // -1.0 convert from SPM convention
+  }
 }
 
-struct triplei { int x; int y; int z; };
-struct triplef { float x; float y; float z; };
-typedef struct triplei tripleint;
-typedef struct triplef triplefloat;
 
 template <class T>
 void get_stats(const volume<int>& labelim, const volume<T>& origim,
-	       std::vector<int>& size,
-	       std::vector<T>& maxvals, std::vector<T>& meanvals,
-	       std::vector<tripleint>& max, std::vector<triplefloat>& cog,
-	       bool invert) 
+	       vector<int>& size,
+	       vector<T>& maxvals, vector<T>& meanvals,
+	       vector<triple<int> >& max, vector<triple<float> >& cog,
+	       bool minv) 
 {
   int labelnum = labelim.max();
   size.resize(labelnum+1,0);
   maxvals.resize(labelnum+1, (T) 0);
   meanvals.resize(labelnum+1,(T) 0);
-  tripleint zero;
+  triple<int> zero;
   zero.x = 0; zero.y = 0; zero.z = 0;
-  triplefloat zerof;
+  triple<float> zerof;
   zerof.x = 0; zerof.y = 0; zerof.z = 0;
   max.resize(labelnum+1,zero);
   cog.resize(labelnum+1,zerof);
-  std::vector<float> sum(labelnum+1,0.0);
+  vector<float> sum(labelnum+1,0.0);
   for (int z=labelim.minz(); z<=labelim.maxz(); z++) {
     for (int y=labelim.miny(); y<=labelim.maxy(); y++) {
       for (int x=labelim.minx(); x<=labelim.maxx(); x++) {
@@ -69,8 +168,8 @@ void get_stats(const volume<int>& labelim, const volume<T>& origim,
 	cog[idx].z+=((float) oxyz)*z;
 	sum[idx]+=(float) oxyz;
 	if ( (size[idx]==1) || 
-	     ( (oxyz>maxvals[idx]) && (!invert) ) || 
-	     ( (oxyz<maxvals[idx]) && (invert) ) ) 
+	     ( (oxyz>maxvals[idx]) && (!minv) ) || 
+	     ( (oxyz<maxvals[idx]) && (minv) ) ) 
 	  {
 	    maxvals[idx] = oxyz;
 	    max[idx].x = x;
@@ -93,17 +192,17 @@ void get_stats(const volume<int>& labelim, const volume<T>& origim,
 }
 
 
-std::vector<int> get_sortindex(const std::vector<int>& vals)
+vector<int> get_sortindex(const vector<int>& vals)
 {
   // return the mapping of old indices to new indices in the
   //   new *ascending* sort of vals
   int length=vals.size();
-  std::vector<pair<int, int> > sortlist(length);
+  vector<pair<int, int> > sortlist(length);
   for (int n=0; n<length; n++) {
     sortlist[n] = pair<int, int>(vals[n],n);
   }
   sort(sortlist.begin(),sortlist.end());  // O(N.log(N))
-  std::vector<int> idx(length);
+  vector<int> idx(length);
   for (int n=0; n<length; n++) {
     idx[n] = sortlist[n].second;
   }
@@ -111,9 +210,9 @@ std::vector<int> get_sortindex(const std::vector<int>& vals)
 }
 
 
-void get_sizeorder(const std::vector<int>& size, std::vector<int>& sizeorder) 
+void get_sizeorder(const vector<int>& size, vector<int>& sizeorder) 
 {
-  std::vector<int> sizecopy(size), idx;
+  vector<int> sizecopy(size), idx;
   idx = get_sortindex(sizecopy);
 
   // second part of pair is now the prior-index of the sorted values
@@ -125,174 +224,203 @@ void get_sizeorder(const std::vector<int>& size, std::vector<int>& sizeorder)
 }
 
 
-template <class T>
+template <class T, class S>
 void relabel_image(const volume<int>& labelim, volume<T>& relabelim,
-		   const std::vector<T>& newlabels)
+		   const vector<S>& newlabels)
 {
   copyconvert(labelim,relabelim);
   for (int z=relabelim.minz(); z<=relabelim.maxz(); z++) {
     for (int y=relabelim.miny(); y<=relabelim.maxy(); y++) {
       for (int x=relabelim.minx(); x<=relabelim.maxx(); x++) {
-	relabelim(x,y,z) = newlabels[labelim(x,y,z)];
+	relabelim(x,y,z) = (T) newlabels[labelim(x,y,z)];
       }
     }
   }
 }
 
 
+template <class T, class S>
+void print_results(const vector<int>& idx, 
+		   const vector<int>& size,
+		   const vector<int>& pthreshindex,
+		   const vector<float>& pvals,
+		   const vector<T>& maxvals,
+		   const vector<triple<S> >& maxpos,
+		   const vector<triple<float> >& cog,
+		   const vector<T>& copemaxval,
+		   const vector<triple<S> >& copemaxpos,
+		   const vector<T>& copemean,
+		   const volume<T>& zvol, const volume<T>& cope)
+{
+  int length=size.size();
+
+  vector<triple<float> > fmaxpos, fcopemaxpos, fcog;
+  copyconvert(maxpos,fmaxpos);
+  copyconvert(cog,fcog);
+  copyconvert(copemaxpos,fcopemaxpos);
+  if (medxcoords.value()) {
+    vox2medxcoord(fmaxpos,zvol);
+    vox2medxcoord(fcog,zvol);
+    if (!copename.unset()) vox2medxcoord(fcopemaxpos,cope);
+  }
+  cerr << "Prior to addoriginoffset" << endl;
+  addoriginoffset(fmaxpos,zvol);
+  addoriginoffset(fcog,zvol);
+  if (!copename.unset()) addoriginoffset(fcopemaxpos,zvol);
+  cerr << "Post addoriginoffset" << endl;
+  if (mm.value()) {
+    vox2mmcoord(fmaxpos,zvol);
+    vox2mmcoord(fcog,zvol);
+    if (!copename.unset()) vox2mmcoord(fcopemaxpos,cope);
+  }
+
+  cerr << "Setting up table" << endl;
+  string units=" (vox)";
+  if (!mm.unset()) units=" (mm)";
+  string tablehead;
+  tablehead = "Cluster Index\tVoxels";
+  if (!pthresh.unset()) tablehead += "\tP\t-log10(P)";
+  tablehead += "\tMax Z\tX" + units + "\tY" + units + "\tZ" + units
+    + "\tCOG X" + units + "\tCOG Y" + units + "\tCOG Z" + units;
+  if (!copename.unset()) {
+    tablehead+= "\tMax COPE\tX" + units + "\tY" + units + "\tZ" + units 
+                 + "\tMean COPE";
+  }
+
+  cout << tablehead << endl;
+
+  for (int n=length-1; n>=1; n--) {
+    int index=idx[n];
+    int k = size[index];
+    float p = pvals[index];
+    if (pthreshindex[index]>0) {
+      float mlog10p;
+      if (p>1e-38) mlog10p=-log10(p);
+      else mlog10p=42;
+      cout << setprecision(3) << pthreshindex[index] << "\t" << k << "\t"; 
+      if (!pthresh.unset()) { cout << p << "\t" << mlog10p << "\t"; }
+      cout << maxvals[index] << "\t" 
+	   << fmaxpos[index].x << "\t" << fmaxpos[index].y << "\t" 
+	   << fmaxpos[index].z << "\t"
+	   << fcog[index].x << "\t" << fcog[index].y << "\t" 
+	   << fcog[index].z;
+      if (!copename.unset()) { 
+	cout   << "\t" << copemaxval[index] << "\t"
+	       << fcopemaxpos[index].x << "\t" << fcopemaxpos[index].y << "\t" 
+	       << fcopemaxpos[index].z << "\t" << copemean[index];
+	}
+      cout << endl;
+    }
+  }
+}
+
+
+
 template <class T>
 int fmrib_main(int argc, char *argv[])
 {
   volume<int> labelim;
-  float p_thresh = atof(argv[5]);
-  int outform = atoi(argv[6]);
-  float dLh = atof(argv[7]);
-  int V = atoi(argv[8]);
-  std::vector<int> size;
-  std::vector<tripleint> max;
-  std::vector<triplefloat> cog;
-  std::vector<T> maxvals, meanvals;
-  float th;
-  bool invert=false, fractionalthresh=false;
-  int numconnected=26;
-  {
-    // read in the volume
-    volume<T> vin;
-    read_volume(vin,argv[1]);
-    
-    // parse the options
-    for (int n=9; n<argc; n++) {
-      string opt=argv[n];
-      if (opt=="-f") { fractionalthresh=true; }
-      else if (opt=="-6") { numconnected=6; }
-      else if (opt=="-18") { numconnected=18; }
-      else if (opt=="-26") { numconnected=26; }
-      else { 
-	cerr << "Unrecognised option "<<opt<<endl<<endl;
-	print_usage(argv); 
-	return -1;
+  vector<int> size;
+  vector<triple<int> > maxpos;
+  vector<triple<float> > cog;
+  vector<T> maxvals, meanvals;
+  float th = thresh.value();
+
+  // read in the volume
+  volume<T> zvol, mask, cope;
+  read_volume(zvol,inputname.value());
+  
+  if ( fractional.value() ) {
+    float frac = th;
+    th = frac*(zvol.robustmax() - zvol.robustmin()) + zvol.robustmin();
+  }
+  mask = zvol;
+  mask.threshold((T) th);
+  if (minv.value()) { mask = ((T) 1) - mask; }
+  
+  // Find the connected components
+  labelim = connected_components(mask,numconnected.value());
+  
+  // process according to the output format
+  get_stats(labelim,zvol,size,maxvals,meanvals,maxpos,cog,minv.value());
+
+
+  // process the cope image if entered
+  vector<int> copesize;
+  vector<triple<int> > copemaxpos;
+  vector<triple<float> > copecog;
+  vector<T> copemaxval, copemean;
+  if (!copename.unset()) {
+    read_volume(cope,copename.value());
+    get_stats(labelim,cope,copesize,copemaxval,copemean,copemaxpos,copecog,
+	      minv.value());
+  }
+
+  // re-threshold for p
+  int length = size.size();
+  size[0] = 0;  // force background to be ordered last
+  vector<int> pthreshsize;
+  vector<float> pvals(length);
+  pthreshsize = size;
+  int nozeroclust=0;
+  if (!pthresh.unset()) {
+    Infer infer(dLh.value(), th, voxvol.value());
+    for (int n=1; n<length; n++) {
+      int k = size[n];
+      pvals[n] = infer(k);
+      if (pvals[n]>pthresh.value()) {
+	pthreshsize[n] = 0;
+	nozeroclust++;
       }
     }
+  }
 
-    // threshold the volume
-    th = atof(argv[4]);
-    if (th<0.0) { 
-      invert=true;
-      th=-th;
-    }
+  // get sorted index (will revert to cluster size if no pthresh)
+  vector<int> idx;
+  idx = get_sortindex(pthreshsize);
 
-    if ( fractionalthresh ) {
-      float frac = th;
-      th = frac*(vin.robustmax() - vin.robustmin()) + vin.robustmin();
+  vector<int> threshidx(length);
+  for (int n=length-1; n>=1; n--) {
+    int index=idx[n];
+    if (pthreshsize[index]>0) {
+      threshidx[index] = n - nozeroclust;
     }
-    vin.threshold((T) th);
-    if (invert) { vin = ((T) 1) - vin; }
-    
-    // Find the connected components
-    labelim = connected_components(vin,numconnected);
-    
-    // process according to the output format
-    read_volume(vin,argv[1]);  // re-read the original volume
-    get_stats(labelim,vin,size,maxvals,meanvals,max,cog,invert);
   }
 
-  // cout << "CLUSTERS " << size.size()-1 << endl;
+  // print table
+  if (!no_table.value()) {
+    print_results(idx, size, threshidx, pvals, maxvals, maxpos, cog,
+		  copemaxval, copemaxpos, copemean, zvol, cope);
+  }
 
-  // save the results
-  if (outform==0) {
+  // save relevant volumes
+  if (!outindex.unset()) {
     volume<int> relabeledim;
-    size[0] = 0;
-    relabel_image(labelim,relabeledim,size);
-    save_volume(relabeledim,argv[2]);
-  } else if (outform==1) {
+    relabel_image(labelim,relabeledim,threshidx);
+    save_volume(relabeledim,outindex.value());
+  }
+  if (!outsize.unset()) {
+    volume<int> relabeledim;
+    relabel_image(labelim,relabeledim,pthreshsize);
+    save_volume(relabeledim,outsize.value());
+  }
+  if (!outmax.unset()) {
     volume<int> relabeledim;
-    std::vector<int> sizeorder;
-    size[0] = 0;
-    get_sizeorder(size,sizeorder);
-    relabel_image(labelim,relabeledim,sizeorder);
-    save_volume(relabeledim,argv[2]);
-  } else if (outform==2) {
-    volume<T> relabeledim;
     relabel_image(labelim,relabeledim,maxvals);
-    save_volume(relabeledim,argv[2]);
-  } else if (outform==3) {
-    volume<T> relabeledim;
+    save_volume(relabeledim,outmax.value());
+  }
+  if (!outmean.unset()) {
+    volume<int> relabeledim;
     relabel_image(labelim,relabeledim,meanvals);
-    save_volume(relabeledim,argv[2]);
-  } else if ( (outform==4) || (outform==5) ) {
-    Infer infer(dLh, th, V);
-    std::vector<int> idx;
-    int length = size.size();
-
-    // process the COPE image
-    volume<T> vin;
-    read_volume(vin,argv[3]);
-    std::vector<int> copesize;
-    std::vector<tripleint> copemax;
-    std::vector<triplefloat> copecog;
-    std::vector<T> copemaxval, copemean;
-    get_stats(labelim,vin,copesize,copemaxval,copemean,copemax,copecog,invert);
-
-    // re-threshold for p
-    size[0] = 0;
-    std::vector<int> pthreshsize;
-    pthreshsize = size;
-    int nozeroclust=0;
-    for (int n=1; n<length; n++) {
-      int k = size[n];
-      float p = infer(k);
-      if (p>p_thresh) {
-	pthreshsize[n] = 0;
-	nozeroclust++;
-      }
-    }
-    // cout << "No of sub-threshold p-clusters = " << nozeroclust << endl;
-    idx = get_sortindex(pthreshsize);
-
-    // Print results
-    cout << "Cluster Index\tVoxels\tP\t-log10(P)\tMax Z\tX\tY\tZ" 
-	 << "\tCOG X\tCOG Y\tCOG Z\tMax COPE\tX\tY\tZ\tMean COPE" << endl;
-    std::vector<int> pthreshindex(length,0);
-    for (int n=length-1; n>=1; n--) {
-      int index=idx[n];
-      int k = size[index];
-      float p = infer(k);
-      if (pthreshsize[index]>0) {
-	pthreshindex[index] = n - nozeroclust;
-	float mlog10p;
-	if (p>1e-38) mlog10p=-log10(p);
-	else mlog10p=42;
-	cout << setprecision(3) << pthreshindex[index] << "\t" << k << "\t" 
-	     << p << "\t" << mlog10p << "\t" 
-	     << maxvals[index] << "\t" 
-	     << max[index].x << "\t" << max[index].y << "\t" << max[index].z << "\t"
-	     << cog[index].x << "\t" << cog[index].y << "\t" << cog[index].z << "\t"
-	     << copemaxval[index] << "\t"
-	     << copemax[index].x << "\t" << copemax[index].y << "\t" 
-	     << copemax[index].z << "\t" << copemean[index] << endl;
-      }
-    }
-      
-    if (outform==4) {
-      // Threshold the input volume st it is 0 for all non-clusters
-      //   and maintains the same values otherwise
-      labelim.threshold(1);
-      volume<T> lcopy;
-      read_volume(vin,argv[1]);
-      copyconvert(labelim,lcopy);
-      vin *= lcopy;
-      save_volume(vin,argv[2]);
-    } else { // outform==5
-      volume<int> relabeledim;
-      relabel_image(labelim,relabeledim,pthreshindex);
-      save_volume(relabeledim,argv[2]);
-    }
-  } else if (outform==6) {
-    save_volume(labelim,argv[2]);
-  } else {
-    cerr << "Unrecognised output format" << endl;
-    print_usage(argv);
-    return -1;
+    save_volume(relabeledim,outmean.value());
+  }
+  if (!outthresh.unset()) {
+    // Threshold the input volume st it is 0 for all non-clusters
+    //   and maintains the same values otherwise
+    volume<T> lcopy;
+    relabel_image(labelim,lcopy,threshidx);
+    lcopy.threshold(0.5);
+    save_volume(zvol * lcopy,outthresh.value());
   }
 
   return 0;
@@ -304,11 +432,46 @@ int main(int argc,char *argv[])
 {
 
   Tracer tr("main");
-  if (argc<5) { 
-    print_usage(argv);
-    return -1; 
-  }
-  
-  return call_fmrib_main(dtype(argv[1]),argc,argv);
+
+  OptionParser options(title, examples);
+
+  options.add(inputname);
+  options.add(thresh);
+  options.add(outindex);
+  options.add(outthresh);
+  options.add(outsize);
+  options.add(outmax);
+  options.add(outmean);
+  options.add(pthresh);
+  options.add(copename);
+  options.add(voxvol);
+  options.add(dLh);
+  options.add(fractional);
+  options.add(numconnected);
+  options.add(medxcoords);
+  options.add(mm);
+  options.add(minv);
+  options.add(no_table);
+  options.add(verbose);
+  options.add(help);
+
+  options.parse_command_line(argc, argv);
+
+  if ( (help.value()) || (!options.check_compulsory_arguments(true)) )
+    {
+      options.usage();
+      exit(EXIT_FAILURE);
+    }
+
+  if ( (!pthresh.unset()) && (dLh.unset() || voxvol.unset()) ) 
+    {
+      options.usage();
+      cerr << endl 
+	   << "Both --dlh and --volume MUST be set if --pthresh is used!" 
+	   << endl;
+      exit(EXIT_FAILURE);
+    }
+    
+  return call_fmrib_main(dtype(inputname.value()),argc,argv);
 }