GLM cleanup

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

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