0

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.

Community
  • 1
  • 1
FaCoffee
  • 7,609
  • 28
  • 99
  • 174

1 Answers1

2

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 [https://stackoverflow.com/a/19141711/774971 ]:

i, j = where(a > (1.0 - 5*np.finfo(a.dtype).eps)*the_peak)
Community
  • 1
  • 1
dmitry_romanov
  • 5,146
  • 1
  • 33
  • 36
  • Well, I don't know. There might be multiple peaks. What do you suggest in that case? – FaCoffee Aug 05 '16 at 15:10
  • Also, my question is not about the value of my max, since I normalize the array. It is about its coordinates in the grid. So your last line, I don't need it, as my max is always 1.0. – FaCoffee Aug 05 '16 at 15:13
  • 2
    I did multiple peaks in new version. My point about `max` is that you might get `0.99999999999` instead of `1.00000000000` and you will miss your peak completely. It happens in floating point arithmetic, especially when you do heavy math or modeling [http://floating-point-gui.de/errors/comparison ]. – dmitry_romanov Aug 05 '16 at 15:18