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

@@ -64,8 +64,7 @@ public:
                   ExtrapolationType et=Zeros,
                   unsigned int order=3,
                   bool copy_low_order=true,
-                  Utilities::NoOfThreads
-                  nthr=Utilities::NoOfThreads(1),
+                  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)
@@ -84,6 +83,22 @@ public:
   Splinterpolator& operator=(const Splinterpolator& src)
   { if(_own_coef) delete [] _coef; assign(src); return(*this); }
 
+  // Copy the spline coefficients into dest.
+  void Copy(std::vector<T>& dest)
+  {
+    unsigned int N = 1;
+    for (auto d : _dim) {
+      N *= d;
+    }
+    dest.resize(N);
+
+    auto coefs = coef_ptr();
+
+    for (unsigned int i = 0; i < N; i++) {
+       dest[i] = coefs[i];
+    }
+  }
+
   // Set new data in Splinterpolator.
   void Set(const T *data_or_coefs,
            const std::vector<unsigned int>& dim,
@@ -797,18 +812,10 @@ unsigned int Splinterpolator<T>::get_start_indicies(const double *coord, int *si
 {
   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=0; i<_ndim; i++) {
+    sinds[i] = std::ceil(coord[i]) - ni/2;
   }
+
   for (unsigned int i=_ndim; i<5; i++) sinds[i] = 0;
 
   return(ni);
@@ -1491,7 +1498,7 @@ void Splinterpolator<T>::deconv_along(unsigned int dim)
     ss *= _dim[i];
   }
 
-  if (_nthr==1) { // If we are to run single-threaded
+  if (_nthr<=1) { // If we are to run single-threaded
     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++) {
@@ -1631,11 +1638,17 @@ double Splinterpolator<T>::SplineColumn::init_fwd_sweep(double z, ExtrapolationT
     double z2i=z;
     for (unsigned int i=1; i<n; i++, ptr--, z2i*=z) iv += z2i * *ptr;
   }
-  else {
+  else if (et == Mirror) {
     double *ptr=&_col[1];
     double z2i=z;
     for (unsigned int i=1; i<n; i++, ptr++, z2i*=z) iv += z2i * *ptr;
   }
+  // et == Constant || et == Zeros
+  else {
+    double *ptr=&_col[0];
+    double z2i=z;
+    for (unsigned int i=0; i<n; i++, ptr++, z2i*=z) iv += z2i * *ptr;
+  }
   return(iv);
 }
 
@@ -1662,9 +1675,14 @@ double Splinterpolator<T>::SplineColumn::init_bwd_sweep(double z, double lv, Ext
     }
     iv /= (z2i-1.0);
   }
-  else {
+  else if (et == Mirror) {
     iv = -z/(1.0-z*z) * (2.0*_col[_sz-1] - lv);
   }
+  // et == Constant || et == Zeros
+  else  {
+    iv = z / (z - 1) * _col[_sz-1];
+  }
+
   return(iv);
 }