Alternativa a `numpy.tile` per la maschera periodica

Question

Ho un'immagine, memorizzata in una serie numerica di uint8 s, di forma (planes, rows, cols). Ho bisogno di confrontarlo con i valori memorizzati in una maschera, anche di uint8 s, di forma (mask_rows, mask_cols). Mentre l'immagine può essere molto grande, la maschera è in genere piccola, normalmente (256, 256) e deve essere affiancata su image. Per semplificare il codice, facciamo finta che sia rows = 100 * mask_rows e cols = 100 * mask_cols.

Il modo in cui sto attualmente la manipolazione di questo thresholding è qualcosa di simile:

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

La più grande gamma posso trattare in questo modo prima di ottenere uno schiaffo in faccia con un MemoryError è un po 'più grande (3, 11100, 11100). Il modo in cui l'ho capito, facendo le cose in questo modo ho fino a tre matrici ginormous che coesistono in memoria: image, il mask piastrellato e il mio ritorno out. Ma la maschera piastrellata è la stessa piccola matrice copiata più di 10.000 volte. Quindi, se potessi risparmiare quella memoria, userei solo 2/3 della memoria, e dovrei essere in grado di elaborare le immagini 3/2 più grandi, quindi delle dimensioni attorno allo (3, 13600, 13600). Questo è, tra l'altro, coerente con quello che ho se faccio il thresholding in posizione con

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

mio (non) tentativo di sfruttare la natura periodica di mask per elaborare le matrici più grandi è stato quello di indice mask con periodica array lineari:

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] 

per piccoli array lo fa lo stesso risultato come l'altra, anche se con una sorta di 20x rallentamento (!!!), ma non riesce a svolgere terribilmente per le taglie più grandi. Invece di un MemoryError alla fine si blocca anche python, anche per valori che l'altro metodo gestisce senza problemi.

Quello che penso sta accadendo è che NumPy è in realtà costruendo la matrice (planes, rows, cols) all'indice mask, è quindi non solo non c'è risparmio di memoria, ma dal momento che è un array di int32 s, in realtà è in corso quattro volte più spazio per store ...

Qualche idea su come procedere? Per risparmiare la fatica, seguito troverete un codice sandbox di giocare con:

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

vorrei quasi vi ho indicato una rolling window tipo di trucco, ma per questa cosa semplice non sovrapposte, rimodellare normale lo fa solo anche. (Le dà nuova forma qui sono al sicuro, sarà NumPy mai fare una copia per loro)

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) 

E poiché tempi sono tutto (non proprio, ma ...):

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