ENH: Allow `Splinterpolator` instances to be created from an existing set of spline coefficients, and extend extrapolation options

• The Splinterpolator class can now accept a pointer to an array of pre-calculated coefficients - in this case, a copy of the pointer is used (the data is not copied), and pre-filtering (coefficient calculation) is skipped
• New SoftZeros extrapolation method - out of bounds coefficients will be set to zero, but out-of-bounds coordinates will still be constructed via interpolation. This differs from the existing Zeros extrapolation method, where out-of-bounds coordinates are forced to zero.
• Adjustment to the coefficient initialisation method used when Constant extrapolation is requested; previously, a mirroring boundary condition was used.
• New Splinterpolator::Copy function, which allows the spline coefficients to be copied/extracted.

splinterpolator.h

@@ -45,35 +45,67 @@ class Splinterpolator
 {
 public:
   // Constructors
-  Splinterpolator() : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(1) {}
-  Splinterpolator(const T *data, const std::vector<unsigned int>& dim, const std::vector<ExtrapolationType>& et, unsigned int order=3, bool copy_low_order=true, Utilities::NoOfThreads nthr=Utilities::NoOfThreads(1), double prec=1e-8) : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(nthr._n)
+  Splinterpolator()
+    : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(1) {}
+  Splinterpolator(const T *data_or_coefs,
+                  const std::vector<unsigned int>& dim,
+                  const std::vector<ExtrapolationType>& et,
+                  unsigned int order=3,
+                  bool copy_low_order=true,
+                  Utilities::NoOfThreads nthr=Utilities::NoOfThreads(1),
+                  double prec=1e-8,
+                  bool data_are_coefs=false)
+    : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(nthr._n)
   {
-    common_construction(data,dim,order,prec,et,copy_low_order);
+    common_construction(data_or_coefs,dim,order,prec,et,copy_low_order,data_are_coefs);
   }
-  Splinterpolator(const T *data, const std::vector<unsigned int>& dim, ExtrapolationType et=Zeros, unsigned int order=3, bool copy_low_order=true, Utilities::NoOfThreads nthr=Utilities::NoOfThreads(1), double prec=1e-8) : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(nthr._n)
+  Splinterpolator(const T *data_or_coefs,
+                  const std::vector<unsigned int>& dim,
+                  ExtrapolationType et=Zeros,
+                  unsigned int order=3,
+                  bool copy_low_order=true,
+                  Utilities::NoOfThreads
+                  nthr=Utilities::NoOfThreads(1),
+                  double prec=1e-8,
+                  bool data_are_coefs=false)
+    : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0), _nthr(nthr._n)
   {
     std::vector<ExtrapolationType>   ett(dim.size(),et);
-    common_construction(data,dim,order,prec,ett,copy_low_order);
+    common_construction(data_or_coefs,dim,order,prec,ett,copy_low_order,data_are_coefs);
   }
   // Copy construction. May be removed in future
-  Splinterpolator(const Splinterpolator<T>& src) : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0) { assign(src); }
+  Splinterpolator(const Splinterpolator<T>& src)
+    : _valid(false), _own_coef(false), _coef(0), _cptr(0), _ndim(0) { assign(src); }
 
   // Destructor
   ~Splinterpolator() { if(_own_coef) delete [] _coef; }
 
   // Assignment. May be removed in future
-  Splinterpolator& operator=(const Splinterpolator& src) { if(_own_coef) delete [] _coef; assign(src); return(*this); }
+  Splinterpolator& operator=(const Splinterpolator& src)
+  { if(_own_coef) delete [] _coef; assign(src); return(*this); }
 
   // 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, bool copy_low_order=true, double prec=1e-8)
+  void Set(const T *data_or_coefs,
+           const std::vector<unsigned int>& dim,
+           const std::vector<ExtrapolationType>& et,
+           unsigned int order=3,
+           bool copy_low_order=true,
+           double prec=1e-8,
+           bool data_are_coefs=false)
   {
     if (_own_coef) delete [] _coef;
-    common_construction(data,dim,order,prec,et,copy_low_order);
+    common_construction(data_or_coefs,dim,order,prec,et,copy_low_order,data_are_coefs);
   }
-  void Set(const T *data, const std::vector<unsigned int>& dim, ExtrapolationType et, unsigned int order=3, bool copy_low_order=true, double prec=1e-8)
+  void Set(const T *data_or_coefs,
+           const std::vector<unsigned int>& dim,
+           ExtrapolationType et,
+           unsigned int order=3,
+           bool copy_low_order=true,
+           double prec=1e-8,
+           bool data_are_coefs=false)
   {
     std::vector<ExtrapolationType>  vet(dim.size(),Zeros);
-    Set(data,dim,vet,order,copy_low_order,prec);
+    Set(data_or_coefs,dim,vet,order,copy_low_order,prec,data_are_coefs);
   }
 
   // Return interpolated value
@@ -204,12 +236,18 @@ private:
   //
   // 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, bool copy);
+  void common_construction(const T *data_or_coefs,
+                           const std::vector<unsigned int>& dim,
+                           unsigned int order,
+                           double prec,
+                           const std::vector<ExtrapolationType>& et,
+                           bool copy,
+                           bool data_are_coefs);
   void assign(const Splinterpolator<T>& src);
-  bool calc_coef(const T *data, bool copy);
+  bool calc_coef(const T *data_or_coefs, bool copy, bool data_are_coefs);
   void deconv_along(unsigned int dim);
-  void deconv_along_mt_helper(unsigned int dim, unsigned int mdim, unsigned int mstep, unsigned int offset, unsigned int step, 
-			      const std::vector<unsigned int>& rdim, const std::vector<unsigned int>& rstep);						
+  void deconv_along_mt_helper(unsigned int dim, unsigned int mdim, unsigned int mstep, unsigned int offset, unsigned int step,
+			      const std::vector<unsigned int>& rdim, const std::vector<unsigned int>& rstep);
   T coef(int *indx) const;
   const T* coef_ptr() const {if (_own_coef) return(_coef); else return(_cptr); }
   unsigned int indx2indx(int indx, unsigned int d) const;
@@ -1329,14 +1367,21 @@ unsigned int Splinterpolator<T>::n_nonzero(const unsigned int *vec) const
 /////////////////////////////////////////////////////////////////////
 
 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, bool copy)
+void Splinterpolator<T>::common_construction(
+  const T *data_or_coefs,
+  const std::vector<unsigned int>& dim,
+  unsigned int order,
+  double prec,
+  const std::vector<ExtrapolationType>& et,
+  bool copy,
+  bool data_are_coefs)
 {
   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");
+  if (!data_or_coefs) throw SplinterpolatorException("common_construction: zero data pointer");
 
   _order = order;
   _prec = prec;
@@ -1345,8 +1390,7 @@ void Splinterpolator<T>::common_construction(const T *data, const std::vector<un
   _ndim = dim.size();
   for (unsigned int i=0; i<5; i++) _dim[i]  = (i < dim.size()) ? dim[i] : 1;
 
-  _own_coef = calc_coef(data,copy);
-
+  _own_coef = calc_coef(data_or_coefs,copy,data_are_coefs);
   _valid = true;
 }
 
@@ -1385,16 +1429,16 @@ void Splinterpolator<T>::assign(const Splinterpolator<T>& src)
 /////////////////////////////////////////////////////////////////////
 
 template<class T>
-bool Splinterpolator<T>::calc_coef(const T *data, bool copy)
+bool Splinterpolator<T>::calc_coef(const T *data_or_coefs, bool copy, bool data_are_coefs)
 {
-  if (_order < 2 && !copy) { _cptr = data; return(false); }
+  if (_order < 2 && !copy) { _cptr = data_or_coefs; return(false); }
 
   // Allocate memory and put the original data into _coef
   //
   unsigned int ts=1;
   for (unsigned int i=0; i<_dim.size(); i++) ts *= _dim[i];
   _coef = new T[ts];
-  memcpy(_coef,data,ts*sizeof(T));
+  memcpy(_coef,data_or_coefs,ts*sizeof(T));
 
   if (_order < 2) return(true);  // If nearest neighbour or linear, that's all we need
 
@@ -1464,8 +1508,8 @@ void Splinterpolator<T>::deconv_along(unsigned int dim)
 
 template<class T>
 void Splinterpolator<T>::deconv_along_mt_helper(unsigned int                     dim,
-						unsigned int                     mdim, 
-						unsigned int                     mstep, 
+						unsigned int                     mdim,
+						unsigned int                     mstep,
 						unsigned int                     offset, // Offset into parallel dimension
 						unsigned int                     step,   // Step size in parallel dimension
 						const std::vector<unsigned int>& rdim,