Sliding a subarray along a 2d array

Question

This is something I've been struggling with for a couple of weeks. The algorithm is the following:

1. Select a subarray as an array of rows and columns from a larger array
2. Compute the median of the subarray
3. Replace cells in subarray with median value
4. Move the subarray to the right by its own length
5. Repeat to end of array
6. Move subarray down by its own height
7. Repeat

I've got steps 1 to 3 as follows:

import numpy as np
w1 = np.arange(100).reshape(10,10)
side = 3
patch = w1[0:side, 0:side]

i, j = patch.shape
for j in range(side):
    for i in range(side):
        patch[i,j] = np.median(patch)

Eventually, I'll be using a 901x877 array from an image but I'm just trying to get a hold of this simple task first. How can I slide the array along and then down with a loop?

Answer 1

Here are a few "code smells" I see. Start with the range(side) since this number is set to 3 then you are going to have a result of [0,1,2]. Is that what you really want?

you set i,j = patch.size then immediately over write these values, in your for loops.

Finally, you're recalculating median every loop.

Ok, here's what I'd do.

1. figure out how many patches you'll need in both width and height. and multiply those by the size of the side.
2. slice your array (matrix) up into those pieces.
3. assign the patch to the median.

---

import numpy as np                                                                                                                                                                                         
w1 = np.arange(100).reshape(10,10)                                                                                                                                                                         
side = 3                                                                                                                                                                                                   
w, h = w1.shape                                                                                                                                                                                            
width_index   = np.array(range(w//side)) * side                                                                                                                                                             
height_index  = np.array(range(h//side)) * side                                                                                                                                                             

def assign_patch(patch, median, side):                                                                                                                                                                     
    """Break this loop out to prevent 4 nested 'for' loops"""                                                                                                                                              
    for j in range(side):                                                                                                                                                                                  
        for i in range(side):                                                                                                                                                                              
            patch[i,j] = median                                                                                                                                                                            
    return patch                                                                                                                                                                                           

for width in width_index:                                                                                                                                                                                  
    for height in height_index:                                                                                                                                                                            
        patch  = w1[width:width+side, height:height+side]                                                                                                                                                  
        median = np.median(patch)                                                                                                                                                                          
        assign_patch(patch, median, side)                                                                           

print w1        

Answer 2

You can use scikit-image's view_as_blocks and NumPy broadcasting to vectorize the operation:

import numpy as np
import skimage

w1 = np.arange(144).reshape(12,12)
print(w1)
# [[  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  81  82  83]
#  [ 84  85  86  87  88  89  90  91  92  93  94  95]
#  [ 96  97  98  99 100 101 102 103 104 105 106 107]
#  [108 109 110 111 112 113 114 115 116 117 118 119]
#  [120 121 122 123 124 125 126 127 128 129 130 131]
#  [132 133 134 135 136 137 138 139 140 141 142 143]]

side = 3
w2 = skimage.util.view_as_blocks(w1, (side, side))
w2[...] = np.median(w2, axis=(-2, -1))[:, :, None, None]
print(w1)
# [[ 13  13  13  16  16  16  19  19  19  22  22  22]
#  [ 13  13  13  16  16  16  19  19  19  22  22  22]
#  [ 13  13  13  16  16  16  19  19  19  22  22  22]
#  [ 49  49  49  52  52  52  55  55  55  58  58  58]
#  [ 49  49  49  52  52  52  55  55  55  58  58  58]
#  [ 49  49  49  52  52  52  55  55  55  58  58  58]
#  [ 85  85  85  88  88  88  91  91  91  94  94  94]
#  [ 85  85  85  88  88  88  91  91  91  94  94  94]
#  [ 85  85  85  88  88  88  91  91  91  94  94  94]
#  [121 121 121 124 124 124 127 127 127 130 130 130]
#  [121 121 121 124 124 124 127 127 127 130 130 130]
#  [121 121 121 124 124 124 127 127 127 130 130 130]]

Note that I had to change the size of your array to 12x12 so that all of your tiles of 3x3 actually fit in there.