Find a 3x3 sliding window over an image

Question

I have an image.

I want to obtain a 3x3 window (neighbouring pixels) for every pixel in the image.

I have this Python code:

for x in range(2,r-1,1):
    for y in range(2,c-1,1):
        mask5=numpy.array([cv.Get2D(copy_img,x-1,y-1),cv.Get2D(copy_img,x-1,y),cv.Get2D(copy_img,x-1,y+1),cv.Get2D(copy_img,x,y-1),cv.Get2D(copy_img,x,y),cv.Get2D(copy_img,x,y+1),cv.Get2D(copy_img,x+1,y-1),cv.Get2D(copy_img,x+1,y),cv.Get2D(copy_img,x+1,y+1)])
        cent=[cv.Get2D(copy_img,x,y)]

mask5 is the 3x3 window. cent is the center pixel.

Is there a more efficient way to do this - i.e. using maps, iterators - anything but the two nested loops I've used?

Aalii

What is your intention? Probably you want to perform a convolution? Tell us what youa re going to do with mask5, then we can help you better, cheers. – Differentiate

@fraxel: After I get the window, I need to sort the pixels in the 3x3 window by intensity, and create another (one-dimensional) sliding window over this list, and based on a complex condition on the means of the pixels in these (1-d) slides, modify the original center pixel accordingly. – Aalii

Sounds like you might be trying to do adaptive threshold? Take a look at that, it may do what you want. – Differentiate

@frexel: Not really. Is there a way to simply optimize the loops above, as shown - using maps or iterators - anything to avoid using for loops? – Aalii

What exactly that loop does ? Is it just copying a 3x3 matrix centered at (x,y) to mask y ? And do you want to do it for all the pixels? Or just one pixel at (x,y)? – Blasting

oh well.. take 3*3 slices of the array instead of getting 9 points. Use your image as a numpy array and replace your main line with: copy_img[x:x+3,y:y+3] – Differentiate

@AbidRahmanK: The loop just copies a 3x3 matrix centered at (x,y) to the variable. This is just dummy sample code. There will be more code added there later. I obviously want to do this for all pixels - hence the loop! – Aalii

@Differentiate : Yeah, that is i thought too.. – Blasting

@fraxel: That's fine, but doesn't do away with my double for loop. – Aalii

@VelvetGhost : check my answer. Is this what you wanted? – Blasting

Without a full understanding of what you are doing, thats the best I can do! :) – Differentiate

This is an old question, but since it still comes up in google I'll note there is support for this in numpy now: numpy.org/devdocs/reference/generated/… – Felicitasfelicitate

Answer 1

This can be done faster, by reshaping and swapping axes, and then repeating over all kernel elements, like this:

im = np.arange(81).reshape(9,9)
print np.swapaxes(im.reshape(3,3,3,-1),1,2)

This gives you an array of 3*3 tiles which tessalates across the surface:

[[[[ 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]]]

 [[[27 28 29]   [[30 31 32]   [[33 34 35]
   [36 37 38]    [39 40 41]    [42 43 44]
   [45 46 47]]   [48 49 50]]   [51 52 53]]]

 [[[54 55 56]   [[57 58 59]   [[60 61 62]
   [63 64 65]    [66 67 68]    [69 70 71]
   [72 73 74]]   [75 76 77]]   [78 79 80]]]]

To get the overlapping tiles we need to repeat this 8 further times, but 'wrapping' the array, by using a combination of vstack and column_stack. Note that the right and bottom tile arrays wrap around (which may or may not be what you want, depending on how you are treating edge conditions):

im =  np.vstack((im[1:],im[0]))
im =  np.column_stack((im[:,1:],im[:,0]))
print np.swapaxes(im.reshape(3,3,3,-1),1,2)

#Output:
[[[[10 11 12]   [[13 14 15]   [[16 17  9]
   [19 20 21]    [22 23 24]    [25 26 18]
   [28 29 30]]   [31 32 33]]   [34 35 27]]]

 [[[37 38 39]   [[40 41 42]   [[43 44 36]
   [46 47 48]    [49 50 51]    [52 53 45]
   [55 56 57]]   [58 59 60]]   [61 62 54]]]

 [[[64 65 66]   [[67 68 69]   [[70 71 63]
   [73 74 75]    [76 77 78]    [79 80 72]
   [ 1  2  3]]   [ 4  5  6]]   [ 7  8  0]]]]

Doing it this way you wind up with 9 sets of arrays, so you then need to zip them back together. This, and all the reshaping generalises to this (for arrays where the dimensions are divisible by 3):

def new(im):
    rows,cols = im.shape
    final = np.zeros((rows, cols, 3, 3))
    for x in (0,1,2):
        for y in (0,1,2):
            im1 = np.vstack((im[x:],im[:x]))
            im1 = np.column_stack((im1[:,y:],im1[:,:y]))
            final[x::3,y::3] = np.swapaxes(im1.reshape(rows/3,3,cols/3,-1),1,2)
    return final

Comparing this new function to looping through all the slices (below), using timeit, its about 4 times faster, for a 300*300 array.

def old(im):
    rows,cols = im.shape
    s = []
    for x in xrange(1,rows):
        for y in xrange(1,cols):
            s.append(im[x-1:x+2,y-1:y+2])
    return s

Differentiate

Thanks. I managed to reduce all of this to precisely one short line of code - unfortunately, in MATLAB (im2col function). Does this function have a direct equivalent in Python? – Aalii

@VelvetGhost im2col is just a function that wraps an implementation like mine or the one above (albeit with a bit more polish). The MATLAB implementation doesn't magically remove code. If you can use your own function, so much the better as you're not reliant on an expensive toolbox for it to work (such is the joy of MATLAB). Also, don't assume that because it's in a toolbox, it's fast or efficient. – Mcglothlin

@HenryGomersall: Thanks! Yes, I do know that MATLAB doesn't magically remove code. I just found it easier and quicker when writing code under a deadline. Thanks very much though for your Python implementation. It'll help me learn more about Python when I have a bit of time. – Aalii

Here's a question: how could we modify this method to work for arbitrary (instead of multiple-of-3) image dimensions? – Sennit

@Magsol, I've created some cython code for "windowizing" 3D image volumes, which works for any window size. You could modify this easily to work for 2D images. – Candicandia

Answer 2

I think the following does what you are after. The loop is only over the 9 elements. I'm sure there is a way of vectorizing it, but it's probably not worth the effort.

import numpy

im = numpy.random.randint(0,50,(5,7))

# idx_2d contains the indices of each position in the array
idx_2d = numpy.mgrid[0:im.shape[0],0:im.shape[1]]

# We break that into 2 sub arrays
x_idx = idx_2d[1]
y_idx = idx_2d[0]

# The mask is used to ignore the edge values (or indeed any values).
mask = numpy.ones(im.shape, dtype='bool')
mask[0, :] = False
mask[:, 0] = False
mask[im.shape[0] - 1, :] = False
mask[:, im.shape[1] - 1] = False

# We create and fill an array that contains the lookup for every
# possible 3x3 array.
idx_array = numpy.zeros((im[mask].size, 3, 3), dtype='int64')

# Compute the flattened indices for each position in the 3x3 grid
for n in range(0, 3):
    for m in range(0, 3):
        # Compute the flattened indices for each position in the 
        # 3x3 grid
        idx = (x_idx + (n-1)) + (y_idx  + (m-1)) * im.shape[1]

        # mask it, and write it to the big array
        idx_array[:, m, n] = idx[mask]


# sub_images contains every valid 3x3 sub image
sub_images = im.ravel()[idx_array]

# Finally, we can flatten and sort each sub array quickly
sorted_sub_images = numpy.sort(sub_images.reshape((idx[mask].size, 9)))

Mcglothlin

Answer 3

Try the following code as matlab function im2col(...)

import numpy as np

def im2col(Im, block, style='sliding'):
    """block = (patchsize, patchsize)
        first do sliding
    """
    bx, by = block
    Imx, Imy = Im.shape
    Imcol = []
    for j in range(0, Imy):
        for i in range(0, Imx):
            if (i+bx <= Imx) and (j+by <= Imy):
                Imcol.append(Im[i:i+bx, j:j+by].T.reshape(bx*by))
            else:
                break
    return np.asarray(Imcol).T

if __name__ == '__main__':
    Im = np.reshape(range(6*6), (6,6))
    patchsize = 3
    print Im
    out =  im2col(Im, (patchsize, patchsize))
    print out
    print out.shape
    print len(out)

Rothschild