Added executables for Geoff

Makefile

@@ -2,21 +2,34 @@ include $(FSLCONFDIR)/default.mk
 
 PROJNAME = siena
 
-USRINCFLAGS = -I${INC_ZLIB}
-USRLDFLAGS = -L${LIB_ZLIB}
+USRINCFLAGS = -I${INC_NEWMAT} -I${INC_ZLIB}
+USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_ZLIB}
 
 LIBS = -lss_32R -lfslio -lniftiio -lznz -lm -lz
+LIBCC = -lnewimage -lmiscmaths -lfslio -lniftiio -lznz -lnewmat -lutils -lm -lz
 
-XFILES = siena_diff
+SN_OBJS = surface_norm.o
+CS_OBJS = cylsamp.o
+
+XFILES = siena_diff surface_norm cylsamp
 
 SCRIPTS = siena siena_flirt siena_cal sienax siena_flow2tal
 
 OPTFLAGS=
 
-all:	siena_diff docscripts
+MJOPTFLAGS=-O3
+
+all:	${XFILES} docscripts
 
 siena_diff: siena_diff.c
 	$(CC) $(CFLAGS) -DFDT="float" -o siena_diff siena_diff.c $(LDFLAGS) $(LIBS)
 
+surface_norm: ${SN_OBJS}
+	${CXX} ${CXXFLAGS} ${MJOPTFLAGS} ${LDFLAGS} -o $@ ${SN_OBJS} ${LIBCC}
+
+cylsamp: ${CS_OBJS}
+	${CXX} ${CXXFLAGS} ${MJOPTFLAGS} ${LDFLAGS} -o $@ ${CS_OBJS} ${LIBCC}
+
 docscripts:
 	./makedoc
+

cylsamp.cc

+/*  cylsamp.cc
+
+    Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 2004 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+// Performs cylindrical sampling over a surface.  That is, averaging
+//  all values within a cylindrical region centred at each surface
+//  point on a mask.
+
+
+#define _GNU_SOURCE 1
+#define POSIX_SOURCE 1
+
+#include "newimage/newimageall.h"
+#include "miscmaths/miscmaths.h"
+#include "utils/options.h"
+
+using namespace MISCMATHS;
+using namespace NEWIMAGE;
+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="cylsamp (Version 1.0)\nCopyright(c) 2004, University of Oxford (Mark Jenkinson)";
+string examples="cylsamp [options] -i <input mask image> -s <smoothing in mm> -o <output surface normal image>";
+
+// 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<float>  radius(string("-r"),5.0,
+		      string("radius of cylinder in mm (default is 5.0)"),
+		      false, requires_argument);
+Option<float>  height(string("-h"),10.0,
+		      string("height of cylinder in mm (default is 10.0)"),
+		      false, requires_argument);
+Option<string> maskname(string("-m"), string(""),
+		      string("input mask image filename"),
+		      true, requires_argument);
+Option<string> normname(string("-n"), string(""),
+		      string("input surface normal image filename"),
+		      true, requires_argument);
+Option<string> edgemaskname(string("-e"), string(""),
+		      string("input edge mask image filename"),
+		      true, requires_argument);
+Option<string> flowname(string("-f"), string(""),
+		      string("input flow image filename"),
+		      true, requires_argument);
+Option<string> outname(string("-o"), string(""),
+		       string("output surface normal filename"),
+		       true, requires_argument);
+int nonoptarg;
+
+////////////////////////////////////////////////////////////////////////////
+
+// Local functions
+
+int do_work(int argc, char* argv[]) 
+{
+  volume<float> vflow, vmask, vedgemask, vsamp;
+  volume4D<float> snorm;
+  read_volume(vflow,flowname.value());
+  read_volume(vmask,maskname.value());
+  read_volume(vedgemask,edgemaskname.value());
+  read_volume4D(snorm,normname.value());
+
+  // force input to be a binary mask
+  if (verbose.value()) print_info(vmask,"vmask");
+  vsamp = vmask;
+
+  float r=radius.value();
+  float r2 = r*r;
+  float h=height.value();
+  int len=(int) (sqrt(r*r + h*h)) + 1;
+
+  bool atedge;
+  if (verbose.value()) { cerr << "Performing Cylindrical Sampling" << endl; }
+  for (int z=vmask.minz(); z<=vmask.maxz(); z++) {
+    for (int y=vmask.miny(); y<=vmask.maxy(); y++) {
+      for (int x=vmask.minx(); x<=vmask.maxx(); x++) {
+	atedge = false;
+	// check to see if it is an edge point
+	if ( (vmask(x,y,z)>0.5) ) {
+	  if (vmask(x,y,z-1)<0.5) atedge=true;
+	  else { 
+	    if (vmask(x,y-1,z)<0.5) atedge=true;
+	    else { 
+	      if (vmask(x-1,y,z)<0.5) atedge=true;
+	      else {
+		if (vmask(x+1,y,z)<0.5) atedge=true;
+		else {
+		  if (vmask(x,y+1,z)<0.5) atedge=true;
+		  else {
+		    if (vmask(x,y,z+1)<0.5) atedge=true;
+		  }
+		}
+	      }
+	    }
+	  }
+	}
+	if (atedge) {
+	  // OK, now do the cylindrical sampling
+	  float tot=0.0;
+	  int num=0;
+	  for (int z1=Max(z-len,vmask.minz()); z1<=Min(z+len,vmask.maxz()); z1++) {
+	    for (int y1=Max(y-len,vmask.miny()); y1<=Min(y+len,vmask.maxy()); y1++) {
+	      for (int x1=Max(x-len,vmask.minx()); x1<=Min(x+len,vmask.maxx()); x1++) {
+		if (vedgemask(x1,y1,z1)>0.5) {
+		  float yy = snorm(x1,y1,z1,0) * (x1 - x) 
+		    + snorm(x1,y1,z1,1) * (y1 - y)  
+		    + snorm(x1,y1,z1,2) * (z1 - z);
+		  if (fabs(yy)<=h) {
+		    float xx2 = ( Sqr(x1-x) + Sqr(y1-y) + Sqr(z1-z) ) - Sqr(yy);
+		    if (xx2 <= r2) {
+		      // inside cylinder
+		      num++;
+		      tot += vflow(x1,y1,z1);
+		    }
+		  }
+		}
+	      }
+	    }
+	  }
+	  vsamp(x,y,z) = tot / Max(num,1);  // prevent divide by zero
+	} else {
+	  vsamp(x,y,z) = 0.0;
+	}
+      }
+    }
+    if (verbose.value()) { cerr << "."; }
+  }
+  if (verbose.value()) { cerr << endl; }
+
+  // save the results
+  save_volume(vsamp,outname.value());
+
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+int main(int argc,char *argv[])
+{
+
+  Tracer tr("main");
+  OptionParser options(title, examples);
+
+  try {
+    // must include all wanted options here (the order determines how
+    //  the help message is printed)
+    options.add(flowname);
+    options.add(edgemaskname);
+    options.add(maskname);
+    options.add(normname);
+    options.add(outname);
+    options.add(radius);
+    options.add(height);
+    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);
+      }
+    
+  }  catch(X_OptionError& e) {
+    options.usage();
+    cerr << endl << e.what() << endl;
+    exit(EXIT_FAILURE);
+  } catch(std::exception &e) {
+    cerr << e.what() << endl;
+  } 
+
+  // Call the local functions
+
+  return do_work(argc,argv);
+}
+

surface_norm.cc

+/*  surface_norm.cc
+
+    Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 2004 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+// Calculates the surface normals for a mask, using a smoothed
+//  gradient calculation (all non-surface points get zero ouput)
+
+#define _GNU_SOURCE 1
+#define POSIX_SOURCE 1
+
+#include "newimage/newimageall.h"
+#include "miscmaths/miscmaths.h"
+#include "utils/options.h"
+
+using namespace MISCMATHS;
+using namespace NEWIMAGE;
+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="surface_norm (Version 1.0)\nCopyright(c) 2004, University of Oxford (Mark Jenkinson)";
+string examples="surface_norm [options] -i <input mask image> -s <smoothing in mm> -o <output surface normal image>";
+
+// 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<float>  smoothing(string("-s"),3.0,
+			 string("smoothing in mm (default is 3.0)"),
+			 false, requires_argument);
+Option<string> inname(string("-i"), string(""),
+		      string("input mask filename"),
+		      true, requires_argument);
+Option<string> outname(string("-o"), string(""),
+		       string("output surface normal filename"),
+		       true, requires_argument);
+int nonoptarg;
+
+////////////////////////////////////////////////////////////////////////////
+
+// Local functions
+
+int do_work(int argc, char* argv[]) 
+{
+  volume<float> vin;
+  volume4D<float> snorm;
+  read_volume(vin,inname.value());
+
+  // force input to be a binary mask
+  vin.binarise(0.5);
+  if (verbose.value()) print_info(vin,"vin");
+
+  // smooth the mask volume and calculate the gradient
+  if (verbose.value()) { cerr << "Performing smoothing" << endl; }
+  { 
+    volume<float> smoothv;
+    smoothv = smooth(vin,smoothing.value());
+    if (verbose.value()) print_info(smoothv,"smoothv");
+    if (verbose.value()) { cerr << "Calculating the gradient" << endl; }
+    gradient(smoothv,snorm);
+  }
+
+  if (verbose.value()) print_info(snorm,"snorm");
+
+  bool atedge;
+  if (verbose.value()) { cerr << "Normalising gradients" << endl; }
+  // zero all non-surface values and normalise vectors at the surface
+  for (int z=vin.minz(); z<=vin.maxz(); z++) {
+    for (int y=vin.miny(); y<=vin.maxy(); y++) {
+      for (int x=vin.minx(); x<=vin.maxx(); x++) {
+	atedge = false;
+	if ( (vin(x,y,z)>0.5) ) {
+	  if (vin(x,y,z-1)<0.5) atedge=true;
+	  else { 
+	    if (vin(x,y-1,z)<0.5) atedge=true;
+	    else { 
+	      if (vin(x-1,y,z)<0.5) atedge=true;
+	      else {
+		if (vin(x+1,y,z)<0.5) atedge=true;
+		else {
+		  if (vin(x,y+1,z)<0.5) atedge=true;
+		  else {
+		    if (vin(x,y,z+1)<0.5) atedge=true;
+		  }
+		}
+	      }
+	    }
+	  }
+	}
+	if (atedge) {
+	  float norm;
+	  norm = 1.0 / sqrt(Sqr(snorm(x,y,z,0)) + Sqr(snorm(x,y,z,1)) 
+			    + Sqr(snorm(x,y,z,2)) );
+	  snorm(x,y,z,0) *= norm;
+	  snorm(x,y,z,1) *= norm;
+	  snorm(x,y,z,2) *= norm;
+	} else {
+	  snorm(x,y,z,0) = 0.0;
+	  snorm(x,y,z,1) = 0.0;
+	  snorm(x,y,z,2) = 0.0;
+	}
+      }
+    }
+    if (verbose.value()) { cerr << "."; }
+  }
+  if (verbose.value()) { cerr << endl; }
+
+  // save the results
+  save_volume4D(snorm,outname.value());
+
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+int main(int argc,char *argv[])
+{
+
+  Tracer tr("main");
+  OptionParser options(title, examples);
+
+  try {
+    // must include all wanted options here (the order determines how
+    //  the help message is printed)
+    options.add(inname);
+    options.add(outname);
+    options.add(smoothing);
+    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);
+      }
+    
+  }  catch(X_OptionError& e) {
+    options.usage();
+    cerr << endl << e.what() << endl;
+    exit(EXIT_FAILURE);
+  } catch(std::exception &e) {
+    cerr << e.what() << endl;
+  } 
+
+  // Call the local functions
+
+  return do_work(argc,argv);
+}
+