pythonimagenumpyscipy

Get coordinates of local maxima in 2D array above certain value


from PIL import Image
import numpy as np
from scipy.ndimage.filters import maximum_filter
import pylab

# the picture (256 * 256 pixels) contains bright spots of which I wanna get positions
# problem: data has high background around value 900 - 1000

im = Image.open('slice0000.png')
data = np.array(im)

# as far as I understand, data == maximum_filter gives True-value for pixels
# being the brightest in their neighborhood (here 10 * 10 pixels)

maxima = (data == maximum_filter(data,10))
# How can I get only maxima, outstanding the background a certain value, let's say 500 ?

I'm afraid I don't really understand the scipy.ndimage.filters.maximum_filter() function. Is there a way to obtain pixel-coordinates only within the spots and not within the background?

https://i.sstatic.net/RImHW.png (16-bit grayscale picture, 256*256 pixels)


Solution

  • import numpy as np
    import scipy
    import scipy.ndimage as ndimage
    import matplotlib.pyplot as plt
    
    fname = '/tmp/slice0000.png'
    neighborhood_size = 5
    threshold = 1500
    
    data = scipy.misc.imread(fname)
    
    data_max = ndimage.maximum_filter(data, neighborhood_size) # apply maximum filter with a size of neighborhood_size
    maxima = (data == data_max) # boolean mask: local maximum within neighborhood_size
    data_min = ndimage.minimum_filter(data, neighborhood_size) # apply minimum filter with a size of neighborhood_size
    diff = ((data_max - data_min) > threshold) # boolean mask where the difference of the filters exceeds threshold
    maxima[~diff] = False # remove the local maxima which do not satisfy the minimum difference in neighborhood
    
    labeled, num_objects = ndimage.label(maxima) # label connected components on maxima binary array (boolean mask)
    slices = ndimage.find_objects(labeled) # slices are 2d rect
    x, y = [], []
    for dy,dx in slices:
        x_center = (dx.start + dx.stop - 1)/2
        x.append(x_center)
        y_center = (dy.start + dy.stop - 1)/2    
        y.append(y_center)
    
    plt.imshow(data)
    plt.savefig('/tmp/data.png', bbox_inches = 'tight')
    
    plt.autoscale(False)
    plt.plot(x,y, 'ro')
    plt.savefig('/tmp/result.png', bbox_inches = 'tight')
    

    Given data.png:

    enter image description here

    the above program yields result.png with threshold = 1500. Lower the threshold to pick up more local maxima:

    enter image description here

    References: