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.

Community
  • 1
  • 1
Jiby
  • 1,865
  • 1
  • 13
  • 22

5 Answers5

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
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
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]]]
Community
  • 1
  • 1
Ed Smith
  • 12,716
  • 2
  • 43
  • 55
  • 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 04 '15 at 12:03
  • Great, glad it's efficient too. I only found scipy shift recently when answering another question (http://stackoverflow.com/questions/30399534/shift-elements-in-a-numpy-array/30534478#30534478). 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 04 '15 at 14:39
  • The `scipy` shift uses spline interpolation. So it is calculating new values, not just moving the existing ones around. – hpaulj Jun 04 '15 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. – Warren Weckesser Jun 04 '15 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 '21 at 08: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

john k
  • 6,268
  • 4
  • 55
  • 59
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

YXD
  • 31,741
  • 15
  • 75
  • 115
  • 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 04 '15 at 11:21
  • Yes, this method would copy once and then the slicing would just be a view of that data. – YXD Jun 04 '15 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 04 '15 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]]])
Lee
  • 29,398
  • 28
  • 117
  • 170