Unexpected behaviour of footprint in skimage.feature.peak_local_max

[ElpadoCan]

Description:

I'm not sure this is unexpected because maybe I don't fully understand how footprint is used in peak_local_max but I thought that only one peak can be detected within the same footprint, or at least that is what I understood.

However, it doesn't seem to work that way. I attached below a minimal example with synthetic data. The distance between the two peaks is set at 8 pixels and the footprint is a circle with diameter 20. The two peaks are detected even if they clearly lie within the same footprint.

What does footprint do? Thank you very much!

PS: I know I can use the min_distance parameter, but I'm running peak detection in 3D and my volume is anisotropic

Way to reproduce:

import numpy as np

import skimage.draw
import skimage.feature

import matplotlib.pyplot as plt

def gauss_2D_peak(xx, yy, xc, yc, sigma):
    xx = xx - xc
    yy = yy - yc
    gauss_x = np.exp(-(xx**2)/(2*(sigma**2)))
    gauss_y = np.exp(-(yy**2)/(2*(sigma**2)))
    gauss_yx = gauss_x*gauss_y
    return gauss_yx

Y = 128
X = 128

yy, xx = np.ogrid[0:Y, 0:X]

sigma = 3
yc, xc = 64, 64


dist = 8
xc1 = round(xc - dist/2)
xc2 = round(xc + dist/2)

img = (
    gauss_2D_peak(xx, yy, xc1, yc, sigma)
    + gauss_2D_peak(xx, yy, xc2, yc, sigma)
)

min_dist = 20
radius_footprint = min_dist/2
wh = int(np.ceil(radius_footprint))
y, x = np.ogrid[-wh:wh+1, -wh:wh+1]
footprint = (x**2 + y**2)/(radius_footprint**2) <= 1

peaks = skimage.feature.peak_local_max(img, footprint=footprint)

print(peaks)

fig, ax = plt.subplots(1, 3)

ax[0].imshow(img)
ax[0].plot(peaks[:, 1], peaks[:, 0], 'r.')

ax[1].plot(img.max(axis=0))
ax[2].imshow(footprint)

plt.show()

Version information:

3.10.10 (tags/v3.10.10:aad5f6a, Feb  7 2023, 17:20:36) [MSC v.1929 64 bit (AMD64)]
Windows-10-10.0.19045-SP0
scikit-image version: 0.22.0
numpy version: 1.26.1

[stefanv]

It sounds like what you need is for min_distance to accept a value along each axis.

The reason why the footprint picks up both peaks is because they are exactly the same height. If, instead, the one was slightly higher than the other, it would be rejected.

I'm not sure what the best way would be to fairly break this tie, otherwise.

Here's code for differing peak heights, that also plots the footprint circle on the input:

import numpy as np

import skimage as ski

import matplotlib.pyplot as plt

def gauss_2D_peak(xx, yy, xc, yc, sigma):
    xx = xx - xc
    yy = yy - yc
    gauss_x = np.exp(-(xx**2)/(2*(sigma**2)))
    gauss_y = np.exp(-(yy**2)/(2*(sigma**2)))
    gauss_yx = gauss_x*gauss_y
    return gauss_yx

Y = 128
X = 128

yy, xx = np.ogrid[0:Y, 0:X]

sigma = 3
yc, xc = 64, 64


dist = 8
xc1 = round(xc - dist/2)
xc2 = round(xc + dist/2)

img = (
    gauss_2D_peak(xx, yy, xc1, yc, sigma)
    + 1.1 * gauss_2D_peak(xx, yy, xc2, yc, sigma)
)

min_dist = 40
radius_footprint = min_dist/2
wh = int(np.ceil(radius_footprint))
y, x = np.ogrid[-wh:wh+1, -wh:wh+1]
footprint = (x**2 + y**2)/(radius_footprint**2) <= 1


peaks = ski.feature.peak_local_max(img, footprint=footprint)

print(peaks)

fig, ax = plt.subplots(1, 3)

ax[0].imshow(img)
ax[0].plot(peaks[:, 1], peaks[:, 0], 'r.')

ax[1].plot(img.max(axis=0))
ax[2].imshow(footprint)

for coord in peaks[:, ::-1]:
    circle = plt.Circle(coord, radius_footprint, edgecolor='red', fill=False)
    ax[0].add_patch(circle)

plt.show()

[ElpadoCan]

Thank you very much for the speedy reply!

I did not know that, very interesting. I don't know how feasible it is, but should an arbitrary peak be returned like with min_distance parameter?

It sounds like what you need is for min_distance to accept a value along each axis.

That is not possible at the moment right?

[stefanv]

No, unfortunately not. But feel free to file a feature request for either of the above! (Or, even better, if you have a moment to glance at the code to see how hard it would be to add: https://github.com/scikit-image/scikit-image/blob/main/skimage/feature/peak.py#L115).

[ElpadoCan]

Perfect, I was already looking into the code to check how this could be implemented. I think allowing one distance per axis would be cool. Thank you very much!

[stefanv]

Great, thanks @ElpadoCan! Please file an issue for that, whether or not you intend to work on it yourself.

[rfezzani]

I think that what you observe is the expected behavior!

As the doc describes, footprint defines the local region within which to search for peaks at every point in image, and min_distance the minimal allowed distance separating peaks.

min_distance is the main parameter to control the number of peaks and in the code snippet you provided, it is set to its default value 1.

I modified the example you provided to test different values of min_distance:

import numpy as np

import skimage.draw
import skimage.feature

import matplotlib.pyplot as plt


def gauss_2D_peak(xx, yy, xc, yc, sigma):
    xx = xx - xc
    yy = yy - yc
    gauss_x = np.exp(-(xx**2)/(2*(sigma**2)))
    gauss_y = np.exp(-(yy**2)/(2*(sigma**2)))
    gauss_yx = gauss_x*gauss_y
    return gauss_yx


Y = 128
X = 128

yy, xx = np.ogrid[0:Y, 0:X]

sigma = 3
yc, xc = 64, 64


dist = 8
xc1 = round(xc - dist/2)
xc2 = round(xc + dist/2)

img = (
    gauss_2D_peak(xx, yy, xc1, yc, sigma)
    + gauss_2D_peak(xx, yy, xc2, yc, sigma)
)

min_dist = 20
radius_footprint = min_dist/2
wh = int(np.ceil(radius_footprint))
y, x = np.ogrid[-wh:wh+1, -wh:wh+1]
footprint = (x**2 + y**2)/(radius_footprint**2) <= 1

fig, ax_list = plt.subplots(1, 6, figsize=(12, 3))

for d, ax in zip([1, 5, 10, 20], ax_list):
    peaks = skimage.feature.peak_local_max(img, min_distance=d,
                                           footprint=footprint)

    print(f"min_distance = {d :2d}: {len(peaks)} peak(s) detected")

    ax.imshow(img)
    ax.plot(peaks[:, 1], peaks[:, 0], 'r.')

    ax.set_title(f"min_distance={d}")
    ax.set_axis_off()


ax_list[-2].plot(img.max(axis=0))
ax_list[-2].set_title("profile")
ax_list[-1].imshow(footprint)
ax_list[-1].set_title("footprint")
ax_list[-1].set_axis_off()

fig.tight_layout()

plt.show()

To obtain

min_distance = 1: 2 peak(s) detected
min_distance = 5: 2 peak(s) detected
min_distance = 10: 1 peak(s) detected
min_distance = 20: 1 peak(s) detected

[ElpadoCan]

Hi @rfezzani,

Thanks for looking into this! And yes, I agree, I know that min_distance is the right parameter to use, but my images are anisotropic and I would ideally need one distance per axis of the input image. I opened #7318 to discuss this implementation.

Regarding the footprint instead, I admit that I don't fully understand how it works. If you see this tutorial on watershed, the footprint is set to a square with 3 pixels side and I think they do it to avoid detecting multiple peaks. That means that footprint should have similar effects to min_distance, right?

You can also see it in @stefanv comment that if you have peaks with different heights, then we get only one peak with a large enough footprint (despite min_distance=1) .

My expectation was that only one peak can be detected within each footprint, but that does not seem to be the case. If you have more insights, please share them with me because I'm very curious, thanks!

Cheers,
Francesco

[rfezzani]

Concerning the anisotropy of your data, I think that the approach proposed by @jni #7318 (comment) is the simplest way to implement this feature.

In the watershed tutorial, the definition of footprint is useless as it is set to its default value (please see the code).

In the solution proposed by @stefanv, if you scale up one of the two peaks, it becomes the local maximum within the region defined by footprint around the other peak, and thus becomes the only detected peak. But I don't know how this can be implemented with real world data :/.

[ElpadoCan]

Concerning the anisotropy of your data, I think that the approach proposed by @jni #7318 (comment) is the simplest way to implement this feature.

I agree, it is quite elegant actually, but does it ensure that if multiple peaks are found within min_distance it returns only the one with the largest height?

In the watershed tutorial, the definition of footprint is useless as it is set to its default value (please see the code).

Oh yeah, you are right, thanks for clarifying this!

In the solution proposed by @stefanv, if you scale up one of the two peaks, it becomes the local maximum within the region defined by footprint around the other peak, and thus becomes the only detected peak. But I don't know how this can be implemented with real world data :/.

Ok, then we agree that the problem is when the peaks have the same height :D. In this case, with min_distance, peak_local_max returns only one peak, so I would expect that footprint does the same, no?

[rfezzani]

I agree, it is quite elegant actually, but does it ensure that if multiple peaks are found within min_distance it returns only the one with the largest height?

Yes, hopefully 🙂 ...

Ok, then we agree that the problem is when the peaks have the same height :D. In this case, with min_distance, peak_local_max returns only one peak, so I would expect that footprint does the same, no?

footprint only defines the region where to search for local minima around each point, and min_distance will "control" the number of detected peaks: one may want to search for local maxima within a disk of radius 5 around each point, but impose a minimum distance of 20 between each peak...