2

I have written a function (given below) to find positions of top num_peaks intensity values in an image along with performing a non-maximum suppression to select only local maxima values:

def find_peaks(img, num_peaks, threshold_fraction=0.5, nhood_size=None): """Find locally maximum intensities in an image""" # calculate threshold as a fraction of intensity range in the image threshold = (threshold_fraction * (img.max() - img.min())) + img.min() # determine offsets in each direction that constitute a neighbourhood if nhood_size is None: nhood_size = np.array(img.shape) * 0.02 nhood_offset = (np.around(nhood_size / 2)).astype(int) # create array with labelled fields to allow intensity-specific sorting rows, cols = np.array(np.where(img >= threshold)) values = [] for i, j in zip(rows, cols): values.append((i, j, img[i, j])) dtype = [('row', int), ('col', int), ('intensity', np.float64)] indices = np.array(values, dtype=dtype) # sort in-place in descending order indices[::-1].sort(order='intensity') # Perform suppression for idx_set in indices: intensity = idx_set[2] if not intensity: continue x0 = idx_set[1] - nhood_offset[1] xend = idx_set[1] + nhood_offset[1] y0 = idx_set[0] - nhood_offset[0] yend = idx_set[0] + nhood_offset[0] indices_to_suppress = np.where((indices['col'] >= x0) & (indices['col'] <= xend) & (indices['row'] >= y0) & (indices['row'] <= yend)) if indices_to_suppress: indices['intensity'][indices_to_suppress] = 0 idx_set[2] = intensity # perform second sorting & return the remaining n (num_peaks) most # intense values indices[::-1].sort(order='intensity') if len(indices) <= num_peaks: return np.array([np.array([i, j]) for i, j in zip(indices['row'], indices['col'])]) # or return all of them return np.array([np.array([i, j]) for i, j in zip(indices['row'][:num_peaks], indices['col'][:num_peaks])])

This seems work properly for small images and a large threshold_fraction (less values to suppress around), but has proven to be quite inefficient for my purposes, where I have lower thresholds like 0.1 to 0.2. I have not been able to make this more efficient with my beginner numpy skills.

I would like to know if I could make any changes to this piece of code that might improve its performance. Also, since I'm using numpy and OpenCV, it would be nice to know if there was a library function(s) that could achieve something similar or utilize it somehow to write an efficient peak finder.

Improve this question

2 Answers 2

Reset to default
1

Your algorithm relies on sorting the image values, so its time complexity is O(N log N) in the number N of image pixels, and requires an extra O(N) memory buffer to store the sorted values.

A standard local-maxima detector relies on dilating the image with a small structuring element (the size of which depends on how "noisy" or "peaky" the image is), and then finding the pixels where the the original image and the dilated one have equal values. This is O(N) in time and requires an extra O(N) buffer. Suppression of non local-maxima pixel is automatic, and removal of "too close" local maxima can be easily done by sorting the found local maxima (whose number is normally << N) by image intensity and removing the close ones in sorted order.

So the latter algorithm is theoretically more performant, and likely practically more performant too, for large enough images. Of course YMMV depending on the implementation and image size.

Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

1

For what it's worth, since you're already using numpy, the additional use of skimage (installed with pip install scikit-image) should be a breeze. You can then use skimage.feature.peak_local_max() or skimage.ndimage.maximum_filter() to create a map of the locations you're interested in, and since the objects are natively considered to be numpy arrays you can suppress those regions by subtracting some fraction of the resultant array from your original.

Improve this answer

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.