0

I'm moving from IDL to Python but there is a routine I require which I can not find in python.

I want to detect the peak(s) within a fairly smooth 2D array which has a shape equal to (60,300). I only expect to find one or two peaks (but never more than 4).

Is there a python routine that is similar to IDL's FindPeaks (link below) that I could use?

http://hesperia.gsfc.nasa.gov/ssw/smei/ucsd/gen/idl/toolbox/findpeaks.pro

Other answers on stackoverflow concerning peak detection in 2D arrays/images seem to be designed to return multiple (e.g., >20) peaks (links below). Can any of these be modified for smoothly varying arrays?

see: Peak detection in a 2D array, Find peak of 2d histogram, Peak detection in a noisy 2d array

(in fact, what I actually want is to detect minimas but there seems to be more questions online regarding maxima/peaks. When I used FindPeaks in IDL I simply inverted my data.)

shclim
  • 115
  • 2
  • 7

2 Answers2

0
numPeaks = 2
my_points = [[9, 30], [11, 61], [29, 48], [30, 97], [44, 33], [64, 76], [83, 45], [84, 22], [90, 28], [93, 81]]
sorted(my_points,reverse=True,key=lambda (x,y):y)[:numPeaks]

might do it (without having looked at the other function)

there is no function like the one you linked to.

this method above explores the entire space and sorts it by Y height giving you various maxima(note this is slightly different than a peak). it may not be an acceptable solution for large data sets ... to find actual peaks (also know as local maxima) an easy algorithm is random hill climbing this will return all the peaks (that is a point that is higher than both its neighboring points)

you can then sort only those peaks by their Y component

Joran Beasley
  • 110,522
  • 12
  • 160
  • 179
  • Thanks Joran. However, I would like the routine to return one peak when (by eye) there is only one peak in the 2D array, and two when there is clearly two etc. i.e., I do not want to specify "numPeaks" – shclim Jan 27 '15 at 17:02
  • 2
    is there some threshold for "by eye" ? eg cant be more than 0.1 away? or cant be more than 1% of total range? or? are these very large arrays where full exploration is not possible? or? – Joran Beasley Jan 27 '15 at 17:04
  • When using FindPeaks in IDL I set flat=1.0 – shclim Jan 27 '15 at 17:10
  • can you think of a way to use this thing that I have provided to return a variable number of peaks depending on some threshold? this may not be an acceptable solution if your problem space is very large ... in which case you need to find the maxima differently (thats what a peak is) – Joran Beasley Jan 27 '15 at 17:12
  • I have just posted an answer to my original question. However, I'm sure there are ways to improve on it. To answer Joran Beasley's comment above, a threshold equal to "np.mean(image)-np.std(image)" seems to do the trick. – shclim Feb 12 '15 at 17:04
0

After some trial and error, and adapting code from Peak detection in a 2D array, I have come up with the following solution. NB: threshold = np.mean(image)-np.std(image)

Please feel free to post any improvements, I'm still kinda new to python

def detect_minima(image, threshold):

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

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

    # we create the mask of the background
    background = (image > threshold)

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

    # we obtain the final mask, containing only minimas, 
    # by removing the background from the local_min mask
    detected_minimas = local_min - eroded_background

    # Label minimas and extract each minima
    labels, nlabels = label(detected_minimas , structure=neighborhood)
    dim1 = np.zeros(nlabels, dtype=np.int)
    dim2 = np.zeros(nlabels, dtype=np.int)
    val  = np.zeros(nlabels)
    for l in range(0,nlabels):
        dim1[l], dim2[l] = minimum_position(image, labels=labels, index=l+1)
        val[l] = image[dim1[l],dim2[l]]

    # ensure points are below threshold
    lt_threshold = val < threshold 
    dim1 = dim1[lt_threshold]
    dim2 = dim2[lt_threshold]
    val  = val[lt_threshold]

    # remove points which neighbour any background or missing points.
    image[background] = threshold
    keep_asl = []

    for l in range(0,len(dim1)):

        bounds1a, bounds1b = dim1[l]-1, dim1[l]+2
        bounds2a, bounds2b = dim2[l]-1, dim2[l]+2   

        if (bounds1a < 0):                  bounds1a = 0
        if (bounds1b >= image.shape[0]):    bounds1a = image.shape[0]-1
        if (bounds2a < 0):                  bounds2a = 0
        if (bounds2b >= image.shape[1]):    bounds2a = image.shape[1]-1

        neighbour_vals = image[ bounds1a:bounds1b, bounds2a:bounds2b ]
        neighbour_max = np.max(neighbour_vals)

        if (type(image.data) == np.ma.core.MaskedArray):
            if ( (neighbour_max < threshold) & (image.fill_value not in neighbour_vals) ):
                    keep_asl.append(l)
        else:
            if (neighbour_max < threshold): keep_asl.append(l)      

    dim1 = dim1[keep_asl]
    dim2 = dim2[keep_asl]
    val  = val[keep_asl]

    # sort by value (lowest first)
    sort_ind = np.argsort(val)
    dim1 = dim1[sort_ind]
    dim2 = dim2[sort_ind]
    val  = val[sort_ind]

    return dim1, dim2, val
shclim
  • 115
  • 2
  • 7