weighted mask option

fslmeants.cc

@@ -72,6 +72,10 @@ Option<int> order(string("--order"), 1,
 Option<bool> transpose(string("--transpose"), false,
 		  string("        output results in transpose format (one row per voxel/mean)"),
 		  false, no_argument);
+Option<bool> weightedMask(string("-w"), false,
+		  string("output weighted mean, using mask values as weights, and exit."),
+		  false, no_argument);
+
 
 int nonoptarg;
 
@@ -96,6 +100,7 @@ int main(int argc,char *argv[])
   options.add(transpose);
   options.add(verbose);
   options.add(help);
+  options.add(weightedMask);
   
   nonoptarg = options.parse_command_line(argc, argv);
   
@@ -122,6 +127,32 @@ int main(int argc,char *argv[])
     mask = 1.0;
   }
 
+
+
+  if ( weightedMask.value() ) {
+    Matrix meanwts(vin.tsize(),1);
+    for(int t=0;t<vin.tsize();t++) {    
+      double mean(0);
+      double meanWeight(0);
+      for (int z=mask.minz(); z<=mask.maxz(); z++) 
+	for (int y=mask.miny(); y<=mask.maxy(); y++) 
+	  for (int x=mask.minx(); x<=mask.maxx(); x++) {
+	    mean+=mask(x,y,z)*vin(x,y,z,t);	
+	    meanWeight+=mask(x,y,z);
+	  }
+      meanwts(t+1,1)=(mean/meanWeight);
+    }
+
+    if ( transpose.value() ) 
+      meanwts=meanwts.t();
+    if ( outmat.set() ) 
+      write_ascii_matrix(meanwts,outmat.value());
+    else 
+      cout << meanwts << endl;
+    return(0);
+  }
+
+
   int nlabs=1;
   if (labelname.set() && showall.unset()) {
     read_volume(label,labelname.value());