New functionality in avwmeants

fslmeants.cc

@@ -12,19 +12,53 @@
 
 #include "newimage/newimageall.h"
 #include "miscmaths/miscmaths.h"
+#include "utils/options.h"
 
 using namespace NEWIMAGE;
 using namespace MISCMATHS;
 
 
-void print_usage(const string& progname) {
-  cerr << "Usage: " << progname << " <input 4D volume> <output text file> [<mask volume> [-v]] " 
-       << endl << endl;
-  cerr << "e.g.   " << progname << " filtered_func_data meants.txt my_mask" 
-       << endl;
-  cerr << "       " << progname << " filtered_func_data meants.txt" 
-       << endl;
-}
+using namespace Utilities;
+
+// The two strings below specify the title and example usage that is
+//  printed out as the help or usage message
+
+string title="avwmeants (Version 1.1)\nCopyright(c) 2004, University of Oxford (Mark Jenkinson)\nPrints average timeseries (intensities) to the screen (or saves to a file).\nThe average is taken over all voxels in the mask (or all voxels in the image if no mask is specified).\n";
+string examples="avwmeants -i filtered_func_data -o meants.txt -m my_mask\navwmeants -i filtered_func_data -m my_mask\navwmeants -i filtered_func_data -c 24 19 10";
+
+// Each (global) object below specificies as option and can be accessed
+//  anywhere in this file (since they are global).  The order of the
+//  arguments needed is: name(s) of option, default value, help message,
+//       whether it is compulsory, whether it requires arguments
+// Note that they must also be included in the main() function or they
+//  will not be active.
+
+Option<bool> verbose(string("-v,--verbose"), false, 
+		     string("switch on diagnostic messages"), 
+		     false, no_argument);
+Option<bool> help(string("-h,--help"), false,
+		  string("display this message"),
+		  false, no_argument);
+Option<bool> usemm(string("--usemm"), false,
+		  string("\tuse mm instead of voxel coordinates (for -c option)"),
+		  false, no_argument);
+Option<bool> showall(string("--showall"), false,
+		  string("show all voxel time series (within mask) instead of averaging"),
+		  false, no_argument);
+Option<string> inname(string("-i"), string(""),
+		  string("~<filename>\tinput 4D image"),
+		  true, requires_argument);
+Option<string> maskname(string("-m"), string(""),
+		  string("~<filename>\tinput 4D mask"),
+		  false, requires_argument);
+Option<string> outmat(string("-o"), string(""),
+		  string("~<filename>\toutput text matrix"),
+		  false, requires_argument);
+Option<float> coordval(string("-c"), 0.0,
+		       string("~<x y z>\trequested spatial coordinate (instead of mask)"), 
+		       false, requires_3_arguments);
+
+int nonoptarg;
 
 
 
@@ -32,66 +66,110 @@ int main(int argc,char *argv[])
 {
 
   Tracer tr("main");
-
-  string progname=argv[0];
-  if (argc<3) { 
-    print_usage(progname);
-    return -1; 
-  }
-
-  bool verbose = false;
-  if (argc>=5) {
-    string optarg=argv[4];
-    if (optarg=="-v") { 
-      verbose = true;
-    } else {
-      cerr << "Unrecognised option: " << optarg << endl;
+  OptionParser options(title, examples);
+
+  options.add(inname);
+  options.add(outmat);
+  options.add(maskname);
+  options.add(coordval);
+  options.add(usemm);
+  options.add(showall);
+  options.add(verbose);
+  options.add(help);
+  
+  nonoptarg = options.parse_command_line(argc, argv);
+  
+  // line below stops the program if the help was requested or 
+  //  a compulsory option was not set
+  if ( (help.value()) || (!options.check_compulsory_arguments(true)) )
+    {
+      options.usage();
+      exit(EXIT_FAILURE);
     }
-  }
+  
+
+  // OK, now do the job ...
 
   volume4D<float> vin;
-  read_volume4D(vin,argv[1]);
+  read_volume4D(vin,inname.value());
 
   volume<float> mask;
-  if (argc>=4) {
-    read_volume(mask,argv[3]);
+  if (maskname.set()) {
+    read_volume(mask,maskname.value());
   } else {
     mask = vin[0];
     mask = 1.0;
   }
 
+
+  // override the mask if a coordinate is specified
+  if (coordval.set()) {
+    mask = vin[0];
+    mask = 0.0;
+    float x, y, z;
+    x = coordval.value(0);
+    y = coordval.value(1);
+    z = coordval.value(2);
+    if (usemm.value()) {
+      // convert from mm to voxels
+      x /= vin[0].xdim();
+      y /= vin[0].ydim();
+      z /= vin[0].zdim();
+    }
+    mask(MISCMATHS::round(x),MISCMATHS::round(y),MISCMATHS::round(z)) = 1.0;
+  }
+
   if (!samesize(vin[0],mask)) {
     cerr << "ERROR: Mask and Input volumes have different (x,y,z) size." 
 	 << endl;
     return -1;
   }
 
+  mask.binarise(1e-8);  // arbitrary "0" threshold
+
   int nt = vin.tsize();
-  Matrix meants(nt,1);
+  int nvox = 1;
+  if (showall.value()) {
+    nvox = MISCMATHS::round(mask.sum());
+    nt += 3;
+  }
+  Matrix meants(nt,nvox);
   meants = 0;
   long int num=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++) {
-	if (fabs(mask(x,y,z))>1e-8) {  // arbitrary "0" threshold
-	  meants += vin.voxelts(x,y,z);
+	if (mask(x,y,z)>0.5) {
 	  num++;
-	  if (verbose) {
-	    cout << x << " " << y << " " << z << " " << vin.voxelts(x,y,z).t();
+	  if (showall.value()) {
+	    meants(1,num) = x;
+	    meants(2,num) = y;
+	    meants(3,num) = z;
+	    meants.SubMatrix(4,nt,num,num) = vin.voxelts(x,y,z).t();
+	  } else {
+	    meants += vin.voxelts(x,y,z);
 	  }
 	}
       }
     }
   }
 
-  cout << "Number of voxels used = " << num << endl;
+  if (verbose.value()) {
+    cout << "Number of voxels used = " << num << endl;
+  }
 
-  // normalise for number of valid entries
-  meants /= (float) num;
+  // normalise for number of valid entries if averaging
+  if (!showall.value()) {
+    if (num>0) meants /= (float) num;
+  }
 
   // save the result
-  write_ascii_matrix(meants,argv[2]);
+  if (outmat.set()) {
+    write_ascii_matrix(meants,outmat.value());
+  } else {
+    cout << meants << endl;
+  }
 
   return 0;
 }