diff --git a/splinterpolator.h b/splinterpolator.h
new file mode 100644
index 0000000000000000000000000000000000000000..c41efde431354301bd04721dc492b6fd8bfb741f
--- /dev/null
+++ b/splinterpolator.h
@@ -0,0 +1,201 @@
+//
+// 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