New 4D decorrelation utility

Makefile

@@ -134,3 +134,6 @@ fslcreatehd: fslcreatehd.o
 fslmaths: fslmaths.o
 	${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ fslmaths.o ${LIBS}
 
+fsldecorr4d: fsldecorr4d.o
+	${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ fsldecorr4d.o ${LIBS}
+

fsldecorr4d.cc

+/*  fsldecorr4d.cc
+
+    Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 2007 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+
+// Decorrelates a set of separable 4D components from
+//  a 4D dataset
+
+#include "newimage/newimageall.h"
+#include "miscmaths/miscmaths.h"
+#include "utils/options.h"
+
+using namespace NEWIMAGE;
+using namespace MISCMATHS;
+
+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="fsldecorr4d (Version 1.0)\nCopyright(c) 2004, University of Oxford (Mark Jenkinson)\nRemoves 4D components (by decorrelation) from a 4D dataset.\nEach component needs to be specified by a separate timecourse and spatial map.\nThe spatial maps are input as a single 4D file and the timecourses as a text matrix (with each column being a timecourse with the same ordering as the corresponding spatial maps.\n";
+string examples="fsldecorr4d -s <spatial maps> -t <timecourses> -i <input> -o <output> [-m mask]";
+
+// 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<string> smapname(string("-s"), string(""),
+		  string("~<filename>\tinput set of spatial maps (4D)"),
+		  true, requires_argument);
+Option<string> tcname(string("-t"), string(""),
+		  string("~<filename>\tinput set of timecourses (text matrix)"),
+		  true, requires_argument);
+Option<string> maskname(string("-m"), string(""),
+		  string("~<filename>\tinput 3D mask"),
+		  false, requires_argument);
+Option<string> outname(string("-o"), string(""),
+		  string("~<filename>\toutput 4D dataset"),
+		  true, requires_argument);
+Option<string> inname(string("-i"), string(""),
+		  string("~<filename>\tinput 4D dataset"),
+		  true, requires_argument);
+
+int nonoptarg;
+
+
+
+int do_work(int argc, char *argv[])
+{
+  volume4D<float> vin;
+  volumeinfo vinfo;
+  read_volume4D(vin,inname.value(),vinfo);
+  volume4D<float> smaps;
+  read_volume4D(smaps,smapname.value());
+
+  if (!samesize(vin[0],smaps[0])) {
+    cerr << "ERROR: Spatial maps and Input volumes have different (x,y,z) size." 
+	 << endl;
+    return 1;
+  }
+
+  volume<float> mask;
+  if (maskname.set()) {
+    read_volume(mask,maskname.value());
+  } else {
+    mask = vin[0];
+    mask = 1.0;
+  }
+
+  if (!samesize(vin[0],mask)) {
+    cerr << "ERROR: Mask and Input volumes have different (x,y,z) size." 
+	 << endl;
+    return 2;
+  }
+
+  mask.binarise(1e-8);  // arbitrary "0" threshold
+
+  Matrix tc;
+  tc = read_matrix(tcname.value());
+
+  if (tc.Nrows() != vin.tsize()) {
+    cerr << "ERROR: Different number of timepoints in timecourse file and input volume." << endl;
+    return 3;
+  }
+  if (tc.Ncols() != smaps.tsize()) {
+    cerr << "ERROR: Different number of components for timecourses and spatial maps." << endl;
+    return 4;
+  }
+
+  // set up required matrices
+  int nc = tc.Ncols();
+  int nt = tc.Nrows();
+  Matrix XX(nc,nc);
+  ColumnVector XY(nc);
+  XY=0.0;  XX=0.0;
+
+  // calculate matrix elements (4D correlations)
+  for (int n=1; n<=nc; n++) {
+    for (int t=0; t<nt; t++) {
+      for (int z=smaps[n-1].minz(); z<=smaps[n-1].maxz(); z++) {
+	for (int y=smaps[n-1].miny(); y<=smaps[n-1].maxy(); y++) {
+	  for (int x=smaps[n-1].minx(); x<=smaps[n-1].maxx(); x++) {
+	    XY(n) += smaps(x,y,z,n-1) * tc(t+1,n) * vin(x,y,z,t);
+	  }
+	}
+      }
+    }
+    for (int m=1; m<=n; m++) {
+      volume<float> tmp;
+      tmp = smaps[n-1] * smaps[m-1];
+      XX(m,n) = tmp.sum();
+      double tval=0.0;
+      for (int t=0; t<nt; t++) {
+	tval += tc(t+1,n) * tc(t+1,m);
+      }
+      XX(m,n) *= tval;
+      XX(n,m) = XX(m,n);
+    }
+  }
+
+  // find amplitudes for each component
+  ColumnVector Beta(nc);
+  Beta = pinv(XX)*XY;
+
+  // remove components from input data
+  volume4D<float> vout(vin);
+  for (int n=1; n<=nc; n++) {
+    for (int t=0; t<nt; t++) {
+      for (int z=smaps[n-1].minz(); z<=smaps[n-1].maxz(); z++) {
+	for (int y=smaps[n-1].miny(); y<=smaps[n-1].maxy(); y++) {
+	  for (int x=smaps[n-1].minx(); x<=smaps[n-1].maxx(); x++) {
+	    vout(x,y,z,t) -= Beta(n) * smaps(x,y,z,n-1) * tc(t+1,n);
+	  }
+	}
+      }
+    }
+  }
+
+  // save the result
+  save_volume4D(vout,outname.value(),vinfo);
+
+  return 0;
+}
+
+
+int main(int argc,char *argv[])
+{
+
+  Tracer tr("main");
+  OptionParser options(title, examples);
+
+  options.add(inname);
+  options.add(outname);
+  options.add(smapname);
+  options.add(tcname);
+  options.add(maskname);
+  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 ...
+  return do_work(argc,argv);
+}
+