Local maxima detection#
For detecting local maxima, pixels surrounded by pixels with lower intensity, we can use some functions in scikit-image and clesperanto.
See also
from skimage.feature import peak_local_max
import pyclesperanto_prototype as cle
from skimage.io import imread, imshow
from skimage.filters import gaussian
import matplotlib.pyplot as plt
We start by loading an image and cropping a region for demonstration purposes. We used image set BBBC007v1 image set version 1 (Jones et al., Proc. ICCV Workshop on Computer Vision for Biomedical Image Applications, 2005), available from the Broad Bioimage Benchmark Collection [Ljosa et al., Nature Methods, 2012].
image = imread("../../data/BBBC007_batch/A9 p7d.tif")[-100:, 0:100]
cle.imshow(image)
Preprocessing#
A common preprocessing step before detecting maxima is blurring the image. This makes sense to avoid detecting maxima that are just intensity variations resulting from noise.
preprocessed = gaussian(image, sigma=2, preserve_range=True)
cle.imshow(preprocessed)
peak_local_max#
The peak_local_max function allows detecting maxima which have higher intensity than surrounding pixels and other maxima according to a defined threshold.
coordinates = peak_local_max(preprocessed, threshold_abs=5)
coordinates
array([[23, 85],
[11, 29],
[41, 40],
[88, 34],
[72, 83],
[69, 89],
[31, 72],
[75, 16],
[80, 22],
[ 6, 56]], dtype=int64)
These coordinates can be visualized using matplotlib’s plot function.
cle.imshow(preprocessed, continue_drawing=True)
plt.plot(coordinates[:, 1], coordinates[:, 0], 'r.')
[<matplotlib.lines.Line2D at 0x2309908fbb0>]
If there are too many maxima detected, one can modify the results by changing the sigma parameter of the Gaussian blur above or by changing the threshold passed to the peak_local_max function.
detect_maxima_box#
The function peak_local_max tends to take long time, e.g. when processing large 3D image data. Thus, an alternaive shall be introduced: clesperanto’s detect_maxima_box is an image filter that sets pixels to value 1 in case surrounding pixels have lower intensity. It typically performs fast also on large 3D image data.
local_maxima_image = cle.detect_maxima_box(preprocessed)
local_maxima_image
|
cle._ image
|
Obviously, it results in a binary image. This binary image can be converted to a label image by labeling individual spots with different numbers. From this label image, we can remove maxima detected at image borders, which might be useful in this case.
all_labeled_spots = cle.label_spots(local_maxima_image)
labeled_spots = cle.exclude_labels_on_edges(all_labeled_spots)
labeled_spots
|
cle._ image
|
To visualize these spots on the original image, it might make sense to increase the size of the spots - just for visualization purposes.
label_visualization = cle.dilate_labels(labeled_spots, radius=3)
cle.imshow(preprocessed, continue_drawing=True)
cle.imshow(label_visualization, labels=True, alpha=0.5)
In the lower center of this image we see now a local maximum that has been detected in the background. We can remove those maxima in lower intensity regions by thresholding.
binary_image = cle.threshold_otsu(preprocessed)
binary_image
|
cle._ image
|
We can now exclude labels from the spots image where the intensity in the binary image is not within range [1..1] (i.e. exactly 1).
final_spots = cle.exclude_labels_with_map_values_out_of_range(
binary_image,
labeled_spots,
minimum_value_range=1,
maximum_value_range=1
)
final_spots
|
cle._ image
|
We can then visualize the spots again using the strategy introduced above, but this time on the original image.
label_visualization2 = cle.dilate_labels(final_spots, radius=3)
cle.imshow(image, continue_drawing=True)
cle.imshow(label_visualization2, labels=True, alpha=0.5)
