From 31b8393737173dad42dfc70a51209d6b13e196f3 Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Fri, 27 Jul 2007 21:00:23 +0000
Subject: [PATCH] *** empty log message ***

---
 meldata.cc   | 34 ++++++++++++++++++++++------------
 meldata.h    | 19 ++++++++++++++++---
 melica.cc    |  2 +-
 meloptions.h | 43 ++++++++++++++++++++++++-------------------
 melreport.cc | 24 ++++++++++++------------
 melreport.h  |  2 +-
 test.cc      |  8 ++++++++
 7 files changed, 84 insertions(+), 48 deletions(-)

diff --git a/meldata.cc b/meldata.cc
index 89fffd3..ef4421f 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -203,7 +203,7 @@ namespace Melodic{
       calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
     }  
    
-    if(numfiles < 2){
+    if(numfiles == 1){
       Data = alldat;
       Matrix tmp = Identity(Data.Nrows());
       DWM.push_back(tmp);
@@ -246,20 +246,18 @@ namespace Melodic{
     //initialize Mean
     read_volume(Mean,opts.inputfname.value().at(0),tempInfo);
 
-    //save first image
-    tmpnam(Mean_fname); // generate a tmp name
-    save_volume(Mean,Mean_fname);    
-
     //create mask
     create_mask(Mask);
 
-    // clean /tmp
-    char callRMstr[1000];
-    ostrstream osc(callRMstr,1000);
-    osc  << "rm " << string(Mean_fname) <<"*  " << '\0';
-
-    system(callRMstr);
- 
+		//setup background image
+		if(opts.bgimage.value()>""){
+			read_volume(background,opts.bgimage.value());
+      if(!samesize(Mean,background)){
+        cerr << "ERROR:: background image and data have different dimensions  \n\n";
+        exit(2);
+      }
+		}
+		
     if(!samesize(Mean,Mask)){
       cerr << "ERROR:: mask and data have different dimensions  \n\n";
       exit(2);
@@ -633,6 +631,10 @@ namespace Melodic{
     else{
       if(opts.perf_bet.value() && opts.use_mask.value()){ //use BET
 				message("Create mask ... ");
+    		//save first image
+    		tmpnam(Mean_fname); // generate a tmp name
+    		save_volume(Mean,Mean_fname);    
+
 				// set up all strings
 				string BET_outputfname = string(Mean_fname)+"_brain";
 
@@ -652,6 +654,14 @@ namespace Melodic{
 	
 				// read back the Mask file   
 				read_volume(theMask,Mask_fname);
+				
+				// clean /tmp
+		    char callRMstr[1000];
+		    ostrstream osc(callRMstr,1000);
+		    osc  << "rm " << string(Mean_fname) <<"*  " << '\0';
+
+		    system(callRMstr);
+		   
 				message("done" << endl);
       }  
       else{
diff --git a/meldata.h b/meldata.h
index 4ebec20..623db92 100644
--- a/meldata.h
+++ b/meldata.h
@@ -133,6 +133,14 @@ namespace Melodic{
       inline volume<float>& get_mean() {return Mean;}
       inline void set_mean(volume<float>& Arg) {Mean = Arg;}
    
+      inline volume<float>& get_bg() {
+				if(opts.bgimage.value()>"")
+					return background;
+				else
+					return Mean;
+			}
+      inline void set_bg(volume<float>& Arg) {background = Arg;}
+   
       inline Matrix& get_Data() {return Data;}
       inline void set_Data(Matrix& Arg) {Data = Arg;}
     
@@ -143,13 +151,16 @@ namespace Melodic{
       inline void set_ICstats(Matrix& Arg) {ICstats = Arg;}
      
       inline Matrix& get_EVP() {return EVP;}
-      inline void set_EVP(Matrix& Arg) {EVP = Arg;}
+      inline void set_EVP(Matrix& Arg) {if(EVP.Storage()==0)
+																					EVP = Arg;}
       
       inline Matrix& get_EV() {return EV;}
-      inline void set_EV(Matrix& Arg) {EV = Arg;}
+      inline void set_EV(Matrix& Arg) {if(EV.Storage()==0)
+																				 EV = Arg;}
 
       inline Matrix& get_PPCA() {return PPCA;}
-      inline void set_PPCA(Matrix& Arg) {PPCA = Arg;}
+      inline void set_PPCA(Matrix& Arg) {if(PPCA.Storage()==0)
+																					 PPCA = Arg;}
 
       inline Matrix& get_stdNoisei() {return stdNoisei;}
       inline void set_stdNoisei(Matrix& Arg) {stdNoisei = Arg;}
@@ -211,6 +222,8 @@ namespace Melodic{
       Matrix stdNoisei;
       volume<float> Mask;
       volume<float> Mean;
+      volume<float> background;
+
       Matrix Data;
       Matrix PPCA;
       Matrix jointCC;
diff --git a/melica.cc b/melica.cc
index ee4fb3f..1379340 100644
--- a/melica.cc
+++ b/melica.cc
@@ -115,7 +115,7 @@ namespace Melodic{
 	  		outMsize(" redUMM ", redUMM);
 			}
     } while((itt_ctr2 < 3 || itt_ctr >= opts.maxNumItt.value()) && 
-			(opts.approach.value() ==  string("tica")) && cum_itt < 25*opts.maxNumItt.value());
+			(opts.approach.value() ==  string("tica")) && cum_itt < 100*opts.maxNumItt.value());
 
     if(itt_ctr>=opts.maxNumItt.value()){
       cerr << "  No convergence after " << cum_itt  <<" steps "<<endl;
diff --git a/meloptions.h b/meloptions.h
index 81d7bc5..3ca7a3f 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -75,6 +75,7 @@ class MelodicOptions {
   	Option<string> filter; 
 
   	Option<bool>   genreport;
+		Option<string> bgimage;
   	Option<float>  tr;
   	Option<bool>   logPower;
 
@@ -172,7 +173,7 @@ class MelodicOptions {
    pca_est(string("--dimest"), string("lap"),
 	   string("use specific dim. estimation technique: lap, bic, mdl, aic, mean (default: lap)"), 
 	   false, requires_argument),
-   joined_whiten(string("--separate_whiten"), true,
+   joined_whiten(string("--sep_whiten"), true,
 	   string("switch on separate whitening"), 
 	   false, no_argument),
    numICs(string("-n,--numICs"), -1,
@@ -188,10 +189,10 @@ class MelodicOptions {
 	   string("switch off variance normalisation"), 
 	   false, no_argument),
    pbsc(string("--pbsc"), false,
-	   string("switch off conversion to percent BOLD signal change"), 
+	   string("        switch off conversion to percent BOLD signal change"), 
 	   false, no_argument),
    pspec(string("--pspec"), false,
-	   string("switch on conversion to powerspectra"), 
+	   string("        switch on conversion to powerspectra"), 
 	   false, no_argument),
    segment(string("--covarweight"), string(""),
 	   string("voxel-wise weights for the covariance matrix (e.g. segmentation information)"),
@@ -226,12 +227,15 @@ class MelodicOptions {
    smodename(string("--smode"),  string(""),
 	   string("\tmatrix of session modes for report generation"), 
 	   false, requires_argument),
-   filter(string("-f,--filter"),  string(""),
-	   string("component numbers to remove\n "), 
-	   false, requires_argument),
+	 filter(string("-f,--filter"),  string(""),
+		 string("component numbers to remove\n "), 
+		 false, requires_argument),
    genreport(string("--report"), false,
 	   string("generate Melodic web report"), 
 	   false, no_argument),
+	 bgimage(string("--bgimage"),  string(""),
+		 string("specify background image for report (default: mean image)\n "), 
+		 false, requires_argument),
    tr(string("--tr"),  0.0,
 	   string("\tTR in seconds"), 
 	   false, requires_argument),
@@ -239,25 +243,25 @@ class MelodicOptions {
 	   string("calculate log of power for frequency spectrum\n"), 
 	    false, no_argument),
 	 fn_Tdesign(string("--Tdes"), string(""),
-	   string("design matrix across time-domain"),
+	   string("        design matrix across time-domain"),
 		 false, requires_argument),
 	 fn_Tcon(string("--Tcon"), string(""),
-		 string("t-contrast matrix across time-domain"),
+		 string("        t-contrast matrix across time-domain"),
 		 false, requires_argument),
    fn_TconF(string("--Tconf"), string(""),
-		 string("F-contrast matrix across time-domain"),
+		 string("        F-contrast matrix across time-domain"),
 		 false, requires_argument),
    fn_Sdesign(string("--Sdes"), string(""),
-		 string("design matrix across subject-domain"),
+		 string("        design matrix across subject-domain"),
 		 false, requires_argument),
 	 fn_Scon(string("--Scon"), string(""),
-		 string("t-contrast matrix across subject-domain"),
+		 string("        t-contrast matrix across subject-domain"),
 		 false, requires_argument),	
 	 fn_SconF(string("--Sconf"), string(""),
-		 string("F-contrast matrix across subject-domain"),
-		 false, requires_argument),	
+		 string("        F-contrast matrix across subject-domain"),
+		 false, requires_argument,false),	
    output_all(string("--Oall"),  false,
-	   string("output everything"), 
+	   string("        output everything"), 
 	   false, no_argument),
    output_unmix(string("--Ounmix"),  false,
 	   string("output unmixing matrix"), 
@@ -290,7 +294,7 @@ class MelodicOptions {
 	   string("prints this help message"), 
 	   false, no_argument),
    debug(string("--debug"),  false,
-	   string("switch on debug messages"), 
+	   string("        switch on debug messages"), 
 	   false, no_argument),
    guessfname(string("-g,--guess"), string(""),
 	   string("file name of guess of mixing matrix"), 
@@ -332,16 +336,16 @@ class MelodicOptions {
    options(string(""), 
 	   string(" melodic -i <filename> <options>")+
 	   string("\n \t \t to run melodic")+
-	   string("\n melodic -i <filename> --mix=melodic_mix")+
-	   string(" --filter=\"string of component numbers\"")+
-	   string("\n \t \t to remove estimated ICs from input")+
+//	   string("\n melodic -i <filename> --mix=melodic_mix")+
+//	   string(" --filter=\"string of component numbers\"")+
+//	   string("\n \t \t to remove estimated ICs from input")+
 	   string("\n melodic -i <filename> --ICs=melodic_IC")+
 	   string(" --mix=melodic_mix <options>")+
 	   string("\n \t \t to run Mixture Model based inference on estimated ICs")+
 	   string("\n melodic --help "))
    {
      try {  
-            options.add(logdir);
+      options.add(logdir);
 	    options.add(inputfname);
 	    options.add(outputfname);
 	    options.add(guessfname);
@@ -372,6 +376,7 @@ class MelodicOptions {
 	    options.add(smodename);
 	    options.add(filter);
 	    options.add(genreport);
+			options.add(bgimage);
 	    options.add(tr);
 	    options.add(logPower);
 			options.add(fn_Tdesign);
diff --git a/melreport.cc b/melreport.cc
index 3597b78..a61db40 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -1,6 +1,6 @@
 /*  MELODIC - Multivariate exploratory linear optimized decomposition into 
     independent components
-    
+    melodat.get_bg()
     melreport.cc - report generation
 
     Christian F. Beckmann, FMRIB Image Analysis Group
@@ -79,9 +79,9 @@ namespace Melodic{
 				float map2max = std::min((map2 + binarise(tempVol[0],float(0.0), 
 						  		tempVol[0].max()) * map2.min()).max(),float(-0.01));
 	
-    		newpic.overlay(newvol, melodat.get_mean(), map1, map2, 
-		       		melodat.get_mean().percentile(0.02),
-		       		melodat.get_mean().percentile(0.98),
+    		newpic.overlay(newvol, melodat.get_bg(), map1, map2, 
+		       		melodat.get_bg().percentile(0.02),
+		       		melodat.get_bg().percentile(0.98),
 		       		map1min, map1max, map2max, map2min, 
 		       		0, 0, &melodat.tempInfo);
 				char instr[10000];
@@ -320,9 +320,9 @@ namespace Melodic{
 	    	//  << map2min << " " << map2max << endl;
 	    
 	    	//	    if(map1.max()-map1.min()>0.01)
-	    	newpic.overlay(newvol, melodat.get_mean(), map1, map2, 
-			   	melodat.get_mean().percentile(0.02),
-			   	melodat.get_mean().percentile(0.98),
+	    	newpic.overlay(newvol, melodat.get_bg(), map1, map2, 
+			   	melodat.get_bg().percentile(0.02),
+			   	melodat.get_bg().percentile(0.98),
 			   	map1min, map1max, map2max, map2min, 
 			   	0, 0, &melodat.tempInfo);
 	    
@@ -411,7 +411,7 @@ namespace Melodic{
 	  		//float map2max = (map2 + binarise(tempVol[0],float(0.0), 
 	  		//		tempVol[0].max()) * map2.min()).robustmax();
 	
-	  		newpic.overlay(newvol, melodat.get_mean(), map1, map2, 
+	  		newpic.overlay(newvol, melodat.get_bg(), map1, map2, 
 			 		float(0.0),
 			 		float(0.0),
 			 		float(0.01), map1max, float(-0.01), map2min, 
@@ -448,9 +448,9 @@ namespace Melodic{
 	  		volume<float> newvol; 
 	  		miscpic newpic;
 
-	  		newpic.overlay(newvol, melodat.get_mean(), map, map, 
-			 		melodat.get_mean().percentile(0.02),
-			 		melodat.get_mean().percentile(0.98),
+	  		newpic.overlay(newvol, melodat.get_bg(), map, map, 
+			 		melodat.get_bg().percentile(0.02),
+			 		melodat.get_bg().percentile(0.98),
 			 		float(0.1), float(1.0), float(0.0), float(0.0),
 			 		0, 0, &melodat.tempInfo);
 
@@ -582,7 +582,7 @@ namespace Melodic{
 	  		//float map2max = (map2 + binarise(tempVol[0],float(0.0), 
 	  		//		tempVol[0].max()) * map2.min()).robustmax();
 	
-	  		newpic.overlay(newvol, melodat.get_mean(), map1, map2, 
+	  		newpic.overlay(newvol, melodat.get_bg(), map1, map2, 
 			 		float(0.0),
 			 		float(0.0),
 			 		float(0.01), map1max, float(-0.01), map2min, 
diff --git a/melreport.h b/melreport.h
index 87567f2..68662e9 100644
--- a/melreport.h
+++ b/melreport.h
@@ -42,7 +42,7 @@ namespace Melodic{
 				logger(plogger){
 	  			if( bool(opts.genreport.value()) ){
 	  				const time_t tmptime = time(NULL);
- 
+ 					
 	  				report.makeDir(logger.appendDir("report"),"00index.html");
 	  				report << "<HTML>" << endl << endl 
 		 					<< "<TITLE>MELODIC Report</TITLE>"<< endl <<endl
diff --git a/test.cc b/test.cc
index 7fc130b..d0e0165 100644
--- a/test.cc
+++ b/test.cc
@@ -81,6 +81,14 @@ int do_work(int argc, char* argv[]) {
 	cerr << weights << endl <<endl << Rows << endl;
 
 	cerr << SP(weights,pow(Rows,-1)) << endl;
+	
+	vector<int> vi;
+	for(int ctr=1;ctr<=5;ctr++)
+		vi.push_back(ctr);
+		
+	for(int ctr=0;ctr<5;ctr++)
+		cerr << " vi["<<ctr<<"] = " <<vi.at(ctr)<<endl; 
+	
 	return 0;
 }
 
-- 
GitLab