CF84 CF84 - 4 months ago 35
Python Question

Python: correct location of peaks in a 2D numpy array?

I know this question falls within the category of peak detection and there are answers available, but I think my problem is pretty simplistic and refers to a proof of principle.

Say I generate multiple,

nxn
2D numpy arrays of float values like this, which refer to a regular grid of
nxn
points (discrete domain):

myarray=
array([[ 0.82760819, 0.82113999, 0.811576 , 0.80308705, 0.81231903,
0.82296263, 0.78448964, 0.79308308, 0.82160627, 0.83475755,
0.8580934 , 0.8857617 , 0.89901092, 0.92479025, 0.91840606,
0.91029942, 0.88523943, 0.85798491, 0.84190422, 0.83350783,
0.83520675],
[ 0.84971526, 0.84759644, 0.81429419, 0.79936736, 0.81750327,
0.81874686, 0.80666801, 0.82297348, 0.84980788, 0.85698662,
0.87819988, 0.89572185, 0.89009723, 0.90347858, 0.89703473,
0.90092666, 0.88362073, 0.86711197, 0.84791422, 0.83632138,
0.83685225], ...] #you can generate any nxn array of random values


Now let's normalize them:

peak_value=myarray.max()
norm_array=myarray/peak_value


And proceed to locate the
(x,y)
coordinates of the peak:

from collections import Counter
longit=[] #x values of the peak
latit=[] #y values of the peak
for x in range(myarray.shape[0]):
for y in range(myarray.shape[1]):
if norm_array[x][y]==1:
longit.append(x)
latit.append(y)

x=numpy.array(longit)
y=numpy.array(latit)

c=zip(x,y)
temp=Counter(elem for elem in c) #Counts the number of peaks in each (x,y) point in the 11x11 grid
d=dict(Counter(temp)) #Has the shape d={(x,y): number of peaks}


Now this is just a single realization of the 2D numpy array. Given multiple arrays, the question is:

Is this the correct way to find the (x,y) of the peaks? Is there a more efficient way? Consider that there might be multiple peaks.

Answer

In C/C++ it is considered dangerous to do == with floating point numbers.

If there are no multiple exactly same pikes, you can use numpy.argmax:

a = random.rand(13, 11)
idx_unrolled = argmax(a)
the_peak = a[idx_unrolled/11, idx_unrolled%11]

If you need all pikes, you can get list of i, j indices using numpy.where:

i, j = where(a > 0.99*the_peak)

Use required number of 9 to tune the margin. For single precision floating points it is not than close to 1.

The best way may be something like [http://stackoverflow.com/a/19141711/774971 ]:

i, j = where(a > (1.0 - 5*np.finfo(a.dtype).eps)*the_peak)