在 3d numpy 数组中查找区域的中心坐标

Question

我有一个大型 numpy 3d 数组 (10000, 3, 3)。我想在其中找到每个区域的中心坐标（具有相同编号的集群）。每个子阵列可以有 1、2、3 或 4 个区域。

我的数组的一个子集是：

largearray= array([[[1, 0, 0],
    [0, 0, 2],
    [3, 0, 2]],

   [[0, 0, 4],
    [0, 0, 4],
    [0, 0, 4]],

   [[5, 0, 0],
    [5, 0, 6],
    [0, 6, 6]],

   [[7, 0, 8],
    [0, 0, 0],
    [9, 0,10]]])

我想要的输出是子数组的位置以及代表中心的 x 和 y 坐标：

#output:
array([[ 0., 0., 0.],
[ 0., 1.5, 2.],
[ 0., 2., 0.],
[ 1., 1.,  2.],
[ 2., 0.5,  0.],
[ 2., 1.66666667, 1.66666667],
[ 3., 0., 0.],
[ 3., 0., 2.],
[ 3., 2., 0.],
[ 3., 2., 2.]])

我对其他输出持开放态度，但这样的事情会很棒！

提前致谢！

Answer 1

你在找这个吗？

n_clusters = 10
for i in range(1, n_clusters + 1):
    matches = np.transpose((largearray == i).nonzero())
    print "The center of cluster {} is at {}".format(i, np.mean(matches, axis=0))

The center of cluster 1 is at [ 0.  0.  0.]
The center of cluster 2 is at [ 0.   1.5  2. ]
The center of cluster 3 is at [ 0.  2.  0.]
The center of cluster 4 is at [ 1.  1.  2.]
The center of cluster 5 is at [ 2.   0.5  0. ]
The center of cluster 6 is at [ 2.          1.66666667  1.66666667]
The center of cluster 7 is at [ 3.  0.  0.]
The center of cluster 8 is at [ 3.  0.  2.]
The center of cluster 9 is at [ 3.  2.  0.]
The center of cluster 10 is at [ 3.  2.  2.]

Answer 2

使用numpy_indexed 包中的功能（免责声明：我是它的作者），可以构建一个完全矢量化的解决方案（即没有for循环）：

import numpy_indexed as npi
idx = np.indices(largearray.shape).reshape(largearray.ndim, largearray.size)
label, mean = npi.group_by(largearray, axis=None).mean(idx, axis=1)

对于大量输入，这应该更有效。

请注意，如果标签在每个子数组中不是唯一的（它们似乎在您的示例中，但未明确说明），但您仍想仅取每个子数组的平均值，您可以简单地写成这样：

(label, subarr), mean = npi.group_by((largearray.flatten(), idx[0])).mean(idx[1:], axis=1)

也就是说，通过子数组索引和标签的唯一元组进行分组。

Answer 3

这是一个完全矢量化的版本，仅使用 numpy：

# the list of all the cluster ids
clusters = np.arange(1, n_clusters+1)

# convert to a boolean array, where mask[i] = largearray != clusters[i]
mask = np.rollaxis(clusters != largearray[...,np.newaxis], axis=-1)

# the coordinate of each item in the array
idx = np.indices(largearray.shape)

# broadcast (cluster_num, 1, ...) with (1, coord, ...)
mask, idx = np.broadcast_arrays(mask[:,np.newaxis], idx[np.newaxis,:])

# an array of the indices, with all the ones we don't care about masked out
idx_mask = np.ma.masked_array(idx, mask)

# flatten out the unneeded dimensions and average over them
means = idx_mask.reshape(idx_mask.shape[:2] + (-1,)).mean(axis=-1)

给予：

masked_array(data =
 [[0.0 0.0 0.0]
 [0.0 1.5 2.0]
 [0.0 2.0 0.0]
 [1.0 1.0 2.0]
 [2.0 0.5 0.0]
 [2.0 1.6666666666666667 1.6666666666666667]
 [3.0 0.0 0.0]
 [3.0 0.0 2.0]
 [3.0 2.0 0.0]
 [3.0 2.0 2.0]],
             mask =
 [[False False False]
 [False False False]
 [False False False]
 [False False False]
 [False False False]
 [False False False]
 [False False False]
 [False False False]
 [False False False]
 [False False False]],
       fill_value = 1e+20)

请注意，这也会通过设置掩码来指示哪些集群不存在

Answer 4

在这样棘手的问题上，难点在于保证最少的操作次数。掩蔽或分组技术通常会增加无用的操作。

简单的python方法大致是：

def center(largearray):
    n=largearray.max()+1
    (x,y,z)=largearray.shape
    sm=np.zeros((n,3),np.float64)
    cnt=np.zeros(n,np.float64)
    for i in range(x):
        for j in range(y):
            for k in range(z):
                l=largearray[i,j,k]
                sm[l,0] += i; sm[l,1] += j; sm[l,2] += k
                cnt[l]+=1
    for l in range(n):
        if cnt[l]>0:
            for m in range(3): sm[l,m] /= cnt[l]  
        else:
            for m in range(3): sm[l,m] = -1
    return sm [1:]

这里是给定数组的（还不错的）性能：

In [16]: %timeit center(largearray)
1000 loops, best of 3: 248 µs per loop

幸运的是，使用numba 可以大大加速此类代码：

In [17]: center2=numba.jit(center)

In [18]: %timeit center2(largearray)
100000 loops, best of 3: 3.29 µs per loop

Answer 5

您可能还想查看处理与此相关的问题的numpy-groupies 包。 [免责声明：我是合著者]。它应该比numpy-indexed（另一个答案中提到的包）更快，因为它使用bincount而不是argsort和reduceat。

但是，您在这里的任务很简单，您可以直接使用bincount：

s0, s1, s2 = a.shape

group_counts = np.bincount(a.ravel())

idx = np.broadcast_to(np.arange(s0).reshape([s0, 1, 1]), [s0,s1,s2])
group_sum_0 = np.bincount(a.ravel(), idx.ravel()) 

idx = np.broadcast_to(np.arange(s1).reshape([1, s1, 1]), [s0,s1,s2])
group_sum_1 = np.bincount(a.ravel(), idx.ravel()) 

idx = np.broadcast_to(np.arange(s2).reshape([1, 1, s2]), [s0,s1,s2])
group_sum_2 = np.bincount(a.ravel(), idx.ravel()) 

group_mean = np.vstack((group_sum_0, group_sum_1, group_sum_2)) / group_counts

group_mean.T[1:] # this is the output you show in the question

或者如果你想“作弊”，你可以使用来自 scipy 的ndimage.measurements 中的一个函数。