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

我有一个大型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

python大神给出的解决方案

我不确定如何定义“局部最大值”，或者您正在使用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)

Python sqlite3数据库已锁定 - python

我在Windows上使用Python 3和sqlite3。我正在开发一个使用数据库存储联系人的小型应用程序。我注意到，如果应用程序被强制关闭（通过错误或通过任务管理器结束），则会收到sqlite3错误（sqlite3.OperationalError：数据库已锁定）。我想这是因为在应用程序关闭之前，我没有正确关闭数据库连接。我已经试过了： connectio…

Python pytz时区函数返回的时区为9分钟 - python

由于某些原因，我无法从以下代码中找出原因：>>> from pytz import timezone >>> timezone('America/Chicago') 我得到：<DstTzInfo 'America/Chicago' LMT-1 day, 18:09:00 STD…

用大写字母拆分字符串，但忽略AAA Python Regex - python

我的正则表达式：vendor = "MyNameIsJoe. I'mWorkerInAAAinc." ven = re.split(r'(?<=[a-z])[A-Z]|[A-Z](?=[a-z])', vendor) 以大写字母分割字符串，例如：'我的名字是乔。 I'mWorkerInAAAinc”变成…

如何打印浮点数的全精度[Python] - python

我编写了以下函数，其中传递了x，y的值：def check(x, y): print(type(x)) print(type(y)) print(x) print(y) if x == y: print "Yes" 现在当我打电话check(1.00000000000000001, 1.0000000000000002)它正在打印：<…

Python:检查新文件是否在文件夹中[重复] - python

This question already has answers here: How do I watch a file for changes? （23个答案） 3年前关闭。 我是python的新手，但是我尝试创建一个自动化过程，其中我的代码将侦听目录中的新文件条目。例如，某人可以手动将zip文件复制到一个特定的文件夹中，并且我希望我的代码能够在文件完全…