Fast peak-finding and centroiding in python
Question: I am trying to develop a fast algorithm in python for finding peaks in an
image and then finding the centroid of those peaks. I have written the
following code using the scipy.ndimage.label and ndimage.find_objects for
locating the objects. This seems to be the bottleneck in the code, and it
takes about 7 ms to locate 20 objects in a 500x500 image. I would like to
scale this up to larger (2000x2000) image, but then the time increases to
almost 100 ms. So, I'm wondering if there is a faster option.

Here is the code that I have so far, which works, but is slow. First I
simulate my data using some gaussian peaks. This part is slow, but in practice
I will be using real data, so I don't care too much about speeding that part
up. I would like to be able to find the peaks very quickly.

    
    
    import time
    import numpy as np
    import matplotlib.pyplot as plt
    import scipy.ndimage
    import matplotlib.patches 
    
    plt.figure(figsize=(10,10))
    ax1 = plt.subplot(221)
    ax2 = plt.subplot(222)
    ax3 = plt.subplot(223)
    ax4 = plt.subplot(224)
    
    size        = 500 #width and height of image in pixels
    peak_height = 100 # define the height of the peaks
    num_peaks   = 20
    noise_level = 50
    threshold   = 60
    
    np.random.seed(3)
    
    #set up a simple, blank image (Z)
    x = np.linspace(0,size,size)
    y = np.linspace(0,size,size)
    
    X,Y = np.meshgrid(x,y)
    Z = X*0
    
    #now add some peaks
    def gaussian(X,Y,xo,yo,amp=100,sigmax=4,sigmay=4):
        return amp*np.exp(-(X-xo)**2/(2*sigmax**2) - (Y-yo)**2/(2*sigmay**2))
    
    for xo,yo in size*np.random.rand(num_peaks,2):
        widthx = 5 + np.random.randn(1)
        widthy = 5 + np.random.randn(1)
        Z += gaussian(X,Y,xo,yo,amp=peak_height,sigmax=widthx,sigmay=widthy)
    
    #of course, add some noise:
    Z = Z + scipy.ndimage.gaussian_filter(0.5*noise_level*np.random.rand(size,size),sigma=5)    
    Z = Z + scipy.ndimage.gaussian_filter(0.5*noise_level*np.random.rand(size,size),sigma=1)    
    
    t = time.time() #Start timing the peak-finding algorithm
    
    #Set everything below the threshold to zero:
    Z_thresh = np.copy(Z)
    Z_thresh[Z_thresh<threshold] = 0
    print 'Time after thresholding: %.5f seconds'%(time.time()-t)
    
    #now find the objects
    labeled_image, number_of_objects = scipy.ndimage.label(Z_thresh)
    print 'Time after labeling: %.5f seconds'%(time.time()-t)
    
    peak_slices = scipy.ndimage.find_objects(labeled_image)
    print 'Time after finding objects: %.5f seconds'%(time.time()-t)
    
    def centroid(data):
        h,w = np.shape(data)   
        x = np.arange(0,w)
        y = np.arange(0,h)
    
        X,Y = np.meshgrid(x,y)
    
        cx = np.sum(X*data)/np.sum(data)
        cy = np.sum(Y*data)/np.sum(data)
    
        return cx,cy
    
    centroids = []
    
    for peak_slice in peak_slices:
        dy,dx  = peak_slice
        x,y = dx.start, dy.start
        cx,cy = centroid(Z_thresh[peak_slice])
        centroids.append((x+cx,y+cy))
    
    print 'Total time: %.5f seconds\n'%(time.time()-t)
    
    ###########################################
    #Now make the plots:
    for ax in (ax1,ax2,ax3,ax4): ax.clear()
    ax1.set_title('Original image')
    ax1.imshow(Z,origin='lower')
    
    ax2.set_title('Thresholded image')
    ax2.imshow(Z_thresh,origin='lower')
    
    ax3.set_title('Labeled image')
    ax3.imshow(labeled_image,origin='lower') #display the color-coded regions
    
    for peak_slice in peak_slices:  #Draw some rectangles around the objects
        dy,dx  = peak_slice
        xy     = (dx.start, dy.start)
        width  = (dx.stop - dx.start + 1)
        height = (dy.stop - dy.start + 1)
        rect = matplotlib.patches.Rectangle(xy,width,height,fc='none',ec='red')
        ax3.add_patch(rect,)
    
    ax4.set_title('Centroids on original image')
    ax4.imshow(Z,origin='lower')
    
    for x,y in centroids:
        ax4.plot(x,y,'kx',ms=10)
    
    ax4.set_xlim(0,size)
    ax4.set_ylim(0,size)
    
    plt.tight_layout
    plt.show()
    

The results for size=500: ![enter image description
here](http://i.stack.imgur.com/CNONo.png)

EDIT: If the number of peaks is large (~100) and the size of the image is
small, then the bottleneck is actually the centroiding part. So, perhaps the
speed of this part also needs to be optimized.


Answer: Your method for finding the peaks (simple thresholding) is of course very
sensitive to the choice of threshold: set it too low and you'll "detect"
things that are not peaks; set it too high and you'll miss valid peaks.

There are more robust alternatives, that will detect all the local maxima in
the image intensity regardless of their intensity value. My preferred one is
applying a dilation with a small (5x5 or 7x7) structuring element, then find
the pixels where the original image and its dilated version have the same
value. This works because, by definition, dilation(x, y, E, img) = { max of
img within E centered at pixel (x,y) }, and therefore dilation(x, y, E, img) =
img(x, y) whenever (x,y) is the location of a local maximum at the scale of E.

With a fast implementation of the morphological operators (e.g. the one in
OpenCV) this algorithm is linear in the size of the image in both space and
time (one extra image-sized buffer for the dilated image, and one pass on
both). In a pinch, it can also be implemented on-line without the extra buffer
and a little extra complexity, and it's still linear time.

To further robustify it in the presence of salt-and-pepper or similar noise,
which may introduce many false maxima, you can apply the method twice, with
structuring elements of different size (say, 5x5 and 7x7), then retain only
the stable maxima, where stability can be defined by unchanging position of
the maxima, or by position not changing by more than one pixel, etc.
Additionally, you may want to suppress low nearby maxima when you have reason
to believe they are due to noise. An efficient way to do this is to first
detect all the local maxima as above, sort them descending by height, then go
down the sorted list and keep them if their value in the image has not changed
and, if they are kept, set to zero all the pixels in a (2d+1) x (2d+1)
neighborhood of them, where d is the min distance between nearby maxima that
you are willing to tolerate.

