Getting the maximum in each rolling window of a 2D numpy array

Question

I have a 2D numpy array and I want to get the maximum value contained in each 2d rolling window that starts from left to right, top to bottom, rolling one row or column each time. The most naive method would be iterating through all rolling windows and get the maximum of all values enclosed in this rolling window. I wrote down this method below:

import numpy as np
shape=(1050,300)
window_size=(120,60)
a = np.arange(shape[1]*shape[0]).reshape(shape[1],shape[0])
max_Map=np.full((shape[1]-window_size[1]+1,shape[0]-window_size[0]+1),0,dtype='uint32')

for i in range(shape[1]-window_size[1]+1):
    for j in range(shape[0]-window_size[0]+1):
        window_max=np.max(a[i:i+window_size[1],j:j+window_size[0]])
        max_Map[i][j]=window_max

But this is terribly inefficient, as there are only 2 rows(or 2 column) changed between each sliding but my code doesn't take into account any correlations between 2 consecutive rolling windows. An improvement I could think of is for each sliding window(assuming rolling horizontally) I will calculate the maximum of the left most column and the maximum of the remaining columns and take the max of the 2 values as the current window maximum. And for the next rolling window the maximum will be max of the newly added column and the previous remaining columns...But I still don't think this is optimized...

I will really appreciate it if someone can point me to the right direction,I feel like this should be a well studied problem but I couldn't find solutions anywhere... Thanks in advance!

Answer 1

Approach #1 Using Scipy's 2D max filter -

from scipy.ndimage.filters import maximum_filter as maxf2D

# Store shapes of inputs
N,M = window_size
P,Q = a.shape

# Use 2D max filter and slice out elements not affected by boundary conditions
maxs = maxf2D(a, size=(M,N))
max_Map_Out = maxs[M//2:(M//2)+P-M+1, N//2:(N//2)+Q-N+1]

Approach #2 Using Scikit's 2D sliding window views -

from skimage.util.shape import view_as_windows

N,M = window_size
max_Map_Out = view_as_windows(a, (M,N)).max(axis=(-2,-1))

Note on window size and its use : The original approach has the window sizes aligned in a flipped manner, i.e. the first shape parameter of window_size slides along the second axis, while the second shape parameter decides how the window slides along the first axis. This might not be the case for other problems that do sliding max filtering, where we usually use the first shape parameter for the first axis of the 2D array and similarly for the second shape parameter. So, to solve for those cases, simply use : M,N = window_size and use the rest of the codes as they are.

Runtime test

Approaches -

def org_app(a, window_size):
    shape = a.shape[1], a.shape[0]
    max_Map=np.full((shape[1]-window_size[1]+1,
                     shape[0]-window_size[0]+1),0,dtype=a.dtype)
    for i in range(shape[1]-window_size[1]+1):
        for j in range(shape[0]-window_size[0]+1):
            window_max=np.max(a[i:i+window_size[1],j:j+window_size[0]])
            max_Map[i][j]=window_max
    return max_Map

def maxf2D_app(a, window_size):
    N,M = window_size
    P,Q = a.shape
    maxs = maxf2D(a, size=(M,N))
    return maxs[M//2:(M//2)+P-M+1, N//2:(N//2)+Q-N+1]

def view_window_app(a, window_size):
    N,M = window_size
    return view_as_windows(a, (M,N)).max(axis=(-2,-1))

Timings and verification -

In [573]: # Setup inputs
     ...: shape=(1050,300)
     ...: window_size=(120,60)
     ...: a = np.arange(shape[1]*shape[0]).reshape(shape[1],shape[0])
     ...: 

In [574]: np.allclose(org_app(a, window_size), maxf2D_app(a, window_size))
Out[574]: True

In [575]: np.allclose(org_app(a, window_size), view_window_app(a, window_size))
Out[575]: True

In [576]: %timeit org_app(a, window_size)
1 loops, best of 3: 2.11 s per loop

In [577]: %timeit view_window_app(a, window_size)
1 loops, best of 3: 1.14 s per loop

In [578]: %timeit maxf2D_app(a, window_size)
100 loops, best of 3: 3.09 ms per loop

In [579]: 2110/3.09  # Speedup using Scipy's 2D max filter over original approach
Out[579]: 682.8478964401295