有效地找到标记图像区域的质心答案

【问题标题】：Efficiently find centroid of labelled image regions有效地找到标记图像区域的质心
【发布时间】：2018-04-01 03:06:54
【问题描述】：

我有一个分割图像作为唯一标签 1 ... k 的二维矩阵。例如：

img = 
    [1 1 2 2 2 2 2 3 3]
    [1 1 1 2 2 2 2 3 3]
    [1 1 2 2 2 2 3 3 3]
    [1 4 4 4 2 2 2 2 3]
    [4 4 4 5 5 5 2 3 3]
    [4 4 4 5 5 6 6 6 6]
    [4 4 5 5 5 5 6 6 6]

我正在尝试确定区域质心。也就是说，每个标签，质心的 X、Y 坐标是多少？例如，标签 1 的质心为 (1.25, 0.625)。只需平均行号 ((0 + 0 + 1 + 1 + 1 + 2 + 2 + 3) / 8 = 1.25) 和列号 ((0 + 0 + 0 + 0 + 1 + 1 + 1 + 2) / 8 = 0.625)

我知道如何做到这一点的唯一方法是使用从 1 到 k（或者在我的示例中为 1 到 6）的 for 循环，找到每个标签的点的索引，并平均它们的坐标通过索引图像的网格。

但是，我希望以一种针对 GPU 计算进行优化的方式来执行此操作。因此，for 循环的使用不太理想（对于几百个标签，在一个不错的 GPU 上，每张图像大约需要 1 秒）。我正在使用 PyTorch，但实际上任何 numpy 解决方案都应该足够了。

此任务是否有 GPU 高效的解决方案？

【问题讨论】：

标签： numpy image-processing gpu image-segmentation pytorch

【解决方案1】：

这种计算需要积累，我不知道在 GPU 上的效率如何。这是伪代码中的顺序算法：

int n[k] = 0
int sx[k] = 0
int sy[k] = 0
loop over y:
   loop over x:
      i = img[x,y]
      ++n[i]
      sx[i] += x
      sy[i] += y
for i = 1 to k
    sx[i] /= n[i]
    sy[i] /= n[i]

当然，(sx[i],sy[i]) 是对象i 的质心。

在 CPU 上确实很快，不值得为此将数据发送到 GPU，除非它已经存在。

【讨论】：

是的，这也是我能想到的最佳解决方案。数据已经在 GPU 上，因此我想在那里做：-/。我看看CPU是怎么走的。这肯定和它一样高效——每个像素被触摸一次，但我想知道是否有办法在 GPU 上并行化它......
如果上述算法的朴素 GPU 实现太慢，您可以做的一件事是为每个图像列创建数组 n、sx 和 sy，然后添加它们以后一起。这可能会减少等待对数组值进行原子更新的核心数量。
我不确定，以原子方式更新全局数组中的单个值是否有效？您是否在 GPU 上遇到错误共享问题？

【解决方案2】：

考虑使用scikit-image 或重用他们的code（基于numpy/scipy）。

这是一个演示：

import numpy as np
from skimage import measure
from time import perf_counter as pc

img = np.array([[1, 1, 2, 2, 2, 2, 2, 3, 3],
                [1, 1, 1, 2, 2, 2, 2, 3, 3],
                [1, 1, 2, 2, 2, 2, 3, 3, 3],
                [1, 4, 4, 4, 2, 2, 2, 2, 3],
                [4, 4, 4, 5, 5, 5, 2, 3, 3],
                [4, 4, 4, 5, 5, 6, 6, 6, 6],
                [4, 4, 5, 5, 5, 5, 6, 6, 6]])

# assuming already labels of 1, 2, ... n
times = [pc()]
props = measure.regionprops(img)
times.append(pc())
for i in range(np.unique(img).shape[0]):
    print(props[i].centroid)
    times.append(pc())

print(np.diff(times))

输出：

(1.25, 0.625)
(1.5555555555555556, 4.4444444444444446)
(1.8999999999999999, 7.4000000000000004)
(4.3636363636363633, 1.1818181818181819)
(5.1111111111111107, 3.6666666666666665)
(5.4285714285714288, 6.7142857142857144)
[  9.05569615e-05   8.51235438e-04   2.48126075e-04   2.59294767e-04
   2.42692657e-04   2.00734598e-04   2.34542530e-04]

【讨论】：

【解决方案3】：

一个想法是使用bincount 来累积每个区域的行和列索引，使用输入数组中的数字作为箱，因此有一个矢量化解决方案，就像这样 -

m,n = a.shape
r,c = np.mgrid[:m,:n]
count = np.bincount(a.ravel())
centroid_row = np.bincount(a.ravel(),r.ravel())/count
centroid_col = np.bincount(a.ravel(),c.ravel())/count

示例运行 -

In [77]: a
Out[77]: 
array([[1, 1, 2, 2, 2, 2, 2, 3, 3],
       [1, 1, 1, 2, 2, 2, 2, 3, 3],
       [1, 1, 2, 2, 2, 2, 3, 3, 3],
       [1, 4, 4, 4, 2, 2, 2, 2, 3],
       [4, 4, 4, 5, 5, 5, 2, 3, 3],
       [4, 4, 4, 5, 5, 6, 6, 6, 6],
       [4, 4, 5, 5, 5, 5, 6, 6, 6]])

In [78]: np.c_[centroid_row, centroid_col]
Out[78]: 
array([[  nan,   nan], 
       [ 1.25,  0.62], # centroid for region-1
       [ 1.56,  4.44], # centroid for region-2
       [ 1.9 ,  7.4 ], # centroid for region-3 and so on.
       [ 4.36,  1.18],
       [ 5.11,  3.67],
       [ 5.43,  6.71]])

【讨论】：