Peak Detection in a 2D Array

Peak detection in a 2D array

I detected the peaks using a local maximum filter. Here is the result on your first dataset of 4 paws:

I also ran it on the second dataset of 9 paws and it worked as well.

Here is how you do it:

import numpy as np
from scipy.ndimage.filters import maximum_filter
from scipy.ndimage.morphology import generate_binary_structure, binary_erosion
import matplotlib.pyplot as pp

#for some reason I had to reshape. Numpy ignored the shape header.
paws_data = np.loadtxt("paws.txt").reshape(4,11,14)

#getting a list of images
paws = [p.squeeze() for p in np.vsplit(paws_data,4)]

def detect_peaks(image):
    """
    Takes an image and detect the peaks usingthe local maximum filter.
    Returns a boolean mask of the peaks (i.e. 1 when
    the pixel's value is the neighborhood maximum, 0 otherwise)
    """

    # define an 8-connected neighborhood
    neighborhood = generate_binary_structure(2,2)

    #apply the local maximum filter; all pixel of maximal value 
    #in their neighborhood are set to 1
    local_max = maximum_filter(image, footprint=neighborhood)==image
    #local_max is a mask that contains the peaks we are 
    #looking for, but also the background.
    #In order to isolate the peaks we must remove the background from the mask.

    #we create the mask of the background
    background = (image==0)

    #a little technicality: we must erode the background in order to 
    #successfully subtract it form local_max, otherwise a line will 
    #appear along the background border (artifact of the local maximum filter)
    eroded_background = binary_erosion(background, structure=neighborhood, border_value=1)

    #we obtain the final mask, containing only peaks, 
    #by removing the background from the local_max mask (xor operation)
    detected_peaks = local_max ^ eroded_background

    return detected_peaks

#applying the detection and plotting results
for i, paw in enumerate(paws):
    detected_peaks = detect_peaks(paw)
    pp.subplot(4,2,(2*i+1))
    pp.imshow(paw)
    pp.subplot(4,2,(2*i+2) )
    pp.imshow(detected_peaks)

pp.show()

All you need to do after is use scipy.ndimage.measurements.label on the mask to label all distinct objects. Then you'll be able to play with them individually.

Note that the method works well because the background is not noisy. If it were, you would detect a bunch of other unwanted peaks in the background. Another important factor is the size of the neighborhood. You will need to adjust it if the peak size changes (the should remain roughly proportional).

Peak detection in a noisy 2d array

I'm adding this answer because it's the solution I ended up using. It's a combination of Bi Rico's comment here (May 30 at 18:54) and the answer given in this question: Find peak of 2d histogram.

As it turns out using the peak detection algorithm from this question Peak detection in a 2D array only complicates matters. After applying the Gaussian filter to the image all that needs to be done is to ask for the maximum bin (as Bi Rico pointed out) and then obtain the maximum in coordinates.

So instead of using the detect-peaks function as I did above, I simply add the following code after the Gaussian 2D histogram is obtained:

# Get 2D histogram.
H, xedges, yedges = np.histogram2d(x, y, range=rang, bins=binsxy)
# Get Gaussian filtered 2D histogram.
H1 = gaussian_filter(H, 2, mode='nearest')
# Get center of maximum in bin coordinates.
x_cent_bin, y_cent_bin = np.unravel_index(H1.argmax(), H1.shape)
# Get center in x,y coordinates.
x_cent_coor , y_cent_coord = np.average(xedges[x_cent_bin:x_cent_bin + 2]), np.average(yedges[y_cent_g:y_cent_g + 2])

Peak finding algorithm in 2d-array with complexity O(n)

The second algorithm given in the linked PDF is O(n). A "window" is defined to collectively be the boundary (i.e. all four outer edges), middle column and middle row of the current sub-square. Here's a summary of the algorithm:

1. Find maximum value in current window
2. Return it if it's a peak
3. Otherwise, find the larger neighbor of this maximum and recurse in the corresponding quadrant.

As described in the slides, the time complexity is defined by T(n) = T(n/2) + cn (the T(n/2) term is due to the edge length being halved on each recursive step; the cn term is the time required to find the maximum in the current window). Hence, the complexity is O(n).

The correctness of this algorithm is based on several observations that are listed on one of the slides:

If you enter a quadrant, it contains a peak of the overall array

This is basically a generalization of the same 1D argument. You only enter a quadrant when it contains some element greater than all elements on the border. So, either that element will be a peak, or you can keep "climbing up" until you find a peak somewhere in the quadrant.

Maximum element of window never decreases as we descend in recursion

The next window in the recursion always contains the maximum element of the current window, so this is true.

Peak in visited quadrant is also peak in overall array

This follows from the definition of peak, since it only depends on immediate neighbors.

2D peak finding algorithm in O(n) worst case time?

1. Let's assume that width of the array is bigger than height, otherwise we will split in another direction.
2. Split the array into three parts: central column, left side and right side.
3. Go through the central column and two neighbour columns and look for maximum.
   
   • If it's in the central column - this is our peak
   • If it's in the left side, run this algorithm on subarray left_side + central_column
   • If it's in the right side, run this algorithm on subarray right_side + central_column

Why this works:

For cases where the maximum element is in the central column - obvious. If it's not, we can step from that maximum to increasing elements and will definitely not cross the central row, so a peak will definitely exist in the corresponding half.

Why this is O(n):

step #3 takes less than or equal to max_dimension iterations and max_dimension at least halves on every two algorithm steps. This gives n+n/2+n/4+... which is O(n). Important detail: we split by the maximum direction. For square arrays this means that split directions will be alternating. This is a difference from the last attempt in the PDF you linked to.

A note: I'm not sure if it exactly matches the algorithm in the code you gave, it may or may not be a different approach.

---

Related Topics