BF: Fix bug in code controlling which adjacent spline coefficients to use when interpolating

This MR fixes a very subtle bug in the Splinterpolator class which was causing artifacts when extrapolating into negative coordinate space.

When using third-order cubic b-spline interpolation to interpolate a value at a particular location, the four nearest spline coefficients are combined to generate the interpolated value. For example, if interpolating a value at location -0.25, the spline coefficients at locations [-2, -1, 0, 1] would be used.

If interpolating at an integer location, e.g. 2, the choice of whether to use coefficients at [0, 1, 2, 3], or [1, 2, 3, 4] is arbitrary, but the convention in the Splinterpolator class is to use the former.

The logic which controls which spline coefficients are used is in the get_start_indicies function. While this function produces sensible values for positive coordinates, it was mis-behaving for negative coordinates.

• For coordinates with a fractional component greater than 0.5, such as -3.75, the code was correctly selecting coefficients [-5, -4, -3, -2]. However, for coordinates with a fractional component smaller than 0.5, such as -3.25, the code was incorrectly selecting coefficients [-4, -3, -2, -1].
• For a negative integer coordinate such as -4, the code was selecting coefficients [-5, -4, -3, -2], which would produce the correct result, but which was not following the same convention as for positive coordinates (i.e. it should be [-6, -5, -4, -3]).

splinterpolator.h

@@ -759,18 +759,15 @@ 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;
+  for (unsigned int i=0; i<_ndim; i++) {
+    if (odd(ni)) {
+      sinds[i] = std::floor(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;
+    else {
+      sinds[i] = std::ceil(coord[i]) - ni/2;
     }
   }
+
   for (unsigned int i=_ndim; i<5; i++) sinds[i] = 0;
 
   return(ni);