/*  slicetimer.cc - FMRIB's Slice Timing Utility
    
    Peter Bannister and Stephen Smith, FMRIB Image Analysis Group
    
    Copyright (C) 2001-2003 University of Oxford  */

/*  CCOPYRIGHT */


#include <cmath>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstdlib>
#include <string>
#include <sstream>
#include <vector>

#include "miscmaths/optimise.h"
#include "newmatap.h"
#include "newmatio.h"
#include "newimage/newimageall.h"
#include "utils/options.h"
#include "miscmaths/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 (default is: slices were acquired bottom-up)"), 
		     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<float>  tglobal(string("--tglobal"), 0.5,
			  string("global shift (default is 0.5 = no shift)"),
			  false, requires_argument);
Option<string> ocustom(string("--ocustom"), string(""),
			  string("filename of single-column custom interleave order file (first slice is referred to as 1 not 0)"),
			  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;
  Matrix timings;
  
  int no_volumes, no_slices;
  float repeat_time, offset=0, slice_spacing;

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

  ColumnVector userkernel = sinckernel1D("hanning", 7, 1201);
  // for(int i=1; i<=1201; i++) cout << i << " " << userkernel(i) << endl;

  float recenter = (((float) no_slices)/2 - 0.5)/ no_slices; // only valid for slice-count-based corrections
  
  for (int slice=1; slice<=no_slices; slice++) {

    if (tglobal.set()) {
      offset = 0.5 - tglobal.value();
    } else if (tcustom.set()) {
      offset = 0.5 - timings(slice, 1);
    } else if (ocustom.set()) { 
      int local_count=1;
      while (local_count <= no_slices) {
	if (timings(local_count, 1) == slice) {
	  offset = recenter -(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 = recenter - ( ceil((float)no_slices / 2) + ((slice -1)/ 2)) * (slice_spacing / repeat_time);
      else
	offset = recenter - ((slice -1) / 2) * (slice_spacing / repeat_time);
    } else if (down.value()) {
      offset = recenter - (no_slices - slice) * (slice_spacing / repeat_time);
    } else { 
      offset = recenter - (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++){
	ColumnVector voxeltimeseries = timeseries.voxelts(x_pos,y_pos,slice-1);
	ColumnVector interpseries = voxeltimeseries;
	for (int time_step=1; time_step <= no_volumes; time_step++){
	  // interpseries(time_step) = interpolate_1d(voxeltimeseries, time_step - offset);
	  interpseries(time_step) = kernelinterpolation_1d(voxeltimeseries, time_step - offset, userkernel, 7);
	}
	timeseries.setvoxelts(interpseries,x_pos,y_pos,slice-1);
      }
    
    if (verbose.value())
      cerr << "Slice " << slice << " offset " << offset << 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(tglobal);
    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;
}