2D peak finding with non-maximum suppression using numpy

Question

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.

Francesco Callari · Accepted Answer · 2018-01-19 17:16:14Z

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.

Hans Musgrave · Accepted Answer · 2018-01-04 16:27:17Z

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.

answered Jan 4, 2018 at 16:27

Hans Musgrave

7,2012 gold badges21 silver badges40 bronze badges

Collectives™ on Stack Overflow

2D peak finding with non-maximum suppression using numpy

2 Answers 2

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related