added GLM cleanup

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;

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();
          
     };
 }

melodic.cc

@@ -162,8 +162,6 @@ int main(int argc, char *argv[])
 	  }
 	}		 
 
-
-
 	if( bool(opts.genreport.value()) ){
 	  report.analysistxt();
 	  report.PPCA_rep();

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";  
 
 }
 

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);
   }

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);