diff --git a/meldata.cc b/meldata.cc
index ec8538187966d91ed0ce0992bb08f6c7aabdcd8e..a735d8b0b539a1b03ac78501f61527790d3aa41b 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -48,7 +48,11 @@ namespace Melodic{
exit(2);
}
-
+ // Do we need to clean up FEAT results?
+ if(opts.glmcleanupmode)
+ GLMcleanup_preproc();
+
+
// clean /tmp
ostringstream osc;
osc << "rm " << string(Mean_fname) <<"* ";
@@ -565,6 +569,159 @@ namespace Melodic{
} //void status()
+
+ void MelodicData::GLMcleanup_preproc()
+ {
+ if(fsl_imageexists(opts.glmcleanup.value() + "/stats/res4d")){
+ volume4D<float> RawData;
+ message("File res4d exists " << endl);
+ message(" Reading data file " << opts.glmcleanup.value() + "/stats/res4d" << " ... ");
+ read_volume4D(RawData,string(opts.glmcleanup.value() + "/stats/res4d"),tempInfo);
+ Data = RawData.matrix(Mask);
+ }
+ else{
+ // res4d does not exist, regress out the PE's instead
+ //meanR=mean(Data);
+ //Data=remmean(Data);
+ message("File res4d does not exists " << endl);
+ message("Re-generating this from the PEs " << endl);
+ Matrix PEs, pemat;
+ volume4D<float> tmppe;
+ int pe_ctr=2;
+ read_volume4D(tmppe,opts.glmcleanup.value() + "/stats/pe1");
+ PEs = tmppe.matrix(Mask);
+ while (fsl_imageexists(opts.glmcleanup.value() + "/stats/pe" + num2str(pe_ctr))){
+ read_volume4D(tmppe,opts.glmcleanup.value() + "/stats/pe" + num2str(pe_ctr));
+ pemat = tmppe.matrix(Mask);
+ PEs = PEs & pemat;
+ pe_ctr++;
+ }
+ outMsize(" ... reading PEs :",PEs);
+ Data = remmean(Data).t();
+ PEs = PEs.t();
+ Matrix tcs, tmpinv;
+ tmpinv = pinv(PEs.t()*PEs);
+ tcs = tmpinv * (PEs.t() * Data);
+ tcs = tcs.t();
+ outMsize(" ... created time courses :",tcs);
+ write_ascii_matrix(opts.glmcleanup.value() + "/stats/timecourses_ICAcleanup",tcs);
+ Data = Data - PEs * tcs.t();
+ Data = Data.t();
+
+ //save Data as res4d_new
+ tmppe.setmatrix(Data,Mask);
+ save_volume4D(tmppe,opts.glmcleanup.value() + "/stats/res4d_new");
+ }
+ }
+
+ void MelodicData::GLMcleanup_postproc()
+ {
+ message("Staring GLM cleanup " << endl << endl);
+ message("Creating new PEs, DOF and sigmasquareds ... " << endl);
+ // First, save all the old files
+ volume4D<float> tmpvol, old_dof;
+ Matrix pemaps;
+ volumeinfo tmpvolinfo;
+ int pe_ctr = 1;
+ while(fsl_imageexists(opts.glmcleanup.value() + "/stats/pe" + num2str(pe_ctr)))
+ {
+ read_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/pe" + num2str(pe_ctr),tmpvolinfo);
+ save_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/old_pe" + num2str(pe_ctr),tmpvolinfo);
+ if (pemaps.Storage()==0) pemaps = tmpvol.matrix(Mask);
+ else pemaps = pemaps & tmpvol.matrix(Mask);
+ pe_ctr++;
+ }
+ read_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/sigmasquareds",tmpvolinfo);
+ save_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/old_sigmasquareds",tmpvolinfo);
+ if(fsl_imageexists(opts.glmcleanup.value() + "/stats/dof"))
+ {
+ read_volume4D(old_dof,opts.glmcleanup.value() + "/stats/dof",tmpvolinfo);
+ save_volume4D(old_dof,opts.glmcleanup.value() + "/stats/old_dof",tmpvolinfo);
+ }
+ else
+ {
+ read_volume4D(old_dof,opts.glmcleanup.value() + "/stats/sigmasquareds");
+ Matrix tmpmat;
+ tmpmat = read_ascii_matrix(opts.glmcleanup.value() + "/stats/dof");
+ char callRMstr[1000];
+ ostrstream osc(callRMstr,1000);
+ osc << "mv " << opts.glmcleanup.value() << "/stats/dof " << opts.glmcleanup.value() << "/stats/old_dof " <<'\0';
+ system(callRMstr);
+ old_dof = old_dof * 0 + tmpmat(1,1);
+ }
+
+ // Read in thresholded maps
+ Matrix threshmaps;
+ int ic_ctr = 2;
+ message(" Reading in thresholded IC maps ...");
+ 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);
+ ic_ctr++;
+ }
+ message(" ... done " << endl);
+ outMsize(" threshmaps ",threshmaps);
+ outMsize(" pemaps ",pemaps);
+
+ //read in data
+ message(" Reading in " << opts.inputfname.value() << " ... ");
+ Matrix orig_data;
+ read_volume4D(tmpvol,opts.inputfname.value(),tempInfo);
+ orig_data = tmpvol.matrix(Mask);
+ orig_data = remmean(orig_data);
+ message(" ... done " << endl);
+
+ //re-regress PE maps and thresholded maps against data
+ message(" Calculating new OLS results ..." << endl);
+ Matrix allmaps, alltcs, resi;
+
+ allmaps = pemaps & threshmaps;
+
+ outMsize(" allmaps ",allmaps);
+ outMsize(" orig_data ",orig_data);
+
+ allmaps = pemaps & threshmaps;
+ alltcs = ( orig_data * allmaps.t() ) * pinv (allmaps * allmaps.t() );
+
+ outMsize(" alltcs",alltcs);
+ resi = orig_data - alltcs * allmaps;
+ pemaps = pemaps - pemaps * threshmaps.t() * pinv (threshmaps * threshmaps.t() ) * threshmaps;
+ 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);
+
+ Matrix new_dof;
+ new_dof = old_dof.matrix(Mask);
+ for(int j_ctr=1;j_ctr<=threshmaps.Ncols();j_ctr++)
+ 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);
+ message(" ... done " << endl);
+
+ //save out new GLM results
+ message(" Saving results ... ");
+ tmpvol.setmatrix(sigsq,Mask);
+ save_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/sigmasquareds");
+
+ for(int tmp_ctr=1;tmp_ctr<=pemaps.Nrows(); tmp_ctr++){
+ tmpvol.setmatrix(pemaps.Row(tmp_ctr),Mask);
+ save_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/pe" + num2str(tmp_ctr));
+ }
+ tmpvol.setmatrix(new_dof,Mask);
+ save_volume4D(tmpvol,opts.glmcleanup.value() + "/stats/dof");
+ message(" ... done " << endl);
+ message(endl << " finished GLM cleanup ! " << endl << endl);
+
+ }
+
Matrix MelodicData::calc_FFT(const Matrix& Mat)
{
Matrix res;
diff --git a/meldata.h b/meldata.h
index fdb693e7a748d1571e81489c67f8752ae81009d0..a01b2535f48a4d97f46c7eb9c375544622fbe4d9 100644
--- a/meldata.h
+++ b/meldata.h
@@ -33,6 +33,13 @@ namespace Melodic{
{
after_mm = false;
}
+
+ //destructor
+ ~MelodicData()
+ {
+ if(opts.glmcleanupmode)
+ GLMcleanup_postproc();
+ }
void save();
@@ -185,9 +192,14 @@ namespace Melodic{
char Mean_fname[1000];
+ int glm_pes;
+
void create_mask(volume4D<float> &theData, volume<float> &theMask);
// void remove_comp(Matrix& Mat, const int numIC);
void create_RXweight();
+
+ void GLMcleanup_preproc();
+ void GLMcleanup_postproc();
};
}
diff --git a/melodic.cc b/melodic.cc
index b960a7a555fdd1a0fd95f7a9320f78624509f24d..1a0a3403b6e93a5a336a142ac46c24d907545733 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -162,8 +162,6 @@ int main(int argc, char *argv[])
}
}
-
-
if( bool(opts.genreport.value()) ){
report.analysistxt();
report.PPCA_rep();
diff --git a/melodic.h b/melodic.h
index 039bb815722f11b08eeb64fcb2942af30e0af34d..fd4841120de95af942556dfe08c3d4883fa2bb14 100644
--- a/melodic.h
+++ b/melodic.h
@@ -11,6 +11,7 @@
#ifndef __MELODIC_h
#define __MELODIC_h
+//#include "fslio/fslio.h"
#include<sstream>
@@ -32,7 +33,7 @@
namespace Melodic{
-const string version = "ln(12) beta";
+const string version = "2.7.1";
}
diff --git a/meloptions.cc b/meloptions.cc
index eb9eb401fc0cef52f512782a8798754898fdcd44..2ef911c8445e532b59cbbe61b4db249e02104c81 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -28,10 +28,11 @@ MelodicOptions* MelodicOptions::gopt = NULL;
void MelodicOptions::parse_command_line(int argc, char** argv, Log& logger, const string &p_version)
{
//set version number and some other stuff
- version = p_version;
- filtermode = false;
- explicitnums = false;
- logfname = string("melodic.log");
+ version = p_version;
+ glmcleanupmode = false;
+ filtermode = false;
+ explicitnums = false;
+ logfname = string("melodic.log");
// work out the path to the $FSLDIR/bin directory
if(getenv("FSLDIR")!=0){
@@ -115,7 +116,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
}
if (filter.value().length()>0){
if (filtermix.value().length()<0){
- cerr << "ERROR:: no mixing matrix for filtering (use --mix='filename') \n\n";
+ cerr << "ERROR:: no mixing matrix for filtering specified (use --mix='filename') \n\n";
print_usage(argc,argv);
exit(2);
} else {
@@ -123,6 +124,24 @@ MelodicOptions* MelodicOptions::gopt = NULL;
varnorm.set_T(false);
}
}
+ if (glmcleanup.value().length()>0){
+ glmcleanupmode = true;
+ output_MMstats.set_T(true);
+ perf_mm.set_T(true);
+ maskfname.set_T(glmcleanup.value() + "/mask");
+ // varnorm.set_T(false);
+ logdir.set_T(glmcleanup.value() + "/GLM_cleanup.ica");
+ if (!fsl_imageexists(inputfname.value()) ||
+ !fsl_imageexists(glmcleanup.value() + "/stats/cope1")){
+ cerr << "ERROR:: one or more file in the specified feat directory does not exist \n\n";
+ exit(2);
+ }
+ // if (FslFileExists(string(glmcleanup.value() + "/stats/res4d").c_str())){
+ // message("\n Warning:: res4d exists - using this as input file\n\n");
+
+ // inputfname.set_T(string(glmcleanup.value() + "/stats/res4d"));
+ //}
+ }
if (threshold.value()<=0){
use_mask.set_T(false);
}
diff --git a/meloptions.h b/meloptions.h
index a51566d89eeeebf84c6ef8f77887712f2b457d96..a2e8c6a44798b174ef70d40b5efe283af5b54931 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -38,6 +38,7 @@ class MelodicOptions {
string binpath;
string logfname;
bool filtermode;
+ bool glmcleanupmode;
bool explicitnums;
Option<string> logdir;
@@ -68,6 +69,7 @@ class MelodicOptions {
Option<string> ICsfname;
Option<string> filtermix;
Option<string> filter;
+ Option<string> glmcleanup;
Option<bool> genreport;
Option<float> tr;
@@ -199,6 +201,9 @@ class MelodicOptions {
filter(string("-f,--filter"), string(""),
string("component numbers to remove\n "),
false, requires_argument),
+ glmcleanup(string("--cleanup"), string(""),
+ string("<feat dir>\n "),
+ false, requires_argument),
genreport(string("--report"), false,
string("generate Melodic web report"),
false, no_argument),
@@ -314,6 +319,7 @@ class MelodicOptions {
options.add(ICsfname);
options.add(filtermix);
options.add(filter);
+ options.add(glmcleanup);
options.add(genreport);
options.add(tr);
options.add(logPower);