如何在三维数组中获取局部极大值周围的区域？

Question

我有一个大的三维numpy数组（1024 x 1024 x 1024），需要找到局部极大值周围的区域，使得所有值大于局部最大值50％的相邻点都被递归地聚集在同一区域内。可以使用迭代算法来解决：

• 从最高值到最低值循环遍历局部最大值。
• 如果该局部最大值已经分配给某个区域，则跳过并继续循环遍历下一个局部最大值。
• 如果该局部最大值未分配给任何区域，则将其分配给新区域。将该局部最大值的值设为M。
• 添加所有值>0.5*M的局部最大值的相邻点到该区域中。
• 递归添加所有值>0.5*M的相邻点、相邻点的相邻点......到该区域中。
• 重复执行，直到所有局部最大值都被分配给某个区域。

对于如此巨大的数组，这种算法效率很低，因此我正在寻找一些使用Python库的向量化解决方案。

具体而言，对于像这样的二维数组：

0.0  1.6  0.0  0.0  0.0
0.0  2.0  1.0  0.0  5.0
1.6  3.0  1.0  0.0  4.6
0.0  0.0  0.0  9.0  4.6

将会有两个这样的区域：

                    5.0
                    4.6
               9.0  4.6

并且

     1.6               
     2.0               
1.6  3.0               

直观地说，我正在寻找局部极大值周围的“山脉”，其中“山脉”由轮廓级别定义，该级别不是绝对的而是相对于局部极大值。

我已经尝试使用scipy.ndimage来查找局部最大值，这非常有用。但是我不知道如何获取它们周围的区域。我还查看了标准聚类算法和图像处理技术，例如斑点检测或本地阈值处理，但似乎没有一个能够重现这个问题。

感谢您的建议。

提前致谢，

编辑：感谢taw，下面是一个解决方案。

import math
import numpy as np

def soln(data, maxind):
    test = np.pad(data,(1,1),'constant',constant_values = [-math.inf,-math.inf])
    regionlist = {} # Dictionary with the results
    for region, ind in enumerate(maxind):  # Loop over all maxima
        M = test[ind[0]+1, ind[1]+1, ind[2]+1] # Value of the maximum

        if M == -np.inf:
            continue

        regionlist[region] = set()
        regionlist[region].add(ind)

        test[ind[0]+1, ind[1]+1, ind[2]+1] = -math.inf # All points that are added to the results are set to -infinity
        neighbors = set()
        neighbors.add((ind[0]+1, ind[1]+1, ind[2]+1))
        while len(neighbors)>0: #create region iteratively
            newneighbors = set() # This will contain the new elements in the region
            for i in neighbors: 
                values = test[i[0]-1:i[0]+2, i[1]-1:i[1]+2, i[2]-1:i[2]+2] # Values of neighbours
                valuesmask = values > .5*M # Neighbours that fall in region
                list1 = range(i[0]-2, i[0]+1)
                list2 = range(i[1]-2, i[1]+1)
                list3 = range(i[2]-2, i[2]+1)
                indlist = list(itertools.product(list1, list2, list3)) # 3-D list with coordinates of neighbours
                for count,j in enumerate(valuesmask): 
                    if j:
                       newneighbors.add(indlist[count])

            #update iteration
            newneighbors = newneighbors.difference(neighbors) # Remove all elements that were already iterated over and added to regionlist
            regionlist[region].update(newneighbors) # Add the new elements in the region to regionlist
            neighbors = set((x[0]-1, x[1]-1, x[2]-1) for x in newneighbors) # In the next iteration, iterate only over new elements in the region
            for i in newneighbors:
                test[i[0]+1, i[1]+1, i[2]+1] = -math.inf #set values to -inf after added to region

    return regionlist

- IvFrat

Answer 1

我不确定“局部最大值”是如何定义的，或者您正在使用scipy.ndimage中的哪些函数来获取它们。这里有一个函数，将给出属于每个区域的索引集（返回索引而不是值）。成本看起来像O（将分配给区域的点数）。常数取决于数组的维度。我认为在复杂度方面不可能做得比这更好。此解决方案还适用于2d数组。

import math
import numpy as np
test = np.array([[0, 1.6, 0, 0, 0,], [0, 2, 1,0,5],[1.6,3,1,0,4.6],[0,0,0,9,4.6]])

maxind = [(3,3),(2,1)] #indices of maxima
def soln(data, maxind):
    test = np.pad(data,(1,1),'constant',constant_values = [-math.inf,-math.inf])
    regionlist = {}
    for region,ind in enumerate(maxind):  #all maxima
        regionlist[region] = set()
        regionlist[region].add(ind)
        M = test[ind[0]+1,ind[1]+1]
        test[ind[0]+1,ind[1]+1] = -math.inf
        neighbors = set()
        neighbors.add((ind[0]+1,ind[1]+1))
        while len(neighbors)>0: #create region iteratively
            newneighbors = set()
            for i in neighbors: 
                values = test[i[0]-1:i[0]+2,i[1]-1:i[1]+2]
                valuesmask = values.flatten() > .5*M
                list1 = np.repeat(list(range(i[0]-2,i[0]+1)),3)
                list2 = np.tile(list(range(i[1]-2,i[1]+1)), 3)
                indlist = list(zip(list1,list2))
                for count,j in enumerate(valuesmask): 
                    if j:
                        newneighbors.add(indlist[count])

            #update iteration
            newneighbors = newneighbors.difference(neighbors)
            regionlist[region].update(newneighbors)
            neighbors = newneighbors
            for i in newneighbors:
                test[i[0]+1,i[1]+1] = -math.inf #set values to -inf after added to region

    return regionlist

regionlist = soln(test, maxind)

- Taw

---

谢谢！这个解决方案对我的需求来说相对较快。我正在使用scipy.ndimage.maximum_filter获取局部最大值（大于所有最近邻点的点），就像这个问题中所示。我已经更新了问题，加入了正确的三维扩展和一个小错误修复：我认为，在你的例子中，neighbors包含了data数组中点的坐标，而根据循环它应该包含test数组中的坐标。 - IvFrat