给定如下峰数组：
peaks = [0, 5, 0, 3, 2, 0, 1, 7, 0]
如何创建指示最新峰值的步骤数组，如下所示：
steps = [0, 5, 5, 3, 3, 3, 3, 7, 7]
要求：
- 这将用于大型 3D 图像 (1000**3( 的图像分析，因此需要快速，这意味着没有 for 循环或列表推导......只有 numpy 矢量化。
- 我上面给出的示例是一个线性列表，但这需要同样适用于ND图像。 这意味着沿单个轴进行操作，但允许同时使用多个轴。
注意
我最近问了一个问题，结果证明是一个骗局(用scipy.maximum.accumulate很容易解决(，但我的问题也包含一个可选的"如果会很好"的转折，如上所述。 事实证明，我实际上也需要第二种行为，所以我只重新发布这部分。
这是一个处理 ND 的解决方案，可以检测"宽峰"，如..., 0, 4, 4, 4, 3, ...但不能检测..., 0, 4, 4, 4, 7, ...。
import numpy as np
import operator as op
def keep_peaks(A, axis=-1):
B = np.swapaxes(A, axis, -1)
# take differences between consecutive elements along axis
# pad with -1 at the start and the end
# the most efficient way is to allocate first, because otherwise
# padding would involve reallocation and a copy
# note that in order to avoid that copy we use np.subtract and its
# out kwd
updown = np.empty((*B.shape[:-1], B.shape[-1]+1), B.dtype)
updown[..., 0], updown[..., -1] = -1, -1
np.subtract(B[..., 1:], B[..., :-1], out=updown[..., 1:-1])
# extract indices where the there is a change along axis
chnidx = np.where(updown)
# get the values of the changes
chng = updown[chnidx]
# find indices of indices 1) where we go up and 2) the next change is
# down (note how the padded -1's at the end are useful here)
# also include the beginning of each 1D subarray
pkidx, = np.where((chng[:-1] > 0) & (chng[1:] < 0) | (chnidx[-1][:-1] == 0))
# use indices of indices to retain only peak indices
pkidx = (*map(op.itemgetter(pkidx), chnidx),)
# construct array of changes of the result along axis
# these will be zero everywhere
out = np.zeros_like(A)
aux = out.swapaxes(axis, -1)
# except where there is a new peak
# at these positions we need to put the differences of peak levels
aux[(*map(op.itemgetter(slice(1, None)), pkidx),)] = np.diff(B[pkidx])
# we could ravel the array and do the cumsum on that, but raveling
# a potentially noncontiguous array is expensive
# instead we keep the shape, at the cost of having to replace the
# value at the beginning of each 2D subarray (we do not need the
# "line-jump" difference but the plain 1st value there)
aux[..., 0] = B[..., 0]
# finally, use cumsum to go from differences to plain values
return out.cumsum(axis=axis)
peaks = [0, 5, 0, 3, 2, 0, 1, 7, 0]
print(peaks)
print(keep_peaks(peaks))
# show off axis kwd and broad peak detection
peaks3d = np.kron(np.random.randint(0, 10, (3, 6, 3)), np.ones((1, 2, 1), int))
print(peaks3d.swapaxes(1, 2))
print(keep_peaks(peaks3d, 1).swapaxes(1, 2))
示例运行：
[0, 5, 0, 3, 2, 0, 1, 7, 0]
[0 5 5 3 3 3 3 7 7]
[[[5 5 3 3 1 1 4 4 9 9 7 7]
[2 2 9 9 3 3 4 4 3 3 7 7]
[9 9 0 0 2 2 5 5 7 7 9 9]]
[[1 1 3 3 9 9 3 3 7 7 0 0]
[1 1 1 1 4 4 5 5 0 0 3 3]
[5 5 5 5 8 8 1 1 2 2 7 7]]
[[6 6 3 3 8 8 2 2 3 3 2 2]
[6 6 9 9 3 3 9 9 3 3 9 9]
[1 1 5 5 7 7 2 2 7 7 1 1]]]
[[[5 5 5 5 5 5 5 5 9 9 9 9]
[2 2 9 9 9 9 4 4 4 4 7 7]
[9 9 9 9 9 9 9 9 9 9 9 9]]
[[1 1 1 1 9 9 9 9 7 7 7 7]
[1 1 1 1 1 1 5 5 5 5 3 3]
[5 5 5 5 8 8 8 8 8 8 7 7]]
[[6 6 6 6 8 8 8 8 3 3 3 3]
[6 6 9 9 9 9 9 9 9 9 9 9]
[1 1 1 1 7 7 7 7 7 7 7 7]]]