python - 找到二维直方图的峰值

Question

我制作了一些(x, y)数据的二维直方图，我得到了这样的图像：

我想要一种方法来获取(x, y)存储最大值的点的坐标H。例如，在上图的情况下，它将是具有近似坐标的两个点：(1090, 1040)和(1110, 1090)。

这是我的代码：

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from os import getcwd
from os.path import join, realpath, dirname

# Path to dir where this code exists.
mypath = realpath(join(getcwd(), dirname(__file__)))
myfile = 'datafile.dat'

x, y = np.loadtxt(join(mypath,myfile), usecols=(1, 2), unpack=True)

fig = plt.figure()
ax = fig.add_subplot(111)

xmin, xmax = min(x), max(x)
ymin, ymax = min(y), max(y)

rang = [[xmin, xmax], [ymin, ymax]]

binsxy = [int((xmax - xmin) / 20), int((ymax - ymin) / 20)]

H, xedges, yedges = np.histogram2d(x, y, range=rang, bins=binsxy)

extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
cp = ax.imshow(H.transpose()[::-1], interpolation='nearest', extent=extent, cmap=cm.jet)
fig.colorbar(cp)

plt.show()

---

编辑

我已经尝试了 Marek 和 qarma 发布的解决方案，试图获取箱的坐标而不是它们的索引，如下所示：

# Marek's answer
x_cent, y_cent = unravel_index(H.argmax(), H.shape)
print('Marek')
print(x_cent, y_cent)
print(xedges[x_cent], yedges[y_cent])

# qarma's answer
idx = list(H.flatten()).index(H.max())
x_cent2, y_cent2 = idx / H.shape[1], idx % H.shape[1]
local_maxs = np.argwhere(H == H.max())
print('\nqarma')
print(x_cent2, y_cent2)
print(xedges[x_cent2], yedges[y_cent2])
print(xedges[local_maxs[0,0]], yedges[local_maxs[0,1]], xedges[local_maxs[1,0]], yedges[local_maxs[1,1]])

这导致：

Marek
(53, 50)
(1072.7838144329899, 1005.0837113402063)

qarma
(53, 50)
(1072.7838144329899, 1005.0837113402063)
(1072.7838144329899, 1005.0837113402063, 1092.8257731958763, 1065.3611340206187)

所以最大坐标是相同的，这很好！现在我有一个小问题，因为如果我放大二维图，我会看到全局最大值和局部最大值的坐标都有点偏离中心：

为什么是这样？

Answer 1

这个问题应该对您有所帮助：Python: get the position of the maximum item in a numpy array

您可以使用H.max()获取最大值，然后将其与H并用于numpy.nonzero查找所有最大值的位置： numpy.nonzero(H.max() == H). 这将比仅仅更昂贵，H.argmax()但您将获得所有最大值。

Answer 2

这是找到第一个全局最大值的方法

idx = list(H.flatten()).index(H.max())
x, y = idx / H.shape[1], idx % H.shape[1]

寻找所有最大值的坐标留给读者作为练习......

numpy.argwhere(H == H.max())

编辑

你的代码：

H, xedges, yedges = np.histogram2d(x, y, range=rang, bins=binsxy)

这里H包含直方图箱的直方图值和xedges, yedges边界。请注意，edges数组的大小比H相应维度的大小大一。因此：

for x, y in numpy.argwhere(H == H.max()):
    # center is between x and x+1
    print numpy.average(xedges[x:x + 2]), numpy.average(yedges[y:y + 2])

Answer 3

该库findpeaks可以使用。

pip install findpeaks

我没有看到您的数据，但让我尝试另一个类似的示例：

from findpeaks import findpeaks

# initialize with default parameters. The "denoise" parameter can be of use in your case
fp = findpeaks()
# import 2D example dataset
img = fp.import_example()
# make the fit
fp.fit(img)
# Make plot
fp.plot()

持久性可用于确定峰值的影响。您会看到点 1、2 和 3 显示最强的峰，然后是其余峰。

# Compute persistence
fp.plot_persistence()