*** empty log message ***

fsl_regfilt.cc

@@ -65,7 +65,7 @@ using namespace std;
 		string("        switch on aggressive filtering (full instead of partial regression)"),
 		false, no_argument);
 	Option<bool> perfvn(string("--vn"),FALSE,
-		string("        perfrom variance-normalisation on data"),
+		string("        perform variance-normalisation on data"),
 		false, no_argument);
 	Option<int> help(string("-h,--help"), 0,
 		string("display this help text"),

melgmix.cc

@@ -597,10 +597,12 @@ namespace Melodic{
 
     props = props / sum(props,2).AsScalar();
     add_params(means,vars,props,logprobY,MDL,Evi,true);
-    
+	
     probmap << SP(sum(tmp1.Rows(2,tmp1.Nrows()),1),
 		  pow(sum(tmp1,1),-1));
 
+	dbgmsg("   mu: " << means << "  sig: " << vars << " prop: " << props << endl);
+
    	if(props(1)<0.4 ){
       //set up GMM
       message("    try Gaussian Mixture Model " << endl);
@@ -618,11 +620,11 @@ namespace Melodic{
       IntSize = Dmax / nummix;
       means(1)=mean(data,2).AsScalar(); 
       for (int ctr=2; ctr <= means.Ncols(); ctr++){
-				means(ctr) =  Dmax - (ctr - 1.5) * IntSize; 
+		means(ctr) =  Dmax - (ctr - 1.5) * IntSize; 
       }
       means(2)=means(1)+sqrt(vars(1));
       if(nummix>2)
-				means(3)=means(1)-sqrt(vars(1));
+		means(3)=means(1)-sqrt(vars(1));
       
       fit(string("GMM"));
     }

melgmix.h

@@ -145,6 +145,14 @@ namespace Melodic{
         probmap = tempVol.matrix(Mask);
       }
 
+	  inline Matrix get_params(){
+		Matrix tmp = zeros(3,means.Ncols());
+		tmp.Row(1) = means;
+		tmp.Row(2) = vars;
+		tmp.Row(3) = props;
+		return tmp;
+      }
+
       double datamean;
       double datastdev;
 

melodic.cc

@@ -68,7 +68,7 @@ int main(int argc, char *argv[]){
     //set up report object     
     MelodicReport report(melodat,opts,logger);
    
-    if (opts.filtermode || opts.filtermix.value().length()>0){
+    if (opts.filtermode || opts.filtermix.value().length()>0 || opts.ICsfname.value().length()>0){
       if(opts.filtermode){ // just filter out some noise from a previous run
     		melodat.setup();
 				if(opts.debug.value())
@@ -237,13 +237,18 @@ void mmonly(Log& logger, MelodicOptions& opts,
     message(" done" << endl);
   }
 
-  message("Reading mixing matrix " << opts.filtermix.value() << " ... ");
-  mixMatrix = read_ascii_matrix(opts.filtermix.value());
-  if (mixMatrix.Storage()<=0) {
-    cerr <<" Please specify the mixing matrix correctly" << endl;
-    exit(2);
+  if(opts.filtermix.value().length() > 0){
+    message("Reading mixing matrix " << opts.filtermix.value() << " ... ");
+    mixMatrix = read_ascii_matrix(opts.filtermix.value());
+    if (mixMatrix.Storage()<=0) {
+      cerr <<" Please specify the mixing matrix correctly" << endl;
+      exit(2);
+    }
+    message(" done" << endl);
+  }else{
+	mixMatrix=unifrnd(ICs.Nrows()+1,ICs.Nrows());	
   }
-  message(" done" << endl);
+
 
   if(opts.smodename.value().length() > 0){
     message("Reading matrix of subject modes: " << opts.smodename.value());
@@ -300,7 +305,8 @@ Matrix mmall(Log& logger, MelodicOptions& opts,MelodicData& melodat, MelodicRepo
     message("  IC map " << ctr << " ... "<< endl;);
     
     Matrix ICmap;
-		dbgmsg(" stdNoisei max : "<< melodat.get_stdNoisei().Maximum() <<" "<< melodat.get_stdNoisei().Minimum() << endl);
+    if(melodat.get_stdNoisei().Storage()>0)
+ 		dbgmsg(" stdNoisei max : "<< melodat.get_stdNoisei().Maximum() <<" "<< melodat.get_stdNoisei().Minimum() << endl);
 
     if(opts.varnorm.value()&&melodat.get_stdNoisei().Storage()>0){
       ICmap = SP(melodat.get_IC().Row(ctr),diagvals(ctr)*melodat.get_stdNoisei());
@@ -322,9 +328,12 @@ Matrix mmall(Log& logger, MelodicOptions& opts,MelodicData& melodat, MelodicRepo
     mixmod.fit("GGM");
 
     if(opts.output_MMstats.value()){
-      message("   saving probability map:");
+      message("   saving probability map:  ");
       melodat.save4D(mixmod.get_probmap(),
 		     string("stats/probmap_")+num2str(ctr));
+%%	  message("   saving mixture model fit:");
+%%	  melodat.saveascii(mixmod.get_params(),
+%%			 string("stats/MMstats_")+num2str(ctr));  
     }
 
     //re-scale spatial maps to mean 0 for nht
@@ -415,7 +424,7 @@ Matrix mmall(Log& logger, MelodicOptions& opts,MelodicData& melodat, MelodicRepo
     }
   }
   
-  if(!opts.filtermode&&opts.filtermix.value().length()==0){
+  if(!opts.filtermode&&opts.ICsfname.value().length()==0){
     //now safe new data
     //    bool what = opts.verbose.value();
     //opts.verbose.set_T(false);