From b8eec68ade0f1edb0d6ad54761f10010acf3a5bb Mon Sep 17 00:00:00 2001
From: Mark Jenkinson <mark@fmrib.ox.ac.uk>
Date: Fri, 5 Nov 2004 17:54:24 +0000
Subject: [PATCH] Added executables for Geoff
---
Makefile | 21 +++++-
cylsamp.cc | 196 ++++++++++++++++++++++++++++++++++++++++++++++++
surface_norm.cc | 167 +++++++++++++++++++++++++++++++++++++++++
3 files changed, 380 insertions(+), 4 deletions(-)
create mode 100644 cylsamp.cc
create mode 100644 surface_norm.cc
diff --git a/Makefile b/Makefile
index 922a516..3a63bd4 100644
--- a/Makefile
+++ b/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
+
diff --git a/cylsamp.cc b/cylsamp.cc
new file mode 100644
index 0000000..f4acc23
--- /dev/null
+++ b/cylsamp.cc
@@ -0,0 +1,196 @@
+/* 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);
+}
+
diff --git a/surface_norm.cc b/surface_norm.cc
new file mode 100644
index 0000000..bdb1aba
--- /dev/null
+++ b/surface_norm.cc
@@ -0,0 +1,167 @@
+/* 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);
+}
+
--
GitLab