替代'numpy.tile`用于周期性掩模

Question

我有一个图像,存储在一个uint8形状为n的numpy数组中(planes, rows, cols).我需要将它与存储在掩模(也是uint8s)形状中的值进行比较(mask_rows, mask_cols).虽然图像可能非常大,但是掩模通常很小,通常(256, 256)是平铺的image.为了简化代码,让我们假装rows = 100 * mask_rows和cols = 100 * mask_cols.

我目前处理这个阈值的方式是这样的:

out = image >= np.tile(mask, (image.shape[0], 100, 100))

我可以用这种方式处理的最大阵列,然后用MemoryError比a大一点的方式拍打脸部(3, 11100, 11100).我想的方式,这样做的方式我在内存中共存最多三个巨大的数组:image平铺mask,我的回归out.但是,平铺的掩码是相同的小数组,复制超过10,000次.因此,如果我可以节省内存,我只会使用2/3的内存,并且应该能够处理大于3/2的图像(3, 13600, 13600).顺便说一句,这与我得到的阈值一致

np.greater_equal(image, (image.shape[0], 100, 100), out=image)

我(尝试失败)尝试利用mask处理更大数组的周期性,已经mask使用周期性线性数组进行索引:

mask = mask[None, ...]
rows = np.tile(np.arange(mask.shape[1], (100,))).reshape(1, -1, 1)
cols = np.tile(np.arange(mask.shape[2], (100,))).reshape(1, 1, -1)
out = image >= mask[:, rows, cols]

对于小型阵列,它确实产生与另一个阵列相同的结果,尽管有20倍的减速(!!!),但是对于较大的尺寸来说却非常糟糕.而不是MemoryError它最终崩溃python,即使对于其他方法处理没有问题的值.

我认为正在发生的是numpy实际上是构造(planes, rows, cols)索引的数组mask,因此不仅没有内存保存,而且因为它是一个int32s 数组,它实际上需要四倍的存储空间...

关于如何解决这个问题的任何想法？为了免除麻烦,请在下面找到一些沙盒代码:

import numpy as np

def halftone_1(image, mask) :
    return np.greater_equal(image, np.tile(mask, (image.shape[0], 100, 100)))

def halftone_2(image, mask) :
    mask = mask[None, ...]
    rows = np.tile(np.arange(mask.shape[1]),
                   (100,)).reshape(1, -1, 1)
    cols = np.tile(np.arange(mask.shape[2]),
                   (100,)).reshape(1, 1, -1)
    return np.greater_equal(image, mask[:, rows, cols])

rows, cols, planes = 6000, 6000, 3
image = np.random.randint(-2**31, 2**31 - 1, size=(planes * rows * cols // 4))
image = image.view(dtype='uint8').reshape(planes, rows, cols)
mask = np.random.randint(256,
                         size=(1, rows // 100, cols // 100)).astype('uint8')

#np.all(halftone_1(image, mask) == halftone_2(image, mask))
#halftone_1(image, mask)
#halftone_2(image, mask)

import timeit
print timeit.timeit('halftone_1(image, mask)',
                    'from __main__ import halftone_1, image, mask',
                    number=1)
print timeit.timeit('halftone_2(image, mask)',
                    'from __main__ import halftone_2, image, mask',
                    number=1)

Answer 1

我几乎已经指出了一种滚动窗口类型的技巧,但对于这个简单的非重叠的东西,正常的重塑也是如此.(这里的重塑是安全的,numpy 永远不会为他们复制)

def halftone_reshape(image, mask):
    # you can make up a nicer reshape code maybe, it is a bit ugly. The
    # rolling window code can do this too (but much more general then reshape).
    new_shape = np.array(zip(image.shape, mask.shape))
    new_shape[:,0] /= new_shape[:,1]
    reshaped_image = image.reshape(new_shape.ravel())

    reshaped_mask = mask[None,:,None,:,None,:]

    # and now they just broadcast:
    result_funny_shaped = reshaped_image >= reshaped_mask

    # And you can just reshape it back:
    return result_funny_shaped.reshape(image.shape)

因为时间是一切(不是真的,但是......):

In [172]: %timeit halftone_reshape(image, mask)
1 loops, best of 3: 280 ms per loop

In [173]: %timeit halftone_1(image, mask)
1 loops, best of 3: 354 ms per loop

In [174]: %timeit halftone_2(image, mask)
1 loops, best of 3: 3.1 s per loop