7

When using scipy.ndimage.interpolation.shift to shift a numpy data array along one axis with periodic boundary treatment (mode = 'wrap'), I get an unexpected behavior. The routine tries to force the first pixel (index 0) to be identical to the last one (index N-1) instead of the "last plus one (index N)".

Minimal example:

# module import
import numpy as np
from scipy.ndimage.interpolation import shift
import matplotlib.pyplot as plt

# print scipy.__version__
# 0.18.1

a = range(10)

plt.figure(figsize=(16,12))

for i, shift_pix in enumerate(range(10)):
    # shift the data via spline interpolation
    b = shift(a, shift=shift_pix, mode='wrap')

    # plotting the data
    plt.subplot(5,2,i+1)
    plt.plot(a, marker='o', label='data')
    plt.plot(np.roll(a, shift_pix), marker='o', label='data, roll')
    plt.plot(b, marker='o',label='shifted data')
    if i == 0:
        plt.legend(loc=4,fontsize=12)
    plt.ylim(-1,10)
    ax = plt.gca()
    ax.text(0.10,0.80,'shift %d pix' % i, transform=ax.transAxes)

Blue line: data before the shift
Green line: expected shift behavior
Red line: actual shift output of scipy.ndimage.interpolation.shift

Is there some error in how I call the function or how I understand its behavior with mode = 'wrap'? The current results are in contrast to the mode parameter description from the related scipy tutorial page and from another StackOverflow post. Is there an off-by-one-error in the code?

Scipy version used is 0.18.1, distributed in anaconda-2.2.0

enter image description here

Drise
  • 4,310
  • 5
  • 41
  • 66
bproxauf
  • 1,076
  • 12
  • 23
  • Scipy version used is 0.18.1, distributed in anaconda-2.2.0 – bproxauf Mar 14 '18 at 14:33
  • don't comment to add clarifications, edit instead (note: already did that for ya) – Jean-François Fabre Mar 14 '18 at 14:42
  • It looks like `shift()` performs the boundary condition after the data has been shifted. So with a shift_pix of 1, it transforms `a` into `[0, 1, 2, 3, 4, 5, 6, 7, 8]`, and *then* does a wrap, bringing `8` to the front. I don't know why it does this though. – Arthur Dent Mar 14 '18 at 16:04
  • @ArthurDent: I'm not sure this is actually the case. Notice that the output values are `[8, 1, 2, 3, 4, 5, 6, 7, 8, 9]`. If I understand correctly what you suggested the wrap would rather give `[8, 0, 1, 2, 3, 4, 5, 6, 7, 8]`. – bproxauf Mar 15 '18 at 07:37
  • @bproxauf Is `[8, 0, 1, 2, 3, 4, 5, 6, 7, 8]` not what you're getting? The second element of the shifted data in your second plot is equal to 0, not 1. When I do `b = shift(a, shift=1, mode='wrap')`, it yields `[8, 0, 1, 2, 3, 4, 5, 6, 7, 8]`. Also, I spent more time than I would like to admit looking at the source code for the shift function, but it's in C and a little too complicated for me to decipher. From what I could tell, it looks like it *should* perform the wrap first, but I couldn't be certain. If you would like to look at it let me know. – Arthur Dent Mar 15 '18 at 16:38
  • @ArthurDent: You are completely right, I misread my own plot axis. Thus what you suggest seems very plausible. As far as I can see from the plot the shifted data arrays never have the value 9, *except* for a shift of 0 or 9 pixels. I am not very solid in C code myself. I already had a brief look, but I could not figure out what is going on either. – bproxauf Mar 16 '18 at 10:09
  • Does anyone in this thread know what accuracy this function shifts data by? It doesn't appear you can specify the factor to zoom by, so I assume it is hard coded. I did a quick search through the C code, but my knowledge in C is pretty limited. – Will.Evo Jun 18 '19 at 16:38

2 Answers2

1

It seems that the behaviour you have observed is intentional.

The cause of the problem lies in the C function map_coordinate which translates the coordinates after shift to ones before shift:

map_coordinate(double in, npy_intp len, int mode)

The function is used as the subroutine in NI_ZoomShift that does the actual shift. Its interesting part looks like this:

enter image description here

Example. Lets see how the output for output = shift(np.arange(10), shift=4, mode='wrap') (from the question) is computed.

NI_ZoomShift computes edge values output[0] and output[9] in some special way, so lets take a look at computation of output[1] (a bit simplified):

# input  =         [0,1,2,3,4,5,6,7,8,9]
# output = [ ,?, , , , , , , , ]          '?' == computed position
# shift  = 4
output_index = 1

in  = output_index - shift    # -3
sz  = 10 - 1                  # 9
in += sz * ((-5 / 9) + 1)
#  +=  9 * ((     0) + 1) == 9
# in == 6

return input[in]  # 6

It is clear that sz = len - 1 is responsible for the behaviour you have observed. It was changed from sz = len in a suggestively named commit dating back to 2007: Fix off-by-on errors in ndimage boundary routines. Update tests.

I don't know why such change was introduced. One of the possible explanations that come to my mind is as follows:

Function 'shift' uses splines for interpolation. A knot vector of an uniform spline on interval [0, k] is simply [0,1,2,...,k]. When we say that the spline should wrap, it is natural to require equality on values for knots 0 and k, so that many copies of the spline could be glued together, forming a periodic function:

0--1--2--3-...-k              0--1--2--3-...-k              0--1-- ...
               0--1--2--3-...-k              0--1--2--3-...-k      ...

Maybe shift just treats its input as a list of values for spline's knots?

Radek
  • 846
  • 4
  • 17
1

It is worth noting that this behavior appears to be a bug, as noted in this SciPy issue: https://github.com/scipy/scipy/issues/2640

The issue appears to effect every extrapolation mode in scipy.ndimage other than mode='mirror'.

shoyer
  • 9,165
  • 1
  • 37
  • 55