Fast grid search for height "h" defining local maxima

I have a series (n = 45) of 3D volumes of sparse fluorophores, taken over time, for which I would like to extract the (x, y, z) positions of. This temporal data has varying dynamic range - some images have very low contrast and still others are very high.

The number of spots I expect to detect should be approx ~ 1000 such spots per time point, so it’s a matter of finding the right height “h” (for each of the 45 3-D volumes) in skimage’s skimage.extrema.h_maxima to call such local maximas. Given that the # of spots I expect to detect should be constant, is there a fast and robust grid search that allows me to find the best height “h” for each volume? That each of my imaging volumes has very different dynamic range - despite histogram matching - makes setting one standard height very difficult! Would appreciate any advice!

Hi @b2jia welcome to the forum! h_maxima is a great function but for your use case it’s probably very expensive to call it with several values of h and adjust to find the right number of maxima - would it be an option at all? You could also try to measure some relevant statistics of the image to estimate which h you need (for example based on percentiles of the image np.percentile).

Another idea could be to use a different function to find local maxima. Have you tried using local_maxima or peak_local_max instead of h_maxima? This would return a large number of local maxima which you could sort by decreasing mean intensity, or local contrast, etc… and use a threshold on the sorted list to have always the same number of extrema.

1 Like

Hi Emmanuelle,

Thanks for your response! skimage.morphology.peak_local_max is really nifty, thank you for the recommendation. Not sure why it is folded into morphology, while many other peak calling functions are part of skimage.exposure.

One thing I noticed is that with peak_local_max one can provide a series of masks, on top of which we can constrain the maximum # of peaks called in each mask area. It’s unclear to me how the parameter threshold_rel translates to each mask - does it set the threshold as threshold_rel * np.max(total_img) or threshold_rel * np.max(only_masked_area)? I went to look at the source code but it is a bit unclear.

Hum, probably you are using an “old” version of scikit-image because in the last stable release (0.17) peak_local_max is found in skimage.feature ( I looked at the function signature and I did not found a parameter for passing a mask, can you please explain? I agree that the parameter num_peaks for a maximum number of peaks should be really convenient in your case.

1 Like

By mask, I was referring to the labels parameter for peak_local_max, which also allows me to specify the number of peaks called per labeled area.

Thanks again for the referral to peak_local_max, I have since updated my version. It is too bad that we can’t estimate directly the peak statistics from peak_local_max. I ended up just using a different function!