splinterpolator.h

//  
//  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("Splinterpolator::" + m_msg).c_str();
  }

  ~SplinterpolatorException() throw() {}
};

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//
// Class Splinterpolator:
//
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

template<class T>
class Splinterpolator
{
public:
  // Constructors
  Splinterpolator() : _valid(false) {}
  Splinterpolator(const T *data, const std::vector<unsigned int>& dim, const std::vector<ExtrapolationType>& et, unsigned int order=3, double prec=1e-8)
  {
    common_construction(data,dim,order,prec,et);
  }
  Splinterpolator(const T *data, const std::vector<unsigned int>& dim, ExtrapolationType et=Zeros, unsigned int order=3, double prec=1e-8)
  {
    std::vector<ExtrapolationType>   ett(dim.size(),et);
    common_construction(data,dim,order,prec,ett);
  }

  // Destructor
  ~Splinterpolator() { if(_valid) delete [] _coef; }

  // Set new data in Splinterpolator.
  void Set(const T *data, const std::vector<unsigned int>& dim, const std::vector<ExtrapolationType>& et, unsigned int order=3)
  {
    if (_valid) delete [] _coef;
    common_construction(data,dim,order,et);
  }
  void Set(const T *data, const std::vector<unsigned int>& dim, ExtrapolationType et, unsigned int order=3)
  {
    std::vector<ExtrapolationType>  vet(dim.size(),Zeros);
    Set(data,dim,vet,order);
  }

  // Return interpolated value
  T operator()(const std::vector<float>&  coord) const;
  T operator()(double x, double y=0, double z=0, double t=0) const
  {
    double coord[5] = {x,y,z,t,0.0};
    return(value_at(coord));
  }

  //
  // The "useful" functionality pretty much ends here.
  // Remaining functions are mainly for debugging/diagnostics.
  //
  unsigned int NDim() const { return(_ndim); }
  const std::vector<unsigned int>& Size() const { return(_dim); }
  unsigned int Size(unsigned int dim) const { if (dim > 4) return(0); else return(_dim[dim]);}
  T Coef(unsigned int x, unsigned int y=0, unsigned int z=0) const
  {
    std::vector<unsigned int> indx(3,0);
    indx[0] = x; indx[1] = y; indx[2] = z;
    return(Coef(indx));
  }
  T Coef(std::vector<unsigned int> indx) const;
  NEWMAT::ReturnMatrix CoefAsNewmatMatrix() const;
  NEWMAT::ReturnMatrix KernelAsNewmatMatrix(double sp=0.1, unsigned int deriv=0) const;
  
  //
  // 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
    void Get(const T  *dp)
    {
      for (unsigned int i=0; i<_sz; i++, dp+=_step) _col[i] = static_cast<double>(*dp);
    }
    // Insert column into volume
    void Set(T *dp) const
    {
      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
    void Deconv(unsigned int order, ExtrapolationType et, double prec);

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

    unsigned int get_poles(unsigned int order, double *z, unsigned int *sf) const;
    double init_bwd_sweep(double z, double lv, ExtrapolationType et, double prec) const;
    double init_fwd_sweep(double z, ExtrapolationType et, double prec) const;

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

  
private:
  bool                            _valid;         // Decides if neccessary information has been set or not
  T                               *_coef;         // Volume of spline coefficients
  unsigned int                    _order;         // Order of splines
  unsigned int                    _ndim;          // # of non-singleton dimensions
  double                          _prec;          // Precision when dealing with edges
  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, double prec, const std::vector<ExtrapolationType>& et);
  void calc_coef(const T *data);
  void deconv_along(unsigned int dim);
  T coef(int *indx) const;
  double value_at(const double *coord) const;
  unsigned int get_start_indicies(const double *coord, int *sinds) const;
  unsigned int get_wgts(const double *coord, const int *sinds, double **wgts) const;
  double get_wgt(double x) const;
  double get_dwgt(double x) const;
  std::pair<double,double> range() const;
  bool odd(unsigned int i) const {return(static_cast<bool>(i%2));}
  bool even(unsigned int i) const {return(!odd(i));}

  //
  // Disallowed member functions
  //
  Splinterpolator(const Splinterpolator& s);             // Don't allow copy-construction
  Splinterpolator& operator=(const Splinterpolator& s);  // Don't allow assignment
};

/////////////////////////////////////////////////////////////////////
//
// Here starts public member functions for Splinterpolator
//
/////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////
//
// Returns the value of the coefficient given by indx (zero-offset)
//
/////////////////////////////////////////////////////////////////////

template<class T>
T Splinterpolator<T>::Coef(std::vector<unsigned int> indx) const
{
  if (!indx.size()) throw SplinterpolatorException("Coef: indx has zeros dimensions");
  if (indx.size() > 5) throw SplinterpolatorException("Coef: indx has more than 5 dimensions");
  for (unsigned int i=0; i<indx.size(); i++) if (indx[i] >= _dim[i]) throw SplinterpolatorException("Coef: indx out of range");

  unsigned int lindx=indx[indx.size()-1];
  for (int i=indx.size()-2; i>=0; i--) lindx = _dim[i]*lindx + indx[i];
  return(_coef[lindx]);
}

/////////////////////////////////////////////////////////////////////
//
// Returns the values of all coefficients as a Newmat matrix. If 
// _ndim==1 it will return a row-vector, if _ndim==2 it will return
// a matrix, if _ndim==3 it will return a tiled matrix where the n
// first rows (where n is the number of rows in one slice) pertain to 
// the first slice, the next n rows to the second slice, etc. And 
// correspondingly for 4- and 5-D.
//
/////////////////////////////////////////////////////////////////////

template<class T>
NEWMAT::ReturnMatrix Splinterpolator<T>::CoefAsNewmatMatrix() const
{
  NEWMAT::Matrix   mat(_dim[1]*_dim[2]*_dim[3]*_dim[4],_dim[0]);
  std::vector<unsigned int>  cindx(5,0);

  unsigned int r=0;
  for (cindx[4]=0; cindx[4]<_dim[4]; cindx[4]++) {
    for (cindx[3]=0; cindx[3]<_dim[3]; cindx[3]++) {
      for (cindx[2]=0; cindx[2]<_dim[2]; cindx[2]++) {
	for (cindx[1]=0; cindx[1]<_dim[1]; cindx[1]++, r++) {
	  for (cindx[0]=0; cindx[0]<_dim[0]; cindx[0]++) {
            mat.element(r,cindx[0]) = Coef(cindx);  
	  }
	}
      }
    }
  }
  mat.Release();
  return(mat);
}

/////////////////////////////////////////////////////////////////////
//
// Return the kernel matrix to verify its correctness.
//
/////////////////////////////////////////////////////////////////////

template<class T>
NEWMAT::ReturnMatrix Splinterpolator<T>::KernelAsNewmatMatrix(double sp,                // Distance (in ksp) between points 
                                                              unsigned int deriv) const // Derivative (only 0/1 implemented).
{
  if (deriv > 1) throw SplinterpolatorException("KernelAsNewmatMatrix: only 1st derivatives implemented");

  std::pair<double,double>           rng = range();

  unsigned int i=0;
  for (double x=rng.first; x<=rng.second; x+=sp, i++) ; // Intentional
  NEWMAT::Matrix  kernel(i,2);
  for (double x=rng.first, i=0; x<=rng.second; x+=sp, i++) {
    kernel.element(i,0) = x;
    kernel.element(i,1) = (deriv) ? get_dwgt(x) : get_wgt(x);
  }
  kernel.Release();
  return(kernel);  
}
/////////////////////////////////////////////////////////////////////
//
// Here starts public member functions for SplineColumn
//
/////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////
//
// This function implements the forward and backwards sweep
// as defined by equation 2.5 in Unsers paper:
//
// B-spline signal processing. II. Efficiency design and applications
//
/////////////////////////////////////////////////////////////////////

template<class T>
void Splinterpolator<T>::SplineColumn::Deconv(unsigned int order, ExtrapolationType et, double prec)
{
  double         z[3] = {0.0, 0.0, 0.0};     // Poles
  unsigned int   np = 0;                     // # of poles
  unsigned int   sf;                         // Scale-factor
  np = get_poles(order,z,&sf);

  for (unsigned int p=0; p<np; p++) {
    _col[0] = init_fwd_sweep(z[p],et,prec);
    double lv = _col[_sz-1];
    // Forward sweep
    double *ptr=&_col[1];
    for (unsigned int i=1; i<_sz; i++, ptr++) *ptr += z[p] * *(ptr-1);
    _col[_sz-1] = init_bwd_sweep(z[p],lv,et,prec);
    // Backward sweep
    ptr = &_col[_sz-2];
    for (int i=_sz-2; i>=0; i--, ptr--) *ptr = z[p]*(*(ptr+1) - *ptr);
  }
  double *ptr=_col;
  for (unsigned int i=0; i<_sz; i++, ptr++) *ptr *= sf; 
}

/////////////////////////////////////////////////////////////////////
//
// Here starts private member functions for Splinterpolator
//
/////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////
//
// Returns the interpolated value at location given by coord. 
// coord must be a pointer to an array of indicies with _ndim
// values.
//
/////////////////////////////////////////////////////////////////////

template<class T>
double Splinterpolator<T>::value_at(const double *coord) const
{
  double       iwgt[8], jwgt[8], kwgt[8], lwgt[8], mwgt[8];
  double       *wgts[] = {iwgt, jwgt, kwgt, lwgt, mwgt};
  int          inds[5];
  unsigned int ni = 0;

  ni = get_start_indicies(coord,inds);
  get_wgts(coord,inds,wgts);

  double val=0.0;
  for (int m=0, me=(_ndim>4)?ni:1; m<me; m++) { 
    for (int l=0, le=(_ndim>3)?ni:1; l<le; l++) { 
      for (int k=0, ke=(_ndim>2)?ni:1; k<ke; k++) { 
        double wgt1 = wgts[4][m]*wgts[3][l]*wgts[2][k];
        for (int j=0, je=(_ndim>1)?ni:1; j<je; j++) { 
          double wgt2 = wgt1*wgts[1][j];
          for (int i=0; i<static_cast<int>(ni); i++) {
            int cindx[] = {inds[0]+i,inds[1]+j,inds[2]+k,inds[3]+l,inds[4]+m};
            val += coef(cindx)*wgts[0][i]*wgt2;
	  }
	}
      }
    }
  }
  return(val);
}

/////////////////////////////////////////////////////////////////////
//
// Returns (in sinds) the indicies of the first coefficient in all
// _ndim directions with a non-zero weight for the location given
// by coord. The caller is responsible for coord and sinds being
// valid pointers to arrays of 5 values.
// The return-value gives the total # of non-zero weights.
//
/////////////////////////////////////////////////////////////////////

template<class T>
unsigned int Splinterpolator<T>::get_start_indicies(const double *coord, int *sinds) const
{
  unsigned int ni = _order+1;
  
  if (odd(ni)) {
    for (unsigned int i=0; i<_ndim; i++) {
      sinds[i] = static_cast<int>(coord[i]+0.5) - ni/2;
    }
  }
  else {
    for (unsigned int i=0; i<_ndim; i++) {
      int ix = static_cast<int>(coord[i]+0.5);
      if (ix < coord[i]) sinds[i] = ix - (ni-1)/2;
      else sinds[i] = ix -ni/2;
    }
  }
  for (unsigned int i=_ndim; i<5; i++) sinds[i] = 0;

  return(ni);
}

/////////////////////////////////////////////////////////////////////
//
// Returns (in wgts) the weights for the coefficients given by sinds 
// for the location given by coord. 
//
/////////////////////////////////////////////////////////////////////

template<class T>
unsigned int Splinterpolator<T>::get_wgts(const double *coord, const int *sinds, double **wgts) const
{
  unsigned int ni = _order+1;

  for (unsigned int dim=0; dim<_ndim; dim++) {
    for (unsigned int i=0; i<ni; i++) {
      wgts[dim][i] = get_wgt(coord[dim]-(sinds[dim]+i));
    }
  }
  for (unsigned int dim=_ndim; dim<5; dim++) wgts[dim][0] = 1.0;

  return(ni);
}

/////////////////////////////////////////////////////////////////////
//
// Returns the weight for a spline at coordinate x, where x is relative
// to the centre of the spline.
//
/////////////////////////////////////////////////////////////////////

template<class T>
double Splinterpolator<T>::get_wgt(double x) const
{
  double val = 0.0;
  double ax = abs(x);  // Kernels all symmetric

  switch (_order) {
  case 0:
    if (ax < 0.5) val = 1.0;
    break;
  case 1:
    if (ax < 1) val = 1-ax;;
    break;
  case 2:
    if (ax < 0.5) val = 0.75-ax*ax;
    else if (ax < 1.5) val = 0.5*(1.5-ax)*(1.5-ax);
    break;
  case 3:
    if (ax < 1) val = 2.0/3.0 + 0.5*ax*ax*(ax-2);
    else if (ax < 2) { ax = 2-ax; val = (1.0/6.0)*(ax*ax*ax); }
    break;
  case 4:
    if (ax < 0.5) { ax *= ax; val = (115.0/192.0) + ax*((2.0*ax-5.0)/8.0); }
    else if (ax < 1.5) val = (55.0/96.0) + ax*(ax*(ax*((5.0-ax)/6.0) - 1.25) + 5.0/24.0);
    else if (ax < 2.5) { ax -= 2.5; ax *= ax; val = (1.0/24.0)*ax*ax; }
    break;
  case 5:
    if (ax < 1) { double xx = ax*ax; val = 0.55 + xx*(xx*((3.0-ax)/12.0) - 0.5); }
    else if (ax < 2) val = 0.425 + ax*(ax*(ax*(ax*((ax-9.0)/24.0) + 1.25) - 1.75) + 0.625);
    else if (ax < 3) { ax = 3-ax; double xx = ax*ax; val = (1.0/120.0)*ax*xx*xx; }
    break;
  case 6:
    if (ax < 0.5) { ax *= ax; val = (5887.0/11520.0) + ax*(ax*((21.0-4.0*ax)/144.0) -77.0/192.0); }
    else if (ax < 1.5) val = 7861.0/15360.0 + ax*(ax*(ax*(ax*(ax*((ax - 7.0)/48.0) + 0.328125) - 35.0/288.0) - 91.0/256.0) -7.0/768.0);
    else if (ax < 2.5) val = 1379.0/7680.0 + ax*(ax*(ax*(ax*(ax*((14.0-ax)/120.0) - 0.65625) + 133.0/72.0) - 2.5703125) + 1267.0/960.0);
    else if (ax < 3.5) { ax -= 3.5; ax *= ax*ax; val = (1.0/720.0) * ax*ax; }
    break;
  case 7:
    if (ax < 1) { double xx = ax*ax; val = 151.0/315.0 + xx*(xx*(xx*((ax-4.0)/144.0) + 1.0/9.0) - 1.0/3.0); }
    else if (ax < 2) val = 103.0/210.0 + ax*(ax*(ax*(ax*(ax*(ax*((12.0-ax)/240.0) -7.0/30.0) + 0.5) - 7.0/18.0) - 0.1) -7.0/90.0);
    else if (ax < 3) val = ax*(ax*(ax*(ax*(ax*(ax*((ax-20.0)/720.0) + 7.0/30.0) - 19.0/18.0) + 49.0/18.0) - 23.0/6.0) + 217.0/90.0) - 139.0/630.0;
    else if (ax < 4) { ax = 4-ax; double xxx=ax*ax*ax; val = (1.0/5040.0)*ax*xxx*xxx; }
    break;
  default:
    throw SplinterpolatorException("get_wgt: invalid order spline");
    break;
  }

  return(val);
}

/////////////////////////////////////////////////////////////////////
//
// Returns the weight for the first derivative of a spline at 
// coordinate x, where x is relative to the centre of the spline.
//
/////////////////////////////////////////////////////////////////////

template<class T>
double Splinterpolator<T>::get_dwgt(double x) const
{
  double val = 0.0;
  double ax = abs(x);             // Kernels all anti-symmetric
  int    sign = (ax) ? x/ax : 1;  // Arbitrary choice for when x=0

  switch (_order) {
  case 0:
    throw SplinterpolatorException("get_dwgt: invalid order spline");
    break;
  case 1:
    if (ax < 1) val = sign * -1;
    break;
  case 2:
    if (ax < 0.5) val = sign * -2.0*ax;
    else if (ax < 1.5) val = sign * (-1.5 + ax);
    break;
  case 3:
    if (ax < 1) val = sign * (1.5*ax*ax - 2.0*ax);
    else if (ax < 2) { ax = 2-ax; val = sign * -0.5*ax*ax; }
    break;
  case 4:
    if (ax < 0.5) val = sign * (ax*ax*ax - 1.25*ax);
    else if (ax < 1.5) val = sign * (5.0/24.0 - ax*(2.5 - ax*(2.5 - (2.0/3.0)*ax)));
    else if (ax < 2.5) { ax -= 2.5; val = sign * (1.0/6.0)*ax*ax*ax; }
    break;
  case 5:
    if (ax < 1) val = sign * ax*(ax*(ax*(1-(5.0/12.0)*ax)) - 1);
    else if (ax < 2) val = sign * (0.625 - ax*(3.5 - ax*(3.75 - ax*(1.5 - (5.0/24.0)*ax))));
    else if (ax < 3) { ax -= 3; ax = ax*ax; val = sign * (-1.0/24.0)*ax*ax; }
    break;
  case 6:
    break;
  case 7:
    break;
  default:
    throw SplinterpolatorException("get_dwgt: invalid order spline");
    break;
  }

  return(val);
}

template<class T>
std::pair<double,double> Splinterpolator<T>::range() const
{
  std::pair<double,double>  rng(0.0,0.0);
  rng.second = static_cast<double>(_order+1.0)/2.0;
  rng.first = - rng.second;
  return(rng);
}
/////////////////////////////////////////////////////////////////////
//
// Returns the value of the coefficient indexed by indx. Unlike the
// public Coef() this routine allows indexing outside the valid
// volume, returning values that are dependent on the extrapolation
// model when these are encountered.
//
// N.B. May change value of input index N.B.
//
/////////////////////////////////////////////////////////////////////

template<class T>
T Splinterpolator<T>::coef(int *indx) const
{
  // First fix any outside-volume indicies
  for (unsigned int i=0; i<_ndim; i++) {
    if (indx[i] < 0) {
      switch (_et[i]) {
      case Zeros:
	return(static_cast<T>(0));
	break;
      case Constant:
	indx[i] = 0;
	break;
      case Mirror:
	indx[i] = 1-indx[i];
	break;
      case Periodic:
	indx[i] = _dim[i]+indx[i];
	break;
      default:
	break;
      }
    }
    else if (indx[i] >= static_cast<int>(_dim[i])) {
      switch (_et[i]) {
      case Zeros:
	return(static_cast<T>(0));
	break;
      case Constant:
	indx[i] = _dim[i]-1;
	break;
      case Mirror:
	indx[i] = 2*_dim[i]-indx[i]-1;
	break;
      case Periodic:
	indx[i] = indx[i]-_dim[i];
	break;
      default:
	break;
      }
    }
  }
  // Now make linear index
  unsigned int lindx=indx[_ndim-1];
  for (int i=_ndim-2; i>=0; i--) lindx = _dim[i]*lindx + indx[i];
  return(_coef[lindx]);
}

/////////////////////////////////////////////////////////////////////
//
// Takes care of the "common" tasks when constructing a 
// Splinterpolator object. Called by constructors and by .Set()
//
/////////////////////////////////////////////////////////////////////

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

  calc_coef(data);

  _valid = true;
}

/////////////////////////////////////////////////////////////////////
//
// Performs deconvolution, converting signal to spline coefficients.
//
/////////////////////////////////////////////////////////////////////

template<class T>
void Splinterpolator<T>::calc_coef(const T *data)
{
  // Put the original data into _coef
  //
  unsigned int ts=1;
  for (unsigned int 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;
}

/////////////////////////////////////////////////////////////////////
//
// Performs deconvolution along one of the dimensions, visiting
// all points along the other dimensions.
//
/////////////////////////////////////////////////////////////////////

template<class T>
void Splinterpolator<T>::deconv_along(unsigned int dim)
{
  // 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, ss=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  col(mdim,mstep);          // Column helps us do the job
  
  for (unsigned int l=0; l<rdim[3]; l++) {
    for (unsigned int k=0; k<rdim[2]; k++) {
      for (unsigned int j=0; j<rdim[1]; j++) {
        T *dp = _coef + l*rstep[3] + k*rstep[2] + j*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],_prec);    // Deconvolve it
          col.Set(dp);                          // Put back the deconvolved column
	}
      }
    }
  }
  return;
}

/////////////////////////////////////////////////////////////////////
//
// Here starts private member functions for SplineColumn
//
/////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////
//
// This function returns the poles and scale-factors for splines
// of order 2--7. The values correspond to those found in
// table 1 in Unsers 1993 paper: 
// B-spline signal processing. II. Efficiency design and applications
// 
// The actual values have been taken from
// http://bigwww.epfl.ch/thevenaz/interpolation/coeff.c
//
/////////////////////////////////////////////////////////////////////

template<class T>
unsigned int Splinterpolator<T>::SplineColumn::get_poles(unsigned int order, double *z, unsigned int *sf) const
{
  unsigned int   np = 0;                     // # of poles

  switch (order) {
  case 2:
    np = 1;
    z[0] = 2.0*sqrt(2.0) - 3.0;
    *sf = 8;
    break;
  case 3:
    np = 1;
    z[0] = sqrt(3.0) - 2.0;
    *sf = 6;
    break;
  case 4:
    np = 2;
    z[0] = sqrt(664.0 - sqrt(438976.0)) + sqrt(304.0) - 19.0;
    z[1] = sqrt(664.0 + sqrt(438976.0)) - sqrt(304.0) - 19.0;
    *sf = 384;
    break;
  case 5:
    np = 2;
    z[0] = sqrt(135.0 / 2.0 - sqrt(17745.0 / 4.0)) + sqrt(105.0 / 4.0) - 13.0 / 2.0;
    z[1] = sqrt(135.0 / 2.0 + sqrt(17745.0 / 4.0)) - sqrt(105.0 / 4.0) - 13.0 / 2.0;
    *sf = 120;
    break;
  case 6:
    np = 3;
    z[0] = -0.48829458930304475513011803888378906211227916123938;
    z[1] = -0.081679271076237512597937765737059080653379610398148;
    z[2] = -0.0014141518083258177510872439765585925278641690553467;
    *sf = 46080;
    break;
  case 7:
    np = 3;
    z[0] = -0.53528043079643816554240378168164607183392315234269;
    z[1] = -0.12255461519232669051527226435935734360548654942730;
    z[2] = -0.0091486948096082769285930216516478534156925639545994;
    *sf = 5040;
    break;
  default:
    throw SplinterpolatorException("SplineColumn::get_poles: invalid order of spline");
    break;
  }
  return(np);
}

/////////////////////////////////////////////////////////////////////
//
// Initialises the first value for the forward sweep. The initialisation
// will always be an approximation (this is where the "infinite" in IIR
// breaks down) so the value will be calculated to a predefined precision.
//
/////////////////////////////////////////////////////////////////////

template<class T>
double Splinterpolator<T>::SplineColumn::init_fwd_sweep(double z, ExtrapolationType et, double prec) const
{
  //
  // Move logs away from here after debugging
  //
  unsigned int n = static_cast<unsigned int>((log(prec)/log(abs(z))) + 1.5);
  n = (n > _sz) ? _sz : n;

  double iv = _col[0];
  if (et == Mirror) {
    double *ptr=&_col[1];
    double z2i=z;
    for (unsigned int i=1; i<n; i++, ptr++, z2i*=z) iv += z2i * *ptr;
  }
  else {
    double *ptr=&_col[_sz-1];
    double z2i=z;
    for (unsigned int i=1; i<n; i++, ptr--, z2i*=z) iv += z2i * *ptr;
  }
  return(iv); 
}

/////////////////////////////////////////////////////////////////////
//
// Initialises the first value for the backward sweep. The initialisation
// will always be an approximation (this is where the "infinite" in IIR
// breaks down) so the value will be calculated to a predefined precision.
//
/////////////////////////////////////////////////////////////////////

template<class T>
double Splinterpolator<T>::SplineColumn::init_bwd_sweep(double z, double lv, ExtrapolationType et, double prec) const
{
  double iv = 0.0;
  if (et == Mirror) {
    iv = -z/(1.0-z*z) * (2.0*_col[_sz-1] - lv);
  }
  else {
    unsigned int n = static_cast<unsigned int>((log(prec)/log(abs(z))) + 1.5);
    n = (n > _sz) ? _sz : n;
    iv = z * _col[_sz-1];
    double z2i = z*z;
    double *ptr=_col;
    for (unsigned int i=1; i<n; i++, ptr++, z2i*=z) {
      iv += z2i * *ptr;
    }
    iv /= (z2i-1.0);
  }
  return(iv);
}
} // End namespace SPLINTERPOLATOR

#endif // End #ifndef splinterpolator.h