通过使用 NumPy/SciPy 检测向量的局部最大值来提取直方图模式

Posted 2023-03-12

【中文标题】通过使用 NumPy/SciPy 检测向量的局部最大值来提取直方图模式

【英文标题】：Extract histogram modes by detecting the local maxima of a vector with NumPy/SciPy

【发布时间】：2017-08-14 13:04:30

【问题描述】：

NumPy/SciPy` 有没有办法在提取局部最大值时只保留直方图模式（在下图中显示为蓝点）？：

这些最大值是使用 scipy.signal.argrelmax 提取的，但我只需要获取两个 modes 值并忽略检测到的其余最大值：

# calculate dB positive image
img_db = 10 * np.log10(img)
img_db_pos = img_db + abs(np.min(img_db))
data = img_db_pos.flatten() + 1

# data histogram
n, bins = np.histogram(data, 100, normed=True)

# trim data
x = np.linspace(np.min(data), np.max(data), num=100)

# find index of minimum between two modes
ind_max = argrelmax(n)
x_max = x[ind_max]
y_max = n[ind_max]

# plot
plt.hist(data, bins=100, normed=True, color='y')
plt.scatter(x_max, y_max, color='b')
plt.show()

注意：

我已经设法使用这个Smoothing filter 来获得与直方图匹配的曲线（但我没有曲线方程）。

【问题讨论】：

区分最大值是什么意思？

@kazemakase 我的意思是只保留这两种模式并摆脱检测到的其他最大值。我可能应该编辑我的问题。

【参考方案1】：

这可以通过适当调整函数argrelmax的参数order来实现，即通过调整每边有多少点被认为是检测局部最大值。

我已使用此代码创建模拟数据（您可以使用不同变量的值来确定它们对生成的直方图的影响）：

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import argrelextrema, argrelmax

m = 50
s = 10
samples = 50000
peak_start = 30
peak_width = 10
peak_gain = 4

np.random.seed(3)
data = np.random.normal(loc=m, scale=s, size=samples)
bell, edges = np.histogram(data, bins=np.arange(2*(m + 1) - .5), normed=True)
x = np.int_(.5*(edges[:-1] + edges[1:]))

bell[peak_start + np.arange(peak_width)] *= np.linspace(1, peak_gain, peak_width)

plt.bar(x, bell)
plt.show()

如下图，仔细选择order的值很重要。事实上，如果order 太小，您可能会检测到嘈杂的局部最大值，而如果order 太大，您可能无法检测到某些模式。

In [185]: argrelmax(bell, order=1)
Out[185]: (array([ 3,  5,  7, 12, 14, 39, 47, 51, 86, 90], dtype=int64),)

In [186]: argrelmax(bell, order=2)
Out[186]: (array([14, 39, 47, 51, 90], dtype=int64),)

In [187]: argrelmax(bell, order=3)
Out[187]: (array([39, 47, 51], dtype=int64),)

In [188]: argrelmax(bell, order=4)
Out[188]: (array([39, 51], dtype=int64),)

In [189]: argrelmax(bell, order=5)
Out[189]: (array([39, 51], dtype=int64),)   

In [190]: argrelmax(bell, order=11)
Out[190]: (array([39, 51], dtype=int64),)

In [191]: argrelmax(bell, order=12)
Out[191]: (array([39], dtype=int64),)

这些结果很大程度上取决于直方图的形状（如果您只更改用于生成数据的参数之一，order 的有效值范围可能会有所不同）。为了使模式检测更加稳健，我建议您将平滑直方图传递给argrelmax，而不是原始直方图。

【讨论】：

order=4 为我的直方图做了它（如上面的问题所示）。

【参考方案2】：

我猜，你想在y_max 中找到第二大数字。希望这个例子对你有所帮助：

np.random.seed(4)  # for reproducibility
data = np.zeros(0)
for i in xrange(10):
    data = np.hstack(( data, np.random.normal(i, 0.25, 100*i) ))

# data histogram
n, bins = np.histogram(data, 100, normed=True)

# trim data
x = np.linspace(np.min(data), np.max(data), num=100)

# find index of minimum between two modes
ind_max = argrelmax(n)
x_max = x[ind_max]
y_max = n[ind_max]

# find first and second max values in y_max
index_first_max = np.argmax(y_max)
maximum_y = y_max[index_first_max]
second_max_y = max(n for n in y_max if n!=maximum_y)
index_second_max = np.where(y_max == second_max_y)

# plot
plt.hist(data, bins=100, normed=True, color='y')
plt.scatter(x_max, y_max, color='b')
plt.scatter(x_max[index_first_max], y_max[index_first_max], color='r')
plt.scatter(x_max[index_second_max], y_max[index_second_max], color='g')
plt.show()