Search code examples
pythonimage-processingscikit-image

skimage peak local max finds multiple spots in close proximity due to image impurities

I have an image that looks like this, with some larger impurities/overexposed spots. It generally doesn't matter if they're detected, as measurements are time resolved, so they'll be removed later on.

However, I'm interested in as many as the small dots as possible - as fast as possible. skimage.feature.peak_local_max does a really good job, and is very easy to use on different data, because there's no need to play around much with intensity scaling.

enter image description here

The problem is though, that large spots for some reason give very strong positives.

import skimage.io
import skimage.feature
import skimage.morphology
from matplotlib.collections import PatchCollection
import matplotlib.pyplot as plt


def plotRoi(spots, img_ax, color, radius):
    patches = []
    for spot in spots:
        y, x = spot
        c = plt.Circle((x, y), radius)
        patches.append(c)
    img_ax.add_collection(PatchCollection(patches, facecolors = "None", edgecolors = color, alpha = 0.3, linewidths = 1))


img = skimage.io.imread("/Path/to/img.png")
img = img[:,:,0]

fig, ax = plt.subplots()

ax.imshow(img, cmap = "Greys")
spots = skimage.feature.peak_local_max(img, min_distance = 0, exclude_border = True, num_peaks = 2000)
plotRoi(spots, ax, "red", radius = 10)

plt.show()

enter image description here

And searching for thousands of spots in some images lead to a large number of local maxima being pretty much on top of each other. Is there a way to avoid this, e.g. by applying a filter on image loading, as I would prefer not to move to a slower type of peak fitting?

Solution

  • The problem is that there are peaks in regions of pixels with exactly the same value. One way to solve this is to merge the peaks to the centre of mass these regions.

    Below I recreate the problem and solve it as stated.

    import numpy as np
import matplotlib.pyplot as plt
from skimage.feature import peak_local_max
from scipy.ndimage.measurements import center_of_mass, label

# Generate test data with two peaks, one of which consists of two pixels of equal value
image = np.zeros((11,11),dtype=np.uint8)
image[5,3] = 128
image[5,2] = 255
image[5,7:9] = 255
image[6,8] = 128

# Finding peaks normally; results in three peaks
peaks = peak_local_max(image)

# Find peaks and merge equal regions; results in two peaks
is_peak = peak_local_max(image, indices=False) # outputs bool image
labels = label(is_peak)[0]
merged_peaks = center_of_mass(is_peak, labels, range(1, np.max(labels)+1))
merged_peaks = np.array(merged_peaks)

# Visualize the results
fig,(ax1,ax2)=plt.subplots(1,2)
ax1.imshow(image.T,cmap='gray')
ax1.plot(peaks[:,0],peaks[:,1],'ro')

ax2.imshow(image.T,cmap='gray')
ax2.plot(merged_peaks[:,0],merged_peaks[:,1],'ro')

    Result:

    merged_peaks