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

Question

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

Answer 1

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

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

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

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

Answer 2

Ответ Fothering heltonbiker, получивший его быстро, будет иметь важное значение. Как и повторная сборка ваших «плиток». Если вы не можете загрузить большой массив в память, вам нужно будет загрузить его как файл с 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. 

теперь process..can продолжить или новый сценарий

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 Извините, я не могу дать больше кода, у меня есть 2-й tiler/reassembler, который я сделал для gis, но это не на этом компьютере. Это может даже не помочь, потому что ваш tiler не вернет фактические плитки, а списки фрагментов, возможно, несколько списков того, где можно взять фрагмент (на большом массиве), где положить его в результаты (массив больших результатов) и где нарезать результат (результат конвейера среза, захваченного из большого массива) среза ... итерации по спискам срезов, и обработка будет простой. Но сделать вашу функцию среза будет сложно.

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 

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

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