Python: Find highest value in neighborhood in a 3D array (stack of images)

Question

I use a template matching algorithm to detect instances of a given object in a target image. The algorithm produces a heatmap looking like this:

So far I have used skimage.feature.peak_local_max to find the local peaks in the heatmap which worked well.

Now I want to run the template matching for different rotations of the template which will give me several heatmaps. Like these:

I have stacked these images on top of each other to create a 3D array and now I want to find the highest values in the local neighborhood. The neighborhood is for example defined as 15 pixels wide and high (but includes all rotations).

The 3D array is of shape (height, width, rotation). The z-axis corresponds to the rotation of the template.

EDIT: I have found a solution myself (with the helpful comments from markuscosinus) and will post this here in case anyone has the same problem:

I looked at the implementation of scipy.ndimage.filters.maximum_filter (which can be found here. I twitched the code slightly to work with my data shape and used footprint = np.ones((size_neighborhood, size_neighborhood, nr_rotations)).

Answer 1

One way you could achieve this is the following:

1. Find local maxima for each rotation using skimage.feature.peak_local_max
2. For every maximum m at (x,y,z) found in step 1 check if there were found others in the neighbourhood. Easiest case would be to check the neighbourhood cube with corners (x-n,y-n,y-n) and (x+n,y+n,y+n), where n is the neighbourhood size. If there are some in the neighbourhood, only keep the largest.

I know, this is not very pleasing in terms of efficiency, but in case you do not have too many or too large pictures it should get the job done :)

Answer 2

I have found a solution myself (with the helpful comments from markuscosinus) and will post this here in case anyone has the same problem:

I looked at the implementation of scipy.ndimage.filters.maximum_filter (which can be found here. I twitched the code slightly to work with my data shape and used footprint = np.ones((size_neighborhood, size_neighborhood, nr_rotations)).