Как сделать эффективную для памяти двумерную свертку на больших массивах

Question

У меня есть проблема, когда мне нужно свернуть один очень большой 2D-массив (файл на диске) с меньшим массивом, который помещается в памяти. scipy.signal.fftconvolve замечательно, когда массивы умещаются в памяти, но не помогают, когда нет. Есть ли какой-то другой разумный подход, кроме зацикливания всех точек в каждом массиве для вычисления свертки вручную? Я не очень хорошо с математикой, но мне интересно, если fftconvolve может быть разделен на части и объединен с небольшим перекрытием? Что-то другое?

Answer 1

Я могу предложить вам два разных подхода (хотя я не рискну предоставить пример кода, надеюсь, вы не будете против его выяснения):

1) Использование numpy.memmap; " Файлы с отображением в памяти используются для доступа к небольшим сегментам больших файлов на диске без чтения всего файла в память. (...) Объект memmap можно использовать везде, где принят ndarray ".

2) Разделите большой массив на плитки, выполните свертку с mode='full' и наложить результаты. Для каждой плитки вы получите "рамку" вокруг плитки с такой же шириной ядра свертки.

Можно объединить оба подхода (например, прочитать плитки из файла memmapped и наложить результаты в другой файл memmapped, который является результатом).

Answer 2

Дальнейшее реагирование хелтонбайкеров и их быстрое выполнение будут важны. Как будет повторная сборка ваших "плиток". Если вы не можете загрузить свой большой массив в память, вам нужно загрузить его как файл memmapped... но чтобы иметь его как файл memmapped, вам нужно сначала создать его.

Вот грубый псевдокод.

#you need to know how big your data matrix is

#todo: assign nrows and ncols based on your data.

fp = np.memmap('your_new_memmap_filename.dat', dtype='float32', mode='w+', shape=(nrows, ncols))#create file with appropriate dimensions

data = open('yourdatafile.txt', 'r')
i = 0
for line in data:
    arr = map(float, line.strip().split(','))   #comma delimited? header row?
    fp[i, :] = arr
    i += 1

del fp  #see documentation...del will flush the output and close the file.

сейчас обрабатывать.. можно продолжить или новый скрипт

convolve_matrix = somenumpyarray
fp_result = np.memmap('your_results_memmap_filename.dat', dtype='float32', mode='w+', shape=(nrows, ncols))

#open big file read only
fp = np.memmap('your_new_memmap_filename.dat', dtype='float32', mode='r', shape=(nrows, ncols))

chunksize = 10000 #?
for i in range(int(nrows/chunksize) - 1):  #don't forget the remainder at the end
    chunk = fp[i * chunksize: (i + 1) * chunksize, :]
    res = fftconvolve(chunk, convolve_matrix)
    fp_result[i * chunksize: (i + 1) * chunksize, :] = res

#deal with remainder

del fp_result
del fp

обратите внимание, что этот псевдокод не перекрывается, и вам необходимо заполнить некоторые пробелы. Кроме того, как только вы получите работу с плиткой, убедитесь, что вы используете Joblib и обрабатываете плитки параллельно. https://pythonhosted.org/joblib/parallel.html Извините, я не могу дать больше кода, у меня есть двумерный сборщик / сборщик, который я сделал для ГИС, но его нет на этом компьютере. Это может даже не сильно помочь, потому что ваш тайлер не будет возвращать фактические плитки, но списки срезов, возможно, несколько списков того, где взять фрагмент (в большом массиве), где поместить его в результаты (большой массив результатов) и где вырезать. результат (убедительный результат среза, извлеченного из большого массива) среза... перебирает списки срезов, и обработка будет легкой. Но сделать вашу функцию среза будет сложно.

for source_slice, result_slice, mini_slice in zip(source_slice, result_slice, mini_slice):
    matrix2convolve = big_fp[source_slice[0]:source_slice[1], :]
    convolve_result = fftconvolve(matrix2convolve, convolve_matrix)
    big_result_fp[result_slice[0]:result_slice[1], :] = convolve_result[mini_slice[0]:mini_slice[1], :]

Answer 3

Звучит наивно, но на 100% верно, если вы переместите пиксель результатов.

Свертка, как это обычно делается, неверна. Фактическая операция не требует памяти вообще. Вы можете просто прочитать файл, выполнить свертку и записать его обратно в то же место.

Недостаток алгоритма в том, что он выполняется, заключается в том, что он требует наличия центрального пикселя. Если вы предпочитаете размещать пиксель результатов в верхнем левом углу, выполните свертку. Вы можете просто прочитать свертку как операцию сканирования. Поскольку ни одна из точек вниз или вправо не была изменена этой операцией, ваш результат верный.

Как только эта совершенно произвольная и бессмысленная зависимость от предыдущих пикселей будет нарушена, вы можете делать классные вещи, такие как объединение ядер сверток и выполнение операции в пределах того же объема памяти (или файла), из которого он в данный момент читается. http://godsnotwheregodsnot.blogspot.com/2015/02/combining-convolution-kernels.html

Смотрите источник здесь http://pastebin.com/bk0A2Z5D

Причина, по которой пиксели находятся в центре, заключается в том, что они не двигаются. Но это того не стоит. Вы действительно можете переместить их обратно в нужное место в той же памяти, если вы были настроены на то, чтобы вокруг результата было одинаковое количество мусора. На самом деле, если вы выполнили следующую свертку с результатом в правом нижнем углу и перебрали массив в обратном направлении, вы бы получили результаты обратно там, где они более или менее начинались, или с ядром:

0,0
0,1

Короче говоря, единственная причина вашей проблемы в том, что кто-то некоторое время назад решил, что должен быть центральный пиксель. Когда вы отбрасываете эту идиотскую идею, вы всегда знаете, где находится пиксель результатов, вы можете выполнять забавные операции, такие как прекомпиляция матриц свертки, выполняя, таким образом, все свертки одновременно с комбинированными ядрами и, для вас, выполнять операцию исключительно путем чтения с диска. Подумайте только, вы могли бы стать первым человеком, получившим размытую копию этого изображения НАСА Галактики Андромеды.

Если кто-нибудь знает историю того, кто первым принял решение о создании центрального пикселя, я бы хотел знать. Потому что без него свёртка - это только операция пиксельной развертки.