From 1108bd892e83a1fa6ea340e98587abaddd7cc10f Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Sun, 12 Feb 2012 14:59:06 -0500 Subject: [PATCH] ENH: move 'util.array_views' to 'util.shape' --- skimage/util/__init__.py | 1 - skimage/util/array_views.py | 231 ------------------ skimage/util/shape.py | 221 +++++++++++++++++ .../{test_array_views.py => test_shape.py} | 62 ++--- 4 files changed, 252 insertions(+), 263 deletions(-) delete mode 100644 skimage/util/array_views.py create mode 100644 skimage/util/shape.py rename skimage/util/tests/{test_array_views.py => test_shape.py} (66%) diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py index d01470bd..980e3880 100644 --- a/skimage/util/__init__.py +++ b/skimage/util/__init__.py @@ -1,2 +1 @@ from .dtype import * -from .array_views import * diff --git a/skimage/util/array_views.py b/skimage/util/array_views.py deleted file mode 100644 index be9072da..00000000 --- a/skimage/util/array_views.py +++ /dev/null @@ -1,231 +0,0 @@ -# Authors: Nicolas Poilvert -# Nicolas Pinto -# License: BSD 3-clause - -__all__ = ['block_view', 'rolling_view'] - -import numpy as np -from numpy.lib.stride_tricks import as_strided as ast - - -def block_view(arr, block): - """Offers a view on array 'arr' which allows one to easily pick a 'block' - and reason within that block when manipulating the array indices. - - Parameters - ---------- - arr: ndarray - input array from which we want to obtain a block view - - block: tuple - each element in the tuple represents the number of input array elements - to include in a block along the corresponding direction - - Returns - ------- - block view on input array. - - Examples - -------- - >>> import numpy as np - >>> from skimage.util.array_views import block_view - >>> A = np.arange(4*4).reshape(4,4) - >>> A - array([[ 0, 1, 2, 3], - [ 4, 5, 6, 7], - [ 8, 9, 10, 11], - [12, 13, 14, 15]]) - >>> B = block_view(A, block=(2,2)) - >>> B[0, 1] - array([[2, 3], - [6, 7]]) - >>> B[1, 0, 1, 1] - 13 - >>> A = np.arange(4*4*6).reshape(4,4,6) - >>> A - array([[[ 0, 1, 2, 3, 4, 5], - [ 6, 7, 8, 9, 10, 11], - [12, 13, 14, 15, 16, 17], - [18, 19, 20, 21, 22, 23]], - - [[24, 25, 26, 27, 28, 29], - [30, 31, 32, 33, 34, 35], - [36, 37, 38, 39, 40, 41], - [42, 43, 44, 45, 46, 47]], - - [[48, 49, 50, 51, 52, 53], - [54, 55, 56, 57, 58, 59], - [60, 61, 62, 63, 64, 65], - [66, 67, 68, 69, 70, 71]], - - [[72, 73, 74, 75, 76, 77], - [78, 79, 80, 81, 82, 83], - [84, 85, 86, 87, 88, 89], - [90, 91, 92, 93, 94, 95]]]) - >>> B = block_view(A, block=(1,2,2)) - >>> B.shape - (4, 2, 3, 1, 2, 2) - >>> B[2:, 0, 2] - array([[[[52, 53], - [58, 59]]], - - - [[[76, 77], - [82, 83]]]]) - """ - - # -- otherwise we make sure the user gave a - # tuple - if not isinstance(block, tuple): - raise TypeError('block needs to be a tuple') - - # -- basic invalid values for 'block' - block_shape = np.array(block).astype(np.int) - if (block_shape <= 0).any(): - raise ValueError('non strictly positive block shape given') - if block_shape.size > arr.ndim: - raise ValueError('block ndim larger than input array ndim') - if block_shape.size < arr.ndim: - raise ValueError('block ndim smaller than input array ndim') - - arr = np.ascontiguousarray(arr) - - # -- checking that the block view is compatible - # with the shape of the input array - A_shape = np.array(arr.shape).astype(np.int) - if (A_shape % block_shape).sum() != 0: - raise ValueError('block shape not compatible with input array') - - # -- actually building the block view - rng = range(len(block)) - shape = ( - tuple([arr.shape[i] / block[i] for i in rng]) - + block - ) - strides = ( - tuple([arr.strides[i] * block[i] for i in rng]) - + arr.strides - ) - - return ast(arr, shape=shape, strides=strides) - - -def rolling_view(arr, window): - """This function offers a 'rolling view' for any N-dimensional array. The - 'window' defines the shape of the elementary N-dimensional orthotope (better - know as hyperrectangle [1]) of the view. - - Parameters - ---------- - arr: ndarray object - N-dimensional input array - - window: N-tuple - tuple of size N that gives the shape of the elementary window - - Returns - ------- - a rolling view on the input array - - Notes - ----- - One should be very careful with rolling views when it comes to memory usage. - Indeed, although a 'view' has the same memory footprint as its base array, - the actual array that emerges when this 'view' is used in a computation is - generally a (much) larger array than the original, especially for - 2-dimensional arrays and above. - - For example, let us consider a 3 dimensional array of size (100, 100, 100) - of ``float64``. This array takes about 8*100**3 Bytes for storage which is - just 8 MB. If one decides to build a rolling view on this array with a - window of (3, 3, 3) the hypothetical size of the rolling view (if one was to - reshape the view for example) would be 8*(100-3+1)**3*3**3 which is about - 203 MB! The scaling becomes even worse as the dimension of the input array - becomes larger. - - References - ---------- - .. [1] http://en.wikipedia.org/wiki/Hyperrectangle - - Examples - -------- - >>> import numpy as np - >>> from skimage.util.array_views import rolling_view - >>> A = np.arange(10) - >>> A - array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) - >>> window = (3,) - >>> B = rolling_view(A, window) - >>> B.shape - (8, 3) - >>> B - array([[0, 1, 2], - [1, 2, 3], - [2, 3, 4], - [3, 4, 5], - [4, 5, 6], - [5, 6, 7], - [6, 7, 8], - [7, 8, 9]]) - >>> A = np.arange(5*4).reshape(5, 4) - >>> A - array([[ 0, 1, 2, 3], - [ 4, 5, 6, 7], - [ 8, 9, 10, 11], - [12, 13, 14, 15], - [16, 17, 18, 19]]) - >>> window = (4, 3) - >>> B = rolling_view(A, window) - >>> B.shape - (2, 2, 4, 3) - >>> B - array([[[[ 0, 1, 2], - [ 4, 5, 6], - [ 8, 9, 10], - [12, 13, 14]], - - [[ 1, 2, 3], - [ 5, 6, 7], - [ 9, 10, 11], - [13, 14, 15]]], - - - [[[ 4, 5, 6], - [ 8, 9, 10], - [12, 13, 14], - [16, 17, 18]], - - [[ 5, 6, 7], - [ 9, 10, 11], - [13, 14, 15], - [17, 18, 19]]]]) - """ - - # -- basic requirements on inputs - if not isinstance(arr, np.ndarray): - raise TypeError('the input should be an ndarray object') - if not isinstance(window, tuple): - raise TypeError('the window shape should be a tuple') - if not (len(window) == arr.ndim): - raise ValueError('array dimension and window length dont match') - - arr = np.ascontiguousarray(arr) - - # -- defining some variables - arr_shape = np.array(arr.shape) - window_shape = np.array(window, dtype=arr_shape.dtype) - - # -- compatibility checks - if ((arr_shape - window_shape) < 0).any(): - raise ValueError("'window_shape' is too large") - - if ((window_shape - 1) < 0).any(): - raise ValueError("'window_shape' is too small") - - # -- shape of output 'rolling view' array - out_shape = tuple(arr_shape - window_shape + 1) + window - - # -- strides of output 'rolling view' array - out_strides = arr.strides + arr.strides - - return ast(arr, shape=out_shape, strides=out_strides) diff --git a/skimage/util/shape.py b/skimage/util/shape.py new file mode 100644 index 00000000..22156e06 --- /dev/null +++ b/skimage/util/shape.py @@ -0,0 +1,221 @@ +# Authors: Nicolas Poilvert +# Nicolas Pinto +# License: BSD 3-clause + +__all__ = ['view_as_blocks', 'view_as_windows'] + +import numpy as np +from numpy.lib.stride_tricks import as_strided + + +def view_as_blocks(arr_in, block_shape): + """Block view of the input n-dimensionaly array (using re-striding). + + Parameters + ---------- + arr: ndarray + The n-dimensional input array. + + block_shape: tuple + The shape of the block. + + Returns + ------- + arr_out: ndarray + Block view of the input array. + + Examples + -------- + >>> import numpy as np + >>> from skimage.util.shape import view_as_blocks + >>> A = np.arange(4*4).reshape(4,4) + >>> A + array([[ 0, 1, 2, 3], + [ 4, 5, 6, 7], + [ 8, 9, 10, 11], + [12, 13, 14, 15]]) + >>> B = view_as_blocks(A, block_shape=(2, 2)) + >>> B[0, 1] + array([[2, 3], + [6, 7]]) + >>> B[1, 0, 1, 1] + 13 + >>> A = np.arange(4*4*6).reshape(4,4,6) + >>> A + array([[[ 0, 1, 2, 3, 4, 5], + [ 6, 7, 8, 9, 10, 11], + [12, 13, 14, 15, 16, 17], + [18, 19, 20, 21, 22, 23]], + + [[24, 25, 26, 27, 28, 29], + [30, 31, 32, 33, 34, 35], + [36, 37, 38, 39, 40, 41], + [42, 43, 44, 45, 46, 47]], + + [[48, 49, 50, 51, 52, 53], + [54, 55, 56, 57, 58, 59], + [60, 61, 62, 63, 64, 65], + [66, 67, 68, 69, 70, 71]], + + [[72, 73, 74, 75, 76, 77], + [78, 79, 80, 81, 82, 83], + [84, 85, 86, 87, 88, 89], + [90, 91, 92, 93, 94, 95]]]) + >>> B = view_as_blocks(A, block_shape=(1, 2, 2)) + >>> B.shape + (4, 2, 3, 1, 2, 2) + >>> B[2:, 0, 2] + array([[[[52, 53], + [58, 59]]], + + + [[[76, 77], + [82, 83]]]]) + """ + + # -- basic checks on arguments + if not isinstance(block_shape, tuple): + raise TypeError('block needs to be a tuple') + + block_shape = np.array(block_shape) + if (block_shape <= 0).any(): + raise ValueError("'block_shape' elements must be strictly positive") + + if block_shape.size != arr_in.ndim: + raise ValueError("'block_shape' must have the same length " + "as 'arr_in.shape'") + + arr_shape = np.array(arr_in.shape) + if (arr_shape % block_shape).sum() != 0: + raise ValueError("'block_shape' is not compatible with 'arr_in'") + + # -- restride the array to build the block view + arr_in = np.ascontiguousarray(arr_in) + + new_shape = tuple(arr_shape / block_shape) + tuple(block_shape) + new_strides = tuple(arr_in.strides * block_shape) + arr_in.strides + + arr_out = as_strided(arr_in, shape=new_shape, strides=new_strides) + + return arr_out + + +def view_as_windows(arr_in, window_shape): + """Rolling window view of the input n-dimensionaly array (using + re-striding). + + + Parameters + ---------- + arr_in: ndarray + The n-dimensional input array. + + window_shape: tuple + Defines the shape of the elementary n-dimensional orthotope + (better know as hyperrectangle [1]) of the rolling window view. + + Returns + ------- + arr_out: ndarray + (rolling) window view of the input array. + + Notes + ----- + One should be very careful with rolling views when it comes to + memory usage. Indeed, although a 'view' has the same memory + footprint as its base array, the actual array that emerges when this + 'view' is used in a computation is generally a (much) larger array + than the original, especially for 2-dimensional arrays and above. + + For example, let us consider a 3 dimensional array of size (100, + 100, 100) of ``float64``. This array takes about 8*100**3 Bytes for + storage which is just 8 MB. If one decides to build a rolling view + on this array with a window of (3, 3, 3) the hypothetical size of + the rolling view (if one was to reshape the view for example) would + be 8*(100-3+1)**3*3**3 which is about 203 MB! The scaling becomes + even worse as the dimension of the input array becomes larger. + + References + ---------- + .. [1] http://en.wikipedia.org/wiki/Hyperrectangle + + Examples + -------- + >>> import numpy as np + >>> from skimage.util.shape import view_as_windows + >>> A = np.arange(10) + >>> A + array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) + >>> window_shape = (3,) + >>> B = view_as_windows(A, window_shape) + >>> B.shape + (8, 3) + >>> B + array([[0, 1, 2], + [1, 2, 3], + [2, 3, 4], + [3, 4, 5], + [4, 5, 6], + [5, 6, 7], + [6, 7, 8], + [7, 8, 9]]) + >>> A = np.arange(5*4).reshape(5, 4) + >>> A + array([[ 0, 1, 2, 3], + [ 4, 5, 6, 7], + [ 8, 9, 10, 11], + [12, 13, 14, 15], + [16, 17, 18, 19]]) + >>> window_shape = (4, 3) + >>> B = view_as_windows(A, window_shape) + >>> B.shape + (2, 2, 4, 3) + >>> B + array([[[[ 0, 1, 2], + [ 4, 5, 6], + [ 8, 9, 10], + [12, 13, 14]], + + [[ 1, 2, 3], + [ 5, 6, 7], + [ 9, 10, 11], + [13, 14, 15]]], + + + [[[ 4, 5, 6], + [ 8, 9, 10], + [12, 13, 14], + [16, 17, 18]], + + [[ 5, 6, 7], + [ 9, 10, 11], + [13, 14, 15], + [17, 18, 19]]]]) + """ + + # -- basic checks on arguments + if not isinstance(arr_in, np.ndarray): + raise TypeError("'arr_in' must be a numpy ndarray") + if not isinstance(window_shape, tuple): + raise TypeError("'window_shape' must be a tuple") + if not (len(window_shape) == arr_in.ndim): + raise ValueError("'window_shape' is incompatible with 'arr_in.shape'") + + arr_shape = np.array(arr_in.shape) + window_shape = np.array(window_shape, dtype=arr_shape.dtype) + + if ((arr_shape - window_shape) < 0).any(): + raise ValueError("'window_shape' is too large") + + if ((window_shape - 1) < 0).any(): + raise ValueError("'window_shape' is too small") + + # -- build rolling window view + arr_in = np.ascontiguousarray(arr_in) + + new_shape = tuple(arr_shape - window_shape + 1) + tuple(window_shape) + new_strides = arr_in.strides + arr_in.strides + + arr_out = as_strided(arr_in, shape=new_shape, strides=new_strides) + + return arr_out diff --git a/skimage/util/tests/test_array_views.py b/skimage/util/tests/test_shape.py similarity index 66% rename from skimage/util/tests/test_array_views.py rename to skimage/util/tests/test_shape.py index c0f97b80..127487f1 100644 --- a/skimage/util/tests/test_array_views.py +++ b/skimage/util/tests/test_shape.py @@ -1,65 +1,65 @@ import numpy as np from nose.tools import raises from numpy.testing import assert_equal -from skimage.util.array_views import block_view, rolling_view +from skimage.util.shape import view_as_blocks, view_as_windows @raises(TypeError) -def test_block_view_block_not_a_tuple(): +def test_view_as_blocks_block_not_a_tuple(): A = np.arange(10) - block_view(A, [5]) + view_as_blocks(A, [5]) @raises(ValueError) -def test_block_view_negative_shape(): +def test_view_as_blocks_negative_shape(): A = np.arange(10) - block_view(A, (-2,)) + view_as_blocks(A, (-2,)) @raises(ValueError) -def test_block_view_block_too_large(): +def test_view_as_blocks_block_too_large(): A = np.arange(10) - block_view(A, (11,)) + view_as_blocks(A, (11,)) @raises(ValueError) -def test_block_view_wrong_block_dimension(): +def test_view_as_blocks_wrong_block_dimension(): A = np.arange(10) - block_view(A, (2,2)) + view_as_blocks(A, (2,2)) @raises(ValueError) -def test_block_view_1D_array_wrong_block_shape(): +def test_view_as_blocks_1D_array_wrong_block_shape(): A = np.arange(10) - block_view(A, (3,)) + view_as_blocks(A, (3,)) -def test_block_view_1D_array(): +def test_view_as_blocks_1D_array(): A = np.arange(10) - B = block_view(A, (5,)) + B = view_as_blocks(A, (5,)) assert_equal(B, np.array([[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]])) -def test_block_view_2D_array(): +def test_view_as_blocks_2D_array(): A = np.arange(4*4).reshape(4,4) - B = block_view(A, (2,2)) + B = view_as_blocks(A, (2,2)) assert_equal(B[0,1], np.array([[2, 3], [6, 7]])) assert_equal(B[1, 0, 1, 1], 13) -def test_block_view_3D_array(): +def test_view_as_blocks_3D_array(): A = np.arange(4*4*6).reshape(4,4,6) - B = block_view(A, (1,2,2)) + B = view_as_blocks(A, (1,2,2)) assert_equal(B.shape, (4, 2, 3, 1, 2, 2)) assert_equal(B[2:, 0, 2], np.array([[[[52, 53], [58, 59]]], @@ -68,45 +68,45 @@ def test_block_view_3D_array(): @raises(TypeError) -def test_rolling_view_input_not_array(): +def test_view_as_windows_input_not_array(): A = [1, 2, 3, 4, 5] - rolling_view(A, (2,)) + view_as_windows(A, (2,)) @raises(TypeError) -def test_rolling_view_window_not_tuple(): +def test_view_as_windows_window_not_tuple(): A = np.arange(10) - rolling_view(A, [2]) + view_as_windows(A, [2]) @raises(ValueError) -def test_rolling_view_wrong_window_dimension(): +def test_view_as_windows_wrong_window_dimension(): A = np.arange(10) - rolling_view(A, (2,2)) + view_as_windows(A, (2,2)) @raises(ValueError) -def test_rolling_view_negative_window_length(): +def test_view_as_windows_negative_window_length(): A = np.arange(10) - rolling_view(A, (-1,)) + view_as_windows(A, (-1,)) @raises(ValueError) -def test_rolling_view_window_too_large(): +def test_view_as_windows_window_too_large(): A = np.arange(10) - rolling_view(A, (11,)) + view_as_windows(A, (11,)) -def test_rolling_view_1D(): +def test_view_as_windows_1D(): A = np.arange(10) window_shape = (3,) - B = rolling_view(A, window_shape) + B = view_as_windows(A, window_shape) assert_equal(B, np.array([[0, 1, 2], [1, 2, 3], [2, 3, 4], @@ -117,11 +117,11 @@ def test_rolling_view_1D(): [7, 8, 9]])) -def test_rolling_view_2D(): +def test_view_as_windows_2D(): A = np.arange(5*4).reshape(5, 4) window_shape = (4, 3) - B = rolling_view(A, window_shape) + B = view_as_windows(A, window_shape) assert_equal(B.shape, (2, 2, 4, 3)) assert_equal(B, np.array([[[[ 0, 1, 2], [ 4, 5, 6],