new option to fsl_sbca

fsl_sbca.cc

@@ -37,7 +37,7 @@ using namespace std;
 		string("output file base name"),
 		true, requires_argument);
   Option<string> fnseed(string("-s,--seed"), string(""),
-		string("seed voxel coordinate or file name of seed mask (3D file)"),
+		string("seed voxel coordinate or file name of seed mask (3D/4D file)"),
 		true, requires_argument);
   Option<string> fntarget(string("-t,--target"), string(""),
 		string("file name of target mask(s) (3D or 4D file)"),
@@ -48,6 +48,9 @@ using namespace std;
   Option<string> fnconf(string("--conf"), string(""),
 		string("        file name (or comma-separated list of file name) for confound ascii txt files"),
 		false, requires_argument);
+  Option<string> fnseeddata(string("--seeddata"), string(""),
+		string("file name of 4D data file for the seed"),
+		false, requires_argument);		
   Option<bool> map_bin(string("--bin"), false,
 		string("        binarise spatial maps prior to calculation of time courses"),
 		false, no_argument);
@@ -177,9 +180,12 @@ void create_mask(string what){
 void create_seeds(string what){
   if(debug.value()) 
 	cerr << "DBG: in create_seeds" << endl;
-  
+
+
+  volume4D<float> tmp_vol;
   if(fsl_imageexists(what)){
-	read_volume(maskS,what);
+	read_volume4D(tmp_vol,what);
+	maskS = tmp_vol[0];
     if(!samesize(orig_data[0],maskS)){
 	  cerr << "ERROR: Seed mask image does not match input image" << endl;
       exit(1);
@@ -187,17 +193,36 @@ void create_seeds(string what){
   }  
   else
     create_mask(what);	
-  
-  volume4D<float> tmp_mask;
-  tmp_mask.addvolume(maskS);
-  maskS.binarise(1e-8);
 
-  Matrix scales = tmp_mask.matrix(maskS);
-  seeds = remmean(orig_data.matrix(maskS),1);	
-  if(!map_bin.value())
-    seeds = SP(seeds, ones(seeds.Nrows(),1) * scales);
+  if(tmp_vol.tsize() > 1 && tmp_vol.tsize() == orig_data.tsize()){
+	maskS *= tmp_vol[0] / tmp_vol.tsize();
+	for(int ctr=1; ctr < tmp_vol.tsize(); ctr++)
+	  maskS += tmp_vol[ctr] * tmp_vol[ctr] / tmp_vol.tsize(); 
+    maskS.binarise(1e-8);
+	seeds = remmean(tmp_vol.matrix(maskS),1);
+  }
+  else{    
+    volume4D<float> tmp_mask;
+    tmp_mask.addvolume(maskS);
+    maskS.binarise(1e-8);
+
+    if(fnseeddata.value()>"" && fsl_imageexists(fnseeddata.value())){
+      volume4D<float> seed_data;
+	  if(verbose.value())
+		cout << " Reading input data for seeds " << fnseeddata.value() << endl;
+      read_volume4D(seed_data,fnseeddata.value());
+      seeds = remmean(seed_data.matrix(maskS),1);
+    }else{
+	  seeds = remmean(orig_data.matrix(maskS),1);	
+      if(!map_bin.value()){
+    	Matrix scales = tmp_mask.matrix(maskS);	
+        seeds = SP(seeds, ones(seeds.Nrows(),1) * scales);
+      }
+    }
+  }
 	
   voxels = seeds.Ncols();
+
   if(debug.value()){
 	cerr << "DBG: " << voxels << " voxels" << endl;
 	cerr << "DBG: seeds matrix is " << seeds.Nrows() << " x " << seeds.Ncols() << endl;
@@ -487,6 +512,7 @@ int main(int argc,char *argv[]){
 		options.add(fntarget);
 		options.add(regress_only);
 		options.add(fnconf);
+		options.add(fnseeddata);
 		options.add(map_bin);
 		options.add(tc_mean);
 		options.add(tc_order);