13

I need to shift a 3D array by a 3D vector of displacement for an algorithm. As of now I'm using this (admitedly very ugly) method :

shiftedArray = np.roll(np.roll(np.roll(arrayToShift, shift[0], axis=0)
                                     , shift[1], axis=1),
                             shift[2], axis=2)  

Which works, but means I'm calling 3 rolls ! (58% of my algorithm time is spent in these, according to my profiling)

From the docs of Numpy.roll:

Parameters:
shift : int

axis : int, optional

No mention of array-like in parameter ... So I can't have a multidimensional rolling ?

I thought I could just call a this kind of function (sounds like a Numpy thing to do) :

np.roll(arrayToShift,3DshiftVector,axis=(0,1,2))

Maybe with a flattened version of my array reshaped ? but then how do I compute the shift vector ? and is this shift really the same ?

I'm surprised to find no easy solution for this, as I thought this would be a pretty common thing to do (okay, not that common, but ...)

So how do we --relatively-- efficiently shift a ndarray by a N-Dimensional vector ?


Note: This question was asked in 2015, back when numpy's roll method did not support this feature.

5 Answers 5

11

In theory, using scipy.ndimage.interpolation.shift as described by @Ed Smith should work, but because of a bug (https://github.com/scipy/scipy/issues/1323), it didn't give a result that is equivalent to multiple calls of np.roll.


UPDATE: "Multi-roll" capability was added to numpy.roll in numpy version 1.12.0. Here's a two-dimensional example, in which the first axis is rolled one position and the second axis is rolled three positions:

In [7]: x = np.arange(20).reshape(4,5)

In [8]: x
Out[8]: 
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19]])

In [9]: numpy.roll(x, [1, 3], axis=(0, 1))
Out[9]: 
array([[17, 18, 19, 15, 16],
       [ 2,  3,  4,  0,  1],
       [ 7,  8,  9,  5,  6],
       [12, 13, 14, 10, 11]])

This makes the code below obsolete. I'll leave it there for posterity.


The code below defines a function I call multiroll that does what you want. Here's an example in which it is applied to an array with shape (500, 500, 500):

In [64]: x = np.random.randn(500, 500, 500)

In [65]: shift = [10, 15, 20]

Use multiple calls to np.roll to generate the expected result:

In [66]: yroll3 = np.roll(np.roll(np.roll(x, shift[0], axis=0), shift[1], axis=1), shift[2], axis=2)

Generate the shifted array using multiroll:

In [67]: ymulti = multiroll(x, shift)

Verify that we got the expected result:

In [68]: np.all(yroll3 == ymulti)
Out[68]: True

For an array this size, making three calls to np.roll is almost three times slower than a call to multiroll:

In [69]: %timeit yroll3 = np.roll(np.roll(np.roll(x, shift[0], axis=0), shift[1], axis=1), shift[2], axis=2)
1 loops, best of 3: 1.34 s per loop

In [70]: %timeit ymulti = multiroll(x, shift)
1 loops, best of 3: 474 ms per loop

Here's the definition of multiroll:

from itertools import product
import numpy as np


def multiroll(x, shift, axis=None):
    """Roll an array along each axis.

    Parameters
    ----------
    x : array_like
        Array to be rolled.
    shift : sequence of int
        Number of indices by which to shift each axis.
    axis : sequence of int, optional
        The axes to be rolled.  If not given, all axes is assumed, and
        len(shift) must equal the number of dimensions of x.

    Returns
    -------
    y : numpy array, with the same type and size as x
        The rolled array.

    Notes
    -----
    The length of x along each axis must be positive.  The function
    does not handle arrays that have axes with length 0.

    See Also
    --------
    numpy.roll

    Example
    -------
    Here's a two-dimensional array:

    >>> x = np.arange(20).reshape(4,5)
    >>> x 
    array([[ 0,  1,  2,  3,  4],
           [ 5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14],
           [15, 16, 17, 18, 19]])

    Roll the first axis one step and the second axis three steps:

    >>> multiroll(x, [1, 3])
    array([[17, 18, 19, 15, 16],
           [ 2,  3,  4,  0,  1],
           [ 7,  8,  9,  5,  6],
           [12, 13, 14, 10, 11]])

    That's equivalent to:

    >>> np.roll(np.roll(x, 1, axis=0), 3, axis=1)
    array([[17, 18, 19, 15, 16],
           [ 2,  3,  4,  0,  1],
           [ 7,  8,  9,  5,  6],
           [12, 13, 14, 10, 11]])

    Not all the axes must be rolled.  The following uses
    the `axis` argument to roll just the second axis:

    >>> multiroll(x, [2], axis=[1])
    array([[ 3,  4,  0,  1,  2],
           [ 8,  9,  5,  6,  7],
           [13, 14, 10, 11, 12],
           [18, 19, 15, 16, 17]])

    which is equivalent to:

    >>> np.roll(x, 2, axis=1)
    array([[ 3,  4,  0,  1,  2],
           [ 8,  9,  5,  6,  7],
           [13, 14, 10, 11, 12],
           [18, 19, 15, 16, 17]])

    """
    x = np.asarray(x)
    if axis is None:
        if len(shift) != x.ndim:
            raise ValueError("The array has %d axes, but len(shift) is only "
                             "%d. When 'axis' is not given, a shift must be "
                             "provided for all axes." % (x.ndim, len(shift)))
        axis = range(x.ndim)
    else:
        # axis does not have to contain all the axes.  Here we append the
        # missing axes to axis, and for each missing axis, append 0 to shift.
        missing_axes = set(range(x.ndim)) - set(axis)
        num_missing = len(missing_axes)
        axis = tuple(axis) + tuple(missing_axes)
        shift = tuple(shift) + (0,)*num_missing

    # Use mod to convert all shifts to be values between 0 and the length
    # of the corresponding axis.
    shift = [s % x.shape[ax] for s, ax in zip(shift, axis)]

    # Reorder the values in shift to correspond to axes 0, 1, ..., x.ndim-1.
    shift = np.take(shift, np.argsort(axis))

    # Create the output array, and copy the shifted blocks from x to y.
    y = np.empty_like(x)
    src_slices = [(slice(n-shft, n), slice(0, n-shft))
                  for shft, n in zip(shift, x.shape)]
    dst_slices = [(slice(0, shft), slice(shft, n))
                  for shft, n in zip(shift, x.shape)]
    src_blks = product(*src_slices)
    dst_blks = product(*dst_slices)
    for src_blk, dst_blk in zip(src_blks, dst_blks):
        y[dst_blk] = x[src_blk]

    return y
3

I think scipy.ndimage.interpolation.shift will do what you want, from the docs

shift : float or sequence, optional

The shift along the axes. If a float, shift is the same for each axis. If a sequence, shift should contain one value for each axis.

Which means you can do the following,

from scipy.ndimage.interpolation import shift
import numpy as np

arrayToShift = np.reshape([i for i in range(27)],(3,3,3))

print('Before shift')
print(arrayToShift)

shiftVector = (1,2,3)
shiftedarray = shift(arrayToShift,shift=shiftVector,mode='wrap')

print('After shift')
print(shiftedarray)

Which yields,

Before shift
[[[ 0  1  2]
  [ 3  4  5]
  [ 6  7  8]]

 [[ 9 10 11]
  [12 13 14]
  [15 16 17]]

 [[18 19 20]
  [21 22 23]
  [24 25 26]]]
After shift
[[[16 17 16]
  [13 14 13]
  [10 11 10]]

 [[ 7  8  7]
  [ 4  5  4]
  [ 1  2  1]]

 [[16 17 16]
  [13 14 13]
  [10 11 10]]]
5
  • Nice ! I guess I was obsessed with numpy's roll and googled around that. Only when I wrote this SO question did I state it as a "shift", which would have led me straight to this function ! Bit of an edit though, since it uses splines for non-round shifts (but I use purely integers), I had to specify order=0 (otherwise this method took 100 times longer than my old method ! =p)
    – Jiby
    Jun 4, 2015 at 12:03
  • Great, glad it's efficient too. I only found scipy shift recently when answering another question (stackoverflow.com/questions/30399534/…). Fortran has an intrinsic shift functionality so I guess scipy uses this which is why it can be much faster than numpy. Guess the order=0 avoids splines...
    – Ed Smith
    Jun 4, 2015 at 14:39
  • The scipy shift uses spline interpolation. So it is calculating new values, not just moving the existing ones around.
    – hpaulj
    Jun 4, 2015 at 14:58
  • 2
    This is not equivalent to multiple uses of np.roll. For example, note that your input contains 0, but your output does not. Jun 4, 2015 at 22:15
  • @WarrenWeckesser the discrepancy in the results is due to the inappropriate mode='wrap'. mode='grid-wrap' is needed to replicate the behavior of np.roll
    – Bastian
    Jun 25, 2021 at 8:35
1

np.roll takes in multiple dimensions. Just do

np.roll(arrayToShift, (shift[0], shift[1], shift[2]), axis=(0,1,2))

It's not very smart, so you have to specify the axis

0

I believe roll is slow because the rolled array can't be expressed as a view of the original data as a slice or reshape operation can. So data is copied every time. For background see: https://scipy-lectures.github.io/advanced/advanced_numpy/#life-of-ndarray

What may be worth trying it to first pad your array (with 'wrap' mode) and then use slices on the padded array to get shiftedArray: http://docs.scipy.org/doc/numpy/reference/generated/numpy.pad.html

3
  • True, and I should have stated that I started with the assumption that a (single) copy of the array was okay to do, but rolling thrice in a row would indeed mean 2 redundant copies, that I want to avoid.
    – Jiby
    Jun 4, 2015 at 11:21
  • Yes, this method would copy once and then the slicing would just be a view of that data.
    – YXD
    Jun 4, 2015 at 11:23
  • np.roll uses arange to construct indexes tuple, and then a take. So it's doing the slower advanced indexing. In theory you could construct your own 3d indexing tuple, and apply it just once. But that could be a lot of work.
    – hpaulj
    Jun 4, 2015 at 14:51
0

take in wrap mode can be used and I think that it does not change the array in memory.

Here is an implementation using @EdSmith's inputs:

arrayToShift = np.reshape([i for i in range(27)],(3,3,3))
shiftVector = np.array((1,2,3))
ind = 3-shiftVector    
np.take(np.take(np.take(arrayToShift,range(ind[0],ind[0]+3),axis=0,mode='wrap'),range(ind[1],ind[1]+3),axis=1,mode='wrap'),range(ind[2],ind[2]+3),axis=2,mode='wrap')

which gives the same as the OP's:

np.roll(np.roll(np.roll(arrayToShift, shift[0], axis=0) , shift[1], axis=1),shift[2], axis=2) 

gives:

array([[[21, 22, 23],
        [24, 25, 26],
        [18, 19, 20]],

       [[ 3,  4,  5],
        [ 6,  7,  8],
        [ 0,  1,  2]],

       [[12, 13, 14],
        [15, 16, 17],
        [ 9, 10, 11]]])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.