python - 如何获得 3-D 阵列中局部最大值周围的区域？

Question

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

循环局部最大值，从最大值到最小值。

如果局部最大值已经分配给一个区域，continue循环下一个局部最大值。

如果局部最大值未分配给任何区域，则将其分配给新区域。令该局部最大值的值为 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 ，这对于首先找到局部最大值非常有用。但我不知道如何获得它们周围的区域。我还研究了标准聚类算法和图像处理技术，例如 Blob 检测或局部阈值处理，但似乎都没有重现这个问题。

任何建议表示赞赏。

提前致谢，

编辑:多亏了 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)