3

I have a large image A and a smaller image B, both expressed as 2-D numpy arrays. I want to use A as the canvas, and write translated copies of B all over it, packed in a hexagonal arrangement. The part I can't get my head around is how to handle it such that the image wraps both vertically and horizontally—essentially what I want is regular tessellation of a (padded, as necessary) sub-image onto a torus.

I've seen the discussion of numpy.take and numpy.roll at wrapping around slices in Python / numpy and that shows me how to access and return a copy of a wrapped slice of an array, but I want to assign to that—i.e., for arbitrary integers rowOffset and columnOffset I want to do the equivalent of:

  A = numpy.zeros((5,11), int)
  B = numpy.array([[1,2,3,4,5,6,7]]) * numpy.array([[10,100,1000]]).T
  # OK, we wouldn't be able to fit more than one or two copies of B into A, but they demonstrate the wrapped placement problem

  wrappedRowIndices = ( numpy.arange(B.shape[0]) + rowOffset ) % A.shape[0]
  wrappedColumnIndices = ( numpy.arange(B.shape[1]) + columnOffset ) % A.shape[1]
  A[ wrappedRowIndices,  : ][ :, wrappedColumnIndices ] = B

I see from a comment on the question, and from a moment's reflection on the way numpy arrays are represented, that there's no way a wrapped slice can be returned as a view in the way this demands.

Is there (Y) a way of assigning to wrapped slices of an array in this way, or (X) an existing utility for performing the kind of tessellation I'm trying to achieve?

Community
  • 1
  • 1
jez
  • 14,867
  • 5
  • 37
  • 64
  • how will it be hexagonal? isn't `B` a square? – maxymoo Oct 24 '16 at 21:47
  • I meant that the centres of the copies of `B` will be distributed across `A` at the corners of packed hexagons, not that `B` itself is somehow hexagonal. I guess the easiest-to-handle variant is the one where the dimensions of `B` are smaller than the hexagon so that there is no overlap between copies. – jez Oct 24 '16 at 21:56
  • I guess if you have to use views from `strides`, you could append cols and rows into `A` and then index into it instead, memory permitting of course. To create the padded version : `A1 = np.column_stack((A,A[:,:B.shape[1]-1]))`, `A2 = np.row_stack((A1,A1[:B.shape[0]-1]))`. – Divakar Oct 24 '16 at 22:02

3 Answers3

3

np.put is a 1d equivalent to np.take:

In [1270]: A=np.arange(10)
In [1271]: np.take(A,[8,9,10,11],mode='wrapped')
Out[1271]: array([8, 9, 0, 1])
In [1272]: np.put(A,[8,9,10,11],[10,11,12,13],mode='wrapped')
In [1273]: A
Out[1273]: array([12, 13,  2,  3,  4,  5,  6,  7, 10, 11])
In [1274]: np.take(A,[8,9,10,11],mode='wrapped')
Out[1274]: array([10, 11, 12, 13])

Its docs suggest np.place and np.putmask (and np.copyto). I haven't used those much, but it might be possible to construct a mask, and rearrangement of B that would do the copy.

=================

Here's an experiment with place:

In [1313]: A=np.arange(24).reshape(4,6)
In [1314]: mask=np.zeros(A.shape,bool)
In [1315]: mask[:3,:4]=True
In [1316]: B=-np.arange(12).reshape(3,4)

So I have mask the same size as A, with a 'hole' the size of B.

I can roll both the mask and B, and place the values in A in a wrapped fashion.

In [1317]: np.place(A, np.roll(mask,-2,0), np.roll(B,1,0).flat)
In [1318]: A
Out[1318]: 
array([[ -8,  -9, -10, -11,   4,   5],
       [  6,   7,   8,   9,  10,  11],
       [  0,  -1,  -2,  -3,  16,  17],
       [ -4,  -5,  -6,  -7,  22,  23]])

And with 2d rolls

In [1332]: m=np.roll(np.roll(mask,-2,0),-1,1)
In [1333]: m
Out[1333]: 
array([[ True,  True,  True, False, False,  True],
       [False, False, False, False, False, False],
       [ True,  True,  True, False, False,  True],
       [ True,  True,  True, False, False,  True]], dtype=bool)
In [1334]: b=np.roll(np.roll(B,1,0),-1,1)
In [1335]: b
Out[1335]: 
array([[ -9, -10, -11,  -8],
       [ -1,  -2,  -3,   0],
       [ -5,  -6,  -7,  -4]])
In [1336]: A=np.zeros((4,6),int)
In [1337]: np.place(A, m, b.flat)
In [1338]: A
Out[1338]: 
array([[ -9, -10, -11,   0,   0,  -8],
       [  0,   0,   0,   0,   0,   0],
       [ -1,  -2,  -3,   0,   0,   0],
       [ -5,  -6,  -7,   0,   0,  -4]])
hpaulj
  • 221,503
  • 14
  • 230
  • 353
2

Your current code breaks down into a __getitem__ and a __setitem__. The __getitem__ does not return a view, as you noted, so the __setitem__ just ends up modifying the copy.

You need to do the whole thing in one __setitem__ (i.e. one set of brackets):

A[wrappedRowIndices[:,np.newaxis], wrappedColumnIndices] = B

Due to broadcasting, this is equivalent to:

A[wrappedRowIndices[:,np.newaxis], wrappedColumnIndices[np.newaxis,:]] = B

When indexing with more than one array, the rules is that:

# ... here is NOT the python Ellipsis!
y = x[a, b, c, ...]
y[i, j, ..] = x[a[i,j,...], b[i,j,...], ...]

There's actually a builtin for this, np.ix_():

A[np.ix_(wrappedRowIndices, wrappedColumnIndices)] = B

Generalizing to ND, you get:

def place_wrapped(canvas, brush, position):
    assert canvas.ndim == brush.ndim == len(position)
    ind = np.ix_(*(
        (np.arange(b_dim) + shift) % c_dim
        for b_dim, c_dim, shift in zip(brush.shape, canvas.shape, position)
    ))
    canvas[ind] = brush
Eric
  • 95,302
  • 53
  • 242
  • 374
  • This answer taught me one valuable thing—namely that you **can** after all take rectangular chunks of an array with only one pair of square brackets—but the assignment you suggest fails with `ValueError: array is not broadcastable to correct shape` – jez Oct 25 '16 at 01:06
  • If the number of elements is right you might fix things with a `flat` or two. – hpaulj Oct 25 '16 at 01:17
  • I would have thought the rule you stated would require `[np.newaxis,:]` rather than `[:, np.newaxis]` for the second indexing array. But both variants work (for *getting* a copy of the chunk) and both fail with the same `ValueError` (for setting the chunk content). Appending a `.flat` to the RHS makes no difference. Appending a `.flat` to the LHS gets rid of the exception but now it seems we're working on a `__getitem__` copy again: `A` remains unchanged. – jez Oct 25 '16 at 01:59
  • @jez: see my update. My slice was the transpose of what it should be. You could have spotted this by comparing `A[...].shape` to `B.shape` – Eric Oct 25 '16 at 11:06
1

Here's a function that works to solve my problem Y based on hpaulj's answer. This is probably not the most efficient way to solve X, however, because of all the heavy lifting done in numpy.roll.

import numpy

def place_wrapped(canvas, brush, position):
    mask = numpy.zeros(canvas.shape, bool)
    mask[[slice(extent) for extent in brush.shape]] = True
    for axis, shift in enumerate(position):
        canvas_extent = canvas.shape[axis]
        brush_extent = brush.shape[axis]
        shift %= canvas_extent
        if shift:
            mask = numpy.roll(mask, shift, axis=axis)
            nwrapped = shift + brush_extent - canvas_extent
            if nwrapped > 0: brush = numpy.roll(brush, nwrapped, axis=axis)
    numpy.place(canvas, mask, brush.flat)


A = numpy.zeros((5,11), int)
B = numpy.array([[1,2,3,4,5,6,7]]) * numpy.array([[10,100,1000]]).T
print(B)

rowOffset = 3
columnOffset = 7
place_wrapped(A, B, (rowOffset, columnOffset))
print(A)

Eric's corrected approach also worked, and I imagine it must be more efficient because it does not have to make copies of any data:

def place_wrapped2(canvas, brush, position):
    ind = [
        ((numpy.arange(brush.shape[axis]) + shift) % canvas.shape[axis]).reshape([
            extent if i == axis else 1
            for i, extent in enumerate(brush.shape)
        ])
        for axis, shift in enumerate(position)
    ]
    canvas[ind] = brush

A *= 0 # reset
place_wrapped2(A, B, (rowOffset, columnOffset))
print(A)
jez
  • 14,867
  • 5
  • 37
  • 64