From 97b20b576a4268aff38161dc5fde326eff5796e7 Mon Sep 17 00:00:00 2001
From: Matthew Webster <matthew.webster@ndcn.ox.ac.uk>
Date: Mon, 20 Jan 2020 16:35:09 +0000
Subject: [PATCH] Changes as requested by MC

---
 fslselectvols.cc | 70 +++++++++++++++++++++++-------------------------
 1 file changed, 33 insertions(+), 37 deletions(-)

diff --git a/fslselectvols.cc b/fslselectvols.cc
index caa1742..9577272 100644
--- a/fslselectvols.cc
+++ b/fslselectvols.cc
@@ -16,7 +16,7 @@ using namespace Utilities;
 
 
 int parse_filterstring(vector<int>& comps,const string& vols){
-  int ctr=0;    
+  int ctr=0;
   char *p;
   char t[1024];
   const char *discard = ", [];{(})abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ~!@#$%^&*_-=+|\':><./?";
@@ -31,7 +31,7 @@ int parse_filterstring(vector<int>& comps,const string& vols){
     p=strtok(NULL,discard);
     if(p){
       ctr = atoi(p);
-      if(ctr>0){
+      if(ctr>=0){
 	comps.push_back(ctr);
       }
     }
@@ -46,7 +46,7 @@ bool file_exists (const std::string& name) {
     } else {
         f.close();
         return false;
-    }   
+    }
 }
 
 int get_vols(Matrix& id_vols,const string& vols){
@@ -105,60 +105,56 @@ int main(int argc,char *argv[]){
     if ( (help.value()) || (!options.check_compulsory_arguments(true)) ){
       options.usage();
       exit(EXIT_FAILURE);
-    }
-    else{
-      
-      volume4D<float> data;
-      volume4D<float> subdata;
+    } else {
+      volume4D<float> data, subdata;
       read_volume4D(data,ifilename.value());
 
       Matrix list;
       get_vols(list,vols.value());
 
-      bool meancalc=calc_mean.value(), varcalc=calc_var.value();
+      bool meancalc( calc_mean.value() ), varcalc( calc_var.value() );
 
       if(!meancalc && !varcalc)
-	subdata.reinitialize(data.xsize(),data.ysize(),data.zsize(),list.Ncols());
+	      subdata.reinitialize(data.xsize(),data.ysize(),data.zsize(),list.Ncols());
       else
-	subdata.reinitialize(data.xsize(),data.ysize(),data.zsize(),1);
+	      subdata.reinitialize(data.xsize(),data.ysize(),data.zsize(),1);
       subdata=0;
       copybasicproperties(data,subdata);
-      
-      float v,v2;
-      for(int z=0;z<data.zsize();z++){
-	for(int y=0;y<data.ysize();y++){
-	  for(int x=0;x<data.xsize();x++){
-	    v=0;v2=0;
-	    for(int i=1;i<=list.Ncols();i++){	  
-	      v += data[list(1,i)](x,y,z);
-	      v2 += data[list(1,i)](x,y,z)*data[list(1,i)](x,y,z);
-	      if(!meancalc && !varcalc)
-		subdata[i-1](x,y,z)= data[list(1,i)](x,y,z);
-	    }
-	    if(meancalc)
-	      subdata(x,y,z,0)= v/float(list.Ncols());
-	    if(varcalc)
-	      subdata(x,y,z,0)= v2/float(list.Ncols()) -v*v/float(list.Ncols()*list.Ncols());
-	  }
-	}
-      }
-      
+      float nImages(list.Ncols()); //float to avoid risk of integer math
+      for(int z=0;z<data.zsize();z++)
+	      for(int y=0;y<data.ysize();y++)
+	        for(int x=0;x<data.xsize();x++) {
+	          double v(0),v2(0);
+	          for(int i=1;i<=list.Ncols();i++) {
+              int index( list(1,i) );
+              double value( data(x,y,z,index) );
+              if( varcalc || meancalc ) {
+	              v += value;
+	              v2 += value * value;
+              } else
+		            subdata(x,y,z,i-1) = value;
+	          }
+	          if(meancalc)
+	            subdata(x,y,z,0)= v/nImages;
+	          if(varcalc)
+	            subdata(x,y,z,0)=(nImages/(nImages-1)) * (v2/nImages - v*v/(nImages*nImages));
+	        }
       save_volume4D(subdata,ofilename.value());
       return 0;
-
     }
-  }catch(X_OptionError& e) {
+  } catch(X_OptionError& e) {
     options.usage();
     cerr << endl << e.what() << endl;
     exit(EXIT_FAILURE);
-  }catch(std::exception &e) {
+  }
+  catch(std::exception &e) {
     cerr << e.what() << endl;
-  } 
+  }
   // either mean or variance but not both
   if(calc_mean.value() && calc_var.value()){
     cerr<<" Cannot output mean *and* variance. Please choose one or the other but not both"<<endl;
     exit(1);
   }
-  
-  
+
+
 }
-- 
GitLab