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;
}