//  
//  splinterpolator.h
//
// Jesper Andersson, FMRIB Image Analysis Group
//
// Copyright (C) 2008 University of Oxford 
//
//     CCOPYRIGHT
//

#ifndef splinterpolator_h
#define splinterpolator_h

#include <vector>
#include <string>
#include <cmath>
#include "newmat.h"
#include "miscmaths/miscmaths.h"

namespace SPLINTERPOLATOR {

enum ExtrapolationType {Zeros, Constant, Mirror, Periodic};

class SplinterpolatorException: public std::exception
{
private:
  std::string m_msg;
public:
  SplinterpolatorException(const std::string& msg) throw(): m_msg(msg) {}

  virtual const char *what() const throw() {
    return string("FslSplines::" + m_msg).c_str();
  }

  ~SplinterpolatorException() throw() {}
};

template<class T>
class Splinterpolator
{
public:
  // Constructors
  Splinterpolator() : _valid(false) {}
  Splinterpolator(const T *data, const std::vector<unsigned int>& dim, unsigned int order=3, std::vector<ExtrapolationType> et(3,Zeros))
  {
    common_construction(data,dim,order,et);
  }
  Splinterpolator(const T *data, const std::vector<unsigned int>& dim, unsigned int order, ExtrapolationType et)
  {
    std::vector<ExtrapolationType>   ett(dim.size(),et);
    common_construction(data,dim,order,ett);
  }
  // Destructor
  ~Splinterpolator() { if(_valid) delete [] _coef; }

  //
  // Here we declare nested helper-class SplineColumn
  //
  class SplineColumn
  {
  public:
    // Constructor
    SplineColumn(unsigned int sz, unsigned int step) : _sz(sz), _step(step) { _col = new double[_sz]; }
    // Destructor
    ~SplineColumn() { delete [] _col; }
    // Extract a column from a volume
    Get(const T  *dp)
    {
      for (unsigned int i=0; i<_sz; i++, dp+=_step) _col[i] = static_cast<double>(*dp);
    }
    // Insert column into volume
    Set(T *dp)
    {
      T   test = 1.5;
      if (test == 1) {  // If T is not float or double
        for (unsigned int i=0; i<_sz; i++, dp+=_step) *dp = static_cast<T>(_col[i] + 0.5);   // Round to nearest integer
      }
      else {
        for (unsigned int i=0; i<_sz; i++, dp+=_step) *dp = static_cast<T>(_col[i]);
      }
    }
    // Deconvolve column
    Deconv(unsigned int order, ExtrapolationType et) {};

  private:
    unsigned int  _sz;
    unsigned int  _step;
    double        *_col;

    SplineColumn(const SplineColumn& sc);              // Don't allow copy-construction
    SplineColumn& operator=(const SplineColumn& sc);   // Dont allow assignment
  };

  void Set(const T *data, const std::vector<unsigned int>& dim, unsigned int order=3, std::vector<ExtrapolationType> et(3,Zeros))
  {
    if (_valid) delete [] _coef;
    common_construction(data,dim,order,et);
  }
  
private:
  bool                            _valid;         // Decides if neccessary information has been set or not
  T                               *_coef;         // Volume of spline coefficients
  unsigned int                    _order;         // Order of splines
  std::vector<unsigned int>       _dim;           // Dimensions of data
  std::vector<ExtrapolationType>  _et;            // How to do extrapolation
  //
  // Private helper-functions
  //
  void common_construction(const T *data, const std::vector<unsigned int>& dim, unsigned int order, const std::vector<ExtrapolationType>& et);
  void calc_coef(const T *data);
  void deconv_along(unsigned int dim);
  //
  // Disallowed member functions
  //
  Splinterpolator(const Splinterpolator& s);             // Don't allow copy-construction
  Splinterpolator& operator=(const Splinterpolator& s);  // Don't allow assignment
};

template<class T>
void Splinterpolator<T>::common_construction(const T *data, const std::vector<unsigned int>& dim, unsigned int order, const std::vector<ExtrapolationType>& et)
{
  if (!dim.size()) throw SplinterpolatorException("Splinterpolator::common_construction: data has zeros dimensions");
  if (!dim.size() > 5) throw SplinterpolatorException("Splinterpolator::common_construction: data cannot have more than 5 dimensions");
  if (dim.size() != et.size()) throw SplinterpolatorException("Splinterpolator::common_construction: dim and et must have the same size");
  for (unsigned int i=0; i<dim.size(); i++) if (!dim[i]) throw SplinterpolatorException("Splinterpolator::common_construction: data cannot have zeros size in any direction");
  if (order > 7) throw SplinterpolatorException("Splinterpolator::common_construction: spline order must be lesst than 7");
  if (!data) throw SplinterpolatorException("Splinterpolator::common_construction: zero data pointer");
  
  for (unsigned int ts=1, i=0; i<dim.size(); i++) ts *= dim[i];
  _coef = new T[ts];
  _order = order;
  _dim.resize(5); = dim;
  _et = et;

  calc_coef(data);

  _valid = true;
}

template<class T>
void Splinterpolator<T>::calc_coef(const T *data)
{
  // Put the original data into _coef
  //
  for (unsigned int ts=1, i=0; i<dim.size(); i++) ts *= dim[i];
  memcpy(_coef,data,ts*sizeof(T));

  if (_order < 2) return;  // If nearest neighbour or linear, that's all we need 

  // Loop over all non-singleton dimensions and deconvolve along them
  //
  std::vector<unsigned int>  tdim(dim.size()-1,0); 
  for (unsigned int cdir=0; cdir<dim.size(); cdir++) {
    if (dim[cdir] > 1) deconv_along(cdir)
  }

  return;
}

template<class T>
void Splinterpolator<T>::deconv_along(unsigned int dim);
{
  T  *dp = _coef;

  // Set up to reflect "missing" dimension
  //
  std::vector<unsigned int>   rdim(4,1);     // Sizes along remaining dimensions
  std::vector<unsigned int>   rstep(4,1);    // Step-sizes (in "volume") of remaining dimensions
  unsigned int                mdim = 1;      // Size along "missing" dimension
  unsigned int                mstep = 1;     // Step-size along "missing" dimension
  for (unsigned int i=0, j=0, s=1; i<5; i++) {
    if (i == dim) {  // If it is our "missing" dimension
      mdim = _dim(i);
      mstep = ss;
    }
    else {
      rdim[j] = _dim[i];
      rstep[j++] = ss;
    }
    ss *= _dim[i];
  }

  SplineColumn<T>  col(mdim,mstep);           // Column helps us do the job

  for (unsigned int l=0; l<rdim[3]; l++, dp+=rstep[3]) {
    for (unsigned int k=0; k<rdim[2]; k++, dp+=rstep[2]) {
      for (unsigned int j=0; j<rdim[1]; j++, dp+=rstep[1]) {
	for (unsigned int i=0; i<rdim[0]; i++, dp+=rstep[0]) {
	  col.Get(dp);     // Extract a column from the volume
          col.Deconv(_order,_et[dim]);    // Deconvolve it
          col.Set(dp);     // Put back the deconvolved column
	}
      }
    }
  }
  return;
}

} // End namespace SPLINTERPOLATOR

#endif // End #ifndef splinterpolator.h