new version with cumulative masking

fslstats.cc

@@ -144,7 +144,6 @@ int generate_masks(const volume4D<float> &mask, volume4D<float> &masknz, const v
   return 0;
 }
 
-
 int fmrib_main_float(int argc, char* argv[]) 
 {
 
@@ -153,15 +152,10 @@ int fmrib_main_float(int argc, char* argv[])
   cout.setf(ios::left, ios::adjustfield); 
   cout.precision(6);  
 
-  volume4D<float> vin, vol, mask, masknz;
+  volume4D<float> vol, mask, masknz;
   read_volume4D(vol,argv[1]);
 
   float lthr=0, uthr=0;  // these initial values are not used
-  if (masks_used) {
-    vin = vol;
-    generate_masks(mask,masknz,vin,lthr,uthr);
-    vol = vin * mask;
-  }
 
   int narg=2;
   string sarg;
@@ -216,12 +210,10 @@ int fmrib_main_float(int argc, char* argv[])
       cout << MISCMATHS::round(coord(1)) << " " << 
 	MISCMATHS::round(coord(2)) << " " << MISCMATHS::round(coord(3)) << " ";
     } else if (sarg=="-w") {
-      if (!masks_used) { 
-	if (vin.nvoxels()<1) { vin = vol; }
-	masks_used=true;
-	generate_masks(mask,masknz,vin,lthr,uthr); 
-	vol = vin * mask; 
-      }
+	if (masknz.nvoxels()<1) { //Need to generate non-zeromask 
+	  generate_masks(mask,masknz,vol,lthr,uthr); 
+	  vol*=mask; 
+	}
       int xmin=masknz.maxx(),xmax=masknz.minx(),ymin=masknz.maxy(),ymax=masknz.miny(),zmin=masknz.maxz(),zmax=masknz.minz(),tmin=masknz.maxt(),tmax=masknz.mint();
       
       for(int t=masknz.mint();t<=masknz.maxt();t++) {
@@ -260,13 +252,11 @@ int fmrib_main_float(int argc, char* argv[])
       if (zmin>zmax) { int tmp=zmax;  zmax=zmin;  zmin=tmp; }
       // now output nifti coords
       cout << xmin << " " << 1+xmax-xmin << " " << ymin << " " << 1+ymax-ymin << " " << zmin << " " << 1+zmax-zmin << " " << tmin << " " << 1+tmax-tmin << " ";
-    } else if (sarg=="-e") {
-      if (!masks_used) { 
-	if (vin.nvoxels()<1) { vin = vol; }
-	masks_used=true;
-	generate_masks(mask,masknz,vin,lthr,uthr); 
-	vol = vin * mask; 
-      }
+      } else if (sarg=="-e") {
+	if (mask.nvoxels()<1) {
+	  generate_masks(mask,masknz,vol,lthr,uthr); 
+	  vol*=mask; 
+	}
       ColumnVector hist;
       int nbins=1000;
       double entropy=0;
@@ -279,7 +269,7 @@ int fmrib_main_float(int argc, char* argv[])
       }
       entropy /= log((double) nbins);
       cout << entropy << " ";
-    } else if (sarg=="-E") { 
+      } else if (sarg=="-E") { 
       ColumnVector hist;
       int nbins=1000;
       double entropy=0;
@@ -313,19 +303,11 @@ int fmrib_main_float(int argc, char* argv[])
    	  mask.addvolume(mask[mask.maxt()]);
         }
       }
-      if (!masks_used) {
-	masks_used=true;
-	vin = vol;
-      }
-      float th= 0.5;
-      if (th!=0) {
-        mask.binarise(th);
-      } else {
-        mask.binarise(1);
-      }
-      generate_masks(mask,masknz,vin);
-        if (mask.tsize()!=1) vol=vin*mask; 
-	else vol=vin*mask[0];
+      
+      mask.binarise(0.5);
+      generate_masks(mask,masknz,vol);
+        if (mask.tsize()!=1) vol*=mask; 
+	else vol*=mask[0];
     } else if (sarg=="-l") {
       narg++;
       if (narg<argc) {
@@ -335,12 +317,8 @@ int fmrib_main_float(int argc, char* argv[])
 	exit(2);
       }
       lthr_used=true;
-      if (!masks_used) {
-	masks_used=true;
-	vin = vol;
-      }
-      generate_masks(mask,masknz,vin,lthr,uthr);
-      vol = vin * mask;
+      generate_masks(mask,masknz,vol,lthr,uthr);
+      vol*=mask;
     } else if (sarg=="-u") {
       narg++;
       if (narg<argc) {
@@ -350,12 +328,8 @@ int fmrib_main_float(int argc, char* argv[])
 	exit(2);
       }
       uthr_used=true;
-      if (!masks_used) {
-	masks_used=true;
-	vin = vol;
-      }
-      generate_masks(mask,masknz,vin,lthr,uthr);
-      vol = vin * mask;
+      generate_masks(mask,masknz,vol,lthr,uthr);
+      vol*=mask;
     } else if (sarg=="-a") {
       vol = abs(vol);
     } else if (sarg=="-v") {
@@ -436,11 +410,9 @@ int fmrib_main_float(int argc, char* argv[])
     	cerr << "Percentile must be between 0 and 100" << endl;
     	exit(1);
       }
-      if (!masks_used) {
-	if (vin.nvoxels()<1) { vin = vol; }
-	masks_used=true;
-	generate_masks(mask,masknz,vin,lthr,uthr); 
-	vol = vin * mask; 
+      if (mask.nvoxels()<1) {
+	generate_masks(mask,masknz,vol,lthr,uthr); 
+	vol*=mask; 
       }
       cout << vol.percentile((float) n/100.0,masknz) << " ";
     } else if (sarg=="-h") {