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

Question

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

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

# 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()

注意:

我已设法使用此平滑滤镜来获得与直方图匹配的曲线(但我没有曲线方程).

Answer 1

我猜，你想在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()

Answer 2

这可以通过适当地调整order函数的参数来实现argrelmax,即通过调整每侧上的多少点来考虑局部最大值.

我已经使用此代码创建模拟数据(您可以使用不同变量的值来计算它们对生成的直方图的影响):

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原始直方图.