Initial revision

slicetimer.cc

+/*  slicetimer.cc - FMRIB's Slice Timing Utility
+    
+    Peter Bannister and Stephen Smith, FMRIB Image Analysis Group
+    
+    Copyright (C) 2001-2002 University of Oxford  */
+
+/*  CCOPYRIGHT */
+
+
+#include <math.h>
+#include <iostream.h>
+#include <iomanip>
+#include <fstream.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string>
+#include <strstream>
+#include <unistd.h>
+#include <vector.h>
+
+#include "optimise.h"
+#include "newmatap.h"
+#include "newmatio.h"
+#include "newimageall.h"
+#include "options.h"
+#include "kernel.h"
+
+using namespace MISCMATHS;
+using namespace NEWMAT;
+using namespace NEWIMAGE;
+using namespace Utilities;
+
+string title="slicetimer (Version 1.8)\nFMRIB's Interpolation for Slice Timing\nCopyright(c) 2001-2002, University of Oxford";
+string examples="slicetimer -i <timeseries> [-o <corrected_timeseries>] [options]\n";
+
+Option<bool> help(string("-h,--help"), false,
+		  string("display this message"),
+		  false, no_argument);
+Option<bool> verbose(string("-v,--verbose"), false, 
+		     string("switch on diagnostic messages"), 
+		     false, no_argument);
+Option<bool> odd(string("--odd"), false, 
+		     string("use interleaved acquisition"), 
+		     false, no_argument);
+Option<bool> down(string("--down"), false, 
+		     string("reverse slice indexing"), 
+		     false, no_argument);
+Option<string> inputname(string("-i,--in"), string(""),
+			 string("filename of input timeseries"),
+			 true, requires_argument);
+Option<string> outputname(string("-o,--out"), string(""),
+			  string("filename of output timeseries"),
+			  false, requires_argument);
+Option<string> tcustom(string("--tcustom"), string(""),
+			  string("filename of single-column custom interleave timing file"),
+			  false, requires_argument);
+Option<string> ocustom(string("--ocustom"), string(""),
+			  string("filename of single-column custom interleave order file"),
+			  false, requires_argument);
+Option<int> direction(string("-d,--direction"), 3,
+		      string("direction of slice acquisition (x=1,y=2,z=3) - default is z"),
+		      false, requires_argument);
+Option<float> repeat(string("-r,--repeat"), 3.0,
+		      string("Specify TR of data - default is 3s"),
+		      false, requires_argument);
+
+int do_slice_correction()
+{
+  volume4D<float> timeseries;
+  volumeinfo in_info;
+  ColumnVector voxeltimeseries, interpseries;
+  Matrix timings;
+  
+  int no_volumes, no_slices;
+  float repeat_time, offset, slice_spacing, recenter;
+
+  if (inputname.set()) {
+    if (verbose.value()) { cout << "Reading input volume" << endl; }
+    read_volume4D(timeseries,inputname.value(),in_info);
+    if (!outputname.set())
+      outputname.set_value(inputname.value() + "_st");
+  } else if (outputname.set()) {
+    cerr << "Must specify an input volume (-i or --in) to generate corrected data." 
+	 << endl;
+    return -1;
+  }
+
+  no_slices = timeseries.zsize();
+  no_volumes = timeseries.tsize();
+  repeat_time = timeseries.tdim();
+  if (repeat_time ==0){
+    cerr << "Zero TR in file header - fixing ... " ;
+    repeat_time = repeat.value();
+  }
+  cerr << "TR = " << repeat_time << endl;
+  slice_spacing = repeat_time / no_slices;
+
+  if (direction.value() == 1) timeseries. swapdimensions(3,2,1); // Flip z and x
+  if (direction.value() == 2) timeseries. swapdimensions(1,3,2); // Flip z and y
+
+  if (tcustom.set()) {
+    timings.ReSize(1, timeseries.zsize());
+    timings = read_ascii_matrix(tcustom.value(), timeseries.zsize(), 1);
+  }  else if (ocustom.set()) {
+    timings.ReSize(1, timeseries.zsize());
+    timings = read_ascii_matrix(ocustom.value(), timeseries.zsize(), 1);
+  }
+
+  offset = 0.0;
+
+  ColumnVector userkernel = sinckernel1D("hanning", 7, 1201);
+  
+  //  if (down.value())
+  //  recenter = (((float) no_slices)/2 + 0.5)/ no_slices;
+  //else
+  recenter = (((float) no_slices)/2 - 0.5)/ no_slices;
+  
+  for (int slice=1; slice<=no_slices; slice++) {
+    if (tcustom.set()) {
+      offset = -timings(slice, 1);
+    } else if (ocustom.set()) { 
+      int local_count=1;
+      while (local_count <= no_slices) {
+	if (timings(local_count, 1) == slice) {
+	  offset = -(local_count -1)* (slice_spacing / repeat_time);
+	  local_count = no_slices + 1;
+	} else 
+	  local_count++;
+      }
+    } else if (odd.value()) { // acquisition order: 1,3,5, ..., 2,4,6 ...
+      if ((slice % 2) == 0) // even
+	offset = -( ceil((float)no_slices / 2) + ((slice -1)/ 2)) * (slice_spacing / repeat_time);
+      else
+	offset = -((slice -1) / 2) * (slice_spacing / repeat_time);
+    } else if (down.value()) {
+      offset = -(no_slices - slice)* (slice_spacing / repeat_time);
+    } else { 
+      offset = -(slice -1)* (slice_spacing / repeat_time);
+    }
+    
+    for (int x_pos = 0; x_pos < timeseries. xsize(); x_pos++)
+      for (int y_pos = 0; y_pos < timeseries. ysize(); y_pos++){
+	voxeltimeseries = timeseries.voxelts(x_pos,y_pos,slice-1);
+	interpseries = voxeltimeseries;
+	for (int time_step=1; time_step <= no_volumes; time_step++)
+	  interpseries(time_step) = interpolate_1d(voxeltimeseries, time_step + offset + recenter);
+	timeseries.setvoxelts(interpseries,x_pos,y_pos,slice-1);
+      }
+    
+    if (verbose.value())
+      cerr << "Slice / Offset (s) " << slice << " / " << offset *repeat_time * -1 << endl;
+  }
+  
+  if (direction.value() == 1) timeseries. swapdimensions(3,2,1); // reverse Flip z and x
+  if (direction.value() == 2) timeseries. swapdimensions(1,3,2); // reverse Flip z and y
+
+  save_volume4D(timeseries, outputname.value());
+  return 0;
+}
+
+int main (int argc,char** argv)
+{
+  Tracer tr("main");
+  
+  OptionParser options(title, examples);
+
+  try {
+    options.add(inputname);
+    options.add(outputname);
+    options.add(help);
+    options.add(verbose);
+    options.add(down);
+    options.add(repeat);
+    options.add(direction);
+    options.add(odd);
+    options.add(tcustom);
+    options.add(ocustom);
+
+    options.parse_command_line(argc, argv);
+
+    if ( (help.value()) || (!options.check_compulsory_arguments(true)) )
+      {
+	options.usage();
+	exit(EXIT_FAILURE);
+      }
+        
+    if ( inputname.unset()) 
+      {
+	options.usage();
+	cerr << endl 
+	     << "--in or -i MUST be used." 
+	     << endl;
+	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;
+  } 
+
+  int retval = do_slice_correction();
+  
+  if (retval!=0) {
+    cerr << endl << endl << "Error detected: try -h for help" << endl;
+  }
+
+  return retval;
+}
+