在numpy数组中快速找到局部最大梯度值的方法？

Question

我有一个二维数组，我想找出所有局部最大值的索引。也就是说，给定一个索引 (i, j)，它的最大梯度是指它周围8个邻居值中，变化最大的绝对值：

Index: (i, j)

Neighbors:
(i-1,j+1)  (i,j+1)  (i+1,j+1)
(i-1,j)    [index]    (i+1,j)
(i-1,j-1)  (i,j-1)  (i+1,j-1)

Neighbor angles:
315           0            45
270        [index]         90
225          180          135

MaxGradient(i,j) = Max(|Val(index) - Val(neighbor)|)

如果这个索引的最大梯度至少和它周围邻居的最大梯度一样大，那么我们就说这个索引是 局部最大值。

这个算法的输出应该是一个二维数组，里面包含元组，或者是一个三维数组。对于原始数组中的每个索引，输出数组会包含一个值，表示这个索引是否是局部最大值，如果是的话，还会给出梯度的角度。

我最开始的实现方法是遍历数组两次，第一次计算最大梯度（存储在一个临时数组中），然后第二次遍历这个临时数组来确定局部最大值的索引。每次我都是通过循环，逐个查看每个索引。

有没有更高效的方法可以在numpy中实现这个呢？

Answer 1

考虑这8个相对索引：

X1 X2 X3
X4 X  X5
X6 X7 X8

你可以对每个像素X计算差值，比如D1=Val(X)-Val(X1)，D2=Val(X)-Val(X2)，D3=Val(X)-Val(X3)，D4=Val(X)-Val(X4)。其实你不需要计算其他的差值，因为它们是前四个的镜像。为了计算这些差值，你可以在图像的周围加上一行和一列的零，然后进行减法运算。

Answer 2

正如Cyborg所提到的，完成你的计算只需要计算四个差值（注意，如果这是在均匀网格上进行空间梯度计算，实际上对角线和反对角线的计算应该有一个1/sqrt(2)的因子）。如果我理解你的问题没错，使用numpy的实现可以像这样：

A=np.random.random(100).reshape(10,10)

# Padded copy of A
B=np.empty((12,12))
B[1:-1,1:-1]=A
B[0,1:-1]=A[0,:]   
B[-1,1:-1]=A[-1,:]
B[1:-1,0]=A[:,0]
B[1:-1,-1]=A[:,-1]
B[0,0]=A[1,1]
B[-1,-1]=A[-1,-1]
B[-1,0]=A[-1,0]
B[0,1]=A[0,1]

# Compute 4 absolute differences
D1=np.abs(B[1:,1:-1]-B[:-1,1:-1]) # first dimension
D2=np.abs(B[1:-1,1:]-B[1:-1,:-1]) # second dimension
D3=np.abs(B[1:,1:]-B[:-1,:-1]) # Diagonal
D4=np.abs(B[1:,:-1]-B[:-1,1:]) # Antidiagonal

# Compute maxima in each direction
M1=np.maximum(D1[1:,:],D1[:-1,:])
M2=np.maximum(D2[:,1:],D2[:,:-1])
M3=np.maximum(D3[1:,1:],D3[:-1,:-1])
M4=np.maximum(D4[1:,:-1],D4[:-1,1:])

# Compute local maximum for each entry
M=np.max(np.dstack([M1,M2,M3,M4]),axis=2)

这样你就能得到输入A在M中四个方向上的最大差值。类似的思路可以用来标记局部最大值，最终可以得到类似这样的结果：

T=np.where((M==np.max(np.dstack([Ma,Mb,Mc,Md,Me,Mf,Mg,Mh]),axis=2)))

这将给你一个数组，里面包含了M中局部最大值的坐标。