Skip to content
Snippets Groups Projects
Commit 49c1e621 authored by Mark Jenkinson's avatar Mark Jenkinson
Browse files

Added transpose option and tidied up indenting of the code

parent 9de16523
No related branches found
No related tags found
No related merge requests found
...@@ -66,6 +66,10 @@ Option<bool> bin_mask(string("--no_bin"), true, ...@@ -66,6 +66,10 @@ Option<bool> bin_mask(string("--no_bin"), true,
Option<int> order(string("--order"), 1, Option<int> order(string("--order"), 1,
string(" select number of Eigenvariates (default 1)"), string(" select number of Eigenvariates (default 1)"),
false, requires_argument); 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; int nonoptarg;
...@@ -85,6 +89,7 @@ int main(int argc,char *argv[]) ...@@ -85,6 +89,7 @@ int main(int argc,char *argv[])
options.add(eig); options.add(eig);
options.add(order); options.add(order);
options.add(bin_mask); options.add(bin_mask);
options.add(transpose);
options.add(verbose); options.add(verbose);
options.add(help); options.add(help);
...@@ -145,34 +150,36 @@ int main(int argc,char *argv[]) ...@@ -145,34 +150,36 @@ int main(int argc,char *argv[])
mask.binarise(1e-8); // arbitrary "0" threshold mask.binarise(1e-8); // arbitrary "0" threshold
if(eig.value()){ if(eig.value()){
Matrix dat, evecs, scales; Matrix dat, evecs, scales;
scales = mask_nonbin.matrix(mask); scales = mask_nonbin.matrix(mask);
dat = vin.matrix(mask); dat = vin.matrix(mask);
if(!bin_mask.value()) if(!bin_mask.value())
dat = SP (dat, ones(dat.Nrows(),1) * scales); dat = SP (dat, ones(dat.Nrows(),1) * scales);
dat = remmean(dat,1); dat = remmean(dat,1);
if (verbose.value()) { if (verbose.value()) {
cout << "Number of voxels used = " << dat.Ncols() << endl; cout << "Number of voxels used = " << dat.Ncols() << endl;
} }
SymmetricMatrix Corr; SymmetricMatrix Corr;
Corr << dat * dat.t() / dat.Ncols(); Corr << dat * dat.t() / dat.Ncols();
DiagonalMatrix tmpD; DiagonalMatrix tmpD;
EigenValues(Corr,tmpD,evecs); EigenValues(Corr,tmpD,evecs);
evecs = fliplr(evecs.Columns(evecs.Ncols()-order.value()+1 , evecs.Ncols())) * std::sqrt(dat.Nrows()); evecs = fliplr(evecs.Columns(evecs.Ncols()-order.value()+1 , evecs.Ncols())) * std::sqrt(dat.Nrows());
Matrix res2 = mean(dat,2); Matrix res2 = mean(dat,2);
res2 = res2.Column(1).t() * evecs.Column(1); res2 = res2.Column(1).t() * evecs.Column(1);
if((float)res2.AsScalar()<0) if((float)res2.AsScalar()<0)
evecs = -1.0 * evecs; evecs = -1.0 * evecs;
if (outmat.set()) { if (transpose.value()) { evecs=evecs.t(); }
if (outmat.set()) {
write_ascii_matrix(evecs,outmat.value()); write_ascii_matrix(evecs,outmat.value());
} else { } else {
cout << evecs << endl; cout << evecs << endl;
} }
if (transpose.value()) { evecs=evecs.t(); }
} }
else{ else{
int nt = vin.tsize(); int nt = vin.tsize();
...@@ -184,43 +191,45 @@ int main(int argc,char *argv[]) ...@@ -184,43 +191,45 @@ int main(int argc,char *argv[])
Matrix meants(nt,nvox); Matrix meants(nt,nvox);
meants = 0; meants = 0;
long int num=0; long int num=0;
for (int z=mask.minz(); z<=mask.maxz(); z++) { for (int z=mask.minz(); z<=mask.maxz(); z++) {
for (int y=mask.miny(); y<=mask.maxy(); y++) { for (int y=mask.miny(); y<=mask.maxy(); y++) {
for (int x=mask.minx(); x<=mask.maxx(); x++) { for (int x=mask.minx(); x<=mask.maxx(); x++) {
if (mask(x,y,z)>0.5) { if (mask(x,y,z)>0.5) {
num++; num++;
if (showall.value()) { if (showall.value()) {
ColumnVector v(4); ColumnVector v(4);
v << x << y << z << 1.0; v << x << y << z << 1.0;
v = (vin[0].niftivox2newimagevox_mat()).i() * v; v = (vin[0].niftivox2newimagevox_mat()).i() * v;
meants(1,num) = v(1); meants(1,num) = v(1);
meants(2,num) = v(2); meants(2,num) = v(2);
meants(3,num) = v(3); meants(3,num) = v(3);
meants.SubMatrix(4,nt,num,num) = vin.voxelts(x,y,z); meants.SubMatrix(4,nt,num,num) = vin.voxelts(x,y,z);
} else { } else {
meants += vin.voxelts(x,y,z); meants += vin.voxelts(x,y,z);
} }
} }
} }
} }
} }
if (verbose.value()) { if (verbose.value()) {
cout << "Number of voxels used = " << num << endl; cout << "Number of voxels used = " << num << endl;
} }
// normalise for number of valid entries if averaging // normalise for number of valid entries if averaging
if (!showall.value()) { if (!showall.value()) {
if (num>0) meants /= (float) num; if (num>0) meants /= (float) num;
} }
// save the result // save the result
if (transpose.value()) { meants=meants.t(); }
if (outmat.set()) { if (outmat.set()) {
write_ascii_matrix(meants,outmat.value()); write_ascii_matrix(meants,outmat.value());
} else { } else {
cout << meants << endl; cout << meants << endl;
} }
if (transpose.value()) { meants=meants.t(); }
} }
return 0; return 0;
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment