Classes for spline interpolation in up to 5 dimesnsions

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("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