From ed7b47a5cbb2369c51dde101fd8ced83636e805a Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Mon, 15 May 2006 09:09:30 +0000
Subject: [PATCH] GLM cleanup
---
meldata.cc | 38 +++++++++++++++++++++++++++-----------
meloptions.h | 5 +++++
2 files changed, 32 insertions(+), 11 deletions(-)
diff --git a/meldata.cc b/meldata.cc
index 362313f..3b86f30 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -652,19 +652,35 @@ namespace Melodic{
// Read in thresholded maps
Matrix threshmaps;
- int ic_ctr = 2;
+ int ic_ctr = 1;
message(" Reading in thresholded IC maps ...");
- read_volume4D(tmpvol,logger.appendDir("stats/thresh_zstat1"));
- threshmaps = tmpvol.matrix(Mask);
+ // read_volume4D(tmpvol,logger.appendDir("stats/thresh_zstat1"));
+ // threshmaps = tmpvol.matrix(Mask);
while(fsl_imageexists(logger.appendDir("stats/thresh_zstat" + num2str(ic_ctr))))
{
read_volume4D(tmpvol,logger.appendDir("stats/thresh_zstat" + num2str(ic_ctr)));
- threshmaps = threshmaps & tmpvol.matrix(Mask);
+
+ Matrix newmap;
+ newmap = tmpvol.matrix(Mask);
+ Matrix corrcf;
+ corrcf = pemaps.t() | newmap.t();
+ corrcf = corrcoef(corrcf);
+ corrcf = corrcf.SubMatrix(corrcf.Nrows(),corrcf.Nrows(),1,pemaps.Nrows());
+ message("IC map " << ic_ctr << ": " << corrcf);
+ if(opts.glmPEreg.value()==0.0 || opts.glmPEreg.value() > max(abs(corrcf)).AsScalar()){
+ message(" removing IC map " << endl);
+ if(threshmaps.Storage() == 0)
+ threshmaps = newmap;
+ else
+ threshmaps = threshmaps & newmap;
+ }
+ else
+ message(" not removing IC map " << endl);
ic_ctr++;
}
message(" ... done " << endl);
- outMsize(" threshmaps ",threshmaps);
- outMsize(" pemaps ",pemaps);
+ // outMsize(" threshmaps ",threshmaps);
+ //outMsize(" pemaps ",pemaps);
//read in data
message(" Reading in " << opts.inputfname.value() << " ... ");
@@ -681,21 +697,21 @@ namespace Melodic{
allmaps = pemaps & threshmaps;
outMsize(" allmaps ",allmaps);
- outMsize(" orig_data ",orig_data);
+ //outMsize(" orig_data ",orig_data);
allmaps = pemaps & threshmaps;
alltcs = ( orig_data * allmaps.t() ) * pinv (allmaps * allmaps.t() );
- outMsize(" alltcs",alltcs);
+ //outMsize(" alltcs",alltcs);
resi = orig_data - alltcs * allmaps;
pemaps = pemaps - pemaps * threshmaps.t() * pinv (threshmaps * threshmaps.t() ) * threshmaps;
- outMsize(" ... new PEs :",pemaps);
+ //outMsize(" ... new PEs :",pemaps);
Matrix R;
R = Identity(orig_data.Nrows()) - alltcs * pinv(alltcs.t() * alltcs) * alltcs.t();
Matrix sigsq;
sigsq = sum(SP(resi,resi)) / R.Trace();
- outMsize(" ... new sigmasquareds :",sigsq);
+ //outMsize(" ... new sigmasquareds :",sigsq);
Matrix new_dof;
new_dof = old_dof.matrix(Mask);
@@ -703,7 +719,7 @@ namespace Melodic{
for(int i_ctr=1;i_ctr<=threshmaps.Nrows();i_ctr++)
if(abs(threshmaps(i_ctr,j_ctr)) > 0)
new_dof(1,j_ctr)--;
- outMsize(" ... new DOFs :", new_dof);
+ //outMsize(" ... new DOFs :", new_dof);
message(" ... done " << endl);
//save out new GLM results
diff --git a/meloptions.h b/meloptions.h
index a2e8c6a..80e8ea6 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -70,6 +70,7 @@ class MelodicOptions {
Option<string> filtermix;
Option<string> filter;
Option<string> glmcleanup;
+ Option<float> glmPEreg;
Option<bool> genreport;
Option<float> tr;
@@ -204,6 +205,9 @@ class MelodicOptions {
glmcleanup(string("--cleanup"), string(""),
string("<feat dir>\n "),
false, requires_argument),
+ glmPEreg(string("--PEreg"), 0.0,
+ string("spatial correlation threshold between ICA and PE maps"),
+ false, requires_argument),
genreport(string("--report"), false,
string("generate Melodic web report"),
false, no_argument),
@@ -320,6 +324,7 @@ class MelodicOptions {
options.add(filtermix);
options.add(filter);
options.add(glmcleanup);
+ options.add(glmPEreg);
options.add(genreport);
options.add(tr);
options.add(logPower);
--
GitLab