From 220bc12d731fe68c52d2895145581d9f4ffe0bdc Mon Sep 17 00:00:00 2001
From: Mark Jenkinson <mark@fmrib.ox.ac.uk>
Date: Fri, 17 Nov 2000 20:45:20 +0000
Subject: [PATCH] Fully working version, with all bells and whistles. Some
testing has been done, but more will be done with FEAT4 soon.
---
cluster.cc | 158 +++++++++++++++++++++++++++++++++++++++--------------
1 file changed, 116 insertions(+), 42 deletions(-)
diff --git a/cluster.cc b/cluster.cc
index 90bfc2c..84358b8 100644
--- a/cluster.cc
+++ b/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);
+
}
--
GitLab