使用3x3D数组作为索引的numpy直方图

Question

我有三个3x3的数组，分别代表一张3D RGB图像的红色、绿色和蓝色通道。有没有什么简单的方法可以用numpy来创建这些输入通道的直方图体积呢？

这个操作相当于

""" assume R, G and B are 3D arrays and output is a 3D array filled with zeros """
for x in x_dim:
     for y in y_dim:
          for z in z_dim:
               output[ R[x][y][z] ][ G[x][y][z] ][ B[x][y][z] ] += 1

不过，这段代码在处理大图像时太慢了。numpy能不能提高上面算法的效率呢？

Answer 1

你可以使用 numpy的 histogramdd 来计算一个n维数组的直方图。如果你不想为每个二维切片都生成一个直方图，记得把那个维度的 bins 设置为1。

为了得到整体的直方图，你可以分别计算R、G和B通道的直方图，然后对每个位置取这三个值中的最大值。

Answer 2

假设我们有8位的通道，三个整数(R, G, B)可以看作是一个256进制的数字：R*256**2 + G*256 + B。这样，我们就可以把这三个数组R、G、B合并成一个“颜色值”的数组，然后用np.bincount来生成我们想要的直方图。

import numpy as np

def using_bincount(r,g,b):
    r=r.ravel().astype('int32')
    g=g.ravel().astype('int32')
    b=b.ravel().astype('int32')
    output=np.zeros((base*base*base),dtype='int32')
    result=np.bincount(r*base**2+g*base+b)
    output[:len(result)]+=result
    output=output.reshape((base,base,base))
    return output

def using_histogramdd(r,g,b):
    data = np.vstack((r.flat, g.flat, b.flat)).astype(np.uint8).T
    del(r); del(g); del(b)
    hist, edges = np.histogramdd(
        data, bins=base, range=([0,base],[0,base],[0,base])
        )
    return hist

np.random.seed(0)
n = 200
base = 256
r = np.random.randint(base, size=(n,n,n)).astype(np.uint8)
g = np.random.randint(base, size=(n,n,n)).astype(np.uint8)
b = np.random.randint(base, size=(n,n,n)).astype(np.uint8)

if __name__=='__main__':
    bhist=using_bincount(r,g,b)
    hhist=using_histogramdd(r,g,b)
    assert np.allclose(bhist,hhist)

这些timeit的结果表明，使用_bincount比使用_histogramdd要快，可能是因为_histogramdd是为了处理浮点数和范围的箱子而设计的，而bincount则专门用于计数整数。

% python -mtimeit -s'import test' 'test.using_bincount(test.r,test.g,test.b)'
10 loops, best of 3: 1.07 sec per loop
% python -mtimeit -s'import test' 'test.using_histogramdd(test.r,test.g,test.b)'
10 loops, best of 3: 8.42 sec per loop

Answer 3

你可以使用 numpy.histogramdd 来实现这个功能，不过正如你所说，@jozzas 提出的那种方法是行不通的。你需要做的是把你那三个三维数组“压扁”，然后把它们组合成一个二维数组，尺寸是 (x_dim*y_dim*z_dim, 3)，然后把这个二维数组传给 histogramdd。其实你原来的数据是三维的并不重要，因为在计算直方图时，空间信息并没有用处。

下面是一个使用随机数据的例子：

import numpy 

n = 400  # approximate largest cube size that works on my laptop
# Fill channel cubes with random 8-bit integers
r = numpy.random.randint(256, size=(n,n,n)).astype(numpy.uint8)
g = numpy.random.randint(256, size=(n,n,n)).astype(numpy.uint8)
b = numpy.random.randint(256, size=(n,n,n)).astype(numpy.uint8)

# reorder data into for suitable for histogramming
data = numpy.vstack((r.flat, g.flat, b.flat)).astype(numpy.uint8).T

# Destroy originals to save space
del(r); del(g); del(b)

m = 256                                  # size of 3d histogram cube
hist, edges = numpy.histogramdd(
    data, bins=m, range=((-0.5,255.5),(-0.5,255.5),(-0.5,255.5))
    )

# Check that it worked
assert hist.sum() == n**3, 'Failed to conserve pixels'

需要注意的是，这样做会消耗比你预期的更多内存，因为 histogramdd 似乎在处理时使用的是64位浮点数，即使我们传给它的是8位整数。