14

编辑我在下面保留了我面临的更复杂的问题,但我的问题np.take可以更好地总结如下。假设您有一个img形状数组(planes, rows),还有另一个lut形状数组(planes, 256),并且您想使用它们来创建一个新out的形状数组(planes, rows),其中out[p,j] = lut[p, img[p, j]]. 这可以通过花哨的索引来实现,如下所示:

In [4]: %timeit lut[np.arange(planes).reshape(-1, 1), img]
1000 loops, best of 3: 471 us per loop

但是,如果不是花哨的索引,而是使用 take 和 python 循环,则planes可以大大加快速度:

In [6]: %timeit for _ in (lut[j].take(img[j]) for j in xrange(planes)) : pass
10000 loops, best of 3: 59 us per loop

可以lut并且img以某种方式重新排列,以便在没有 python 循环的情况下进行整个操作,而是使用numpy.take(或替代方法)而不是传统的花式索引来保持速度优势?


原始问题 我有一组要在图像上使用的查找表 (LUT)。保存 LUT 的数组是有形的(planes, 256, n),而图像是有形的(planes, rows, cols)。两者都是,与LUTdtype = 'uint8'的轴相匹配。256这个想法是从LUT的第p-平面运行图像的第-平面通过每个LUT。np

如果我的lutimg是以下:

planes, rows, cols, n = 3, 4000, 4000, 4
lut = np.random.randint(-2**31, 2**31 - 1,
                        size=(planes * 256 * n // 4,)).view('uint8')
lut = lut.reshape(planes, 256, n)
img = np.random.randint(-2**31, 2**31 - 1,
                    size=(planes * rows * cols // 4,)).view('uint8')
img = img.reshape(planes, rows, cols)

使用像这样的精美索引后,我可以实现我的目标

out = lut[np.arange(planes).reshape(-1, 1, 1), img]

这给了我一个 shape 数组(planes, rows, cols, n),其中out[i, :, :, j]保存了穿过LUT 的 -thi平面的-th LUT 的-th 平面......imgji

一切都很好,除了这个:

In [2]: %timeit lut[np.arange(planes).reshape(-1, 1, 1), img]
1 loops, best of 3: 5.65 s per loop

这是完全不可接受的,特别是因为我有以下所有看起来不太好看的替代品,np.take而不是运行得更快:

  1. 单个平面上的单个 LUT 运行速度大约 x70:

    In [2]: %timeit np.take(lut[0, :, 0], img[0])
    10 loops, best of 3: 78.5 ms per loop
    
  2. 运行所有所需组合的 python 循环几乎快 x6 完成:

    In [2]: %timeit for _ in (np.take(lut[j, :, k], img[j]) for j in xrange(planes) for k in xrange(n)) : pass
    1 loops, best of 3: 947 ms per loop
    
  3. 即使在 LUT 和图像中运行所有平面组合,然后丢弃planes**2 - planes不需要的平面组合,也比花哨的索引更快:

    In [2]: %timeit np.take(lut, img, axis=1)[np.arange(planes), np.arange(planes)]
    1 loops, best of 3: 3.79 s per loop
    
  4. 我能想到的最快的组合是一个 python 循环在平面上迭代并更快地完成 x13:

    In [2]: %timeit for _ in (np.take(lut[j], img[j], axis=0) for j in xrange(planes)) : pass
    1 loops, best of 3: 434 ms per loop
    

当然,问题是如果np.take没有任何 python 循环就没有办法做到这一点?理想情况下,需要的任何重塑或调整大小都应该发生在 LUT 上,而不是图像上,但我对你们能想出的任何事情持开放态度......

4

1 回答 1

6

首先,我不得不说我真的很喜欢你的问题。无需重新排列LUTIMG以下解决方案有效:

%timeit a=np.take(lut, img, axis=1)
# 1 loops, best of 3: 1.93s per loop

但是从结果中你必须查询对角线:a[0,0], a[1,1], a[2,2]; 得到你想要的。我试图找到一种仅对对角线元素进行索引的方法,但仍然没有成功。

以下是重新排列LUTand的一些方法IMG:如果索引IMG为 0-255,第 1 平面,256-511 用于第 2 平面,512-767 用于第 3 平面,则以下方法有效,但这会阻止您using 'uint8',这可能是一个大问题......:

lut2 = lut.reshape(-1,4)
%timeit np.take(lut2,img,axis=0)
# 1 loops, best of 3: 716 ms per loop
# or
%timeit np.take(lut2, img.flatten(), axis=0).reshape(3,4000,4000,4)
# 1 loops, best of 3: 709 ms per loop

在我的机器中,您的解决方案仍然是最佳选择,而且非常足够,因为您只需要对角线评估,即 plane1-plane1、plane2-plane2 和 plane3-plane3:

%timeit for _ in (np.take(lut[j], img[j], axis=0) for j in xrange(planes)) : pass
# 1 loops, best of 3: 677 ms per loop

我希望这可以让您对更好的解决方案有所了解。使用 寻找更多选项以及与或flatten()类似的方法会很不错,这似乎很有希望。np.apply_over_axes()np.apply_along_axis()

我使用下面的代码来生成数据:

import numpy as np
num = 4000
planes, rows, cols, n = 3, num, num, 4
lut = np.random.randint(-2**31, 2**31-1,size=(planes*256*n//4,)).view('uint8')
lut = lut.reshape(planes, 256, n)
img = np.random.randint(-2**31, 2**31-1,size=(planes*rows*cols//4,)).view('uint8')
img = img.reshape(planes, rows, cols)
于 2013-05-09T11:03:05.237 回答