From 49c1e62183470e4982ee331cd5e358b926f58582 Mon Sep 17 00:00:00 2001
From: Mark Jenkinson <mark@fmrib.ox.ac.uk>
Date: Tue, 24 Mar 2009 11:35:47 +0000
Subject: [PATCH] Added transpose option and tidied up indenting of the code
---
fslmeants.cc | 87 +++++++++++++++++++++++++++++-----------------------
1 file changed, 48 insertions(+), 39 deletions(-)
diff --git a/fslmeants.cc b/fslmeants.cc
index a777243..35d8ac5 100644
--- a/fslmeants.cc
+++ b/fslmeants.cc
@@ -66,6 +66,10 @@ Option<bool> bin_mask(string("--no_bin"), true,
Option<int> order(string("--order"), 1,
string(" select number of Eigenvariates (default 1)"),
false, requires_argument);
+Option<bool> transpose(string("--transpose"), false,
+ string(" output results in transpose format (one row per voxel/mean)"),
+ false, no_argument);
+
int nonoptarg;
@@ -85,6 +89,7 @@ int main(int argc,char *argv[])
options.add(eig);
options.add(order);
options.add(bin_mask);
+ options.add(transpose);
options.add(verbose);
options.add(help);
@@ -145,34 +150,36 @@ int main(int argc,char *argv[])
mask.binarise(1e-8); // arbitrary "0" threshold
if(eig.value()){
- Matrix dat, evecs, scales;
- scales = mask_nonbin.matrix(mask);
- dat = vin.matrix(mask);
- if(!bin_mask.value())
- dat = SP (dat, ones(dat.Nrows(),1) * scales);
- dat = remmean(dat,1);
-
- if (verbose.value()) {
+ Matrix dat, evecs, scales;
+ scales = mask_nonbin.matrix(mask);
+ dat = vin.matrix(mask);
+ if(!bin_mask.value())
+ dat = SP (dat, ones(dat.Nrows(),1) * scales);
+ dat = remmean(dat,1);
+
+ if (verbose.value()) {
cout << "Number of voxels used = " << dat.Ncols() << endl;
}
-
- SymmetricMatrix Corr;
- Corr << dat * dat.t() / dat.Ncols();
- DiagonalMatrix tmpD;
+
+ SymmetricMatrix Corr;
+ Corr << dat * dat.t() / dat.Ncols();
+ DiagonalMatrix tmpD;
EigenValues(Corr,tmpD,evecs);
- evecs = fliplr(evecs.Columns(evecs.Ncols()-order.value()+1 , evecs.Ncols())) * std::sqrt(dat.Nrows());
-
- Matrix res2 = mean(dat,2);
- res2 = res2.Column(1).t() * evecs.Column(1);
-
- if((float)res2.AsScalar()<0)
- evecs = -1.0 * evecs;
-
- if (outmat.set()) {
+ evecs = fliplr(evecs.Columns(evecs.Ncols()-order.value()+1 , evecs.Ncols())) * std::sqrt(dat.Nrows());
+
+ Matrix res2 = mean(dat,2);
+ res2 = res2.Column(1).t() * evecs.Column(1);
+
+ if((float)res2.AsScalar()<0)
+ evecs = -1.0 * evecs;
+
+ if (transpose.value()) { evecs=evecs.t(); }
+ if (outmat.set()) {
write_ascii_matrix(evecs,outmat.value());
} else {
cout << evecs << endl;
}
+ if (transpose.value()) { evecs=evecs.t(); }
}
else{
int nt = vin.tsize();
@@ -184,43 +191,45 @@ int main(int argc,char *argv[])
Matrix meants(nt,nvox);
meants = 0;
long int num=0;
-
+
for (int z=mask.minz(); z<=mask.maxz(); z++) {
for (int y=mask.miny(); y<=mask.maxy(); y++) {
for (int x=mask.minx(); x<=mask.maxx(); x++) {
- if (mask(x,y,z)>0.5) {
- num++;
- if (showall.value()) {
- ColumnVector v(4);
- v << x << y << z << 1.0;
- v = (vin[0].niftivox2newimagevox_mat()).i() * v;
- meants(1,num) = v(1);
- meants(2,num) = v(2);
- meants(3,num) = v(3);
- meants.SubMatrix(4,nt,num,num) = vin.voxelts(x,y,z);
- } else {
- meants += vin.voxelts(x,y,z);
- }
- }
+ if (mask(x,y,z)>0.5) {
+ num++;
+ if (showall.value()) {
+ ColumnVector v(4);
+ v << x << y << z << 1.0;
+ v = (vin[0].niftivox2newimagevox_mat()).i() * v;
+ meants(1,num) = v(1);
+ meants(2,num) = v(2);
+ meants(3,num) = v(3);
+ meants.SubMatrix(4,nt,num,num) = vin.voxelts(x,y,z);
+ } else {
+ meants += vin.voxelts(x,y,z);
+ }
+ }
}
}
}
-
+
if (verbose.value()) {
cout << "Number of voxels used = " << num << endl;
}
-
+
// normalise for number of valid entries if averaging
if (!showall.value()) {
if (num>0) meants /= (float) num;
}
-
+
// save the result
+ if (transpose.value()) { meants=meants.t(); }
if (outmat.set()) {
write_ascii_matrix(meants,outmat.value());
} else {
cout << meants << endl;
}
+ if (transpose.value()) { meants=meants.t(); }
}
return 0;
}
--
GitLab