如何获取 3D 数组中局部最大值周围的区域？

Question

我有一个大型3D numpy 数组（1024 x 1024 x 1024），我需要找到局部最大值周围的区域，以便所有值大于局部最大值的相邻点（例如，局部最大值的 50%）递归地聚集在同一地区。迭代算法是：

• 从最高值到最低值循环局部最大值。
• 如果局部最大值已分配给某个区域，continue则循环下一个局部最大值。
• 如果局部最大值没有分配给任何区域，则将其分配给新区域。令该局部最大值的值为 M。
• 将值 > 0.5*M 的局部最大值的所有邻居添加到区域中。
• 递归地将所有邻居、邻居的邻居......添加到该区域，其值 > 0.5*M。
• 重复直到所有局部最大值都分配给某个区域。

对于如此巨大的数组，这种算法是令人望而却步的，所以我一直在寻找一些带有 Python 库的矢量化解决方案。

更具体地说，对于 2D 数组，例如

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

Answer 1

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

该解决方案也适用于二维数组。

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)