31
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?

http://i.stack.imgur.com/RImHW.png (16-bit grayscale picture, 256*256 pixels)

Chris
  • 42,054
  • 16
  • 135
  • 151
feinmann
  • 932
  • 1
  • 11
  • 20

3 Answers3

61
import numpy as np
import scipy
import scipy.ndimage as ndimage
import scipy.ndimage.filters as filters
import matplotlib.pyplot as plt

fname = '/tmp/slice0000.png'
neighborhood_size = 5
threshold = 1500

data = scipy.misc.imread(fname)

data_max = filters.maximum_filter(data, neighborhood_size)
maxima = (data == data_max)
data_min = filters.minimum_filter(data, neighborhood_size)
diff = ((data_max - data_min) > threshold)
maxima[diff == 0] = 0

labeled, num_objects = ndimage.label(maxima)
slices = ndimage.find_objects(labeled)
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:

Community
  • 1
  • 1
unutbu
  • 777,569
  • 165
  • 1,697
  • 1,613
  • 2
    hello unutbu, I'm afraid I don't really get your solution, meaning the output. at the moment I managed to kick out all maxima that have absolute value less than let's say 1500. I'm just trying if the outcome is satisfying. – feinmann Feb 02 '12 at 14:32
  • Most likely it is I who does not understand your question. Are you looking for a way to find the `(x,y)` coordinates of the maxima? If so, you can find them using `np.where(maxima)`. – unutbu Feb 02 '12 at 14:43
  • 1
    you're right. but I want to get rid of the local maxima being in the background. like saying: a local maximum is only a local maximum if it stands out from its neighborhood more than a certain value. At the moment I cancel the background by setting all pixels to zero that have a value below 1500, but I am not really satisfied with this. Do you know ImageJ? The 'Find Maxima' function does a pretty good job and I'd like to reproduce this output. To be clear: I want to have the coordinates of the brightest pixels within the bright spots on the picture. – feinmann Feb 03 '12 at 08:47
13
import numpy as np
import scipy
import scipy.ndimage as ndimage
import scipy.ndimage.filters as filters
import matplotlib.pyplot as plt

fname = '/tmp/slice0000.png'
neighborhood_size = 5
threshold = 1500

data = scipy.misc.imread(fname)

data_max = filters.maximum_filter(data, neighborhood_size)
maxima = (data == data_max)
data_min = filters.minimum_filter(data, neighborhood_size)
diff = ((data_max - data_min) > threshold)
maxima[diff == 0] = 0

labeled, num_objects = ndimage.label(maxima)
xy = np.array(ndimage.center_of_mass(data, labeled, range(1, num_objects+1)))

plt.imshow(data)
plt.savefig('/tmp/data.png', bbox_inches = 'tight')

plt.autoscale(False)
plt.plot(xy[:, 1], xy[:, 0], 'ro')
plt.savefig('/tmp/result.png', bbox_inches = 'tight')

The previous entry was super useful to me, but the for loop slowed my application down. I found that ndimage.center_of_mass() does a great and fast job to get the coordinates... hence this suggestion.

unutbu
  • 777,569
  • 165
  • 1,697
  • 1,613
legazier
  • 131
  • 1
  • 2
10

This can now be done with skimage.

from skimage.feature import peak_local_max
xy = peak_local_max(data, min_distance=2,threshold_abs=1500)

On my computer, for a VGA image size it runs about 4x faster than the above solution and also returned a more accurate position in certain cases.

Eyal S.
  • 888
  • 3
  • 14
  • 26
  • Here is a nice example/tutorial in `skimage` for finding 'blobs' in an image: https://scikit-image.org/docs/stable/auto_examples/features_detection/plot_blob.html#sphx-glr-auto-examples-features-detection-plot-blob-py – Stefano Mar 16 '20 at 22:59