From da707a78edfbb8d1089d7d98a53ffbd33a85f8ab Mon Sep 17 00:00:00 2001 From: Nicolas Poilvert Date: Sun, 12 Feb 2012 00:09:49 -0500 Subject: [PATCH] initial push of array view codes. tests for block views --- skimage/util/__init__.py | 1 + skimage/util/array_views.py | 236 +++++++++++++++++++++++++ skimage/util/tests/test_array_views.py | 67 +++++++ skimage/version.py | 2 +- 4 files changed, 305 insertions(+), 1 deletion(-) create mode 100644 skimage/util/array_views.py create mode 100644 skimage/util/tests/test_array_views.py diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py index 6ef65e20..aa8ead55 100644 --- a/skimage/util/__init__.py +++ b/skimage/util/__init__.py @@ -1,2 +1,3 @@ from .dtype import * +from .array_views import * diff --git a/skimage/util/array_views.py b/skimage/util/array_views.py new file mode 100644 index 00000000..26c5a191 --- /dev/null +++ b/skimage/util/array_views.py @@ -0,0 +1,236 @@ +# 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 + >>> 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]]]]) + """ + + # -- if 'block' is None, we simply return the + # original array. + if block == None: + return arr + # -- otherwise we make sure the user gave a + # tuple + if not isinstance(block, tuple): + raise ValueError('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') + + # -- 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_shape): + """ + 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_shape: 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 + -------- + >>> A = np.arange(10) + >>> A + array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) + >>> window_shape = (3,) + >>> B = rolling_view(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 = rolling_view(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 requirements on inputs + assert isinstance(arr, np.ndarray) + assert isinstance(window_shape, tuple) + assert len(window_shape) == arr.ndim + + # -- input array dimension + N = arr.ndim + + # -- compatibility checks + if ((np.array(arr.shape).astype(int) - \ + np.array(window_shape).astype(int)) < 0).any(): + raise ValueError('window shape is too large') + + if ((np.array(window_shape).astype(int) - \ + np.ones(N).astype(int)) < 0).any(): + raise ValueError('window shape is too small') + + # -- shape of output 'rolling view' array + out_shape = tuple([arr.shape[i] - window_shape[i] + 1 + for i in range(N)]) + \ + window_shape + + # -- 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/tests/test_array_views.py b/skimage/util/tests/test_array_views.py new file mode 100644 index 00000000..0d062f59 --- /dev/null +++ b/skimage/util/tests/test_array_views.py @@ -0,0 +1,67 @@ +import numpy as np +from nose.tools import raises +from numpy.testing import assert_equal +from skimage.util.array_views import block_view + + +@raises(ValueError) +def test_block_view_block_not_a_tuple(): + + A = np.arange(10) + block_view(A, [5]) + + +@raises(ValueError) +def test_block_view_negative_shape(): + + A = np.arange(10) + block_view(A, (-2)) + + +@raises(ValueError) +def test_block_view_block_too_large(): + + A = np.arange(10) + block_view(A, (11,)) + + +@raises(ValueError) +def test_block_view_wrong_block_dimension(): + + A = np.arange(10) + block_view(A, (2,2)) + + +@raises(ValueError) +def test_block_view_1D_array_wrong_block_shape(): + + A = np.arange(10) + block_view(A, (3,)) + + +def test_block_view_1D_array(): + + A = np.arange(10) + B = block_view(A, (5,)) + assert_equal(B, np.array([[0, 1, 2, 3, 4], + [5, 6, 7, 8, 9]])) + + +def test_block_view_2D_array(): + + A = np.arange(4*4).reshape(4,4) + B = block_view(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(): + + A = np.arange(4*4*6).reshape(4,4,6) + B = block_view(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]]], + [[[76, 77], + [82, 83]]]])) diff --git a/skimage/version.py b/skimage/version.py index ab15b0be..d9d6b43c 100644 --- a/skimage/version.py +++ b/skimage/version.py @@ -1,2 +1,2 @@ # THIS FILE IS GENERATED FROM THE SKIMAGE SETUP.PY -version='unbuilt-dev' +version='0.5dev'