CODEWITHSUNDEEP
pythonmemory-managementimage-processingtreenumpy
Zehan
Zehan

Reputation: 166

Numpy - creating overlapping 3D subarrays as vectors that's memory efficient

I'm trying to create a list of all the overlapping subarrays of equal size from a larger 3D array (for patch-based segmentation) where each subarray needs to be flattened (as a 1D vector) so I can make use of the ball tree in sklearn.neighbours.BallTree.

So for example, given a 100x100x100 image, if I were to break this down into 5x5x5 overlapping patches (subarrays), I would have 96x96x96 = 884,736 of them.

However I have not found any way of doing so without numpy allocating more memory for each flattened/vectorized subarray. This seems to be because each subarray is not contiguous in memory.

e.g. For the 100x100x100 image, if I want each 5x5x5 patch as a 1D vector (of length 125), numpy decides to allocate a brand new array in memory for all 884,736 of them which then becomes rather large especially if I want to work with more than a single 100x100x100 image!

I would welcome any solutions for overcoming this memory challenge in python/numpy. I was considering creating a subclass of the numpy.ndarray object which stores a pointer to the location of the patch in the bigger image and but returns the data as a 1D numpy array only when called (and this is then deleted again when not used) but I have not come across enough details on subclassing ndarray objects to do so. I will be really disappointed if the only solution is to implement everything in C/C++ instead. I appreciate any help that can be provided, thanks!

Upvotes: 2

Views: 760

Answers (1)

Joe Kington
Joe Kington

Reputation: 284830

Based on your question, you may already be aware of all of this. However, I'm posting this "answer" as more of a discussion of what the problems are, because many people may not be aware of them....

If you're not, though, you can make a 96x96x96x5x5x5 array from a 100x100x100 image that acts as a 5x5x5 moving window without allocating any additional memory.

However, because you can only have one stride per dimension, there's no way to reshape this into a 96x96x96x125 array without making a copy.

At any rate, here's an example (basically taken straight from one of my previous answers):

import numpy as np

def rolling_window_lastaxis(a, window):
    """Directly taken from Erik Rigtorp's post to numpy-discussion.
    <http://www.mail-archive.com/[email protected]/msg29450.html>"""
    if window < 1:
       raise ValueError, "`window` must be at least 1."
    if window > a.shape[-1]:
       raise ValueError, "`window` is too long."
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

def rolling_window(a, window):
    """Takes a numpy array *a* and a sequence of (or single) *window* lengths
    and returns a view of *a* that represents a moving window."""
    if not hasattr(window, '__iter__'):
        return rolling_window_lastaxis(a, window)
    for i, win in enumerate(window):
        if win > 1:
            a = a.swapaxes(i, -1)
            a = rolling_window_lastaxis(a, win)
            a = a.swapaxes(-2, i)
    return a

x = np.zeros((100,100,100), dtype=np.uint8)
y = rolling_window(x, (5,5,5))
print 'Now *y* will be a 96x96x96x5x5x5 array...'
print y.shape
print 'Representing a "rolling window" into *x*...'
y[0,0,0,...] = 1
y[1,1,0,...] = 2
print x[:10,:10,0] # Note that *x* and *y* share the same memory!

This yields:

Now *y* will be a 96x96x96x5x5x5 array...
(96, 96, 96, 5, 5, 5)
Representing a "rolling window" into *x*...
[[1 1 1 1 1 0 0 0 0 0]
 [1 2 2 2 2 2 0 0 0 0]
 [1 2 2 2 2 2 0 0 0 0]
 [1 2 2 2 2 2 0 0 0 0]
 [1 2 2 2 2 2 0 0 0 0]
 [0 2 2 2 2 2 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0]]

However, as you've already noted, we can't reshape this to 96x96x96x125 without creating a copy. y.shape = (96,96,96,-1) will raise an error, and z = y.reshape((96,96,96,-1)) will work, but will return a copy.

(The relevant documentation is in numpy.reshape, if this seems confusing. Basically reshape will avoid copying when possible, and return a copy if its not, whereas setting the shape attribute will raise an error when copying isn't possible.)

However, even if you build a more efficient array container, sklearn.neighbors.BallTree is almost definitely going to make temporary intermediate copies.

You mention that you're doing image segmentation. Why not look into a more efficient algorithm than the rather "brute-force" you seem to be trying? (Or if that's not feasible, give us a few more details as to why... Maybe someone will have a better idea?)

Upvotes: 1

Related Questions