Efficient way to split image into overlapping patches, do math on the patches, then merge the image back together by taking the max of at every pixel

Question

Basically trying to just do what the tital says and I have a working function, but it is running much to slow. The input image to this function will always have the same shape and dimension (2d).

Here's my working but slow code

def z_score_detect():
    data = currentData - bkgData
    tileY =12
    tileX =12
    h,w = data.shape

    DETECTED = np.zeros(data.shape).astype(np.float32)
    for curY in range(tileY//2, h-tileY//2,tileY//2):
        for curX in range(tileX//2, w-tileX//2,tileX//2):
            # Get one of the patches of data
            S = data[curY-tileY//2:curY+tileY//2,curX-tileX//2:curX+tileX//2]
            # comput stdev, median
            rms_bot_75 = rms(np.sort(S,axis=None)[:int((tileX*tileY)*.75)])
            std_dev = np.std(S)
            if std_dev != 0:
                # IF SIGMA NOT ZERO
                S = (S-rms_bot_75)/std_dev
                # Done with processing this patch so update our output image
            else:
                continue
            DETECTED[curY-tileY//2:curY+tileY//2,curX-tileX//2:curX+tileX//2] = np.maximum(S,DETECTED[curY-tileY//2:curY+tileY//2,curX-tileX//2:curX+tileX//2])
                        
    return DETECTED

Answer 1

You can use view_as_windows which uses strides under the hood to tile your image and calculate on the tile. Be aware that this creates a view and any attempt to do calculation that changes the content of it should be done with caution. Here, your if statement calculation is sequential and has to be dealt cautiously (or even sequentially). As far as std and median goes, you can use:

from skimage.util import view_as_windows
#this tiles your image into shape (tileX,tileY) overlapping tiles
tiles = view_as_windows(data, (tileX,tileY))
#std and median of each tile
std_dev = tiles.std(axis=(-2,-1))
median = np.median(tiles, axis=(-2,-1))