Fully working version, with all bells and whistles.

Some testing has been done, but more will be done with FEAT4 soon.

cluster.cc

@@ -78,6 +78,19 @@ Option<string> outmax(string("--omax"), string(""),
 Option<string> outmean(string("--omean"), string(""),
 		       string("filename for output of mean image"),
 		       false, requires_argument);
+Option<string> transformname(string("-x,--xfm"), string(""),
+		       string("filename for transform of input->talairach volume"),
+		       false, requires_argument);
+Option<string> talvolname(string("--talvol"), string(""),
+		       string("filename for talairach volume"),
+		       false, requires_argument);
+
+template <class T>
+void print_dims(const volume<T>& vol)
+{
+  cout << "dims=(" << vol.xdim()*vol.xsign() << ","
+       << vol.ydim()*vol.ysign() << "," << vol.zdim()*vol.zsign() << ")";
+}
 
 
 template <class T>
@@ -139,6 +152,23 @@ void addoriginoffset(vector<triple<T> >& coords, const volume<S>& vol)
 }
 
 
+template <class T, class S>
+void applyvoxelxfm(vector<triple<T> >& coords, const Matrix& mat, 
+	      const volume<S>& source, const volume<S>& dest)
+{
+  Matrix voxmat;
+  voxmat = dest.sampling_mat().i() * mat * source.sampling_mat();
+  ColumnVector vec(4);
+  for (unsigned int n=0; n<coords.size(); n++) {
+    vec << coords[n].x << coords[n].y << coords[n].z << 1.0;
+    vec = voxmat * vec;     // apply voxel xfm
+    coords[n].x = vec(1);
+    coords[n].y = vec(2);
+    coords[n].z = vec(3);
+  }
+}
+
+
 template <class T>
 void get_stats(const volume<int>& labelim, const volume<T>& origim,
 	       vector<int>& size,
@@ -258,18 +288,38 @@ void print_results(const vector<int>& idx,
   copyconvert(maxpos,fmaxpos);
   copyconvert(cog,fcog);
   copyconvert(copemaxpos,fcopemaxpos);
-  if (medxcoords.value()) {
+  volume<T> talvol;
+  if ( !transformname.unset() && !talvolname.unset() ) {
+    read_volume(talvol,talvolname.value());
+    Matrix trans;
+    read_matrix(trans,transformname.value(),zvol);
+    if (verbose.value()) { 
+      cout << "Transformation Matrix filename = "<<transformname.value()<<endl;
+      cout << trans.Nrows() << " " << trans.Ncols() << endl; 
+      cout << "Transformation Matrix = " << endl;
+      cout << trans << endl;
+    }
+    applyvoxelxfm(fmaxpos,trans,zvol,talvol);
+    applyvoxelxfm(fcog,trans,zvol,talvol);
+    if (!copename.unset()) applyvoxelxfm(fcopemaxpos,trans,zvol,talvol);
+  }
+  if ( medxcoords.value() && (!mm.value()) ) {
     vox2medxcoord(fmaxpos,zvol);
     vox2medxcoord(fcog,zvol);
     if (!copename.unset()) vox2medxcoord(fcopemaxpos,cope);
   }
   if (mm.value()) {
-    addoriginoffset(fmaxpos,zvol);
-    addoriginoffset(fcog,zvol);
-    if (!copename.unset()) addoriginoffset(fcopemaxpos,zvol);
-    vox2mmcoord(fmaxpos,zvol);
-    vox2mmcoord(fcog,zvol);
-    if (!copename.unset()) vox2mmcoord(fcopemaxpos,cope);
+    const volume<T> *refvol = &zvol;
+    if ( !transformname.unset() && !talvolname.unset() ) refvol = &talvol;
+    if (verbose.value()) { 
+      cout << "Using origin : " << (refvol->getorigin()).t() << endl; 
+    }
+    addoriginoffset(fmaxpos,*refvol);
+    addoriginoffset(fcog,*refvol);
+    if (!copename.unset()) addoriginoffset(fcopemaxpos,*refvol);
+    vox2mmcoord(fmaxpos,*refvol);
+    vox2mmcoord(fcog,*refvol);
+    if (!copename.unset()) vox2mmcoord(fcopemaxpos,*refvol);  // used cope before
   }
 
   string units=" (vox)";
@@ -326,6 +376,7 @@ int fmrib_main(int argc, char *argv[])
   // read in the volume
   volume<T> zvol, mask, cope;
   read_volume(zvol,inputname.value());
+  if (verbose.value()) {cout<<"Zvol: "; print_dims(zvol); cout<<endl;}
   
   if ( fractional.value() ) {
     float frac = th;
@@ -334,9 +385,11 @@ int fmrib_main(int argc, char *argv[])
   mask = zvol;
   mask.threshold((T) th);
   if (minv.value()) { mask = ((T) 1) - mask; }
+  if (verbose.value()) {cout<<"Mask: "; print_dims(mask); cout<<endl;}
   
   // Find the connected components
   labelim = connected_components(mask,numconnected.value());
+  if (verbose.value()) {cout<<"Labelim: "; print_dims(labelim); cout<<endl;}
   
   // process according to the output format
   get_stats(labelim,zvol,size,maxvals,meanvals,maxpos,cog,minv.value());
@@ -431,43 +484,64 @@ int main(int argc,char *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);
-    }
+  try {
+    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(transformname);
+    options.add(talvolname);
+    options.add(verbose);
+    options.add(help);
+    
+    options.parse_command_line(argc, argv);
 
-  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);
-    }
+    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);
+      }
     
+    if ( ( !transformname.unset() && talvolname.unset() ) ||
+	 ( transformname.unset() && (!talvolname.unset()) ) ) 
+      {
+	options.usage();
+	cerr << endl 
+	     << "Both --xfm and --talvol MUST be set together." 
+	     << endl;
+	exit(EXIT_FAILURE);
+      }
+  }  catch(X_OptionError& e) {
+    options.usage();
+    cerr << endl << e.what() << endl;
+    exit(EXIT_FAILURE);
+  } catch(std::exception &e) {
+    cerr << e.what() << endl;
+  } 
+
   return call_fmrib_main(dtype(inputname.value()),argc,argv);
+
 }