6

我保存了一个 120 GB 的文件（通过二进制格式pickle），其中包含大约 50,000 个（600x600）2d numpy 数组。我需要使用中位数堆叠所有这些数组。最简单的方法是简单地将整个文件作为数组列表读取并使用np.median(arrays, axis=0). 但是，我没有太多的 RAM 可以使用，所以这不是一个好的选择。

因此，我尝试逐个像素地堆叠它们，因为我一次只关注一个像素位置(i, j)，然后逐个读取每个数组，将给定位置的值附加到列表中。一旦保存了所有数组中某个位置的所有值，我使用np.median然后只需将该值保存在一个列表中——最终将具有每个像素位置的中值。最后，我可以将其重塑为 600x600，我就完成了。代码如下。

import pickle
import time
import numpy as np

filename = 'images.dat' #contains my 50,000 2D numpy arrays

def stack_by_pixel(i, j):
    pixels_at_position = []
    with open(filename, 'rb') as f:
        while True:
            try:
                # Gather pixels at a given position
                array = pickle.load(f)
                pixels_at_position.append(array[i][j])
            except EOFError:
                break
    # Stacking at position (median)
    stacked_at_position = np.median(np.array(pixels_at_position))
    return stacked_at_position

# Form whole stacked image
stacked = []
for i in range(600):
    for j in range(600):
        t1 = time.time()
        stacked.append(stack_by_pixel(i, j))
        t2 = time.time()
        print('Done with element %d, %d: %f seconds' % (i, j, (t2-t1)))

stacked_image = np.reshape(stacked, (600,600))

在看到一些时间打印输出后，我意识到这是非常低效的。每次完成一个位置(i, j)大约需要 150 秒左右，这并不奇怪，因为它正在一个一个地读取大约 50,000 个数组。鉴于我的大型阵列中有 360,000(i, j)个位置，预计这需要 22 个月才能完成！显然这是不可行的。但我有点不知所措，因为没有足够的 RAM 可以读取整个文件。或者，也许我可以一次保存数组的所有像素位置（每个位置的单独列表），因为它一个一个地打开它们，但不会在 Python 中保存 360,000 个列表（大约 50,000 个元素长）使用很多RAM 也是？

欢迎任何关于如何在不使用大量 RAM 的情况下显着加快运行速度的建议。谢谢！

4

2 回答 2

2

这是 numpy 的内存映射数组的完美用例。内存映射数组允许您将.npy磁盘上的文件视为将其作为 numpy 数组加载到内存中，而无需实际加载它。这很简单

arr = np.load('filename', mmap_mode='r')

在大多数情况下，您可以将其视为任何其他数组。数组元素仅根据需要加载到内存中。不幸的是，一些快速实验表明它median不能很好地处理内存映射数组*，它似乎仍然一次将大部分数据加载到内存中。所以median(arr, 0)可能行不通。

但是，您仍然可以遍历每个索引并计算中位数，而不会遇到内存问题。

[[np.median([arr[k][i][j] for k in range(50000)]) for i in range(600)] for j in range(600)]

其中 50,000 反映了数组的总数。

如果没有为了提取单个像素而解开每个文件的开销，运行时间应该会更快（大约 360000 次）。

当然，这就留下了创建.npy包含所有数据的文件的问题。可以按如下方式创建文件，

arr = np.lib.format.open_memmap(
    'filename',              # File to store in
    mode='w+',               # Specify to create the file and write to it
    dtype=float32,           # Change this to your data's type
    shape=(50000, 600, 600)  # Shape of resulting array
)

然后，像以前一样加载数据并将其存储到数组中（它只是在后台将其写入磁盘）。

idx = 0
with open(filename, 'rb') as f:
    while True:
        try:
            arr[idx] = pickle.load(f)
            idx += 1
        except EOFError:
            break

让它运行几个小时，然后回到这个答案的开头，看看如何加载它并取中位数。再简单不过了**。

*我刚刚在一个 7GB 的文件上对其进行了测试，取了 5,000,000 个元素的 1,500 个样本的中位数，内存使用量约为 7GB，这表明整个数组可能已加载到内存中。不过，先尝试这种方式并没有什么坏处。如果其他人对 memmapped 数组的中位数有经验，请随时发表评论。

** 如果您相信互联网上的陌生人。

于 2018-08-30T19:16:29.970 回答
1

注意：我使用 Python 2.x，将其移植到 3.x 应该不难。

我的想法很简单——磁盘空间很充足，所以让我们做一些预处理，把那个大的pickle文件变成更容易处理成小块的东西。

准备

为了测试这一点，我编写了一个小脚本来生成一个类似于你的 pickle 文件。我假设您的输入图像是灰度的并且具有 8 位深度，并使用numpy.random.randint.

该脚本将作为我们可以比较预处理和处理阶段的基准。

import numpy as np
import pickle
import time

IMAGE_WIDTH = 600
IMAGE_HEIGHT = 600
FILE_COUNT = 10000

t1 = time.time()

with open('data/raw_data.pickle', 'wb') as f:
    for i in range(FILE_COUNT):
        data = np.random.randint(256, size=IMAGE_WIDTH*IMAGE_HEIGHT, dtype=np.uint8)
        data = data.reshape(IMAGE_HEIGHT, IMAGE_WIDTH)
        pickle.dump(data, f)
        print i,

t2 = time.time()
print '\nDone in %0.3f seconds' % (t2 - t1)

在测试运行中，此脚本在 372 秒内完成，生成约 10 GB 的文件。

预处理

让我们逐行分割输入图像——我们将有 600 个文件，其中文件N包含N来自每个输入图像的行。我们可以使用二进制存储行数据numpy.ndarray.tofile（然后使用 加载这些文件numpy.fromfile）。

import numpy as np
import pickle
import time

# Increase open file limit
# See https://stackoverflow.com/questions/6774724/why-python-has-limit-for-count-of-file-handles
import win32file
win32file._setmaxstdio(1024)

IMAGE_WIDTH = 600
IMAGE_HEIGHT = 600
FILE_COUNT = 10000

t1 = time.time()

outfiles = []
for i in range(IMAGE_HEIGHT):
    outfilename = 'data/row_%03d.dat' % i
    outfiles.append(open(outfilename, 'wb'))


with open('data/raw_data.pickle', 'rb') as f:
    for i in range(FILE_COUNT):
        data = pickle.load(f)
        for j in range(IMAGE_HEIGHT):
            data[j].tofile(outfiles[j])
        print i,

for i in range(IMAGE_HEIGHT):
    outfiles[i].close()

t2 = time.time()
print '\nDone in %0.3f seconds' % (t2 - t1)

在测试运行中，该脚本在 134 秒内完成，生成 600 个文件，每个文件 600 万字节。它使用了约 30MB 或 RAM。

加工

很简单，只需使用 加载每个数组numpy.fromfile，然后使用numpy.median获取每列的中位数，将其还原为单行，并将这些行累积到列表中。

最后，用于numpy.vstack重新组装中间图像。

import numpy as np
import time

IMAGE_WIDTH = 600
IMAGE_HEIGHT = 600

t1 = time.time()

result_rows = []

for i in range(IMAGE_HEIGHT):
    outfilename = 'data/row_%03d.dat' % i
    data = np.fromfile(outfilename, dtype=np.uint8).reshape(-1, IMAGE_WIDTH)
    median_row = np.median(data, axis=0)
    result_rows.append(median_row)
    print i,

result = np.vstack(result_rows)
print result

t2 = time.time()
print '\nDone in %0.3f seconds' % (t2 - t1)

在测试运行中，此脚本在 74 秒内完成。你甚至可以很容易地并行化它，但这似乎不值得。该脚本使用了约 40MB 的 RAM。

鉴于这两个脚本都是线性的，使用的时间也应该线性扩展。对于 50000 张图像，预处理大约需要 11 分钟，最终处理大约需要 6 分钟。这是在 i7-4930K @ 3.4GHz 上，故意使用 32 位 Python。

于 2018-08-30T19:55:58.520 回答